Project Euler - Problem 479
Roots on the Rise
式 1/x = (k/x)2(k+x2) - k x の3つの解(実数か複素数)を ak, bk, ck で表すとしよう.
例えば. k = 5 の場合, {a5, b5, c5} はおよそ {5.727244, -0.363622+2.057397i, -0.363622-2.057397i} となる.
1 ≤ p, k ≤ n となるようなすべての整数 p, k に対し S(n) = Σ (ak+bk)p(bk+ck)p(ck+ak)p としよう.
面白いことに, S(n) は常に整数となる. 例えば, S(4) = 51160.
S(106) modulo 1 000 000 007 を求めよ.
(日本語訳より)
例えば. k = 5 の場合, {a5, b5, c5} はおよそ {5.727244, -0.363622+2.057397i, -0.363622-2.057397i} となる.
1 ≤ p, k ≤ n となるようなすべての整数 p, k に対し S(n) = Σ (ak+bk)p(bk+ck)p(ck+ak)p としよう.
面白いことに, S(n) は常に整数となる. 例えば, S(4) = 51160.
S(106) modulo 1 000 000 007 を求めよ.
(日本語訳より)
ハッピーニューイヤーということで 250 問目,Level 10 になりました めでたい.
プログラミング自体を始める切っ掛けにもなった Problem 1 を C 言語入門みたいなサイト見ながら悪戦苦闘して解いたのが 2013 年の 9 月ということで,間にブランクを挟みながらもプログラミング始めてもう 4 年以上経っているらしい.4 年でこれは先が思いやられるが…
問題の式を整理すると
$ k x^{3} - k^{2} x^{2} + x - k^{3} = 0 $
高校数学でおなじみ解と係数の関係より,
$a_k + b_k + c_k = k$
$a_k b_k + b_k c_k + c_k a_k = \frac{1}{k}$
$a_k b_k c_k = k^{2}$
$a_k b_k + b_k c_k + c_k a_k = \frac{1}{k}$
$a_k b_k c_k = k^{2}$
なので,
$ S(n) $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} (a_k + b_k)^{p} (b_k + c_k)^{p} (c_k + a_k)^{p}$
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} \{ (a_k + b_k)(b_k + c_k)(c_k + a_k)\}^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} \{ (k - a_k)(k - b_k)(k - c_k) \}^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} \{ k^{3} - (a_k + b_k + c_k) k^{2} + (a_k b_k + b_k c_k + c_k a_k) k - a_k b_k c_k \}^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} ( 1 - k^{2} )^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \frac{ (1 - k^{2}) - ( 1 - k^{2} )^{n + 1} }{ 1 - ( 1 - k^{2} ) } $
$ \displaystyle = \sum_{k = 1}^{n} k^{-2} \{ 1 - k^{2} - ( 1 - k^{2} )^{n + 1} \} $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} (a_k + b_k)^{p} (b_k + c_k)^{p} (c_k + a_k)^{p}$
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} \{ (a_k + b_k)(b_k + c_k)(c_k + a_k)\}^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} \{ (k - a_k)(k - b_k)(k - c_k) \}^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} \{ k^{3} - (a_k + b_k + c_k) k^{2} + (a_k b_k + b_k c_k + c_k a_k) k - a_k b_k c_k \}^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \sum_{p = 1}^{n} ( 1 - k^{2} )^{p} $
$ \displaystyle = \sum_{k = 1}^{n} \frac{ (1 - k^{2}) - ( 1 - k^{2} )^{n + 1} }{ 1 - ( 1 - k^{2} ) } $
$ \displaystyle = \sum_{k = 1}^{n} k^{-2} \{ 1 - k^{2} - ( 1 - k^{2} )^{n + 1} \} $
ここまでくればあとは Mod をとればいい.
$ k^{-1} \mod 1 000 000 007$
はモジュラ逆数でなんとかします.これがないと (1 - k2)n + 1の部分を二項展開して計算量が 1012 とかになって死ぬ.
typedef long long lint; // 拡張ユークリッド互除法 (入力 : 整数 a, n, 出力 : GCD(a, n) と xa + yb = GCD(a, n) を満たす整数 x, y) void ExEuclid(lint a, lint b, lint& gcd, lint& x, lint& y) { lint r1, r2, r3, x1, x2, x3, y1, y2, y3, q; r1 = a; x1 = 1; y1 = 0; r2 = b; x2 = 0; y2 = 1; while (r2 != 0) { r3 = r1 % r2; q = r1 / r2; x3 = x1 - q * x2; y3 = y1 - q * y2; r1 = r2; r2 = r3; x1 = x2; x2 = x3; y1 = y2; y2 = y3; } gcd = r1; x = x1; y = y1; if (gcd < 0) { gcd *= -1; x *= -1; y *= -1; }; } // a と m は互いに素が前提 lint invmod(lint a, lint m) { lint gcd, x, q; ExEuclid(a, -m, gcd, x, q); if (gcd != 1) return -1; return x; } // n^p mod m lint modpow(lint n, lint p, lint m) { lint result = 1; while (p > 0) { if (p & 1) result = (result * n) % m; p >>= 1; n = (n * n) % m; } return result; } int main() { clock_t start, end; start = clock(); const lint MOD = 1000000007; const lint MAX = 1000000; lint S = 0; for (lint k = 1; k <= MAX; k++) { lint x = modpow(1 - k, MAX + 1, MOD); lint y = modpow(1 + k, MAX + 1, MOD); lint z = (1 - k * k - x * y) % MOD; lint inv = invmod(k, MOD); z = z * inv % MOD; z = z * inv % MOD; S = (S + z) % MOD; } cout << S << endl; end = clock(); printf("%f sec.\n", double(end - start) / CLOCKS_PER_SEC); return 0; }
1.5 sec くらい.