SRM 287 DIV1 Hard CoinGame

SRM 行列 方程式 文字列

問題 Editorial

問題

2人でゲームをする。
シークエンスの集合sequencesが与えられる。
最初のプレイヤーはsequencesから1つ選ぶ。
相手のプレイヤーは最初のプレイヤーが選んだものとは違うシークエンスをsequencesから1つ選ぶ。
その後ランダムなコインを振り続ける。
コインが表になったか裏になったかの履歴が、どちらかのプレイヤーの選んだシークエンスと一致したならばそのプレイヤーの勝ちである。
両方のプレイヤーが最良にプレイするとき、最初のプレイヤーの勝つ確率を求めよ

  • 2 <= |sequences| <= 50
  • 1 <= |sequences[i]| <= 10
  • sequencesの要素は全て異なる
  • sequencesの要素の長さは全て同じ

解答

まず、どのシークエンスを選ぶかは全通り試す。
ここで、1つのシークエンスがもう1つより早く出る確率が求められればよい。
さて、これはDPで求められる。
パターンへのマッチング状態を状態として、HとTのそれぞれで0.5づつ足し込めばよい。「それぞれのプレイヤーが上がった」という状態も別個に用意する。
ここでマッチング状態はSRM 283 DIV1 Hard SuspiciousStringsの問題と同じようにすれば最大21状態程度でできる。
そしてこのDPを収束するまで回したいところだ。が、これは収束するのが遅すぎて間に合わない。
でもこれは行列の累乗で表現できることを思い出す。SRM 283 DIV1 Hard SuspiciousStringsでも使っている。
2^16程度の累乗で十分収束し、行列の乗算を速めに実装すれば十分(最大0.8s程度)間に合う。
ところで、この行列の乗算で表現できる漸化式の極限は連立方程式の解として表現できる(たぶん)
この連立方程式は以下のように構成する。

  1. xから始まった時の勝つ確率 = 0.5*xの後にHが出た後の状態から始まったときの勝つ確率 + 0.5*xの後にTが出た後の状態から始まったときの勝つ確率
  2. 相手が勝った状態から始まったときの勝つ確率 = 0
  3. 自分が勝った状態から始まったときの勝つ確率 = 1

そしてこれを解き、その空文字列のxが答え。
なんとなくそうであることはわかるけれども、どうしてこのようにできるのかはわかってない。
ガウスの消去法で解いた場合、行列累乗よりまあまあ速い(最大0.4s程度)

コメント

ガウスの消去法は行列累乗と比べて、計算量では(log k)分あるといっても、行列の乗算は定数が非常に軽い操作であるので、実際には半分以下程度の速さにしかならないのか。
ただしガウスの消去法の実装はもっと高速化できてもおかしくない。
しかし、ガウスの消去法が想定解法のものでも行列累乗でいけることは大いにありそうだ(ただし行列と連立方程式がどのように変換できるかはわかってない)。
マッチング状態のこの持ち方と行列累乗は短期間に2問も出ていて驚き。行列を作るコードは、前のコードより綺麗に書けたと思う(前のコードが冗長)
この連立方程式を解くというのは自分にとって初めてだった。数学全然苦手なので、常識なことがわかってないと思う。とりあえずこの形は使えそうだと覚えておく

コード

行列累乗バージョン
#include <string>
#include <vector>
#include <cstring>
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define mset(m,v) memset(m,v,sizeof(m))
using namespace std;
typedef double Num;
struct Matrix {
	Num a[2][22][22];	//vector<vector<T> >より配列のほうが速そう？
	Num (*v)[22], (*w)[22];
	Matrix(): v(a[0]), w(a[1]) { mset(a, 0); }
	inline int height() const { return 22; }
	inline int width() const { return 22; }
	inline Num& at(int i, int j) { return v[i][j]; }
	inline const Num& at(int i, int j) const { return v[i][j]; }
	static Matrix identity(int n) {
		Matrix A;
		rep(i, n) A.at(i, i) = 1;
		return A;
	}
	inline static Matrix identity(const Matrix& A) { return identity(A.height()); }
	Matrix& operator*=(const Matrix& B) {
		int n = height(), m = B.width(), p = B.height();
		rep(i, n) rep(j, m) {
			Num x = 0;
			rep(k, p)
				x += at(i, k) * B.at(k, j);
			w[i][j] = x;
		}
		swap(v, w);
		return *this;
	}
};

struct CoinGame {
	double bestProbability(vector <string> sequences) {
		double r = 0;
		int n = sequences.size(), m = sequences[0].size();
		rep(i, n) {
			double p = 0;
			rep(j, n) if(i != j) {
				static int nextH[22], nextT[22];
				string ss[2] = {sequences[i], sequences[j]};
				mset(nextH, 0); mset(nextT, 0);
				rep(k, 2) rep(l, m) {
					string prefix = ss[k].substr(0, l);
					rer(t, 0, l) rep(f, 2) if(prefix.substr(t, l - t) == ss[f].substr(0, l - t)) {
						int *next = ss[f][l - t] == 'H' ? nextH : nextT;
						if(l - t + 1 == m) next[k*10+l] = 20+f;
						else if(next[k*10+l] == 0) next[k*10+l] = f*10+(l-t+1);
					}
				}
				nextH[20] = nextT[20] = 20;
				nextH[21] = nextT[21] = 21;
				Matrix A;
				rep(k, 22) A.at(k, nextH[k]) += .5, A.at(k, nextT[k]) += .5;
				rep(k, 16) A *= A;
				p = max(p, (double)(A.at(0, 21) / (A.at(0, 20) + A.at(0, 21))));
			}
			r = max(r, 1 - p);
		}
		return r;
	}
};
ガウスの消去法バージョン
#include <string>
#include <vector>
#include <cstring>
#include <cmath>
#define rep(i,n) for(int (i)=0;(i)<(int)(n);++(i))
#define rer(i,l,u) for(int (i)=(int)(l);(i)<=(int)(u);++(i))
#define reu(i,l,u) for(int (i)=(int)(l);(i)<(int)(u);++(i))
#define mset(m,v) memset(m,v,sizeof(m))
using namespace std;
typedef double Num;
const int VecSize = 22;
struct Vec {
	Num a[VecSize];
	Num &operator[](int i) { return a[i]; }
	Vec() {mset(a, 0);}
};
struct Mat {
	Num a[VecSize][VecSize];
	Num *operator[](int i) { return a[i]; }
	Mat() {mset(a, 0);}
};

Vec gaussian_elimination(Mat A, Vec b) {
	const int n = VecSize, m = VecSize;
	int pi = 0, pj = 0;
	while(pi < n && pj < m) {
		reu(i, pi+1, n)
			if(fabs(A[i][pj]) > fabs(A[pi][pj])) {
				rep(k, m) swap(A[i][k], A[pi][k]);
				swap(b[i], b[pi]);
			}
		if(fabs(A[pi][pj]) > 1e-7) {
			Num d = A[pi][pj];
			rep(j, m) A[pi][j] /= d;
			b[pi] /= d;
			reu(i, pi+1, n) {
				Num k = A[i][pj];
				rep(j, m) A[i][j] -= k * A[pi][j];
				b[i] -= k * b[pi];
			}
			pi ++;
		}
		pj ++;
	}
	reu(i, pi, n) if(fabs(b[i]) > 0) return Vec();	//Inconsistent
	if(pi < m || pj < m) return Vec();	//Ambiguous
	for(int j = m-1; j >= 0; j --)
		rep(i, j)
			b[i] -= b[j] * A[i][j];
	return b;
}

struct CoinGame {
	double bestProbability(vector <string> sequences) {
		double r = 0;
		int n = sequences.size(), m = sequences[0].size();
		rep(i, n) {
			double p = 0;
			rep(j, n) if(i != j) {
				static int nextH[22], nextT[22];
				string ss[2] = {sequences[i], sequences[j]};
				mset(nextH, 0); mset(nextT, 0);
				rep(k, 2) rep(l, m) {
					string prefix = ss[k].substr(0, l);
					rer(t, 0, l) rep(f, 2) if(prefix.substr(t, l - t) == ss[f].substr(0, l - t)) {
						int *next = ss[f][l - t] == 'H' ? nextH : nextT;
						if(l - t + 1 == m) next[k*10+l] = 20+f;
						else if(next[k*10+l] == 0) next[k*10+l] = f*10+(l-t+1);
					}
				}
				Mat A;
				rep(k, 22) {
					A[k][k] = 1;
					if(k < m || (10<=k&&k<10+m)) {
						A[k][nextH[k]] -= .5;
						A[k][nextT[k]] -= .5;
					}
				}
				Vec b;
				b[21] = 1;
				Vec x = gaussian_elimination(A, b);
				p = max(p, x[0]);
			}
			r = max(r, 1 - p);
		}
		return r;
	}
};