(17) 大きな整数の高速乗算 - 読み書きプログラミング

例えば、8ビット*8ビットの乗算器を持つマイコンを想像してください。このマイコンで2つのmバイトの整数の掛け算を行うには、普通、 $m^2$ 回、乗算器を使います。
この掛け算を、 $O(m\log m)$ のオーダーの乗算回数で行うアルゴリズムを考えます。

一般化して整数x, yを基数wで表すとすると、

$x=\sum_{j=0}^{m-1}{x_{j}w^{j}}$

$y=\sum_{j=0}^{m-1}{y_{j}w^{j}}$

2のべき乗で、2m以上の最小のものをnとすると、xとyの積zの基数wでの表現は以下になります。

$z_{k}=\sum_{j=0}^{n-1}{x_{j}y_{k-j}}$

これは畳み込み演算ですから、各配列を離散Fourier変換し、内積を取り、逆Fourier変換することで求めることができます。
離散Fourier変換は、要素の数が2のべき乗の時、乗算回数が $O(m\log m)$ で計算する高速Fourier変換アルゴリズムが発見されているので、このオーダーで大きな整数の積も計算できるという訳です。

$555221453977032891233$ と $442186118023559174$ の掛け算を、8ビットで計算することを想定して、このコンセプトが正しく動くことをMaximaで実験してみましょう。

/* 4.4.4m */
load("fft")$

int2listWithBase(i, b) := if i = 0 then [] else cons(mod(i, b), int2listWithBase(floor(i / b), b))$

x : 555221453977032891233$
y : 442186118023559174$

w : 256$
x : int2listWithBase(x, w)$
y : int2listWithBase(y, w)$
lenx : length(x)$
leny : length(y)$
len : max(lenx, leny)$
n : 2 ^ (1 + ceiling(float(log(len) / log(2))))$
for j : 1 thru n - lenx do x: endcons(0, x)$
for j : 1 thru n - leny do y: endcons(0, y)$
x : fft(x)$
y : fft(y)$
x : x * y$
x : inverse_fft(x)$
x : makelist(round(e), e, float(realpart(x)*n))$
z : sum(x[s] * w ^ (s - 1), s, 1, n);