modpowがわからない

Question

###modpowを高速で計算したい

プログラミングの問題でn^m(modP)を計算する時、1回ずつnをm回かけるとTLEになってしまいました。
(Pythonだと勝手にやってくれる関数あるらしい)
調べてみると1回のループで1回かけるのではなく複数回かけるといいらしいというのを知りました。
そりゃそうですが

参照元URLでは、二進数とbit演算で高速に処理していました。

以下が参照元URLのコードです。

ll modpow(ll a, ll n=mo-2) {
	ll r=1;
	while(n) r=r*((n%2)?a:1)%mo,a=a*a%mo,n>>=1;
	return r;
}

いまいちやってることがわからないから出力してみた

とりあえず、それぞれの変数をcoutで各変数を出力してみました。
mはシフトしてるので2進数で出力しました。

ll modpow(ll n, ll m=mod-2) {
	ll r=1;
	while(m){
    cout << "r:" << r << endl;
    cout << "m:" << bitset<16>(m) << endl;
    cout << "n:" << n << endl << endl;
    r=r*((m%2)?n:1)%mod,n=n*n%mod,m>>=1;
  }
	return r;
}

出力結果は以下になります。
modpow(99,99):

わからないこと

m>>=1はシフト演算で右に1シフトしてるのはわかります。(桁を１つ小さくしてる感じ)

(m%2)?n:1)は3項演算子で今見てるbitが1(奇数)ならnを0なら1をrにかけるということだと思います。
n=n*n%modもnにnをかけてmodをとってる。
それぞれの式の意味はわかりますが、全体で何をやってるのかわからない。

分かる方教えていただけるとありがたいです。

Answer 1

この回答において記号^は累乗を意味します。

x^(a + b) = (x ^ a) * (x ^ b)を利用します。

x^(2n) = (x ^ n) ^ 2となり、x^(2n + 1) = x * ( (x ^ n) ^ 2)となります。これを利用すれば、計算を随分と省略できます。

例えば(x ^ 1024)の場合、

x ^ 2 → x ^ 4 →x ^ 8　→x ^ 16 →x ^ 32→x ^ 64→x ^ 128→x ^ 256 →x ^ 512→ x^ 1024

と大分計算回数が節約できます。二で割った余りで条件分岐しているのは端数処理のためです。

Comments

[episteme]

それに加えて: a * b (mod P) = (a mod P) * (b mod P) だから

「複数回のかけ算を済ませて最後に mod P」ではなく
「一回のかけ算ごとに mod P」することで
計算途中の値も常にP以下になる。 (途中で桁あふれを起こさない)

...を利用してますね。

[kyokio]

ありがとうございます。
わかってきました！

ですが、modpowの引数にm=mod-2としているのは逆元を計算にこの関数を使う際にmod-2を打たなくても良くするためと考えていいのでしょうか？

[episteme]

それは書いたヒトに訊いて。