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{align}\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}
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 |
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 |