Competitive Programming (2) Advent Calendar 2019の22日目の記事です。

さて、もう1年前ですがこんなツイートがありました

RT: 「数学」「パソコン」「ゲーム」を好きな順に並べたときにパソコンが最後に来る人はそういうツール覚えるのが本当に億劫で、それでも戦えるからこそ競プロ続けてるんですよね。パソコンが最後に来ない人は途中で学問や開発とかの方が楽しくなって強くなる前にフェードアウトする人が多い印象。

— えびま (@evima0) November 21, 2018

要するに競プロ勢は「数学」と「ゲーム」が好きだからこそここまで競プロにハマるという主張です。実際に自分も競プロを始めたころを思い出すと、問題を解けるかどうか、解けるならどれだけ早く解けるかを競うという部分に面白さを感じましたし、ある程度競プロに慣れてくると、様々なアルゴリズム、そしてそれを組み合わせる数学的思考力が問われる部分に面白さを感じるようになりました。

さて、そんな数学とゲームが大好きな皆様ならハマること間違いなしのものがあります。それが、

バックギャモン（別名：西洋すごろく）

です。バックギャモンは名前は聞いたことがあってもどういうゲームかは実はよく知らないという人も多いのではないかと思いますが、本当に数学的思考力の重要な、囲碁や将棋とはまた一味違った面白さのあるゲームです。丁度本日行われた競プロ忘年会のLTでも話させて頂いたので、どういった面白さがあるのかはその時使った以下のスライドをご覧ください。

スライドの内容をまとめると、バックギャモンは数学要素の多いにあるゲームです。そして運があるので初心者でも上級者にも勝てうる一方で、ちゃんと実力も競えるし測れるという、日々レーティングを上げようと格闘している競プロ勢の皆様ならのめりこむこと間違いなしのゲームとなっています。

皆さん、是非バックギャモンを始めましょう！

問題

「Aに同じ数が含まれている」「Aの数が全てバラバラ」という2つの場合があるのですが、それぞれそうじゃない方だと解けない解法を使って解くというのが面白いなと。正解者数的に1100点の割にはかなり難しそう。

#include<bits/stdc++.h>
using namespace std;
#define pb push_back
#define pf push_front
typedef long long lint;
typedef complex<double> P;
#define mp make_pair
#define fi first
#define se second
typedef pair<int,int> pint;
#define All(s) s.begin(),s.end()
#define rAll(s) s.rbegin(),s.rend()
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
lint dp[25010][420];
lint d3[25010][420][2];
int a[420];
lint zyo[25010],kai[420],ika[420];
lint mo=1000000007;
map<int,int> me;
lint extgcd(lint a, lint b, lint &x, lint &y) {
  lint g = a; x = 1; y = 0;
  if (b != 0) g = extgcd(b, a % b, y, x), y -= (a / b) * x;
  return g;
}
lint invMod(lint a, lint m) {
  lint x, y;
  if (extgcd(a, m, x, y) == 1) return (x + m) % m;return 0;
}
int cal(void)
{
	int n,K,m,x=0,y=0;lint out=0;
	cin>>n>>K>>m;
	zyo[0]=1;rep(i,n+10) zyo[i+1]=(zyo[i]*K)%mo;
	kai[0]=1;rep(i,410) kai[i+1]=(kai[i]*(i+1))%mo;
	rep(i,410) ika[i]=invMod(kai[i],mo);
	
	rep(i,m) cin>>a[i];
	rep(i,m-K+1){
		int f=0;me.clear();
		rep(j,K){
			if(me[a[i+j]]>0) f=1;me[a[i+j]]++;
		}
		if(f<1){
			return zyo[n-m]*(n-m+1)%mo;
		}
	}
	me.clear();
	while(x<m){
		if(me[a[x]]>0) break;
		me[a[x]]++;x++;
	}
	me.clear();
	while(y<m){
		if(me[a[m-1-y]]>0) break;
		me[a[m-1-y]]++;y++;
	}
	
	memset(dp,0,sizeof(dp));
	memset(d3,0,sizeof(d3));
	rep(i,n) REP(j,1,K){
		dp[i+1][j]+=dp[i+1][j-1];
		if(j==K-1) dp[i+1][j]+=zyo[i];
		else dp[i+1][j]+=(dp[i][j+1]-dp[i][j])*(K-j);
		dp[i+1][j]%=mo;
		dp[i+1][j]+=dp[i][j];
		dp[i+1][j]%=mo;
	}
	if(x>=m && y>=m){
		d3[0][0][0]=1;
		rep(i,n+1) rep(j,K+1) rep(l,2){
			if(j>0 && i>0){
				d3[i][j][l]+=d3[i][j-1][l];
			}
			d3[i][j][l]%=mo;d3[i][j][l]+=mo;d3[i][j][l]%=mo;
			d3[i+1][1][l]+=d3[i][j][l];
			d3[i+1][j+1][l]+=mo-d3[i][j][l];
			if(j<K-1){
				d3[i+1][j+1][l]+=d3[i][j][l]*(K-j);
				d3[i+1][j+2][l]-=d3[i][j][l]*(K-j);
			}
			else if(j==K-1){
				d3[i+1][j+1][1]+=d3[i][j][l];
				d3[i+1][j+2][1]-=d3[i][j][l];
			}
		}
		rep(i,n){
			REP(j,m,K+1){
				//既にできてるやつ
				out+=d3[i][j][1]*zyo[n-i];
				out%=mo;
				//これから足すやつ。
				out+=d3[i][j][0]*(dp[n-i][j]-dp[n-i][j-1]);
				out%=mo;
			}
		}
		REP(j,m,K+1){
			out+=d3[n][j][1];out%=mo;
		}
		out*=kai[K-m];out%=mo;out*=ika[K];out%=mo;out+=mo;out%=mo;
		return out;
	}
	out=zyo[n-m]*(n-m+1)%mo;
	rep(i,n-m+1){
		out-=(zyo[i]-dp[i][x]+dp[i][x-1])*(zyo[n-m-i]-dp[n-m-i][y]+dp[n-m-i][y-1]);
		out%=mo;
	}
	out+=mo;out%=mo;
	return out;
}
int main()
{
	cout<<cal()<<endl;
}

問題

解説によるとO(N2^N)でできるらしいけどO(N^2 2^N)でも通った。こっちは逆に2300点という点数の割にはかなり簡単に思える。

#include<bits/stdc++.h>
using namespace std;
#define pb push_back
#define pf push_front
typedef long long lint;
typedef complex<double> P;
#define mp make_pair
#define fi first
#define se second
typedef pair<int,int> pint;
#define All(s) s.begin(),s.end()
#define rAll(s) s.rbegin(),s.rend()
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
int dp[(1<<21)+10][22];
int main()
{
	memset(dp,0,sizeof(dp));
	int n,K;cin>>n>>K;string s;
	rep(i,n+1){
		cin>>s;
		rep(j,(1<<i)){
			if(s[j]=='1') dp[(1<<i)+j][i]++;
		}
	}
	for(int i=n;i>0;i--) rep(j,(1<<i)) rep(l,i+1){
		if(dp[(1<<i)+j][l]<1) continue;
		rep(k,l){
			if(k+1<i && ((j&(1<<k))>>k)==((j&(1<<(k+1)))>>(k+1))) continue;
			int to=(j>>(k+1)),bo=(j&((1<<k)-1));
			dp[(1<<(i-1))+(to<<k)+bo][k]+=dp[(1<<i)+j][l];
		}
	}
	for(int i=n;i>0;i--) rep(j,(1<<i)){
		int sum=0;
		rep(k,i+1) sum+=dp[(1<<i)+j][k];
		if(sum>=K){
			string out="";
			rep(l,i){
				out+=('0'+(j%2));j/=2;
			}
			reverse(All(out));
			cout<<out<<endl;return 0;
		}
	}
	cout<<""<<endl;
}

Competitive Programming (1) Advent Calendar 2018の20日目の記事です。この記事では、最近流行り（？）のConvexHullTrickを少し別の角度から見て行こうと思います。

では、まずこの問題を見て見ましょう

CF189DIV1C Kalila and Dimna in the Logging Industry

問題概要をきわめてざっくり説明すると、長さn(<=10^5)の2つの数列a,bが与えられて、dp[i]=min(j<i)dp[j]+a[i]*b[j]というDPを解けばいいです。ただし、aは狭義単調増加でbは狭義単調減少です。</ppp>

ConvexHullTrickによる解法

この問題は、2次元平面上に直線がどんどん追加されていき、あるx座標においてy座標が最も小さくなる直線を求めるということを繰り返すという問題ということができます。そして、「追加クエリにおける直線の傾きが単調」「最小値クエリにおけるxが単調」の2つを満たしているので、ConvexHullTrick（最近はCHTと略されることが多いらしいですね･･･最初に見たときは何のことか全く分かりませんでした）を用いてO(n)で解くことができます。ConvexHullTrickの詳しいことは例えばこちらの記事をご覧ください。ConvexHullTrickを用いた解答コードは以下のようになります。

#include<bits/stdc++.h>
using namespace std;
typedef long long lint;
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
lint a[100010],b[100010],dp[100010];
int deq[100010];
lint cal(int i,int j){return dp[j]+b[j]*a[i];}
double cro(int i,int j){return (0.0+dp[j]-dp[i])/(b[i]-b[j]);}
int main()
{
	int s=0,t=1,n;
	cin>>n;
	rep(i,n) scanf("%I64d",&a[i]);
	rep(i,n) scanf("%I64d",&b[i]);
	dp[0]=0;deq[0]=0;
	REP(i,1,n){
		while(s+1<t && cal(i,deq[s])>cal(i,deq[s+1])) s++;
		dp[i]=cal(i,deq[s]);
		while(s+1<t && cro(deq[t-1],i)<=cro(deq[t-1],deq[t-2])) t--;
		deq[t++]=i;
	}
	cout<<dp[n-1]<<endl;
}

これ、コードの行数は少なくて極めてシンプルなコードに見えますよね？でも自分はこのコードを書く時に意外と苦労した記憶があります。配列を使ってdequeの中身を操作するのですが、どのタイミングでiteratorの値を動かせばいいのか、中身が十分に入ってなくて値の比較ができないというのはどこで判定すればいいのか･･･などなど、コードにも-1や-2をどこで入れるのかなど細かいことに気を使う印象です。

分割統治によって最小値たちを効率よく求める方法

さて、今の問題はいったん置いといて次の問題について考えてみましょう

n×nの2次元配列cが与えられる。各行ごとに最小値をとる列を全て求めよ。ただし、各行ごとに最小値をとる列の場所は広義単調増加になっている

つまり各行ごとに最小値をとる列を赤く塗るとこのようになっているということですね。

これはもちろん愚直にやればO(n^2)ですが、この性質を利用してより効率的に求めることができます。

まず、真ん中の行の最小値をとる列の場所を調べます。そうすると単調性より、その場所の（図の上での）右上と左下には最小値をとる場所は存在しないことが分かります。

次に、左上と右下のブロックそれぞれにおいて同様に真ん中の行の最小値をとる場所を調べて･･･ということを再帰的に繰り返すと、全体でO(nlogn)で求めることができます。より詳しい解説はこちらの記事をご覧下さい。

これを使うと、例えば2次元DPで、配列dp[i]の各値が配列dp[i-1]の各値から遷移され、かつその遷移元の位置が遷移先の位置に対して広義単調であることが分かる場合（図で表すと下のようになります。これを私は勝手に「遷移が交差しない」状態と呼んでいます。）、dp[i-1]からdp[i]への遷移は愚直にやればO(n^2)のところを上に述べたテクニックを使うことでO(nlogn)で全て行うことができます。

ConvexHullTrickの代わりに分割統治を使う

で、この分割統治法がConvexHullTrickとどういう関係があるのかというと、「追加クエリにおける直線の傾きが単調」「最小値クエリにおけるxが単調」の2つを満たしている場合、クエリに対して最小値をとる直線が追加された順番は単調になるという性質が成り立つからです。実際にdequeに直線が追加されたり捨てられたりする順番を考えれば成り立つのが分かると思います。また今回の問題に関しても、dp[j]+a[i]*b[j]が最小となるjがiに対して広義単調増加になる(あるi=i1に対して最小となるjをj1とした時、i=i2>i1に対して最小となるjがj1より小さいはずがない)のは、aが狭義単調増加でbが狭義単調減少なことから分かります。よって、前の節で述べた「遷移が交差しない」という条件を満たします。

しかし、今回の問題では遷移前の状態というのがまさに遷移によって求まるものなので、前節のような方法で求めることはできません。では駄目かというと、ここでもまた分割統治法を使うことで解決できます。

さらなる分割統治法で最適解を求める

結論から言うと、1次元DPで「遷移が交差しない」という条件が満たされている場合、以下のアルゴリズムを用いれば範囲[0,n)の全ての値の最適解を求めることができます。

１．範囲[0,m)の最適解を求める

２．範囲[m,n)の範囲[0,m)から遷移される部分に関する最適解を求める

３．範囲[m,n)の最適解を求める

ここでm=n/2とします。手順1,3に関しては再帰的に行い、手順2に関しては上に述べた分割統治法を行います。このアルゴリズムの詳しい解説に関してはこちらの記事をご覧下さい。手順2の計算量がO(nlogn)である場合、全体でO(nlog^2n)で全ての値の最適解を求めることができます。このアルゴリズムを用いた最初の問題の解答コードは以下のようになります。

#include<bits/stdc++.h>
using namespace std;
typedef long long lint;
#define REP(i,a,b) for(int i=a;i<b;i++)
#define rep(i,n) REP(i,0,n)
lint inf=1e18;
lint a[100010],b[100010],dp[100010];
lint cal(int i,int j){return dp[j]+b[j]*a[i];}
void cal2(int l1,int h1,int l2,int h2){
	if(l1>=h1) return;
	int mi=(l1+h1)/2,it=l2;
	REP(i,l2,h2){
		if(cal(mi,i)<cal(mi,it)) it=i;
	}
	dp[mi]=min(dp[mi],cal(mi,it));
	cal2(l1,mi,l2,it+1);
	cal2(mi+1,h1,it,h2);
}
void rec(int lo,int hi){
	if(lo+2>hi) return;
	int mi=(lo+hi)/2;
	rec(lo,mi);
	cal2(mi,hi,lo,mi);
	rec(mi,hi);
}
int main()
{
	int s=0,t=1,n;
	cin>>n;
	rep(i,n) scanf("%I64d",&a[i]);
	rep(i,n) scanf("%I64d",&b[i]);
	dp[0]=0;
	REP(i,1,n+10) dp[i]=inf;
	rec(0,n);
	cout<<dp[n-1]<<endl;
}

解答の長さ自体は最初のより長いですが、本当に値を渡して、その範囲の探索を行うだけなので少なくとも個人的にはこちらのコードの方が頭を使わずに書けました。また計算量もO(n)からO(nlog^2n)へと増えているのですが、感覚としてこのアルゴリズムは定数倍が極めて早く、O(nlogn)と同等か下手したらそれより早いイメージがあります。実際、実行時間も最初のコードが78msなのに対しこのコードは108msで済んでいます。

もちろん、この方法はConvexHullTrickが使える場合に限らず、遷移元と遷移先の関係に関して単調性が満たされてさえいれば使えるのですが、そのような単調性を満たす自然な設定にしようとすると大体ConvexHullTrickも使える形になってる場合が多いと思います。なので、単調性を利用した解法を思いついたがそれだとTLEするといった場合は、式変形などを経る事でConvexHullTrickが使える形に帰着できないか考えてみるといいと思います。

最後に

結局話としては想定解よりも計算量の大きい解法を説明しただけで、この話は知らなくても解ける問題というのが減るわけではないと思います。でも実装のアプローチを複数持っておくこと自体は悪くないと思いますし、問題の背景にこういう性質が隠れてるのを知るのもそれはそれで面白いんじゃないかと思います。この記事が少しでも皆さんの競プロ人生の役に立てば幸いです。

参考文献

Dynamic Programming Optimizations

Convex-Hull Trick

Totally Monotone Matrix Searching (SMAWK algorithm)

動的計画法高速化テクニック（オフライン・オンライン変換）

問題概要

英語だけど短いし問題文を読んで