プロジェクトオイラー602.

この記事の所要時間: 511

No.602 Product of Head Counts

問題.

Q. アリスは友達に手伝ってもらって, 正しくない(表裏の出る確率が等しくない)コインを使ってランダムな数を生成する. 彼女らは円形に座って, アリスから順番にコインを投げる. それぞれの人は自分が表を出した回数を覚えておく. アリスが表を出した時点で終了とし, アリスはそれぞれの(他の)人の表を出した回数を乗じることでランダムな数を得る.

例えば, アリスのほかにボブ, チャーリー, ドーンで行って, コインの表裏は THHH-TTTT-THHT-H となったとする(T:裏, H:表). この場合, ボブ, チャーリー, ドーンの表の回数はそれぞれ 2, 2, 1 回なので, ランダムな数は \(2\times 2\times1=4\) となる.

ここで, 友達が \(n\) 人 (アリス自身は含めない), コインの裏が出る確率を \(p\) とするときのランダムな数の期待値を \(e(n, p)\) とする. 一般に, \(n\) を固定すると, \(e(n, p)\) は \(p\) の多項式となる. 例えば,

\begin{align*}
e(3, p) = p^3+4p^2+p
\end{align*}

\(c(n, k)\) を多項式 \(e(n, p)\) における \(p^k\) の係数とする. 例えば

\begin{align*}
c(3, 1) &= 1\\
c(3, 2) &= 4\\
c(3, 3) &= 1\\
c(100, 40) &\equiv 987799437\pmod{10^9+7}
\end{align*}

では,

\begin{align*}
c(10000000, 4000000)\bmod 10^9+7
\end{align*}

を求めよ.

考え方とプログラム例(Python)

問題の設定が複雑なので, 数学的に解析していきます.

まず, アリスが \(m\) 回目にはじめて表を出す確率は, \(p^{m-1}(1-p)\) です.

他の \(n\) 人はそれぞれ \(m-1\) 回コインを投げて, 1 回ごとに表が出る確率は \((1-p)\) で, 1 人が表を出す回数の期待値を計算すると

\begin{align*}
&\,\sum_{j=0}^{m-1} j\times {}_{m-1}C_{j}\,(1-p)^jp^{m-1-j}\\ &= \sum_{j=1}^{m-1} j\cdot\dfrac{(m-1)!}{j!(m-1-j)!}(1-p)^jp^{m-1-j}\\
&= \sum_{j=1}^{m-1} (m-1)\cdot\dfrac{(m-2)!}{(j-1)!(m-1-j)!}(1-p)^jp^{m-1-j}\\
&= (m-1)(1-p)\sum_{j=1}^{m-1} {}_{m-2}C_{j-1}(1-p)^{j-1}p^{m-1-p}\\
&= (m-1)(1-p)
\end{align*}

この部分は上のように丁寧に計算しなくても結果を知っている人は多いでしょう.

では, ランダムな数の期待値 \(e(n, p)\) を計算してみると,

\begin{align*}
e(n, p) &= \sum_{m=1}^\infty p^{m-1}(1-p)\{(m-1)(1-p)\}^n\\
&= (1-p)^{n+1}\sum_{m=1}^\infty m^np^m.
\end{align*}

右辺は無限級数表示になっていますが, \(|p|<1\) のとき \(p\) の \(n\) 次多項式に収束することが知られています. その係数 \(c(n, k)\) は, 右辺の和を \(m=k\) までで止めたものの係数と一致します. すると, \(c(n, k)\) は次のような和の形で書けます.

\begin{align*}
c(n, k) = \sum_{l=0}^{k-1} (-1)^l{}_{n+1}C_{l} (k-l)^n
\end{align*}

後は \(n=10000000, k=4000000\) を代入して計算すれば答えが出ます. コンビネーション(二項係数)の部分は

\begin{align*}
\binom{n+1}{m} = \binom{n+1}{m-1}\cdot\dfrac{n-m+2}{m}
\end{align*}

を使うことで早く計算できます.

ここで, 計算途中の数値が巨大なものになるため, 常に \(M=10^9+7\) で割った余りを計算します. 足し算や掛け算は問題ないのですが, 剰余(mod) で計算していると割り算(上のコンビネーションの計算で必要)を行えません. そこで, 剰余での逆数を考える必要があります.

\(M\) を法としたときの \(a\) の逆数 \(a^{-1}\) は,

\begin{align*}
a\cdot a^{-1}\equiv 1\pmod{M}
\end{align*}

を満たす数のことです.

\(M=pa+q\quad (0\leqq q<a)\) のとき,

\begin{align*}
a^{-1} = M – (a * b^{-1}) \bmod M
\end{align*}

が成り立つので, これを使います.

\((k-l)^n\) の部分は \(n\) の2進数表示を利用して, 掛け算のたびに \(M\) で割った余りを計算しています. また, 逆数は必要な分だけ先に計算して inv[] に用意しておきます.

スポンサーリンク
sub2
sub2
  • このエントリーをはてなブックマークに追加