数学の力

京大生が数学の定理・公式の証明や入試問題の解説をするブログ.

プロジェクトオイラー479. ~Roots on the Rise~

No.479 Roots on the Rise

問題. a_k, b_k, c_kを方程式

\begin{align}
\frac{1}{x} &= \left(\frac{k}{x}\right)^2(k+x^2)-kx
\end{align}

の3つの解とする.

例えば k=5のとき \{a_5, b_5, c_5\}はおよそ

\begin{align}
\{5.727244, -0.363622+2.057397i, -0.363622-2.057397i\}
\end{align}

である.

\begin{align}
S(n) &= \sum_{p=1}^n\sum_{k=1}^n(a_k+b_k)^p(b_k+c_k)^p(c_k+a_k)^p
\end{align}

とおく.面白いことに S(n)は常に整数であり,例えば S(4)=51160となる.

 S(10^6) 1000000007で割ったあまりを求めよ.


一見すると S(n)が整数になることに驚きますが,実は簡単に示すことができます.



考え方とプログラム例.

まず与えられた方程式を変形して分かりやすい形にします.両辺に x^2をかけて,

\begin{align}
x = k^2(k+x^2)-kx^3
\end{align}

移項,整理すると

\begin{align}
kx^3-k^2x^2+x-k^3=0
\end{align}

後で使うので, f(x)=x^3-kx^2+\frac{x}{k}-k^2とおきます.(上の式の左辺を kで割り,モニック(最高次の係数が1)にしたもの)

3次方程式なので,解と係数の関係から

\begin{align}
a_k+b_k+c_k &= k
\end{align}

が得られるので,

\begin{align}
\left\{\begin{array}{ll}
a_k+b_k &= k-c_k\\
b_k+c_k &= k-a_k\\
c_k+a_k &= k-b_k
\end{array}\right.
\end{align}

 S(n)の式に代入すると

\begin{align}
S(n) &= \sum_{p=1}^n \sum_{k=1}^n (k-c_k)^p(k-a_k)^p(k-b_k)^p\\
&= \sum_{p=1}^n \sum_{k=1}^n \{(k-a_k)(k-b_k)(k-c_k)\}^p
\end{align}

さて, a_k, b_k, c_k は3次方程式 f(x)=0の解なので,

\begin{align}
f(x) =  (x-a_k)(x-b_k)(x-c_k)
\end{align}

と書けるので,

\begin{align}
(k-a_k)(k-b_k)(k-c_k) &= f(k)\\
&= k^3-k^3+1-k^2\\
&= 1-k^2
\end{align}

よって,

\begin{align}
S(n) &= \sum_{p=1}^n \sum_{k=1}^n (1-k^2)^p
\end{align}

この式から, S(n) nによらず整数になることが分かります.さて,和を取る順序を変えて先に pについての和を取れば,

\begin{align}
S(n) &= \sum_{k=1}^n (1-k^2)\frac{1-(1-k^2)^n}{1-(1-k^2)}\\
&= \sum_{k=1}^n (1-k^2)\frac{1-(1-k^2)^n}{k^2}
\end{align}

import math
N = 10**6
mod = 10**9+7
binary_mod = []
Mod = mod-2
while Mod:
    binary_mod.append(Mod % 2)
    Mod //= 2
def rev(x):
    ret = 1
    x %= mod
    for i in xrange(len(binary_mod)):
        if binary_mod[i] == 1:
            ret = ret * x % mod
        x = x * x % mod
    return ret
binary_N = []
N_p = N
while N_p:
    binary_N.append(N_p % 2)
    N_p //= 2
def pow_N(x):
    ret = 1
    x %= mod
    for i in xrange(len(binary_N)):
        if binary_N[i] == 1:
            ret = ret * x % mod
        x = x * x % mod
    return ret
Sum = 0
print rev(2)
for k in xrange(2, N+1):
    a = (1 - k*k)
    Sum = (Sum - ((pow_N(a) - 1) * a % mod) * rev(k*k) % mod) % mod 
print Sum

 

プロジェクトオイラー323. ~Bitwise-OR operations on random integers~

No.323 Bitwise-OR operations on random integers

久々のProject Euler は確率・期待値の問題です.高校数学でも使うことのある確率の考え方も出てきます.

問題.

 y_0, y_1, y_2, \ldots をランダムな符号なし32ビット整数とする(つまり, 0\leq y_i<2^{32}, 任意の値が同様に確からしく現れるとする).

 x_iに対して,次の漸化式が与えられている:

  •  x_0=0,
  •  x_i=x_{i-1}|y_{i-1}, (i>0). (  | はビットごとのORをとる操作)
最終的に,ある添字 Nが存在して全ての i\geq Nに対して x_i=2^{32}-1(すべての桁が1であるビットパターン)となることが分かる.
 Nの期待値を求めよ.答えは小数点以下10桁に四捨五入して答えよ.


コンピュータの中では整数は2進法の0と1の列として扱われます. 1\leq n\leq 2^{32}-1を満たす自然数 nは,(左端に0が並ぶことを許すと)32桁の0, 1列で表されます.

  | はビットごと(桁ごと)にOR演算を行う計算です.OR演算は,
\begin{align*}
0 \quad\mathrm{OR}\quad 0 &= 0\\
0 \quad\mathrm{OR}\quad 1 &= 1\\
1 \quad\mathrm{OR}\quad 0 &= 1\\
1 \quad\mathrm{OR}\quad 1 &= 1
\end{align*}

という演算で,1が1個以上あれば答えが1になる演算です.そのため,上の問題は各桁に対してランダムに0または1が与えられて,すべての桁において1が1回以上登場するまでの回数の期待値,と考えることができます.

また, N=kになる確率 \mathrm{P}(N=k)を直接求めるのは難しいため, N\leq kとなる確率 Q(k)=\mathrm{P}(N\leq k)を求めて,

\begin{align*}
\mathrm{P}(N=k) &=  \mathrm{P}(N\leq k) - \mathrm{P}(N \leq k-1)\\
&=  Q(k) - Q(k-1)
\end{align*}

として求めます.この考え方は高校数学などの場合の数・確率の分野でも役に立つので知っておいて損はないでしょう.



解説とプログラム例.(c++)

 x j桁目( j=1, 2, \ldots, 32)が初めて 1になるときの添字を表す確率変数を S_jとする.このとき,
\begin{align*}
\mathrm{P}(S_j=i) &= \left(\frac{1}{2}\right)^i\\
\mathrm{P}(S_j\leq i) &= \sum_{k=1}^i \left(\frac{1}{2}\right)^k\\
&= 1-\left(\frac{1}{2}\right)^i.
\end{align*}
このとき, x_i=2^{32}-1である確率は,
\begin{align*}
\mathrm{P}(x_i=2^{32}-1) &= \mathrm{P}(S_1\leq i, S_2\leq i, \ldots, S_{32}\leq i)\\
&= \mathrm{P}(S_1\leq i)\cdot\mathrm{P}(S_2\leq i)\cdots\mathrm{P}(S_{32}\leq i)\\
&= \left\{1-\left(\frac{1}{2}\right)^i\right\}^{32}.
\end{align*}
よって, x_iではじめて 2^{32}-1になる確率( N=iである確率)は,
\begin{align*}
\mathrm{P}(x_{i-1}\neq 2^{32}-1, x_i=2^{32}-1) &= \left\{1-\left(\frac{1}{2}\right)^i\right\}^{32} - \left\{1-\left(\frac{1}{2}\right)^{i-1}\right\}^{32}\\
&= \sum_{k=0}^{32} \binom{32}{k}(-1)^k\left(\frac{1}{2}\right)^{ik}(1-2^k)\\
&= \sum_{k=1}^{32} \binom{32}{k}(-1)^k\left(\frac{1}{2}\right)^{ik}(1-2^k).
\end{align*}
したがって, Nの期待値は
\begin{align*}
\mathrm{E}[N] &= \sum_{i=1}^{\infty} i \sum_{k=0}^{32} \binom{32}{k}(-1)^k\left(\frac{1}{2}\right)^{ik}(1-2^k)\\
&= \sum_{k=1}^{32} \binom{32}{k}(-1)^k(1-2^k)\sum_{i=1}^\infty i\left(\frac{1}{2}\right)^{ik}
\end{align*}
ここで,

\begin{align*}
T = \sum_{i=1}^\infty i\left(\frac{1}{m}\right)^i
\end{align*}

とおくと,

\begin{align*}
\frac{m-1}{m}T = T - \frac{1}{m}T &=  \left(\frac{1}{m}+\frac{2}{m^2}+\frac{3}{m^3}+\cdots\right)-\left(\frac{1}{m^2}+\frac{2}{m^3}+\frac{3}{m^4}+\cdots\right)\\
&=  \frac{1}{m}+\frac{1}{m^2}+\frac{1}{m^3}+\cdots\\
&=  \frac{1}{m}\cdot\frac{1}{1-\frac{1}{m}}\\
&=  \frac{1}{m-1}
\end{align*}

より

\begin{align*}
T = \sum_{i=1}^\infty i\left(\frac{1}{m}\right)^i &=  \frac{m}{(m-1)^2}
\end{align*}

となる.

これを利用すると,

\begin{align*}
\mathrm{E} [N] &=  \sum_{k=1}^{32} \binom{32}{k}(-1)^k(1-2^k)\frac{2^k}{(2^k-1)^2}\\
&=  \sum_{k=1}^{32} \binom{32}{k}(-1)^{k+1}\frac{2^k}{2^k-1}.
\end{align*}

 k=1, 2, \ldots, 32について各項は容易に計算できます.

以下がc++によるプログラム例です.ちなみに,double型で計算すると誤差で小数点以下10桁までの答えが合わなかったので,long double型で計算しました.

#include<bits/stdc++.h>
using namespace std;
void print(int n, double d){cout << fixed << setprecision(n) << d << endl;}
int code(int i){if(i % 2 == 0) return 1; else return -1;}

int main(){
    const int N = 32;
    long double ans = 0.0;
    long long comb = 1;
    long long P = 1;
    for(int k = 1; k <= N; k++){
        comb = comb * (N - k + 1) / k;
        P *= 1;
        ans = ans + (code(k+1) * comb * P) / (ld)(P-1);
    }
    print(10, ans);
    return 0;
}