[ruby-dev:20071] Re: Roundoff problem with Float and Marshal

なかだです。

At Wed, 23 Apr 2003 16:52:03 +0900,
Yukihiro Matsumoto wrote:
> |> 確かにglibcのstrtodを使うとrubyが使っているものよりも正確な
> |> 値を出します。これは内部的にMPを使って精度を稼いでいるようで
> |> す。速度は重要ではないので、Bignumを使えば良いってことなのか
> |> なあ。
> |
> |Bignum#to_fも誤差を含んでるようです。
> |$ ruby -e 'p 39560152763386141.to_f'
> |3.9560152763386144e+16
> 
> 単純にかけ算すると誤差が累積するんでしょうね。これってどうや
> ると誤差が減らせるんですかねえ。

これは加算の丸めが原因のようです。

> |10進小数を2進で正確に表現できないのはしょうがないので、どう丸め
> |るかという話になるわけですが、それがexponentによって変わってし
> |まうというのが問題ではないかと。
> 
> いや、まあ、その通りなんですが、glibcのstrtodよりも明らかに
> 負けているのはなんか悔しいかも。もっとも向こうは1600行もある
> 巨大な関数なんですが。

こんなところで[ruby-talk:69808]は一応通るようになりました。
Celeron Mendocino, gcc 3.2.1 20021207, glibc 2.3.1 math関連のコ
ンパイルオプションはなしです。


Index: bignum.c
===================================================================
RCS file: /cvs/ruby/src/ruby/bignum.c,v
retrieving revision 1.91
diff -u -2 -p -r1.91 bignum.c
--- bignum.c	21 Apr 2003 09:36:31 -0000	1.91
+++ bignum.c	23 Apr 2003 09:42:13 -0000
@@ -15,4 +15,7 @@
 #include <math.h>
 #include <ctype.h>
+#ifdef HAVE_FLOAT_H
+#include <float.h>
+#endif
 
 VALUE rb_cBignum;
@@ -835,7 +838,29 @@ rb_big2dbl(x)
     long i = RBIGNUM(x)->len;
     BDIGIT *ds = BDIGITS(x);
+#ifdef DBL_MANT_DIG
+    int n = DBL_MANT_DIG - BITSPERDIG;
+#endif
 
+#ifdef DBL_MANT_DIG
+    if (i--) {
+	BDIGIT b = ds[i];
+	d = (double)b;
+	while (!(b & (1<<(BITSPERDIG-1)))) {
+	    ++n;
+	    b <<= 1;
+	}
+    }
+#endif
     while (i--) {
-	d = ds[i] + BIGRAD*d;
+	BDIGIT b = ds[i];
+#ifdef DBL_MANT_DIG
+	if (n < BITSPERDIG) {
+	    b &= ~0 << (BITSPERDIG - n);
+	}
+#endif
+	d = b + BIGRAD * d;
+#ifdef DBL_MANT_DIG
+	if ((n -= BITSPERDIG) <= 0) break;
+#endif
     }
     if (isinf(d)) d = HUGE_VAL;
Index: marshal.c
===================================================================
RCS file: /cvs/ruby/src/ruby/marshal.c,v
retrieving revision 1.88
diff -u -2 -p -r1.88 marshal.c
--- marshal.c	22 Apr 2003 10:08:57 -0000	1.88
+++ marshal.c	23 Apr 2003 09:45:30 -0000
@@ -185,5 +185,5 @@ w_long(x, arg)
 }
 
-#ifdef DBL_MANT_DIG
+#if 0 /* defined DBL_MANT_DIG */
 #define DECIMAL_MANT (53-16)	/* from IEEE754 double precision */
 
@@ -292,11 +292,6 @@ w_float(d, arg)
     }
     else {
-	int len;
-
 	/* xxx: should not use system's sprintf(3) */
 	sprintf(buf, "%.*g", FLOAT_DIG, d);
-	len = strlen(buf);
-	w_bytes(buf, len + save_mantissa(d, buf + len), arg);
-	return;
     }
     w_bytes(buf, strlen(buf), arg);
@@ -1027,7 +1022,5 @@ r_object0(arg, proc)
 	    }
 	    else {
-		char *e;
-		d = strtod(ptr, &e);
-		d = load_mantissa(d, e, RSTRING(str)->len - (e - ptr));
+		d = strtod(ptr, 0);
 	    }
 	    v = rb_float_new(d);
Index: util.c
===================================================================
RCS file: /cvs/ruby/src/ruby/util.c,v
retrieving revision 1.33
diff -u -2 -p -r1.33 util.c
--- util.c	17 Apr 2003 16:49:23 -0000	1.33
+++ util.c	23 Apr 2003 09:47:48 -0000
@@ -675,5 +675,6 @@ static int maxExponent = 511;	/* Largest
 				 * no need to worry about additional digits.
 				 */
-static double powersOf10[] = {	/* Table giving binary powers of 10.  Entry */
+static const double
+powersOf10[] = {		/* Table giving binary powers of 10.  Entry */
     10.0,			/* is 10^2^i.  Used to convert decimal */
     100.0,			/* exponents into floating-point numbers. */
@@ -686,4 +687,16 @@ static double powersOf10[] = {	/* Table 
     1.0e256
 };
+static const double
+npowersOf10[] = {
+    0.1,
+    0.01,
+    1.0e-4,
+    1.0e-8,
+    1.0e-16,
+    1.0e-32,
+    1.0e-64,
+    1.0e-128,
+    1.0e-256
+};
 
 /*
@@ -726,5 +739,6 @@ ruby_strtod(string, endPtr)
 {
     int sign, expSign = FALSE;
-    double fraction, dblExp, *d;
+    double fraction;
+    const double *d;
     register const char *p;
     register int c;
@@ -809,5 +823,5 @@ ruby_strtod(string, endPtr)
     }
     else {
-	int frac1, frac2;
+	unsigned frac1, frac2;
 	frac1 = 0;
 	for ( ; mantSize > 9; mantSize -= 1) {
@@ -875,44 +889,16 @@ ruby_strtod(string, endPtr)
 	    errno = ERANGE;
 	}
-	fracExp = exp;
-	exp += 9;
-	if (exp < 0) {
-	    expSign = TRUE;
-	    exp = -exp;
-	}
-	else {
-	    expSign = FALSE;
-	}
-	dblExp = 1.0;
-	for (d = powersOf10; exp != 0; exp >>= 1, d += 1) {
-	    if (exp & 01) {
-		dblExp *= *d;
-	    }
-	}
-	if (expSign) {
-	    fraction = frac1 / dblExp;
-	}
-	else {
-	    fraction = frac1 * dblExp;
-	}
-	exp = fracExp;
 	if (exp < 0) {
-	    expSign = TRUE;
 	    exp = -exp;
+	    d = npowersOf10;
 	}
 	else {
-	    expSign = FALSE;
+	    d = powersOf10;
 	}
-	dblExp = 1.0;
-	for (d = powersOf10; exp != 0; exp >>= 1, d += 1) {
+	fraction = frac1 * 1e9 + frac2;
+	for (; exp != 0; exp >>= 1, d += 1) {
 	    if (exp & 01) {
-		dblExp *= *d;
+		fraction *= *d;
 	    }
-	}
-	if (expSign) {
-	    fraction += frac2 / dblExp;
-	}
-	else {
-	    fraction += frac2 * dblExp;
 	}
     }


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