CodeIQ MAGAZINECodeIQ MAGAZINE

数学の問題をプログラミングで解こう!「ルート・パワー」問題解説

2015.12.25 Category:CodeIQ問題解説・リーダーボード Tag: ,

  • 5
  • このエントリーをはてなブックマークに追加
riverplus-300x288

「数学の問題をプログラミングで解こう!」シリーズ。

平方根の和のべき乗に関する問題です。難易度としては、やや難しめの★3つとして出題させて頂きました。複雑な手計算やプログラミングが必要となり、苦労された方が多かったかもしれません。
ぜひ本記事で解法をチェックして下さい!
by @riverplus Kawazoe

問題のおさらい

自然数 n に対して、(1 + √2 + √3 + √5)n を展開して整理したものを考えます。

例えば、n = 2 のとき、(1 + √2 + √3 + √5)2 を展開して整理すると、11 + 2√2 + 2√3 + 2√5 + 2√6 + 2√10 + 2√15 となります。

このときの有理数の項の値は 11 です。

(1 + √2 + √3 + √5)n を展開して整理した有理数の項を F(n) と定義します。

例えば、F(2) = 11、F(3) = 31、F(12) = 577862325 です。
また、F(12) を 107 で割った余りは 7862325 となることが確かめられます。

標準入力から、自然数 n(1 ≦ n ≦ 1010)が与えられます。
標準出力に F(n) を 107 で割った余りの値を出力するプログラムを書いてください。

方針

平方根の和の n 乗を展開したものに関する問題です。3, 4 乗ぐらいなら手計算でも何とかできますが、それ以上となるとそうはいきません。さらに本問は n の最大が 1010 と非常に大きな値を扱います。これだけ大きな乗数でもうまく計算できる方法を考えていきましょう。

本問を解くための一つ目のポイントは、一つ手前の状態を考えることです。1 + √2 + √3 + √5 を n 乗したものと、n – 1 乗したものとの関係を見ていきましょう。そのために、(1 + √2 + √3 + √5)n を展開したものの 1, √2, √3, √5, √6, √10, √15, √30 の各係数を、an, bn, cn, dn, en, fn, gn, hn とおきます。

そして、(an-1,・・・, hn-1) と (an,・・・, hn) の間の関係を調べます。ようするに漸化式を導出するということです。(an,・・・, hn) は、 n – 1 のときの結果に 1 + √2 + √3 + √5 を掛け算して展開した時の各項の係数です。たいへん根気のいる作業ですが、次の計算で得ることができます。

このように、8 変数間の漸化式が得られました。単純には、初期値 (a0, b0, c0, d0, e0, f0, g0, h0) = (1, 0, 0, 0, 0, 0, 0, 0) から始めて、漸化式を繰り返し適用していけば、求めるべき an の値が得られます。コードは次のようになります(C++)。

#include <iostream>
using namespace std;
long long D = 10000000;

int main(int argc, char *argv[]) {
    long long n;
    cin >> n;

    long long a = 1, b = 0, c = 0, d = 0, e = 0, f = 0, g = 0, h = 0;
    long long a_new, b_new, c_new, d_new, e_new, f_new, g_new, h_new;
    for (long long i = 1; i <= n; i++) {
        a_new = a + 2 * b + 3 * c + 5 * d;
        b_new = a + b + 3 * e + 5 * f;
        c_new = a + c + 2 * e + 5 * g;
        d_new = a + d + 2 * f + 3 * g;
        e_new = b + c + e + 5 * h;
        f_new = b + d + f + 3 * h;
        g_new = c + d + g + 2 * h;
        h_new = e + f + g + h;
        a = a_new % D;
        b = b_new % D;
        c = c_new % D;
        d = d_new % D;
        e = e_new % D;
        f = f_new % D;
        g = g_new % D;
        h = h_new % D;
    }
    cout << a << endl;
}

ただ、この方法では本問の 1010 という上限では、長い実行時間を要してしまいます。クリアするためにはさらなる高速化が必要です。

行列のべき乗に帰着させる

そこで、まずは上の 8 つの漸化式を行列の形で表してみましょう。次のようになります。

二つの列ベクトルをそれぞれ矢印付きの xnxn-1、右辺の 8×8 行列を A と表します。

この形の漸化式を次々と適用していくと、最終的に次のように変形できます。

結局のところ、この問題は、行列 An 乗(の各要素を 107 で割ったあまり)をいかに高速に計算するかという問題に帰着されるわけです。

行列のべき乗に関しては、効率的な計算方法が知られています。以前の出題の記事(数学の問題をプログラミングで解こう!「エース・ナンバー」問題解説)で、整数のべき剰余の高速な算出方法を紹介しました。このテクニックを使うことが、本問を解くもう一つのポイントです。本質的には、べき乗される整数の部分が行列に置き換わっただけです。n 回の行列の掛け算を行う必要はなく、log n 回程度の回数で算出できます。An が計算できたなら、そこに 初期値 x0 = (1, 0, 0, 0, 0, 0, 0, 0) を掛け算すれば、an, bn, cn, dn, en, fn, gn, hn の各値が得られます。(つまりは An の (1, 1) 成分が求めるべき an の値です。)

参考:冪剰余(Wikipedia)

コード例は以下の通りです(C++14)。

#include <iostream>
#include <vector>
using namespace std;
long long D = 10000000;

class Matrix {
public:
    Matrix(void) {};
    Matrix(int size) {
        v.assign(size, vector(size));
    }
    Matrix(vector> u) {
        v = u;
    }
    ~Matrix(void) {};
    vector& operator [] (const int index) {
        return v[index];
    }
    int GetSize(void) {
        return v.size();
    }
    Matrix Pow(long long n, long long D = 0) {
        if (n == 1) return *this;
        Matrix tmp(v.size());
        if (n == 0) {
            for (int i = 0; i < v.size(); i++) tmp[i][i] = 1;
        }
        else if (n % 2 == 0) {
            tmp = Pow(n / 2, D);
            tmp = tmp.Mul(tmp, D);
        }
        else {
            tmp = Pow(n - 1, D);
            tmp = Mul(tmp, D);
        }
        return tmp;
    }
    Matrix Mul(Matrix &a, long long D) {
        int size = a.GetSize();
        Matrix c(size);
        for (int i = 0; i < size; i++) {
            for (int j = 0; j < size; j++) {
                c[i][j] = 0;
                for (int k = 0; k < size; k++) {
                    c[i][j] = (c[i][j] + v[i][k] * a[k][j])  % D;
                }
            }
        }
        return c;
    }
private:
    vector> v;
};

int main(int argc, char *argv[]) {
    long long n;
    cin >> n;

    Matrix m({
        { 1, 2, 3, 5, 0, 0, 0, 0 },
        { 1, 1, 0, 0, 3, 5, 0, 0 },
        { 1, 0, 1, 0, 2, 0, 5, 0 },
        { 1, 0, 0, 1, 0, 2, 3, 0 },
        { 0, 1, 1, 0, 1, 0, 0, 5 },
        { 0, 1, 0, 1, 0, 1, 0, 3 },
        { 0, 0, 1, 1, 0, 0, 1, 2 },
        { 0, 0, 0, 0, 1, 1, 1, 1 }
    });
    Matrix p = m.Pow(n, D);
    cout << p[0][0] << endl;
}

みんなのコード

本問の挑戦者で、コードを公開して頂いた方のコードを Togetter でまとめました。(随時追加していきます。)こちらもぜひチェックして下さい!

CodeIQ「ルート・パワー」問題 みんなのコード

CodeIQ運営事務局より

Kawazoeさん、ありがとうございました!
現在、Kawazoeさんの最新問題が出題中です。
ぜひ挑戦してみてくださいね!

  • 5
  • このエントリーをはてなブックマークに追加

■この記事を書いた人

avatar

@riverplus Kawazoe

数学・パズル好きのソフトウェアエンジニアです。某社で研究開発をしています。Project Euler の問題開発チームメンバー(最近はあまり活動せず)。twitter: @riverplus Web Site: riverplus.net

■関連記事

【謎解きプログラム】条件に当てはまる文字列は?【正規表現】解答と解説... 【謎解きプログラム】条件に当てはまる文字列は?【正規表現】 本問題は、表題のテーマで、プログラムにちなんだ謎を解くというものでした。 それでは以下、各問題とその解答を見ていきましょう。 問題のオープニング ある日、出社すると、あなたのPCのログイン画面に、謎の挑戦状が表示されていた。 「2...
数学の問題をプログラミングで解こう!「クロッシング・ワード」問題解説... 問題のおさらい クロスワードの盤面では、格子状のマス目に白マスまたは黒マスを配置します。 以下は、縦 3 個×横 4 個のマス目に白マス・黒マスを配置する例です。 白マス・黒マスの配置には次のルールがあります。 黒マスによって白マスの領域が分断されてはならない。 黒マスが縦・横に連続し...
【謎解きプログラム】乱数で発生する数値は?【組み合わせ】解答と解説... 【謎解きプログラム】乱数で発生する数値は?【組み合わせ】 本問題は、表題のテーマで、プログラムにちなんだ謎を解くというものでした。 それでは以下、各問題とその解答を見ていきましょう。 問題のオープニング ある日、出社すると、あなたのPCのログイン画面に、謎の挑戦状が表示されていた。 「24...
「放物線とマス目の関係」問題の解答と解説... table.nabe{ margin-left:30px; } .nabefloat{ float:right; } table.nabe td, table.nabe th{ padding:3px; } table.nabe th{ ...
【謎解きプログラム】データをバイナリで見てみよう【バイナリ】解答と解説... 【謎解きプログラム】データをバイナリで見てみよう【バイナリ】 本問題は、表題のテーマで、プログラムにちなんだ謎を解くというものでした。 それでは以下、各問題とその解答を見ていきましょう。 問題のオープニング ある日、出社すると、あなたのPCのログイン画面に、謎の挑戦状が表示されていた。 「...
数学の問題をプログラミングで解こう!「プライム・ペア」問題解説... 問題のおさらい 自然数 k に対し、1 から k までの自然数のうち k と互いに素なものの個数を F(k) と定義します。 (F(k) はオイラーのφ関数とも呼ばれています。参考:オイラーのφ関数(Wikipedia)) 例えば F(12)=4 です。1 から 12 のうち 12 と互いに素...

今週のPickUPレポート

新着記事

週間ランキング

CodeIQとは

CodeIQ(コードアイキュー)とは、自分の実力を知りたいITエンジニア向けの、実務スキル評価サービスです。

CodeIQご利用にあたって
関連サイト
codeiq

リクルートグループサイトへ