[C#] 繰り返し二乗法で n^p % mod を拘束に求める

繰り返し二乗法とは

Np を高速に計算する方法です。NもPも10⁹くらいあるときに、NをP回かけ合わせると計算量 O(N) かかりますが、繰り返し二乗法を使えば O(logN) で実現できます。

Pが偶数のとき

例えばPが偶数(P=2k)の時、

N2k = Nk * Nk

となるので、N^kが計算できればそれを2乗するだけでよくなります。

具体的な例で、P=8の時には

N8 = N4 * N4
N4 = N2 * N2
N2 = N1 * N1

となり、再帰的に計算してやれば高々3回の計算を行えば8乗が求まります。

Pが奇数の時

Pが奇数(P=2k+1)の時は、

N2k+1 = N2k * N1

となります。1引いた数が偶数なのでこのように分解できます。

具体的な例で、P=7の時には

N7 = N6 * N1
N6 = N3 * N3
N3 = N2 * N1
N2 = N1 * N1

となります。やはりオーダーがlogになり、計算が高速化できることがわかります。

C# 再帰関数による実装

では上記の計算過程を再帰関数に落とし込んで実装してみます。なお計算結果は10⁹+7でMODをとるようにします。

static long Pow(long n, long p, int mod = (int)1e9 + 7)
{
    if (p == 0) return 1;　// nの0乗は1
    else if (p % 2 == 0)
    {
        var x = Pow(n, p / 2, mod);
        return (x * x) % mod;
    }
    else
    {
        return Pow(n, p - 1, mod) * n % mod;
    }
}

judge.u-aizu.ac.jp/onlinejudge/description.jsp?id=NTL_1_B&lang=jp

よさげな問題がAOJにあるので繰り返し二乗法を使って解いてみましょう。

C# 2進数で考える実装

再帰関数を使わずに実装することもできます。まず2進数表現を使ってPが2のべき乗の和でいくつになるかを考えます。

P=13のとき

P = (1101)2
  = 23 + 22 + 20
  = 8 + 4 + 1

となるので

N13 = N8+4+1
    = N8 * N4 * N1

表せます。

つまりにPについて、2進数表記の一桁目から見ていき、1になっている桁kがあれば、N^k を答えにかけ合わせればよいです。

2進数で下の桁から順にループして1になっているかどうかを判定する処理を書くと以下のような関数になります。

static long Pow(long n, long p, int mod = (int)1e9 + 7)
{
    long ans = 1;
    long now = n;
    for(; p > 0; p >>= 1)
    {
        if ((p & 1) == 1)
        {
            ans *= now;
            ans %= mod;
        }
        // nのべき乗を作っておく
        now *= now;
        now %= mod;
    }
    return ans;
}

上で実装した再帰関数と同じく、オーダーは logP になります。