[ruby-dev:20133] Re: Roundoff problem with Float and Marshal
From:
nobu.nakada@...
Date:
2003-05-01 09:39:55 UTC
List:
ruby-dev #20133
なかだです。
At Tue, 29 Apr 2003 01:25:08 +0900,
Nobuyoshi-Nakada wrote:
> 強引にBignumで桁数を増やして精度を確保するというのを試してみま
> したが、あまりにも遅すぎなのでもう少し工夫が要りそうです。
若干改善したものの、やっぱりString#to_fは現在の三倍の遅さです。
Index: bignum.c
===================================================================
RCS file: /pub/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 1 May 2003 09:39:33 -0000
@@ -15,4 +15,7 @@
#include <math.h>
#include <ctype.h>
+#ifdef HAVE_FLOAT_H
+#include <float.h>
+#endif
VALUE rb_cBignum;
@@ -94,4 +97,18 @@ rb_big_2comp(x) /* get 2's complement
}
+static long
+bigshrink(x)
+ VALUE x;
+{
+ long len = RBIGNUM(x)->len;
+ BDIGIT *ds = BDIGITS(x);
+
+ if (len > 0) {
+ while (!ds[--len] && len);
+ RBIGNUM(x)->len = ++len;
+ }
+ return len;
+}
+
static VALUE
bignorm(x)
@@ -99,10 +116,7 @@ bignorm(x)
{
if (!FIXNUM_P(x)) {
- long len = RBIGNUM(x)->len;
+ long len = bigshrink(x);
BDIGIT *ds = BDIGITS(x);
- while (len-- && !ds[len]) ;
- RBIGNUM(x)->len = ++len;
-
if (len*SIZEOF_BDIGITS <= sizeof(VALUE)) {
long num = 0;
@@ -828,6 +842,6 @@ rb_dbl2big(d)
}
-double
-rb_big2dbl(x)
+static double
+big2dbl(x)
VALUE x;
{
@@ -844,4 +858,11 @@ rb_big2dbl(x)
}
+double
+rb_big2dbl(x)
+ VALUE x;
+{
+ return rb_int_scale(x, 0);
+}
+
static VALUE
rb_big_to_f(x)
@@ -1078,6 +1099,6 @@ rb_big_minus(x, y)
}
-VALUE
-rb_big_mul(x, y)
+static VALUE
+bigmul(x, y)
VALUE x, y;
{
@@ -1087,20 +1108,4 @@ rb_big_mul(x, y)
BDIGIT *zds;
- if (FIXNUM_P(x)) x = rb_int2big(FIX2LONG(x));
- switch (TYPE(y)) {
- case T_FIXNUM:
- y = rb_int2big(FIX2LONG(y));
- break;
-
- case T_BIGNUM:
- break;
-
- case T_FLOAT:
- return rb_float_new(rb_big2dbl(x) * RFLOAT(y)->value);
-
- default:
- return rb_num_coerce_bin(x, y);
- }
-
j = RBIGNUM(x)->len + RBIGNUM(y)->len + 1;
z = bignew(j, RBIGNUM(x)->sign==RBIGNUM(y)->sign);
@@ -1122,5 +1127,28 @@ rb_big_mul(x, y)
}
- return bignorm(z);
+ return z;
+}
+
+VALUE
+rb_big_mul(x, y)
+ VALUE x, y;
+{
+ if (FIXNUM_P(x)) x = rb_int2big(FIX2LONG(x));
+ switch (TYPE(y)) {
+ case T_FIXNUM:
+ y = rb_int2big(FIX2LONG(y));
+ break;
+
+ case T_BIGNUM:
+ break;
+
+ case T_FLOAT:
+ return rb_float_new(rb_big2dbl(x) * RFLOAT(y)->value);
+
+ default:
+ return rb_num_coerce_bin(x, y);
+ }
+
+ return bignorm(bigmul(x, y));
}
@@ -1390,4 +1418,21 @@ rb_big_quo(x, y)
}
+static VALUE
+bigpow(x, y)
+ VALUE x;
+ unsigned long y;
+{
+ VALUE z = x;
+
+ while (--y > 0) {
+ while (!(y & 1)) {
+ y >>= 1;
+ bigshrink(x = bigmul(x, x));
+ }
+ bigshrink(z = bigmul(z, x));
+ }
+ return z;
+}
+
VALUE
rb_big_pow(x, y)
@@ -1411,16 +1456,5 @@ rb_big_pow(x, y)
yy = NUM2LONG(y);
if (yy > 0) {
- VALUE z = x;
-
- for (;;) {
- yy -= 1;
- if (yy == 0) break;
- while (yy % 2 == 0) {
- yy /= 2;
- x = rb_big_mul(x, x);
- }
- z = rb_big_mul(z, x);
- }
- return bignorm(z);
+ return bignorm(bigpow(x, yy));
}
d = (double)yy;
@@ -1590,12 +1624,10 @@ rb_big_xor(x, y)
}
-static VALUE rb_big_rshift _((VALUE,VALUE));
-
-VALUE
-rb_big_lshift(x, y)
- VALUE x, y;
+static VALUE
+big_lshift(x, shift)
+ VALUE x;
+ int shift;
{
BDIGIT *xds, *zds;
- int shift = NUM2INT(y);
int s1 = shift/BITSPERDIG;
int s2 = shift%BITSPERDIG;
@@ -1604,5 +1636,4 @@ rb_big_lshift(x, y)
long len, i;
- if (shift < 0) return rb_big_rshift(x, INT2FIX(-shift));
len = RBIGNUM(x)->len;
z = bignew(len+s1+1, RBIGNUM(x)->sign);
@@ -1618,13 +1649,13 @@ rb_big_lshift(x, y)
}
*zds = BIGLO(num);
- return bignorm(z);
+ return z;
}
static VALUE
-rb_big_rshift(x, y)
- VALUE x, y;
+big_rshift(x, shift)
+ VALUE x;
+ int shift;
{
BDIGIT *xds, *zds;
- int shift = NUM2INT(y);
long s1 = shift/BITSPERDIG;
long s2 = shift%BITSPERDIG;
@@ -1633,6 +1664,4 @@ rb_big_rshift(x, y)
long i, j;
- if (shift < 0) return rb_big_lshift(x, INT2FIX(-shift));
-
if (s1 > RBIGNUM(x)->len) {
if (RBIGNUM(x)->sign)
@@ -1660,7 +1689,28 @@ rb_big_rshift(x, y)
get2comp(z, Qfalse);
}
+ return z;
+}
+
+VALUE
+rb_big_lshift(x, y)
+ VALUE x, y;
+{
+ int shift = NUM2INT(y);
+ VALUE z = shift < 0 ? big_rshift(x, -shift) : big_lshift(x, shift);
+
+ return bignorm(z);
+}
+
+VALUE
+rb_big_rshift(x, y)
+ VALUE x, y;
+{
+ int shift = NUM2INT(y);
+ VALUE z = shift < 0 ? big_lshift(x, -shift) : big_rshift(x, shift);
+
return bignorm(z);
}
+
static VALUE
rb_big_aref(x, y)
@@ -1759,4 +1809,110 @@ rb_big_size(big)
{
return LONG2FIX(RBIGNUM(big)->len*SIZEOF_BDIGITS);
+}
+
+unsigned int
+rb_int_bitsize(x)
+ VALUE x;
+{
+ unsigned int n, i;
+ unsigned long v, m;
+
+ if (!FIXNUM_P(x)) {
+ long len = RBIGNUM(x)->len;
+ BDIGIT *ds = BDIGITS(x);
+
+ while (len-- && !(v = ds[len])) ;
+ if (len < 0) return 0;
+ n = len * BITSPERDIG;
+ }
+ else {
+ long l = FIX2LONG(x);
+ v = l < 0 ? -l : l;
+ n = 0;
+ }
+ for (m = ~0UL << (i = SIZEOF_LONG * CHAR_BIT / 2); i; m >>= (i >>= 1)) {
+ if (v & m) {
+ n += i;
+ v >>= i;
+ }
+ }
+ if (v) n++;
+ return n;
+}
+
+#define EXTRA_DIG ((DBL_MANT_DIG / BITSPERDIG + 1) * BITSPERDIG)
+
+static VALUE
+exponent5(z, e)
+ VALUE z;
+ int e;
+{
+ static VALUE exp5s = Qnil;
+ long i;
+
+ if (NIL_P(exp5s)) {
+ VALUE x = rb_int2big(5L);
+
+ exp5s = rb_ary_new2(10);
+ rb_ary_push(exp5s, x);
+ for (i = 0; (1<<++i) < DBL_MAX_10_EXP;) {
+ bigshrink(x = rb_obj_freeze(bigmul(x, x)));
+ rb_ary_push(exp5s, x);
+ }
+ rb_obj_freeze(exp5s);
+ rb_gc_register_address(&exp5s);
+ }
+
+ for (i = 0; e > 0; e >>= 1, i++) {
+ if (e & 1) {
+ bigshrink(z = bigmul(z, RARRAY(exp5s)->ptr[i]));
+ }
+ }
+ return z;
+}
+
+double
+rb_int_scale(self, e)
+ VALUE self;
+ int e;
+{
+ int k;
+ VALUE x, y;
+ double z = 0.0;
+
+ if (FIXNUM_P(self)) {
+ if (self == INT2FIX(0)) return 0.0;
+ self = rb_int2big(FIX2LONG(self));
+ }
+ else if (BIGZEROP(self)) {
+ return 0.0;
+ }
+ x = rb_big_abs(self);
+ y = INT2FIX(1);
+ if (e > 0) {
+ if (e > DBL_MAX_10_EXP) {
+ z = 1.0 / z;
+ if (!RBIGNUM(self)->sign) z = -z;
+ return z;
+ }
+ x = exponent5(x, e);
+ }
+ else if (e < 0) {
+ if (e < DBL_MIN_10_EXP) {
+ z = 0.0;
+ if (!RBIGNUM(self)->sign) z = -z;
+ return z;
+ }
+ y = exponent5(rb_int2big(1), -e);
+ }
+ k = EXTRA_DIG - (int)rb_int_bitsize(x);
+ if (!FIXNUM_P(y)) k += (int)rb_int_bitsize(y);
+ if (k > 0) {
+ x = big_lshift(x, k);
+ e -= k;
+ }
+ if (!FIXNUM_P(y)) bigdivrem(x, y, &x, 0);
+ RBIGNUM(x)->sign = RBIGNUM(self)->sign;
+ return ldexp(big2dbl(x), e);
}
--
--- 僕の前にBugはない。
--- 僕の後ろにBugはできる。
中田 伸悦