From: nobu.nakada@... Date: 2003-05-01T18:39:55+09:00 Subject: [ruby-dev:20133] Re: Roundoff problem with Float and Marshal なかだです。 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 #include +#ifdef HAVE_FLOAT_H +#include +#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はできる。 中田 伸悦