# [1.9:Feature] faster bignum multiplication by karatsuba method

e\$B1sF#\$G\$9!#e(B

Python e\$B\$N%=!<%9\$r;29M\$K\$7\$F!"e(Bbignum e\$B\$N>h;;\$re(B Karatsuba
e\$BK!\$J\$I\$Ge(B
e\$B9bB.2=\$7\$F\$_\$^\$7\$?!#e(B

e\$B%Y%s%A%^!<%/\$9\$k\$H\$3\$s\$J46\$8\$G\$9!#e(B

``````                                    e\$B5le(B    e\$B?7e(B
``````

e\$BG\N(e(B
x = 10 ** 10; 500000.times { x * x } 1.11s 0.99s 112.13%
x = 10 ** 20; 500000.times { x * x } 1.28s 1.25s 102.52%
x = 10 ** 30; 500000.times { x * x } 1.52s 1.36s 111.19%
x = 10 ** 1000; 10000.times { x * x } 2.42s 1.51s 159.97%
x = 31000; 50000.times { x * x } 4.25s 2.65s 160.07%
x = 2
30000; 50000.times { x * x } 4.49s 2.42s 185.65%
50.times { 161000000 } 4.37s 4.21s 103.72%
10.times { 10
100000 } 5.06s 1.98s 255.85%
(12345 ** 12345) * (54321 ** 54321) 8.84s 2.57s 344.45%

e\$B\$a\$C\$?\$K;H\$o\$J\$5\$=\$&\$JBg\$-\$JCM\$GFC\$K8z2L\$,\$“\$j\$^\$9\$,!”>.\$5\$Je(B
e\$BCM\$G\$bCY\$/\$O\$J\$i\$J\$\$\$h\$&\$K9)IW\$7\$?\$D\$b\$j\$G\$9!#e(B

e\$B%Q%C%A\$NBg\$^\$+\$JFbMF\$O0J2<\$N\$H\$*\$j\$G\$9!#e(B

(bigmul1_karatsuba e\$B\$J\$I\$G;H\$&e(B)
• bigmul1_normal(x, y)

• e\$BI.;;e(B (e\$B=>Mh\$HF1\$8e(B)
• bigmul1_balance(x, y)

• y e\$B\$N7e\$re(B x e\$B\$N7e?t\$GJ,3d\$7\$F!"8D!9\$K>h;;\$7\$?7k2L\$rB-\$9e(B
• bigmul1_karatsuba(x, y)

• bigsqr_fast(x)

• big_sparse_p(x)

• sparse e\$B\$+\$I\$&\$+\$rH=CG\$9\$ke(B
(e\$B7e\$re(B 3 e\$B\$D%i%s%@%`Cj=P\$7!"e(B0 e\$B\$N?t\$,e(B 1
e\$B8D0J2<\$J\$ie(B sparse)
• bigmul0(x, y)

• x e\$B\$+e(B y e\$B\$,e(B KARATSUBA_MUL_DIGITS (= 70)
e\$B7e\$h\$j>.\$5\$1\$l\$Pe(B
bigsqr_fast e\$B\$^\$?\$Oe(B bigmul1_normal
• x e\$B\$+e(B y e\$B\$,e(B sparse (e\$B\$[\$H\$s\$I\$N7e\$,e(B 0) e\$B\$G\$be(B
bigsqr_fast e\$B\$^\$?\$Oe(B
bigmul1_normal
• x e\$B\$He(B y e\$B\$N7e?t\$N:9\$,FsG\0J>e\$J\$ie(B bigmul1_balance
• e\$B\$=\$l0J30\$J\$ie(B bigmul1_karatsuba

e\$B0l1~e(B make test && make test-all e\$B\$ODL\$C\$F\$\$\$^\$9e(B (bignum
e\$B\$N%F%9%H\$Oe(B
e\$B=<<B\$7\$F\$J\$\$\$N\$G\$"\$^\$jEv\$F\$K\$J\$i\$J\$\$\$G\$9\$,e(B) e\$B!#e(B

e\$B\$I\$&\$G\$7\$g\$&\$+!#e(B

# Index: bignum.c

— bignum.c (revision 20645)
+++ bignum.c (working copy)
@@ -17,6 +17,7 @@
#ifdef HAVE_IEEEFP_H
#include <ieeefp.h>
#endif
+#include <assert.h>

VALUE rb_cBignum;

@@ -1380,12 +1381,36 @@
return bignorm(z);
}

+static void
+bigsub_core(BDIGIT *xds, long xn, BDIGIT *yds, long yn, BDIGIT *zds,
long zn)
+{

• BDIGIT_DBL_SIGNED num;
• long i;
• for (i = 0, num = 0; i < yn; i++) {
• num += (BDIGIT_DBL_SIGNED)xds[i] - yds[i];
• zds[i] = BIGLO(num);
• num = BIGDN(num);
• }
• while (num && i < xn) {
• num += xds[i];
• zds[i++] = BIGLO(num);
• num = BIGDN(num);
• }
• while (i < xn) {
• zds[i] = xds[i];
• i++;
• }
• assert(i <= zn);
• while (i < zn) {
• zds[i++] = 0;
• }
+}

static VALUE
bigsub(VALUE x, VALUE y)
{
VALUE z = 0;

• BDIGIT *zds;

• BDIGIT_DBL_SIGNED num;
long i = RBIGNUM_LEN(x);

/* if x is larger than y, swap */
@@ -1406,32 +1431,52 @@
}

z = bignew(RBIGNUM_LEN(x), z==0);

• zds = BDIGITS(z);

• bigsub_core(BDIGITS(x), RBIGNUM_LEN(x),
• BDIGITS(y), RBIGNUM_LEN(y),
• BDIGITS(z), RBIGNUM_LEN(z));
• for (i = 0, num = 0; i < RBIGNUM_LEN(y); i++) {
• num += (BDIGIT_DBL_SIGNED)BDIGITS(x)[i] - BDIGITS(y)[i];
• zds[i] = BIGLO(num);
• return z;
+}

+static void
+bigadd_core(BDIGIT *xds, long xn, BDIGIT *yds, long yn, BDIGIT *zds,
long zn)
+{

• BDIGIT_DBL num = 0;
• long i;
• if (xn > yn) {
• BDIGIT *tds;
• tds = xds; xds = yds; yds = tds;
• i = xn; xn = yn; yn = i;
• }
• i = 0;
• while (i < xn) {
• num += (BDIGIT_DBL)xds[i] + yds[i];
• zds[i++] = BIGLO(num);
num = BIGDN(num);
}
• while (num && i < RBIGNUM_LEN(x)) {
• num += BDIGITS(x)[i];
• while (num && i < yn) {
• num += yds[i];
zds[i++] = BIGLO(num);
num = BIGDN(num);
}
• while (i < RBIGNUM_LEN(x)) {
• zds[i] = BDIGITS(x)[i];
• while (i < yn) {
• zds[i] = yds[i];
i++;
}
• return z;
• if (num) zds[i++] = (BDIGIT)num;
• assert(i <= zn);
• while (i < zn) {
• zds[i++] = 0;
• }
}

static VALUE
bigadd(VALUE x, VALUE y, int sign)
{
VALUE z;

• BDIGIT_DBL num;
• long i, len;
• long len;

sign = (sign == RBIGNUM_SIGN(y));
if (RBIGNUM_SIGN(x) != sign) {
@@ -1441,30 +1486,15 @@

if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
len = RBIGNUM_LEN(x) + 1;

• z = x; x = y; y = z;
}
else {
len = RBIGNUM_LEN(y) + 1;
}
z = bignew(len, sign);

• len = RBIGNUM_LEN(x);

• for (i = 0, num = 0; i < len; i++) {

• num += (BDIGIT_DBL)BDIGITS(x)[i] + BDIGITS(y)[i];

• BDIGITS(z)[i] = BIGLO(num);

• num = BIGDN(num);

• }

• len = RBIGNUM_LEN(y);

• while (num && i < len) {

• num += BDIGITS(y)[i];

• BDIGITS(z)[i++] = BIGLO(num);

• num = BIGDN(num);

• }

• while (i < len) {

• BDIGITS(z)[i] = BDIGITS(y)[i];

• i++;

• }

• BDIGITS(z)[i] = (BDIGIT)num;

• BDIGITS(y), RBIGNUM_LEN(y),

• BDIGITS(z), RBIGNUM_LEN(z));

return z;
}
@@ -1519,24 +1549,20 @@
}
}

-static void
-rb_big_stop(void *ptr)
+static long
+big_real_len(VALUE x)
{

• VALUE stop = (VALUE)ptr;
• *stop = Qtrue;
• long i = RBIGNUM_LEN(x);
• while (–i && !BDIGITS(x)[i]);
• return i + 1;
}

-struct big_mul_struct {

• VALUE x, y, z, stop;
-};

static VALUE
-bigmul1(void *ptr)
+bigmul1_normal(VALUE x, VALUE y)
{

• struct big_mul_struct bms = (struct big_mul_struct)ptr;
long i, j;
BDIGIT_DBL n = 0;
• VALUE x = bms->x, y = bms->y, z = bms->z;
• VALUE z = bignew(RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1,
RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
BDIGIT *zds;

j = RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1;
@@ -1544,7 +1570,6 @@
while (j–) zds[j] = 0;
for (i = 0; i < RBIGNUM_LEN(x); i++) {
BDIGIT_DBL dd;

• if (bms->stop) return Qnil;
dd = BDIGITS(x)[i];
if (dd == 0) continue;
n = 0;
@@ -1558,45 +1583,251 @@
zds[i + j] = n;
}
}
return z;
}

+static VALUE bigmul0(VALUE x, VALUE y);
+
+/* balancing multiplication by slicing larger argument */
static VALUE
-rb_big_mul0(VALUE x, VALUE y)
+bigmul1_balance(VALUE x, VALUE y)
{

• struct big_mul_struct bms;
• volatile VALUE z;
• VALUE z, t1, t2;
• long i, xn, yn, r, n;
• switch (TYPE(y)) {
• `````` case T_FIXNUM:
``````
• y = rb_int2big(FIX2LONG(y));
• break;
• xn = RBIGNUM_LEN(x);
• yn = RBIGNUM_LEN(y);
• assert(2 * xn <= yn);
• `````` case T_BIGNUM:
``````
• break;
• z = bignew(xn + yn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
• t1 = bignew(xn, 1);
• `````` case T_FLOAT:
``````
• return DBL2NUM(rb_big2dbl(x) * RFLOAT_VALUE(y));
• for (i = 0; i < xn + yn; i++) BDIGITS(z)[i] = 0;
• `````` default:
``````
• return rb_num_coerce_bin(x, y, ‘*’);
• n = 0;
• while (yn > 0) {
• r = xn > yn ? yn : xn;
• MEMCPY(BDIGITS(t1), BDIGITS(y) + n, BDIGIT, r);
• RBIGNUM_SET_LEN(t1, r);
• t2 = bigmul0(x, t1);
• bigadd_core(BDIGITS(z) + n, RBIGNUM_LEN(z) - n,
• ``````   BDIGITS(t2), big_real_len(t2),
``````
• ``````   BDIGITS(z) + n, RBIGNUM_LEN(z) - n);
``````
• yn -= r;
• n += r;
}
• rb_gc_force_recycle(t1);
• bms.x = x;
• bms.y = y;
• bms.z = bignew(RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1,
RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
• bms.stop = Qfalse;
• return z;
+}
• if (RBIGNUM_LEN(x) + RBIGNUM_LEN(y) > 10000) {
• z = rb_thread_blocking_region(bigmul1, &bms, rb_big_stop, &bms.stop);
+/* split a bignum into high and low bignums */
+static void
+big_split(VALUE v, long n, VALUE *ph, VALUE *pl)
+{
• long hn, ln;
• VALUE h, l;
• ln = RBIGNUM_LEN(v) > n ? n : RBIGNUM_LEN(v);
• hn = RBIGNUM_LEN(v) - ln;
• while (–hn && !BDIGITS(v)[hn + ln]);
• h = bignew(++hn, 1);
• MEMCPY(BDIGITS(h), BDIGITS(v) + ln, BDIGIT, hn);
• while (–ln && !BDIGITS(v)[ln]);
• l = bignew(++ln, 1);
• MEMCPY(BDIGITS(l), BDIGITS(v), BDIGIT, ln);
• *pl = l;
• *ph = h;
+}

+/* multiplication by karatsuba method */
+static VALUE
+bigmul1_karatsuba(VALUE x, VALUE y)
+{

• long i, n, xn, yn, t1n, t2n;
• VALUE xh, xl, yh, yl, z, t1, t2, t3;
• BDIGIT *zds;
• xn = RBIGNUM_LEN(x);
• yn = RBIGNUM_LEN(y);
• n = yn / 2;
• big_split(x, n, &xh, &xl);
• if (x == y) {
• yh = xh; yl = xl;
}
• else big_split(y, n, &yh, &yl);
• /* x = xh * b + xl
• ``````* y = yh * b + yl
``````
• ``````*
``````
• ``````* Karatsuba method:
``````
• ``````*   x * y = z2 * b^2 + z1 * b + z0
``````
• ``````*   where
``````
• ``````*     z2 = xh * yh
``````
• ``````*     z0 = xl * yl
``````
• ``````*     z1 = (xh + xl) * (yh + yl) - x2 - x0
``````
• ``````*
``````
• ``````*  ref: http://en.wikipedia.org/wiki/Karatsuba_algorithm
``````
• ``````*/
``````
• /* allocate a result bignum */
• z = bignew(xn + yn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
• zds = BDIGITS(z);
• for (i = 0; i < xn + yn; i++) zds[i] = 0xdfdfdfdf;
• /* t1 ← xh * yh */
• t1 = bigmul0(xh, yh);
• t1n = big_real_len(t1);
• /* copy t1 into high bytes of the result (z2) */
• MEMCPY(zds + 2 * n, BDIGITS(t1), BDIGIT, t1n);
• for (i = 2 * n + t1n; i < xn + yn; i++) BDIGITS(z)[i] = 0;
• /* t2 ← xl * yl */
• if (BIGZEROP(xl) || BIGZEROP(yl)) {
• `````` t2 = xl;
``````
• t2n = 0;
• }
else {
• z = bigmul1(&bms);
• `````` t2 = bigmul0(xl, yl);
``````
• `````` t2n = big_real_len(t2);
``````

}

• /* copy t2 into low bytes of the result (z0) */

• MEMCPY(zds, BDIGITS(t2), BDIGIT, t2n);

• for (i = t2n; i < 2 * n; i++) BDIGITS(z)[i] = 0;

• /* subtract t2 from middle bytes of the result (z1) */

• i = xn + yn - n;

• bigsub_core(zds + n, i, BDIGITS(t2), t2n, zds + n, i);

• bigsub_core(zds + n, i, BDIGITS(t1), t1n, zds + n, i);

• rb_gc_force_recycle(t1);

• if (t2n) rb_gc_force_recycle(t2);

• /* t1 ← xh + xl */

• t1 = bigadd(xh, xl, 1);

• /* t2 ← yh + yl */

• t2 = (x == y) ? t1 : bigadd(yh, yl, 1);

• /* t3 ← t1 * t2 */

• t3 = bigmul0(t1, t2);

• rb_gc_force_recycle(t1);

• if (t1 != t2) rb_gc_force_recycle(t2);

• /* add t3 to middle bytes of the result (z1) */

• bigadd_core(zds + n, i, BDIGITS(t3), big_real_len(t3), zds + n, i);

• rb_gc_force_recycle(t3);

• return z;
}

+/* efficient squaring (2 times faster than normal multiplication)

• ref: Handbook of Applied Cryptography, Algorithm 14.16
• `````` http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
``````
• */
+static VALUE
+bigsqr_fast(VALUE x)
+{
• long len = RBIGNUM_LEN(x), i, j;
• VALUE z = bignew(2 * len + 1, 1);
• BDIGIT *xds = BDIGITS(x), *zds = BDIGITS(z);
• BDIGIT_DBL c, v, w;
• for (i = 2 * len + 1; i–; ) zds[i] = 0;
• for (i = 0; i < len; i++) {
• v = (BDIGIT_DBL)xds[i];
• if (!v) continue;
• c = (BDIGIT_DBL)zds[i + i] + v * v;
• zds[i + i] = BIGLO(c);
• c = BIGDN(c);
• v *= 2;
• for (j = i + 1; j < len; j++) {
• `````` w = (BDIGIT_DBL)xds[j];
``````
• `````` c += (BDIGIT_DBL)zds[i + j] + BIGLO(v) * w;
``````
• `````` zds[i + j] = BIGLO(c);
``````
• `````` c = BIGDN(c);
``````
• `````` if (BIGDN(v)) c += w;
``````
• }
• if (c) {
• `````` c += (BDIGIT_DBL)zds[i + len];
``````
• `````` zds[i + len] = BIGLO(c);
``````
• `````` c = BIGDN(c);
``````
• }
• if (c) zds[i + len + 1] += c;
• }
• return z;
+}

+#define KARATSUBA_MUL_DIGITS 70
+
+
+/* determine whether a bignum is sparse or not by random sampling */
+static inline VALUE
+big_sparse_p(VALUE x)
+{

• long c = 0, n = RBIGNUM_LEN(x);
• unsigned long rb_rand_internal(unsigned long i);
• if ( BDIGITS(x)[rb_rand_internal(n / 2) + n / 4]) c++;
• if (c <= 1 && BDIGITS(x)[rb_rand_internal(n / 2) + n / 4]) c++;
• if (c <= 1 && BDIGITS(x)[rb_rand_internal(n / 2) + n / 4]) c++;
• return (c <= 1) ? Qtrue : Qfalse;
+}

+#if 0
+static void
+dump_bignum(VALUE x)
+{

• long i;
• printf(“0x0”);
• for (i = RBIGNUM_LEN(x); i–; ) {
• `````` printf("_%08x", BDIGITS(x)[i]);
``````
• }
• puts(“”);
+}
+#endif

+static VALUE
+bigmul0(VALUE x, VALUE y)
+{

• long xn, yn;
• xn = RBIGNUM_LEN(x);
• yn = RBIGNUM_LEN(y);
• /* make sure that y is longer than x */
• if (xn > yn) {
• VALUE t;
• long tn;
• t = x; x = y; y = t;
• tn = xn; xn = yn; yn = tn;
• }
• assert(xn <= yn);
• /* normal multiplication when x is small */
• if (xn < KARATSUBA_MUL_DIGITS) {
• `````` normal:
``````
• if (x == y) return bigsqr_fast(x);
• `````` return bigmul1_normal(x, y);
``````
• }
• /* normal multiplication when x or y is a sparse bignum */
• if (big_sparse_p(x)) goto normal;
• if (big_sparse_p(y)) return bigmul1_normal(y, x);
• /* balance multiplication by slicing y when x is much smaller than
y */
• if (2 * xn <= yn) return bigmul1_balance(x, y);
• /* multiplication by karatsuba method */
• return bigmul1_karatsuba(x, y);
+}

/*

• call-seq:
• ``````big * other  => Numeric
``````

@@ -1607,7 +1838,22 @@
VALUE
rb_big_mul(VALUE x, VALUE y)
{

• return bignorm(rb_big_mul0(x, y));
• switch (TYPE(y)) {
• `````` case T_FIXNUM:
``````
• y = rb_int2big(FIX2LONG(y));
• break;
• `````` case T_BIGNUM:
``````
• break;
• `````` case T_FLOAT:
``````
• return DBL2NUM(rb_big2dbl(x) * RFLOAT_VALUE(y));
• `````` default:
``````
• return rb_num_coerce_bin(x, y, ‘*’);
• }
• return bignorm(bigmul0(x, y));
}

struct big_div_struct {
@@ -1661,6 +1907,13 @@
return Qnil;
}

+static void
+rb_big_stop(void *ptr)
+{

• VALUE stop = (VALUE)ptr;
• *stop = Qtrue;
+}

static VALUE
bigdivrem(VALUE x, VALUE y, VALUE *divp, VALUE *modp)
{
@@ -2037,7 +2290,7 @@
BDIGIT_DBL num;

`````` if (len < 4000 / BITSPERDIG) {
``````
• return bigtrunc(rb_big_mul0(x, x));
• return bigtrunc(bigmul0(x, x));
}

a = bignew(len - k, 1);
@@ -2054,7 +2307,7 @@
}
MEMCPY(BDIGITS(z) + 2 * k, BDIGITS(a2), BDIGIT, RBIGNUM_LEN(a2));
RBIGNUM_SET_LEN(z, len);

• a2 = bigtrunc(rb_big_mul0(a, b));
• a2 = bigtrunc(bigmul0(a, b));
len = RBIGNUM_LEN(a2);
for (i = 0, num = 0; i < len; i++) {
num += (BDIGIT_DBL)BDIGITS(z)[i + k] + ((BDIGIT_DBL)BDIGITS(a2)[i] <<
1);
@@ -2125,7 +2378,7 @@
if (z) z = bigtrunc(bigsqr(z));
• ``````   z = z ? bigtrunc(rb_big_mul0(z, x)) : x;
``````
• ``````   z = z ? bigtrunc(bigmul0(z, x)) : x;
``````
}
}
return bignorm(z);
Index: random.c
===================================================================
— random.c (revision 20645)
+++ random.c (working copy)
@@ -452,6 +452,16 @@
return rb_big_norm((VALUE)val);
}

+unsigned long
+rb_rand_internal(unsigned long i)
+{

• struct MT *mt = &default_mt.mt;
• if (!genrand_initialized(mt)) {
• rand_init(mt, random_seed());
• }
• return limited_rand(mt, i);
+}

/*

• call-seq:
• ``````rand(max=0)    => number
``````

e\$B\$^\$D\$b\$He(B e\$B\$f\$-\$R\$m\$G\$9e(B

In message “Re: [ruby-dev:37392] [1.9:Feature] faster bignum
multiplication by karatsuba method”
on Fri, 12 Dec 2008 04:02:42 +0900, “Yusuke ENDOH” [email protected]
writes:

|Python e\$B\$N%=!<%9\$r;29M\$K\$7\$F!"e(Bbignum e\$B\$N>h;;\$re(B Karatsuba e\$BK!\$J\$I\$Ge(B
|e\$B9bB.2=\$7\$F\$_\$^\$7\$?!#e(B

|e\$B0l1~e(B make test && make test-all e\$B\$ODL\$C\$F\$\$\$^\$9e(B (bignum e\$B\$N%F%9%H\$Oe(B
|e\$B=<<B\$7\$F\$J\$\$\$N\$G\$"\$^\$jEv\$F\$K\$J\$i\$J\$\$\$G\$9\$,e(B) e\$B!#e(B
|
|e\$B\$I\$&\$G\$7\$g\$&\$+!#e(B

trunke\$B\$K%3%_%C%H\$7\$F\$/\$@\$5\$\$!#e(B

e\$B:XF#\$H?=\$7\$^\$9!#e(B

On Fri, 12 Dec 2008 04:02:42 +0900
“Yusuke ENDOH” [email protected] wrote:

e\$B1sF#\$G\$9!#e(B

Python e\$B\$N%=!<%9\$r;29M\$K\$7\$F!"e(Bbignum e\$B\$N>h;;\$re(B Karatsuba e\$BK!\$J\$I\$Ge(B
e\$B9bB.2=\$7\$F\$_\$^\$7\$?!#e(B

e\$B\$*\$)!“\$9\$P\$i\$7\$\$!#\$”\$j\$,\$H\$&\$4\$6\$\$\$^\$9!#e(B

e\$B0JA0e(B[ruby-dev:32629]e\$B\$G=q\$\$\$F\$\$\$i\$C\$7\$c\$C\$?e(BFFTe\$BHG\$HHf\$Y\$F\$O\$I\$&\$J\$N\$G\$7\$g\$&!#e(B
e\$BAG?M9M\$(\$G\$O!"N>J}6&7e\$,==J,\$KBg\$-\$\$>l9g\$N\$h\$j87\$7\$\$>r7o2<\$G\$O!"e(BKaratsubae\$B\$+\$ie(B
e\$B@Z\$jBX\$(\$l\$PB.\$/\$J\$k5\$\$,\$7\$^\$9!#e(B

e\$B\$^\$D\$b\$H\$5\$s\$,!V%/%j%9%^%9\$r2a\$.\$?\$i<h\$j9~\$`!W\$H8@\$C\$F\$=\$N\$^\$^\$G\$"\$k\$H;W\$\$\$^\$9!#e(B

e\$B1sF#\$G\$9!#e(B

2008/12/12 13:02 Tadashi S. [email protected]:

e\$B0JA0e(B[ruby-dev:32629]e\$B\$G=q\$\$\$F\$\$\$i\$C\$7\$c\$C\$?e(BFFTe\$BHG\$HHf\$Y\$F\$O\$I\$&\$J\$N\$G\$7\$g\$&!#e(B
e\$BAG?M9M\$(\$G\$O!"N>J}6&7e\$,==J,\$KBg\$-\$\$>l9g\$N\$h\$j87\$7\$\$>r7o2<\$G\$O!"e(BKaratsubae\$B\$+\$ie(B
e\$B@Z\$jBX\$(\$l\$PB.\$/\$J\$k5\$\$,\$7\$^\$9!#e(B

e\$B3P\$(\$F\$/\$l\$F\$\$\$k?M\$,\$\$\$?\$H\$O!#\$"\$j\$,\$H\$&\$4\$6\$\$\$^\$9!#e(B

e\$B<B\$O;d<+?H\$,\$b\$&\$"\$s\$^\$j3P\$(\$F\$J\$\$\$s\$G\$9\$,!“K\Ev\$K7e\$,Bg\$-\$\$>l9g\$Oe(B
FFT e\$B\$Ne(B
e\$BJ}\$,Aa\$\$\$O\$:\$J\$N\$G!”;n\$7\$F\$\$\$^\$;\$s\$,B?J,\$=\$NJ}K!\$OM-8z\$@\$H;W\$\$\$^\$9!#e(B

e\$B\$?\$@\$“\$N%Q%C%A\$O!“e(Bdouble e\$B7?\$,e(B IEEE754
e\$B\$NG@:EY\$G\$”\$k\$3\$He(B (e\$BFC\$K2>?tIt\$,e(B
52bit e\$B0J>e\$”\$k\$3\$He(B)
e\$B\$r2>Dj\$7\$F\$\$\$k\$?\$a!“\$=\$&\$G\$J\$\$JQ\$J4D6-\$GF0\$+\$J\$\$\$N\$Ge(B
e\$B%@%a!”\$H8@\$&%3%a%s%H\$,\$D\$\$\$?5-21\$,\$"\$j\$^\$9e(B (ML e\$B\$G\$J\$/e(B IRC
e\$B\$Ge(B) e\$B!#e(B

e\$B%3%s%Q%\$%k;~\$Ke(B double e\$B7?\$N2>?tIt\$,e(B 52 bit
e\$B0J>e\$+\$I\$&\$+\$rH=Dj\$G\$-\$l\$P<h\$je(B
e\$B9~\$a\$k\$+\$b\$7\$l\$^\$;\$s\$,!"e(B

• e\$B!VJQ\$J4D6-!W\$G\$A\$c\$s\$HF0\$/\$+%F%9%H\$r\$9\$k\$N\$,Fq\$7\$\$e(B
• Karatsuba e\$B>h;;\$G\$b\$^\$"==J,\$+\$J\$H;W\$&e(B
• e\$B:#\$H\$J\$C\$F\$O\$“\$N%Q%C%A\$K\$”\$^\$j<+?.\$,\$J\$\$e(B
(e\$B8m:9I>2A\$r\$I\$&\$d\$C\$?\$+\$H\$+e(B
e\$BA4\$/3P\$(\$F\$J\$\$e(B)

e\$B\$N\$G!“@5D>\$J\$H\$3\$m!”;d<+?H\$O:n6H\$9\$k5\$\$,\$J\$\$\$G\$9!#e(B

e\$B\$J\$+\$@\$G\$9!#e(B

At Sun, 14 Dec 2008 13:01:52 +0900,
Yusuke ENDOH wrote in [ruby-dev:37435]:

e\$B%3%s%Q%\$%k;~\$Ke(B double e\$B7?\$N2>?tIt\$,e(B 52 bit e\$B0J>e\$+\$I\$&\$+\$rH=Dj\$G\$-\$l\$P<h\$je(B
e\$B9~\$a\$k\$+\$b\$7\$l\$^\$;\$s\$,!"e(B

#include <float.h>
#if DBL_MANT_DIG >= 52

#endif

e\$B\$G!#e(B

e\$B1sF#\$G\$9!#e(B

2008/12/12 9:19 Yukihiro M. [email protected]:

|
|e\$B\$I\$&\$G\$7\$g\$&\$+!#e(B

trunke\$B\$K%3%_%C%H\$7\$F\$/\$@\$5\$\$!#e(B

e\$B\$"\$j\$,\$H\$&\$4\$6\$\$\$^\$9!#e(B
e\$B4JC1\$J%F%9%H\$H\$H\$b\$K%3%_%C%H\$7\$^\$7\$?!#e(B

e\$B\$"\$H!"e(Bbigsqr(x) e\$B\$,\$J\$s\$+CY\$\$e(B (bigmul0(x, x)
e\$B\$7\$?J}\$,B.\$\$e(B) e\$B\$3\$H\$K5\$\$E\$\$\$?\$N\$G!"e(B
e\$B\$b\$&\$A\$g\$C\$HD4\$Y\$F\$5\$i\$K=\$@5\$9\$k\$H;W\$\$\$^\$9!#e(B