数学の力

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

プロジェクトオイラー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