From: Yusuke ENDOH Date: 2007-12-19T00:18:17+09:00 Subject: [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