絶滅

どうでもいい

Project Euler - Problem 479

Roots on the Rise

Problem 479 (日本語訳)

式 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 を求めよ.

(日本語訳より)




f:id:my316g:20180101210119p:plain

ハッピーニューイヤーということで 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}$

なので,

$ 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} \} $

ここまでくればあとは 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 くらい.