背景
繰り返し二乗法を行う場合、たとえ結果が32bit整数に収まるとしても、一度64bit整数を介さなければ、オーバーフローのため正しく計算できなくなる。
しかし、言語によっては64bit整数が提供されていない場合がある。32bitを越えると自動的に浮動小数点型に変換される場合が多い。
剰余にもよるが、大抵は誤差が発生してしまうだろう。
方法
繰り返し二乗法の基本的な 擬似 コードを示す。
var mod=123456789;
var n=1234567;
var p=2345678;
var r=1;
for(;p;p>>=1){
if(p&1)r=r*n%mod;
n=n*n%mod;
}
console.log(r);
JavaScriptの文法を用いて書いている。すなわち、このコードは誤差のため 正しくない 。
これをなんとかするには、このrとnに桁区切りを設ける。
$ factor 123456789 #=> 123456789: 3 3 3607 3803なので、10821と11409を使う。
r=10821a+b, n=11409c+dと書くとすると、r*nは123456789ac+10821ad+11409bc+bdと表せる。
欲しいのは 123456789で割った余り なのだから、1項目はまるまる不要ということになる。
bをくくりだすと、計算すればよいのはb*n + a*d*10821ということになる。
これを用いて正しいコードを記述することができる。
var mod=123456789;
var F1=10821;
var F2=mod/F1; //11409
var n=1234567;
var p=2345678;
var r=1;
for(;p;p>>=1){
if(p&1)r=(r%F1*n + (r/F1^0)*(n%F2)*F1)%mod;
n=(n%F1*n + (n/F1^0)*(n%F2)*F1)%mod;
}
console.log(r);
条件
b*n + a*d*10821という式において、aとdは11409未満、bは10821未満である。
つまりこの式はmod * (10821+11409)未満ということになる。double型の仮数部は52bitなので、「modの因数があるとき(Fとする)、mod * (F+mod/F)が1<<52未満」なら条件を満たせることになる。
そのため、実は上のコードでF1=9としても問題なく動作する。
まとめ
本記事では剰余がある程度大きい因数を持つならば、浮動小数点型を用いて繰り返し二乗法を適用できることを示した。
150811追記
このテクニックは場合によりC等に応用できうる。
HackerRank ProjectEuler+ 97。この問題では剰余が10**12である。よって、C言語で普通に解くとしても64bit整数に収まらない。
ところが、この手法を使うと見事10**12 * 2*10**6は1<<63未満になるので、問題なく解くことができる。
尤もJavaScript/AWKでこの問題が解けるのかはわかりませんが…。
(HackerRankはPythonは実行時間5倍のgraceがあるのでpow()使えば普通に突破できるのは不問で。Rubyでは解けませんでしたので。)