繰り返し自乗法

概要

$N^P \pmod{M}$ の値を効率的に解きます。

使える状況

「 $N^P$ を $M$ で割ったあまり」が欲しいときに、 $P$ があまりにも大きいとその計算だけで時間がかかってしまうことがあります。そのような時には、この繰り返し自乗法を用いることで素早く求めることが出来ます。

説明

$N^P \pmod{M}$ の値を求めたいとき、愚直な方法として次のような手法を考えることが出来ます。

using ll = long long int;
ll PowMod(ll N, ll P, ll M){
    ll ans = 1;
    for(ll i=0; i<P; ++i){
        ans *= N;
        ans %= M;
    }
    return ans;
}

しかし、上記のような方法だと $P$ 回ループしてしまうため、 $P=10^9$ などの場合では非常に時間がかかってしまいます。

そこで使うのが、繰り返し自乗法です。

繰り返し自乗法は、

「 $N$ を一回一回掛けてmodを取るのでは時間がかかるので、まとめて掛けてしまおう！」

という考えをもとにしたアルゴリズムです。

この「まとめて掛ける」というのがどういうことか、具体的に $P=20$ のときについて考えながら説明します。

以下、求める値を $X$ 、 $N^{(i)}:=N^i \pmod{M}$ とし、また、" $\times$ "の計算をした時には $M$ で割った余りを求める操作も一緒に行っているものとします。

ボトムアップ(小さい値から大きい値へと計算していく)な考え方

先ほどの愚直な方法では、次のような計算をしていることになります。
$\displaystyle \begin{eqnarray} X &=& N^{(1)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\\ &=& N^{(2)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\\ &=& N^{(3)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\\ &&\vdots\\ &=& N^{(19)}\times N^{(1)}\\ &=& N^{(20)} \end{eqnarray}$
確かに、 $P$ 回、つまり $20$ 回分の計算をしていることがわかりますね。

ここで、上で計算しているように、
$N^{(2)} = N^{(1)}\times N^{(1)}$
となっていることから、
$\displaystyle \begin{eqnarray} X &=& N^{(1)}\times N^{(1)}\times N^{(1)}\times N^{(1)}\times \cdots \times N^{(1)}\times N^{(1)}\\ &=& (N^{(1)}\times N^{(1)})\times (N^{(1)}\times N^{(1)})\times \cdots \times (N^{(1)}\times N^{(1)})\\ &=& N^{(2)}\times N^{(2)}\times \cdots \times N^{(2)}\\ \end{eqnarray}$

とまとめてしまい、あとは、
$\displaystyle \begin{eqnarray} X &=& N^{(2)}\times N^{(2)}\times N^{(2)}\times N^{(2)}\times \cdots \times N^{(2)}\\ &=& N^{(4)}\times N^{(2)}\times N^{(2)}\times \cdots \times N^{(2)}\\ &=& N^{(6)}\times N^{(2)}\times \cdots \times N^{(2)}\\ &&\vdots\\ &=& N^{(18)}\times N^{(2)}\\ &=& N^{(20)} \end{eqnarray}$

と出来るので、計算回数は最初の $N^{(2)}$ を計算する $1$ 回と、まとめた後の $10$ 回となっているので、見事に半分にすることが出来ていますね。

これをさらに、
$N^{(4)}=N^{(2)}\times N^{(2)}$
とすることにより $N^{(4)}$ を求め、同様に $N^{(8)}$ を求め、…と繰り返すことで、最終的には
$\displaystyle \begin{eqnarray} X &=& N^{(16)} \times N^{(4)}\\ &=& N^{(20)} \end{eqnarray}$

とまで式を短くすることが出来ます。
このとき、次の $5$ 回しか計算をしていないことが分かりますね。

$N^{(2)}$ を求める( $N^{(1)}\times N^{(1)}$ )
$N^{(4)}$ を求める( $N^{(2)}\times N^{(2)}$ )
$N^{(8)}$ を求める( $N^{(4)}\times N^{(4)}$ )
$N^{(16)}$ を求める( $N^{(8)}\times N^{(8)}$ )
$N^{(20)}$ を求める( $N^{(16)}\times N^{(4)}$ )

このように、 $N$ の自乗を繰り返し計算した結果を用いて計算を行っているので、「繰り返し自乗法」と言うんですね。

トップダウン(大きい値から小さい値へと計算していく)な考え方

先ほどの方法とほとんど変わりませんが、実装が(少しだけ)簡単な方法もあります。
おそらく、一般的に「繰り返し自乗法」と言ったときにはこちらの方を指すのではないでしょうか。

この方法は、次のように再帰的な式で書き表すことが出来ます。
$\displaystyle N^P = \left\{ \begin{array}{ll} 1 & (P=0)\\ \left(N^{P/2}\right)^2 & (P:偶数)\\ N \times N^{P-1} & (P:奇数) \end{array} \right.$

この式を用いて、先ほどの $P=20$ の時についてを計算してみましょう。
$\displaystyle \begin{eqnarray} N^{(20)} &=& \left(N^{(10)}\right)^2\\ &=& \left(\left(N^{(5)}\right)^2\right)^2\\ &=& \left(\left(N^{(1)}\times N^{(4)}\right)^2\right)^2\\ &=& \left(\left(N^{(1)}\times \left(N^{(2)}\right)^2\right)^2\right)^2\\ &=& \left(\left(N^{(1)}\times \left(\left(N^{(1)}\right)^2\right)^2\right)^2\right)^2\\ \end{eqnarray}$

よって、先ほどと同様に $5$ 回の計算で $N^{(20)}$ を求めることが出来ています。

コード(C++)

using ll = long long int;
ll RepeatSquaring(ll N, ll P, ll M){
    if(P==0) return 1;
    if(P%2==0){
        ll t = RepeatSquaring(N, P/2, M);
        return t*t % M;
    }
    return N * RepeatSquaring(N, P-1, M);
}

コードの説明

トップダウンな考え方の節で示した式を、そのまま関数に落とし込んでいるだけです。

具体例

$2^{1000000000} \pmod{1000000007}$ の値を求める。

今回は、愚直な方法(関数PowMod_Simple)と繰り返し自乗法(PowMod_RepeatSquaring)でどのくらい性能の差があるかを、

再帰回数(再帰orループするたびにカウント)
処理時間(std::chronoを用いて計測)

の2点について見てみます。

#include <iostream>
#include <algorithm>
#include <chrono>

using ll = long long int;

ll count_ps = 0;
ll PowMod_Simple(ll N, ll P, ll M){
    ll ans = 1;
    for(ll i=0; i<P; ++i){
        ++count_ps;
        ans *= N;
        ans %= M;
    }
    return ans;
}

ll count_prs = 0;
ll PowMod_RepeatSquaring(ll N, ll P, ll M){
    ++count_prs;
    if(P==0) return 1;
    if(P%2==0){
        ll t = PowMod_RepeatSquaring(N, P/2, M);
        return t*t % M;
    }
    return N * PowMod_RepeatSquaring(N, P-1, M);
}

int main(){
    std::ios::sync_with_stdio(false);
    std::cin.tie(0);
    
    //愚直な方法-----------------------------
    std::cout << "愚直な方法\n";
    auto startTime = std::chrono::system_clock::now();
    //処理部------
    std::cout << "2^1000000000 mod 1000000007 = " << PowMod_Simple(2, 1000000000, 1000000007) << "\n";
    std::cout << "再帰回数\t:" << count_ps << "\n";
    //------------
    auto endTime = std::chrono::system_clock::now();
    std::cout << "処理時間:" << std::chrono::duration_cast<std::chrono::nanoseconds>(endTime - startTime).count() << "ns\n";
    //---------------------------------------
    
    std::cout << "------\n";
    
    //愚直な方法-----------------------------
    std::cout << "繰り返し自乗法\n";
    startTime = std::chrono::system_clock::now();
    //処理部------
    std::cout << "2^1000000000 mod 1000000007 = " << PowMod_RepeatSquaring(2, 1000000000, 1000000007) << "\n";
    std::cout << "再帰回数\t:" << count_prs << "\n";
    //------------
    endTime = std::chrono::system_clock::now();
    std::cout << "処理時間:" << std::chrono::duration_cast<std::chrono::nanoseconds>(endTime - startTime).count() << "ns\n";
    //---------------------------------------
    
    return 0;
}

実行結果(一例)

愚直な方法
2^1000000000 mod 1000000007 = 140625001
再帰回数	:1000000000
処理時間	:22322817871ns
------
繰り返し自乗法
2^1000000000 mod 1000000007 = 140625001
再帰回数	:43
処理時間	:2709ns

この結果からも分かるように、再帰回数、処理時間ともに、文字通り桁違いで小さく計算することが出来ていますね。
ちなみに、処理時間の「ns」は「ナノ秒(ナノ $=10^{-9}$ )」を表しているため、

愚直な方法 :22.322817871秒
繰り返し自乗法 : 0.000002709秒

の時間がかかっていることを表しています。

備考

愚直な方法は計算量が $O(P)$ であるのに対し、繰り返し自乗法では $O(\log P)$ にまで抑えられています。

あとがき

繰り返し自乗法は使いどころ(累乗の割った余りを求めたいとき)が分かりやすく、また実装も軽めなので、習得しやすいアルゴリズムだと思います。
一度この関連の問題を解いたら、また同じ場面に出会ったときに「あの再帰のやつだ…！」と思い出すことが出来ると思うので、ATC002Bなどで経験しておくと良いでしょう。

sataniC++

C++、競プロ、数学などについて色々考えたりします。