繰り返し二乗法について

はじめに

この度私のライブラリに冪乗の余りを求める、いわゆるmodPow関数を追加しました。 丁度いい機会ということで、本稿ではmodpowの仕組みを説明し、詳細な実装方法と実装例を示します。 数学をガバっていたらすみません。

modpowとは

以下、$\times{}$は整数の積を表します。

まず、剰余を定める重要な定理を紹介します。

(証明略)

ここで、$q, r$をそれぞれ、「$n$を$m$で割った商」、「$n$を$m$で割った非負最小剰余」と定めます。

次に、$\forall{}x \in{} \mathbb{Z}$と非負整数$n$に対して、冪乗$x^n$を次で定めます。 $$ x^n = \begin{cases} x \times{} x^{n-1} & \text{if $n \geq{} 1$} \\ 1 & \text{if $n = 0$} \end{cases} $$

本稿で紹介するmodPow関数で求める値は、$a^x$をMODで割った非負最小剰余です。

(数式苦手な方へ)これは何を求めているのかをできるだけ形式的に定めているだけなので、あまり気にしなくても良いです。 「$a$ を $x$ 回掛け算して、MODで割った余りを求める。」でOKです。

(注意): 便宜上、$0^0=1$ としています。

計算原理

以下、割る数は$0$を含まない正整数とします。 割る数とはプログラムにおけるMODのことです。 (定義を考えれば負の整数に対しても非負最小剰余はただ一つに定まりますが、実用上正整数を用いることが多いのかな？)

通常、冪乗は比較的簡単に非常に大きな値になります。 例えば、$2^{10}=1024$ ですが、$2^{64}=18446744073709551616$ であり、 上の定義をそのまま計算するだけでは非常に厳しいことがわかります。 そこで、まず次の定理が大事になります。

この定理は、積に関してなら どのタイミングで非負最小剰余に変換しても最後は同じ値になるということを示しています。 また、定理では右辺は両方非負最小剰余に変換していますが、もちろん片方だけで行っても同様の結果が得られます。

これで $a^x$ の大きさに関する問題は解決します。 なぜなら、一回 $a$ を掛け算するたびに非負最小剰余に逐次変換していけば良いので、 $(割る数) \times{} (割る数)$さえ正しく計算できれば $a^x$ の非負最小剰余もまた計算可能になるからです。 プログラム的には、用いる整数の型が $k$ bitであるとき、MODの大きさが $k/2$ bit程度であれば計算できます。

したがって、次のようなプログラムは正しく動作します。

long long modPow (long long a, long long x, const int MOD) {
    // 不正な入力を弾く
    assert(0 <= x);
    assert(1 <= MOD);

    // aを正規化
    a %= MOD; a += MOD; a %= MOD;

    long long res = 1;
    for (long long i = 0; i < x; i++) {
        res *= a;
        res %= MOD;
    }

    return res % MOD;
}

$a$ の正規化というのは、底を最初の段階で非負最小剰余に変換しているということです。 上のコードは積や剰余演算を$O(1)$でできると仮定すれば、全体で$O(x)$で抑えられます。 しかし、(競技プログラミングなど)実用上はxは$10^9$程度であったりすることが多く、このままでは使えません。 そこで、通常は次に説明するような高速化を施します。

まず、以下の定理を導入します。

少々天下り的ですが、$p=2$ を選んで、$x = \sum_{i=0} 2^i c_i ~ (c_i \in{} {0, 1})$ とします。 定理より、非負整数 $x$ はこのような表示を先頭の $0$ の自由を除いて一意に持ちます。 このとき、 $$ \begin{split} a^x &= a^{\sum_{i=0} 2^i c_i} \\ &= a^{c_0 + 2 c_1 + 2^2 c_2 + \dots{}} \\ &= a^{c_0} \times{} a^{2 c_1} \times{} a^{2^2 c_2} \times{} \dots{} \\ &= \prod_{i=0} a^{ 2^i c_i} \\ &= \prod_{i=0} \left( a^{2^i} \right) ^{c_i} \end{split} $$ となります。 (式変形がわからない人はこちらを参考にしてください。)

実は、この形まで変形すると非常に高速に計算できるようになっています！ 先程の定義通りに計算する方法では、$x$ が大きいときに時間がかかるということが問題でした。 上の変形は、視覚的に説明するなら、 $$ a^{100} = a^{64} \times{} a^{32} \times{} a^{4} $$ のように分解しているということです。 定理3は、このような分解をしたときに、右辺の指数が必ず $2^k$ と表せるということを保証しているとも言えます。 この定理のおかげで $a^{2^k}$ の形の数さえ高速に列挙できれば良くなり、結果的に冪乗を求めるのも高速になるということです。

また、$a^{2^k} = a^{2^{k-1}} \times{} a^{2^{k-1}}$ という関係が成立するため、$a^{2^k}$ の形の数は簡単に求められます。 以上より、次のアルゴリズムを得ます。

このアルゴリズムの時間計算量は、手順3から手順6のループ1回あたり $O(1)$ とみなせば、 $x < 2^k$ となった時点で停止するため($k$ 以降の数 $l$ において $c_l = 0$ となるから)、全体で $O(log(x))$ になります。

元々の $O(x)$ のアルゴリズムからかなり改善されました！

実装

実は、$x$ を2進法展開するのはほとんどのプログラミング言語で必要ないです。 というのは、コンピュータは内部的に整数を2進法で表現しているからです。 ビット演算と言われるような機能を持つプログラミング言語ならこの過程を飛ばすことができます。

具体的には、$x$ を2進法展開したときの $c_0$ は、プログラム上ではx&1 で得ることができます。 更に、x >>= 1などで「ビットシフト」をすると、次にx&1をしたときには $c_1$ が得られます。 (ただし、このあたりはプログラミング言語によります。)

負の数が関わると2の補数表現など少しややこしくなりますが、$0\leq{} x$ を仮定しているので問題ありません。

以下にC言語、C++、D言語、python3での実装例を示します。

C言語/C++

#include <assert.h>
// C言語なら<assert.h>
// C++なら<cassert.h>

long long modPow (long long a, long long x, const int MOD) {
    // assertion
    assert(0 <= x);
    assert(1 <= MOD);

    // normalize
    a %= MOD; a += MOD; a %= MOD;

    // calculate
    long long ans = 1L % MOD;
    long long base = a;
    while (x != 0) {
        if ((x&1) != 0) {
            ans *= base; ans %= MOD;
        }
        base = base*base; base %= MOD;
        x >>= 1;
    }

    return ans;
}

D言語

long modPow (long a, long x, const int MOD) {
    // assertion
    assert(0 <= x);
    assert(1 <= MOD);

    // normalize
    a %= MOD; a += MOD; a %= MOD;

    // calculate
    long ans = 1L % MOD;
    long base = a % MOD;
    while (x != 0) {
        if ((x&1) != 0) {
            ans *= base; ans %= MOD;
        }
        base = base*base; base %= MOD;
        x >>= 1;
    }

    return ans;
}

python3

def modPow (a, x, MOD):
    # assertion
    assert 0 <= x, "x must be an integer greater than or equal to 0"
    assert 1 <= MOD, "MOD must be an integer greater or equal to 1"

    # normalize
    a %= MOD

    # calculate
    ans = 1 % MOD
    base = a
    while x != 0:
        if (x&1) != 0:
            ans *= base
            ans %= MOD
        base = base*base
        base %= MOD
        x >>= 1

    return ans

終わりに

まとめるのすごく大変だった...

内容に不備があれば著者のtwitterに連絡していただければ助かります。 どれだけ些細な内容でも大歓迎です。

参考文献

本稿における主張、紹介した定理及びその証明はほとんど以下に依ります。

1. 尾関和彦、 情報技術のための離散型数学入門、 共立出版(2023)

主に第5章、第6章を参考にしました。