[ruby-dev:32629] faster Bignum#*

遠藤と申します。

Bignum#* に FFT を用いた乗算を実装してみました。
10 進で 100000 桁以上あるような巨大な整数の乗算が速くなります。
そんな需要はないでしょうか。

$ cat ../test_fft.rb
require "benchmark"

n = 10**10000
puts Benchmark.measure { 100.times { n * n } }

n = 10**100000
puts Benchmark.measure { 10.times { n * n } }

n = 10**1000000
puts Benchmark.measure { 1.times { n * n } }

# 現状
$ ./ruby.org ../test_fft.rb
  2.590000   0.000000   2.590000 (  2.874859)
 26.040000   0.010000  26.050000 ( 28.927583)
266.540000   0.040000 266.580000 (276.491797)

# FFT 版
$ ./ruby.fft ../test_fft.rb
  2.630000   0.010000   2.640000 (  2.636339)
  2.790000   0.040000   2.830000 (  2.836020)
  6.270000   0.090000   6.360000 (  6.368478)

FFT の実装は教科書通りです。基数も 65536 (=2**(8*sizeof(short))) か
256 (=2**(8*sizeof(char))) の決め打ちです。フーリエ変換とかよく
わかってないので、バグってたらすみません。詳しい方に検証や再実装を
して頂けたら最高です。

# あと誰か、ニュートン法による除算を実装しません？


Index: bignum.c
===================================================================
--- bignum.c	(revision 14308)
+++ bignum.c	(working copy)
@@ -1442,10 +1442,145 @@
     }
 }

+static void
+big_fft(double *ar, double *ai, int al, double theta)
+{
+    int m, mh, i, j, k, irev;
+    double wr, wi, xr, xi;
+
+    i = 0;
+    for (j = 1; j < al - 1; j++) {
+	for (k = al >> 1; k > (i ^= k); k >>= 1);
+	if (j < i) {
+	    xr = ar[j];
+	    xi = ai[j];
+	    ar[j] = ar[i];
+	    ai[j] = ai[i];
+	    ar[i] = xr;
+	    ai[i] = xi;
+	}
+    }
+    for (mh = 1; (m = mh << 1) <= al; mh = m) {
+	irev = 0;
+	for (i = 0; i < al; i += m) {
+	    wr = cos(theta * irev);
+	    wi = sin(theta * irev);
+	    for (k = al >> 2; k > (irev ^= k); k >>= 1);
+	    for (j = i; j < mh + i; j++) {
+		k = j + mh;
+		xr = ar[j] - ar[k];
+		xi = ai[j] - ai[k];
+		ar[j] += ar[k];
+		ai[j] += ai[k];
+		ar[k] = wr * xr - wi * xi;
+		ai[k] = wr * xi + wi * xr;
+	    }
+	}
+    }
+}
+
+static void
+big_digits_to_bases(int base, double *xr, long rlen, BDIGIT *xds, long dlen)
+{
+    long i, j, k, div = 1 << (base * CHAR_BIT);
+    BDIGIT v;
+
+    k = 0;
+    for (i = 0; i < dlen; i++) {
+	v = xds[i];
+	for (j = 0; j < sizeof(BDIGIT) / base; j++) {
+	    xr[k++] = (double) (v % div);
+	    v /= div;
+	}
+    }
+    for (; k < rlen; k++) xr[k] = 0;
+}
+
+static void
+big_bases_to_digits(int base, BDIGIT *zds, long dlen, double *xr)
+{
+    long i, j, k, div = 1 << (base * CHAR_BIT);
+    BDIGIT v, s;
+    double x, xm;
+
+    x = 0;
+    i = j = 0;
+    v = 0;
+    s = 1;
+    for (k = 0;; k++) {
+	x += floor(xr[k] + 0.5);
+	xm = fmod(x, div);
+	v |= ((BDIGIT) xm) * s;
+	x = (x - xm) / div;
+	s *= div;
+	if (++j >= sizeof(BDIGIT) / base) {
+	    zds[i++] = v;
+	    if (i == dlen) return;
+	    j = 0;
+	    v = 0;
+	    s = 1;
+	}
+    }
+}
+
 static VALUE
+big_mul_fft(VALUE x, VALUE y)
+{
+    long i, len, base;
+    VALUE z;
+    double *xr, *xi, *yr, *yi, pi = atan(1.0) * 4.0;
+
+    /* calculate FFT length and decide the base of radix */
+    base = sizeof(short);
+    i = (RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1) * sizeof(BDIGIT) / base;
+    for (len = 1; len < i; len *= 2);
+    if (len > 131072) {
+	len *= sizeof(short) / sizeof(char);
+	base = sizeof(char);
+    }
+
+    /* allocate two arraies of complexes */
+    xr = ALLOC_N(double, len);
+    xi = ALLOC_N(double, len);
+    yr = ALLOC_N(double, len);
+    yi = ALLOC_N(double, len);
+
+    /* initialize the arraies */
+    big_digits_to_bases(base, xr, len, BDIGITS(x), RBIGNUM_LEN(x));
+    big_digits_to_bases(base, yr, len, BDIGITS(y), RBIGNUM_LEN(y));
+    for (i = 0; i < len; i++) xi[i] = yi[i] = 0;
+
+    /* perform FFT */
+    big_fft(xr, xi, len, 2 * pi / len);
+    big_fft(yr, yi, len, 2 * pi / len);
+
+    /* multiply the arraies */
+    for (i = 0; i < len; i++) {
+	double t = xr[i] * yr[i] - xi[i] * yi[i];
+	xi[i]    = xr[i] * yi[i] + xi[i] * yr[i];
+	xr[i]    = t;
+    }
+    free(yr);
+    free(yi);
+
+    /* perform IFFT */
+    big_fft(xr, xi, len, -2 * pi / len);
+    free(xi);
+    for (i = 0; i < len; i++) xr[i] /= len;
+
+    /* make a result bignum */
+    len = RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1;
+    z = bignew(len, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
+    big_bases_to_digits(base, BDIGITS(z), len, xr);
+    free(xr);
+
+    return z;
+}
+
+static VALUE
 rb_big_mul0(VALUE x, VALUE y)
 {
-    long i, j;
+    long i, j, nx, ny, len;
     BDIGIT_DBL n = 0;
     VALUE z;
     BDIGIT *zds;
@@ -1465,6 +1600,22 @@
 	return rb_num_coerce_bin(x, y);
     }

+    /* estimate time */
+    for (i = nx = 0; i < RBIGNUM_LEN(x); i++) if (BDIGITS(x)[i]) nx++;
+    for (i = ny = 0; i < RBIGNUM_LEN(y); i++) if (BDIGITS(y)[i]) ny++;
+    for (len = 1; RBIGNUM_LEN(x) + RBIGNUM_LEN(y) > len; len *= 2);
+
+    /* swap x and y (when y is more sparse than x) */
+    if (ny < nx) {
+	z = x; x = y; y = z; nx = ny;
+    }
+
+    /* multiply by FFT (if FFT seems to be faster) */
+    if ((double) nx * RBIGNUM_LEN(y) > (double) len * log(len) * 30) {
+	return big_mul_fft(x, y);
+    }
+
+    /* multiply in a normal way */
     j = RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1;
     z = bignew(j, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
     zds = BDIGITS(z);

-- 
Yusuke ENDOH <[e-mail]>