2014年06月

FFT乗算

• 乗算したい多倍長配列AとBを用意し、FFT乗算では配列の長さは不変なので、A×Bの長さであるA.length+B.lengthより大きく、2の乗数になるようにAとBの配列を増やして0で埋める。
• AとBをそれぞれフーリエ変換したものをaとbとし、i番目の配列要素についてc[i] = a[i]×b[i]を計算する。ただしaとbは複素数になっているので実部と虚部をそれぞれ計算することになる。
• cを逆フーリエ変換してCを得て、Cの各配列要素は多倍長数の各配列の基数よりも大きくなっている部分があるので、基数で剰余をとって繰り上げ計算をして乗算結果が求められる。

フーリエ変換後の乗算は同じ位置の配列要素同士の乗算で済み、計算コストはAとBの長さの積ではなくAとBの長さの和になるためとても高速である。

FFT累乗+修正バイナリ法

FFT乗算の性質上同じ数同士を乗算し続けるほうが高速である。2乗計算をし続ける分には同じ数同士の乗算となるので最適になる。

例えば2^100を2×2×2×・・・を100回繰り返すよりも、100=64+32+4であることを利用して、2^2を2乗し続けて2^4、2^8、2^16、2^32、2^64と計算したあとに2^64×2^32×2^4 = 2^100とするほうが速い・・・と言われるがまだ効率的ではない。2^64までの計算はお互いに同じ数同士なので最高速度となるが、最後の2^64×2^32×2^4はお互いに違う数で、順番を見極めないと余分な配列を必要としてしまう。

そんな手間をかけるぐらいなら次のように計算するべきだ。2^1→2^2→2^4→2^8→2^7→2^14→2^13→2^26→2^25→2^50→2^100。ゴールが2^100になるように、途中で適宜2で割る計算(コストは長さに等しい)を挟んで、2乗計算のみで目標の累乗を目指すようにする。問題点は累乗する数が大きい場合、多倍長除算のコストはとても大きいので非効率になる。累乗する数に大きさに応じて計算方法を振り分けるとよい。

FFT階乗

階乗計算は馬鹿正直に1〜nまでを順番に乗算するよりもFFT乗算が最高効率で働くように乗算する数を調節するほうがよい。単純な連続乗算の場合、n!は急激に大きくなるのにnは多倍長でさえない小さな数なので、FFT乗算の恩恵は一切受けられない。

詳しくはJavaScriptによるFFT乗算(2)に書かれている。予めFFT乗算をしやすいように、2の乗数個の多倍長配列にnまでの積を分割して乗算しておき、2の乗数個の多倍長配列をFFT乗算で求めていく。

修正バイナリ法に関連したニュートン法の改良

ニュートン法は1回のループで精度が2倍になる。精度100桁の結果が欲しい時、100桁に最も近いループを選ぶことになるので単純なケースでは最初の精度を2とするなら2→4→8→16→32→64→128で計算を終えるが、実はこの計算は累乗の項で書いた方法がそのまま適用でき、適宜精度を削って最後に100になるようにするほうが平均的に速い。100桁の場合128桁で計算が終わるため無駄に計算する桁数が28桁だが、130桁欲しい時は256桁まで求めることになるので効率が悪い。従って2→4→8→7→14→26→25→50→100となるように途中で精度を1桁削るほうがトータルでは無駄がなくてよい。

HugeIntにおけるFFT除算とは、内部にFFT乗算を使ったニュートン法のことである。

b進数変換

多倍長10進数をb進数に変換する操作は、多倍長配列の基数(配列の区切り)を変えることと同じだ。通常多倍長配列はN桁区切りとなっているが、これは10^Nを基数とする配列で、10^N進数と呼んでもよい。10進数ではないb進数に変換するためには、bを基数とする配列を作る。そのためには多倍長配列をbで割った剰余を求めなければいけない。

何も考えなければ、多倍長配列をbで割り続けた剰余を新しい配列にしていくとb進数に変換される。しかし巨大な多倍長配列に対して小さな基数bで除算を行ってもFFTの恩恵を受けられないので工夫する。

小さなbで割るのではなく、割った時の商が除数よりも小さくなるような、大きなb^2^nを選んで割る。具体的にはb^2^(n+1)>多倍長数>b^2^nを選ぶ。できた商と剰余は共にb^2^nより小さく、商と剰余を配列にすればこれはb^2^n進数になっている。次はこれをそれぞれb^2^(n-1)で割っても出てくる商はb^2^(n-1)よりも小さく、商と剰余を並べた配列はb^2^(n-1)進数の多倍長配列になる。これをn回繰り返して、b^2^0進数の配列を作ると、求めたいb進数が得られる。

b進数の10進数変換

b進数を10進数に変換する時は愚直にやると、b進数の各i桁目×b^iを計算して全て乗算する。しかしこれが非効率なのは言うまでもなく、効率的なやり方は上の手順を逆回しにすることになる。

b進数の多倍長配列Bについて、B[2i]+B[2i+1]×bを求めると、配列の数は半分のb^2進数になる。b^2進数の多倍長配列B2に対しても同じようにB2[2i]+B2[2i+1]×b^2を計算してb^4進数を作り、これを配列の数が1つになるまで繰り返すとできあがるb^2^n進数はそのまま多倍長10進数そのものになっている。bの累乗計算はFFTで高速に行うことができ、b^2^k進数の配列要素にb^2^kを掛ける計算もお互いに桁数が等しいのでFFTで高速に求めることができる。

変換デモ

関数概要

encode関数の入力は1バイト長の数値配列で出力はbase64の文字列、decode関数の入力はbase64の文字列で出力は1バイト長の数値配列。入出力が1バイト長の数値配列になっているのは作っている別のプログラム向けの仕様で、文字列を1バイト長の数値配列に変換する方法は下を参照。

base64 = {}
base64.encode = function(input) {
        var base64list = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
        var powlist = [1,2,4,8,16,32,64,128,256];
        var base64_code = String();
        var pi = 0;
        var i = 0;
        var tmp = 0;
        while(i<input.length || pi>6) {
                if(pi<6) {
                        tmp = (tmp<<8) | (input[i++]&255);
                        pi += 8;
                }
                base64_code += base64list.charAt(tmp>>(pi-6));
                tmp &= (powlist[pi-6]-1);
                pi -= 6;
        }
        tmp <<= (6-pi);
        base64_code += base64list.charAt(tmp);
        switch(base64_code.length%4) {
                case 1 : base64_code+="===";break;
                case 2 : base64_code+="==";break;
                case 3 : base64_code+="=";break;
                default: break;
        }
        return base64_code;
}


base64.decode = function(input) {
        var base64list = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/=";
        var powlist = [1,2,4,8,16,32,64,128,256];
        var decode_ary = new Array();
        var i = 0;
        var pi = 0;
        var tmp = 0;
        var j;
        while(i<input.length || pi>8) {
                j = base64list.indexOf(input[i++],0);
                if(j<64) {
                        tmp = (tmp<<6) | (j&255);
                        pi += 6;
                }
                if(pi>=8) {
                        decode_ary.push(tmp>>(pi-8));
                        tmp &= (powlist[pi-8]-1);
                        pi -= 8
                }
        }
        return decode_ary;
}

1バイトの数値配列を使うため、そのままでは日本語などの2バイト文字文字を扱えない。これの処理方法がなかなかWeb上に見つからなかったので書いておく。

入力された文字列input_textはそのままコードにすると2バイトになってしまう。1バイト長の文字列に変換するためにunescape(encodeURIComponent(input_text))を通す。この操作で1文字2バイトの文字はエスケープされ1文字1バイトの3文字に変換される。元に戻す時はdecodeURIComponent(escape(decode_text))と逆に通す。そして文字列とコードのそれぞれの変換はcharCodeAtとfromCharCodeを使う。以下は入力された文字列を1バイト長の数値配列に変換するものと、1バイト長の数値配列を元の文字列に変換するもの。

        var decode_text = unescape(encodeURIComponent(input_text));
        var decode_ary = new Array();
        for(var i=0; i<decode_text.length; i++)
                decode_ary[i] = decode_text.charCodeAt(i);

        var decode_text = String();
        for(var i=0; i<decode_ary.length; i++)
                decode_text += String.fromCharCode(decode_ary[i]);
        var output_text = decodeURIComponent(escape(decode_text));