From: Nobuyoshi Nakada <nobu@...>
Date: 2007-05-02T14:45:28+09:00
Subject: [ruby-dev:30753] Re: 大きな数の大まかな割り算

なかだです。

At Wed, 2 May 2007 08:56:56 +0900,
Yukihiro Matsumoto wrote in [ruby-dev:30743]:
> [ruby-core:11069]で報告があったのですが、rationalの分子分母
> があまりに大きいとFloatに変換できず、本来有限の値なのに結果
> がNaNになってしまうようです。

Bignum#quoで対応するというのも可能かも。

$ ./ruby -e 'num = 18511**365; den = 18250**365; p num.quo(den)'
178.221209978501


Index: bignum.c
===================================================================
--- bignum.c	(revision 12240)
+++ bignum.c	(working copy)
@@ -14,4 +14,5 @@
 
 #include <math.h>
+#include <float.h>
 #include <ctype.h>
 #ifdef HAVE_IEEEFP_H
@@ -853,6 +854,6 @@ rb_dbl2big(double d)
 }
 
-double
-rb_big2dbl(VALUE x)
+static double
+big2dbl(VALUE x)
 {
     double d = 0.0;
@@ -863,9 +864,17 @@ rb_big2dbl(VALUE x)
 	d = ds[i] + BIGRAD*d;
     }
+    if (!RBIGNUM(x)->sign) d = -d;
+    return d;
+}
+
+double
+rb_big2dbl(VALUE x)
+{
+    double d = big2dbl(x);
+
     if (isinf(d)) {
 	rb_warn("Bignum out of Float range");
 	d = HUGE_VAL;
     }
-    if (!RBIGNUM(x)->sign) d = -d;
     return d;
 }
@@ -1497,4 +1506,27 @@ rb_big_divmod(VALUE x, VALUE y)
 }
 
+static int
+bdigbitsize(BDIGIT x)
+{
+    int size = 1;
+    int nb = BITSPERDIG / 2;
+    BDIGIT bits = (~0 << nb);
+
+    if (!x) return 0;
+    while (x > 1) {
+	if (x & bits) {
+	    size += nb;
+	    x >>= nb;
+	}
+	x &= ~bits;
+	nb /= 2;
+	bits >>= nb;
+    }
+
+    return size;
+}
+
+static VALUE rb_big_rshift(VALUE,VALUE);
+
 /*
  *  call-seq:
@@ -1512,7 +1544,35 @@ static VALUE
 rb_big_quo(VALUE x, VALUE y)
 {
-    double dx = rb_big2dbl(x);
+    double dx = big2dbl(x);
     double dy;
 
+    if (isinf(dx)) {
+#define DBL_BIGDIG ((DBL_MANT_DIG + BITSPERDIG) / BITSPERDIG)
+	VALUE z;
+	int ex, ey;
+
+	ex = (RBIGNUM(bigtrunc(x))->len - 1) * BITSPERDIG;
+	ex += bdigbitsize(BDIGITS(x)[RBIGNUM(x)->len - 1]);
+	ex -= 2 * DBL_BIGDIG * BITSPERDIG;
+	if (ex) x = rb_big_rshift(x, INT2FIX(ex));
+
+	switch (TYPE(y)) {
+	  case T_FIXNUM:
+	    y = rb_int2big(FIX2LONG(y));
+	  case T_BIGNUM: {
+	    ey = (RBIGNUM(bigtrunc(y))->len - 1) * BITSPERDIG;
+	    ey += bdigbitsize(BDIGITS(y)[RBIGNUM(y)->len - 1]);
+	    ey -= DBL_BIGDIG * BITSPERDIG;
+	    if (ey) y = rb_big_rshift(y, INT2FIX(ey));
+	  bignum:
+	    bigdivrem(x, y, &z, 0);
+	    return rb_float_new(ldexp(big2dbl(z), ex - ey));
+	  }
+	  case T_FLOAT:
+	    y = dbl2big(ldexp(frexp(RFLOAT(y)->value, &ey), DBL_MANT_DIG));
+	    ey -= DBL_MANT_DIG;
+	    goto bignum;
+	}
+    }
     switch (TYPE(y)) {
       case T_FIXNUM:
@@ -1819,6 +1879,4 @@ rb_big_xor(VALUE xx, VALUE yy)
 }
 
-static VALUE rb_big_rshift(VALUE,VALUE);
-
 /*
  * call-seq:


-- 
--- 僕の前にBugはない。
--- 僕の後ろにBugはできる。
    中田 伸悦