「M-SOLUTIONS プロコンオープン」C問題

M-SOLUTIONS Programming Contest - AtCoder

2019年6月1日開催のAtCoder M-SOLUTIONS プロコンオープンに参加した．結果はA, B 2完．このnoteでは時間内に解けなかったC問題の私なりの理解，解説を紹介する．

C問題 Best-of-(2n-1)

問題はこちら．配点は500点．今までのC問題と比べ配点も高く，難易度も高いように感じた．この問題を正解するには，ゲームが行われる回数の期待値 P/Q を求める方法と，それをmod 10^9 +7 で計算する方法の理解が必要．後者についてはこちらの記事がとてもわかりやすかった．

ゲームが行われる回数の期待値の求め方

まず，引き分けが発生しない(Cが0)場合を考える． ゲーム終了時点での高橋君，青木君の勝数のパターンは(N,0), (N, 1), (N, 2), ... , (N,N-1) または，(0,N),(1,N),...,(N-1,N) が考えられる．これらのパターンは全て排反であるため，それぞれに対して(生じる確率)x(ゲームが行われる回数)を計算し和を取ると期待値になる．例えば (N, i)の時は，A/(A+B) = a, B/(A+B) = b とおいて

(N+i回目に高橋君が勝つ確率)*(1~N+i-1回目の間で高橋君がN-1回勝ち，青木君がi回勝つ確率) * (N+i) 

= a * a^(N-1) * b^i * (N+i-1個からi個選ぶ組み合わせ) * (N+i)

となる．これを地道に足し上げたものが，引き分けが発生しない場合の期待値になる．次に引き分けによるゲーム回数の嵩増しを考える．引き分けが発生しない場合，ゲームの推移は，高: 高橋君の勝ち，青:青木君の勝ちとして

高，青，高，青，高，高，青，高

のようになるが，引き分けがある場合は，上の推移の間に引き分けが割り込んで

分高，青，分高，分青，分分高，分高，分青，分分高

のようになる．それぞれの高か青の前に引き分けが何回か(0回も含む)発生すると考えられる．例えば，一番右では

高 → 分分高

となっている．一度，高or青となってから次に高or青になるまでのゲームの回数(分の回数+1)の期待値は，「1回のゲームで高or青となる確率の逆数」になる．直感的には，例えば，高or青となる確率が1/4，引き分けの確率が3/4の時，平均的に4回ゲームをすれば高or青となりそうだとわかる．(数学的に詳しいことは「幾何分布」を調べるとわかる)．

さて，一度高or青となってから次に高or青になるまでのゲームの回数の期待値は，「1回のゲームで高or青となる確率の逆数」になるので，引き分けがある場合のゲームが行われる回数の期待値は，引き分け無しのときの期待値の「1回のゲームで高or青となる確率の逆数」倍になる．つまり，引き分けも考慮したゲーム回数の期待値は

a * a^(N-1) * b^i * (N+i-1個からi個選ぶ組み合わせ) * (N+i) * 「1回のゲームで高or青となる確率の逆数」

= a * a^(N-1) * b^i * (N+i-1個からi個選ぶ組み合わせ) * (N+i) * 100/(100-C)

= a * a^(N-1) * b^i * (N+i-1)! * (N-1)! * i! * (N+i) * 100/(100-C)

となる．これをC++で実装すると次のようになる．(#includeは省略)

#define rep(i,n) for(int i=0;i<n;i++)
using namespace std;

long long factorial(long long i){ // 階乗
 if(i > 0)
   return i * factorial(i-1);
 else
   return 1;
}

int main() {
 long long N, A, B, C;
 cin >> N >> A >> B >> C;

 // N!をあらかじめ計算しておく
 vector<long long> factmod;
 factmod.push_back(1); // 0! = 1 と定義
 for(int i=1; i<2*N+10; i++) { // (2*N)! まで計算しておけば十分
   long long tmp = (factmod[i-1]*i);
   factmod.push_back(tmp);
 }

 long long res = 0;

 rep(i, N) {   // AがN回勝ち，Bがi回勝ちで終了．
   res += (N+i) * factorial(N-1+i) * pow(A, N) * pow(B, i)/ pow(A+B, N+i) / factorial(N-1) /factorial(i) * 100.0 /  (100.0-C);
 }
 rep(i, N) {   // Aがi回勝ち，BがN回勝ちで終了．
   res += (N+i) * factorial(N-1+i) * pow(B, N) * pow(A, i)/ pow(A+B, N+i) / factorial(N-1) /factorial(i) * 100.0 /  (100.0-C);
 }

 cout << res << endl;
 return 0;
}

注) rep(i,n)はfor文を短く書くために定義しているだけ．ループ 内で繰り返し使用する数(N!など)は，ループの外側であらかじめ計算しておく．特に階乗は事前にまとめて実施しておく．

試しに，入力例1: 1 25 25 50 を入力すると2 が出力される．

あとは，mod 10^9 +7 で計算すればよい．

mod 10^9 +7 で計算する方法

冒頭でも触れたがこれに関しては，Kensuke Otsukiさんの大変詳しい記事をtwitter経由で見つけたのでそちらを参照してほしい．mod p における割り算がややこしいが，上で求めた期待値の計算をひたすらmod p での計算していけばよい．どうやら出力すべき R は P ÷ Q をmod 10^9 +7 で行ったものらしい．C++で実装すると次のようになる．(#includeは省略)

#define rep(i,n) for(int i=0;i<n;i++)
using namespace std;

long long factorial(long long i){ // 階乗
 if(i > 0)
   return i * factorial(i-1);
 else
   return 1;
}

// a^n mod
long long modpow(long long a, long long n, long long mod) {
 long long res = 1;
 while (n>0) {
   if(n & 1) res = res * a % mod;
   a = a * a % mod;
   n >>= 1;
 }
 return res;
}

// a^(-1) mod
long long modinv(long long a, long long mod) {
 return modpow(a, mod-2, mod);
}

int main() {
 long long N, A, B, C;
 const long long mod = 1000000007;
 cin >> N >> A >> B >> C;

 // N!(mod 1000000007)をあらかじめ計算しておく
 vector<long long> factmod;
 factmod.push_back(1); // 0! = 1 と定義
 for(int i=1; i<2*N+10; i++) { // (2*N)! まで計算しておけば十分
   long long tmp = (factmod[i-1]*i) %mod;
   factmod.push_back(tmp);
 }

 // よく使う値を事前に計算
 long long invAplusB = modinv(A+B, mod);
 long long invFactorialNm1 = modinv(factmod[N-1], mod);
 long long invHandmC = modinv(100-C, mod);
 long long mothA = modpow(A, N, mod)%mod * invFactorialNm1%mod * 100 %mod * invHandmC %mod;
 long long mothB = modpow(B, N, mod)%mod * invFactorialNm1%mod * 100 %mod * invHandmC %mod;
 
 long long res = 0;

 rep(i, N) {   // AがN回勝ち，Bがi回勝ちで終了．
   // 通常の期待値
   //res += (N+i) * factorial(N-1+i) * pow(A, N) * pow(B, i)/ pow(A+B, N+i) / factorial(N-1) /factorial(i) * 100.0 /  (100.0-C);
   res += (N+i) * factmod[N-1+i]%mod  * modpow(B, i, mod)%mod * modpow(invAplusB, N+i, mod)%mod  * modinv(factmod[i], mod)%mod * mothA%mod;
   res = res % mod;
 }
 rep(i, N) {   // Aがi回勝ち，BがN回勝ちで終了．
   // 通常の期待値
   //res += (N+i) * factorial(N-1+i) * pow(B, N) * pow(A, i)/ pow(A+B, N+i) / factorial(N-1) /factorial(i) * 100.0 /  (100.0-C);
   res += (N+i) * factmod[N-1+i]%mod * modpow(A, i, mod)%mod * modpow(invAplusB, N+i, mod)%mod * modinv(factmod[i], mod)%mod * mothB % mod;
   res = res % mod;
 }

 cout << res << endl;
 return 0;
}

注) ループ 内で繰り返し使用する数は，ループの外側であらかじめ計算しておく．掛けたり足したりするたびに%mod を行う．特に階乗は事前にまとめて実施しておく．