Ruby  2.0.0p247(2013-06-27revision41674)
bignum.c
Go to the documentation of this file.
00001 /**********************************************************************
00002 
00003   bignum.c -
00004 
00005   $Author: nagachika $
00006   created at: Fri Jun 10 00:48:55 JST 1994
00007 
00008   Copyright (C) 1993-2007 Yukihiro Matsumoto
00009 
00010 **********************************************************************/
00011 
00012 #include "ruby/ruby.h"
00013 #include "ruby/thread.h"
00014 #include "ruby/util.h"
00015 #include "internal.h"
00016 
00017 #ifdef HAVE_STRINGS_H
00018 #include <strings.h>
00019 #endif
00020 #include <math.h>
00021 #include <float.h>
00022 #include <ctype.h>
00023 #ifdef HAVE_IEEEFP_H
00024 #include <ieeefp.h>
00025 #endif
00026 #include <assert.h>
00027 
00028 VALUE rb_cBignum;
00029 
00030 static VALUE big_three = Qnil;
00031 
00032 #if defined __MINGW32__
00033 #define USHORT _USHORT
00034 #endif
00035 
00036 #define BDIGITS(x) (RBIGNUM_DIGITS(x))
00037 #define BITSPERDIG (SIZEOF_BDIGITS*CHAR_BIT)
00038 #define BIGRAD ((BDIGIT_DBL)1 << BITSPERDIG)
00039 #define BIGRAD_HALF ((BDIGIT)(BIGRAD >> 1))
00040 #define DIGSPERLONG (SIZEOF_LONG/SIZEOF_BDIGITS)
00041 #if HAVE_LONG_LONG
00042 # define DIGSPERLL (SIZEOF_LONG_LONG/SIZEOF_BDIGITS)
00043 #endif
00044 #define BIGUP(x) ((BDIGIT_DBL)(x) << BITSPERDIG)
00045 #define BIGDN(x) RSHIFT((x),BITSPERDIG)
00046 #define BIGLO(x) ((BDIGIT)((x) & (BIGRAD-1)))
00047 #define BDIGMAX ((BDIGIT)-1)
00048 
00049 #define BIGZEROP(x) (RBIGNUM_LEN(x) == 0 || \
00050                      (BDIGITS(x)[0] == 0 && \
00051                       (RBIGNUM_LEN(x) == 1 || bigzero_p(x))))
00052 
00053 #define BIGNUM_DEBUG 0
00054 #if BIGNUM_DEBUG
00055 #define ON_DEBUG(x) do { x; } while (0)
00056 static void
00057 dump_bignum(VALUE x)
00058 {
00059     long i;
00060     printf("%c0x0", RBIGNUM_SIGN(x) ? '+' : '-');
00061     for (i = RBIGNUM_LEN(x); i--; ) {
00062         printf("_%08"PRIxBDIGIT, BDIGITS(x)[i]);
00063     }
00064     printf(", len=%lu", RBIGNUM_LEN(x));
00065     puts("");
00066 }
00067 
00068 static VALUE
00069 rb_big_dump(VALUE x)
00070 {
00071     dump_bignum(x);
00072     return x;
00073 }
00074 #else
00075 #define ON_DEBUG(x)
00076 #endif
00077 
00078 static int
00079 bigzero_p(VALUE x)
00080 {
00081     long i;
00082     BDIGIT *ds = BDIGITS(x);
00083 
00084     for (i = RBIGNUM_LEN(x) - 1; 0 <= i; i--) {
00085         if (ds[i]) return 0;
00086     }
00087     return 1;
00088 }
00089 
00090 int
00091 rb_bigzero_p(VALUE x)
00092 {
00093     return BIGZEROP(x);
00094 }
00095 
00096 int
00097 rb_cmpint(VALUE val, VALUE a, VALUE b)
00098 {
00099     if (NIL_P(val)) {
00100         rb_cmperr(a, b);
00101     }
00102     if (FIXNUM_P(val)) {
00103         long l = FIX2LONG(val);
00104         if (l > 0) return 1;
00105         if (l < 0) return -1;
00106         return 0;
00107     }
00108     if (RB_TYPE_P(val, T_BIGNUM)) {
00109         if (BIGZEROP(val)) return 0;
00110         if (RBIGNUM_SIGN(val)) return 1;
00111         return -1;
00112     }
00113     if (RTEST(rb_funcall(val, '>', 1, INT2FIX(0)))) return 1;
00114     if (RTEST(rb_funcall(val, '<', 1, INT2FIX(0)))) return -1;
00115     return 0;
00116 }
00117 
00118 #define RBIGNUM_SET_LEN(b,l) \
00119     ((RBASIC(b)->flags & RBIGNUM_EMBED_FLAG) ? \
00120      (void)(RBASIC(b)->flags = \
00121             (RBASIC(b)->flags & ~RBIGNUM_EMBED_LEN_MASK) | \
00122             ((l) << RBIGNUM_EMBED_LEN_SHIFT)) : \
00123      (void)(RBIGNUM(b)->as.heap.len = (l)))
00124 
00125 static void
00126 rb_big_realloc(VALUE big, long len)
00127 {
00128     BDIGIT *ds;
00129     if (RBASIC(big)->flags & RBIGNUM_EMBED_FLAG) {
00130         if (RBIGNUM_EMBED_LEN_MAX < len) {
00131             ds = ALLOC_N(BDIGIT, len);
00132             MEMCPY(ds, RBIGNUM(big)->as.ary, BDIGIT, RBIGNUM_EMBED_LEN_MAX);
00133             RBIGNUM(big)->as.heap.len = RBIGNUM_LEN(big);
00134             RBIGNUM(big)->as.heap.digits = ds;
00135             RBASIC(big)->flags &= ~RBIGNUM_EMBED_FLAG;
00136         }
00137     }
00138     else {
00139         if (len <= RBIGNUM_EMBED_LEN_MAX) {
00140             ds = RBIGNUM(big)->as.heap.digits;
00141             RBASIC(big)->flags |= RBIGNUM_EMBED_FLAG;
00142             RBIGNUM_SET_LEN(big, len);
00143             if (ds) {
00144                 MEMCPY(RBIGNUM(big)->as.ary, ds, BDIGIT, len);
00145                 xfree(ds);
00146             }
00147         }
00148         else {
00149             if (RBIGNUM_LEN(big) == 0) {
00150                 RBIGNUM(big)->as.heap.digits = ALLOC_N(BDIGIT, len);
00151             }
00152             else {
00153                 REALLOC_N(RBIGNUM(big)->as.heap.digits, BDIGIT, len);
00154             }
00155         }
00156     }
00157 }
00158 
00159 void
00160 rb_big_resize(VALUE big, long len)
00161 {
00162     rb_big_realloc(big, len);
00163     RBIGNUM_SET_LEN(big, len);
00164 }
00165 
00166 static VALUE
00167 bignew_1(VALUE klass, long len, int sign)
00168 {
00169     NEWOBJ_OF(big, struct RBignum, klass, T_BIGNUM);
00170     RBIGNUM_SET_SIGN(big, sign?1:0);
00171     if (len <= RBIGNUM_EMBED_LEN_MAX) {
00172         RBASIC(big)->flags |= RBIGNUM_EMBED_FLAG;
00173         RBIGNUM_SET_LEN(big, len);
00174     }
00175     else {
00176         RBIGNUM(big)->as.heap.digits = ALLOC_N(BDIGIT, len);
00177         RBIGNUM(big)->as.heap.len = len;
00178     }
00179     OBJ_FREEZE(big);
00180     return (VALUE)big;
00181 }
00182 
00183 #define bignew(len,sign) bignew_1(rb_cBignum,(len),(sign))
00184 
00185 VALUE
00186 rb_big_new(long len, int sign)
00187 {
00188     return bignew(len, sign != 0);
00189 }
00190 
00191 VALUE
00192 rb_big_clone(VALUE x)
00193 {
00194     long len = RBIGNUM_LEN(x);
00195     VALUE z = bignew_1(CLASS_OF(x), len, RBIGNUM_SIGN(x));
00196 
00197     MEMCPY(BDIGITS(z), BDIGITS(x), BDIGIT, len);
00198     return z;
00199 }
00200 
00201 /* modify a bignum by 2's complement */
00202 static void
00203 get2comp(VALUE x)
00204 {
00205     long i = RBIGNUM_LEN(x);
00206     BDIGIT *ds = BDIGITS(x);
00207     BDIGIT_DBL num;
00208 
00209     if (!i) return;
00210     while (i--) ds[i] = ~ds[i];
00211     i = 0; num = 1;
00212     do {
00213         num += ds[i];
00214         ds[i++] = BIGLO(num);
00215         num = BIGDN(num);
00216     } while (i < RBIGNUM_LEN(x));
00217     if (num != 0) {
00218         rb_big_resize(x, RBIGNUM_LEN(x)+1);
00219         ds = BDIGITS(x);
00220         ds[RBIGNUM_LEN(x)-1] = 1;
00221     }
00222 }
00223 
00224 void
00225 rb_big_2comp(VALUE x)                   /* get 2's complement */
00226 {
00227     get2comp(x);
00228 }
00229 
00230 static inline VALUE
00231 bigtrunc(VALUE x)
00232 {
00233     long len = RBIGNUM_LEN(x);
00234     BDIGIT *ds = BDIGITS(x);
00235 
00236     if (len == 0) return x;
00237     while (--len && !ds[len]);
00238     if (RBIGNUM_LEN(x) > len+1) {
00239         rb_big_resize(x, len+1);
00240     }
00241     return x;
00242 }
00243 
00244 static inline VALUE
00245 bigfixize(VALUE x)
00246 {
00247     long len = RBIGNUM_LEN(x);
00248     BDIGIT *ds = BDIGITS(x);
00249 
00250     if (len == 0) return INT2FIX(0);
00251     if ((size_t)(len*SIZEOF_BDIGITS) <= sizeof(long)) {
00252         long num = 0;
00253 #if 2*SIZEOF_BDIGITS > SIZEOF_LONG
00254         num = (long)ds[0];
00255 #else
00256         while (len--) {
00257             num = (long)(BIGUP(num) + ds[len]);
00258         }
00259 #endif
00260         if (num >= 0) {
00261             if (RBIGNUM_SIGN(x)) {
00262                 if (POSFIXABLE(num)) return LONG2FIX(num);
00263             }
00264             else {
00265                 if (NEGFIXABLE(-num)) return LONG2FIX(-num);
00266             }
00267         }
00268     }
00269     return x;
00270 }
00271 
00272 static VALUE
00273 bignorm(VALUE x)
00274 {
00275     if (RB_TYPE_P(x, T_BIGNUM)) {
00276         x = bigfixize(bigtrunc(x));
00277     }
00278     return x;
00279 }
00280 
00281 VALUE
00282 rb_big_norm(VALUE x)
00283 {
00284     return bignorm(x);
00285 }
00286 
00287 VALUE
00288 rb_uint2big(VALUE n)
00289 {
00290     BDIGIT_DBL num = n;
00291     long i = 0;
00292     BDIGIT *digits;
00293     VALUE big;
00294 
00295     big = bignew(DIGSPERLONG, 1);
00296     digits = BDIGITS(big);
00297     while (i < DIGSPERLONG) {
00298         digits[i++] = BIGLO(num);
00299         num = BIGDN(num);
00300     }
00301 
00302     i = DIGSPERLONG;
00303     while (--i && !digits[i]) ;
00304     RBIGNUM_SET_LEN(big, i+1);
00305     return big;
00306 }
00307 
00308 VALUE
00309 rb_int2big(SIGNED_VALUE n)
00310 {
00311     long neg = 0;
00312     VALUE u;
00313     VALUE big;
00314 
00315     if (n < 0) {
00316         u = 1 + (VALUE)(-(n + 1)); /* u = -n avoiding overflow */
00317         neg = 1;
00318     }
00319     else {
00320         u = n;
00321     }
00322     big = rb_uint2big(u);
00323     if (neg) {
00324         RBIGNUM_SET_SIGN(big, 0);
00325     }
00326     return big;
00327 }
00328 
00329 VALUE
00330 rb_uint2inum(VALUE n)
00331 {
00332     if (POSFIXABLE(n)) return LONG2FIX(n);
00333     return rb_uint2big(n);
00334 }
00335 
00336 VALUE
00337 rb_int2inum(SIGNED_VALUE n)
00338 {
00339     if (FIXABLE(n)) return LONG2FIX(n);
00340     return rb_int2big(n);
00341 }
00342 
00343 #if SIZEOF_LONG % SIZEOF_BDIGITS != 0
00344 # error unexpected SIZEOF_LONG : SIZEOF_BDIGITS ratio
00345 #endif
00346 
00347 /*
00348  * buf is an array of long integers.
00349  * buf is ordered from least significant word to most significant word.
00350  * buf[0] is the least significant word and
00351  * buf[num_longs-1] is the most significant word.
00352  * This means words in buf is little endian.
00353  * However each word in buf is native endian.
00354  * (buf[i]&1) is the least significant bit and
00355  * (buf[i]&(1<<(SIZEOF_LONG*CHAR_BIT-1))) is the most significant bit
00356  * for each 0 <= i < num_longs.
00357  * So buf is little endian at whole on a little endian machine.
00358  * But buf is mixed endian on a big endian machine.
00359  *
00360  * The buf represents negative integers as two's complement.
00361  * So, the most significant bit of the most significant word,
00362  * (buf[num_longs-1]>>(SIZEOF_LONG*CHAR_BIT-1)),
00363  * is the sign bit: 1 means negative and 0 means zero or positive.
00364  *
00365  * If given size of buf (num_longs) is not enough to represent val,
00366  * higier words (including a sign bit) are ignored.
00367  */
00368 void
00369 rb_big_pack(VALUE val, unsigned long *buf, long num_longs)
00370 {
00371     val = rb_to_int(val);
00372     if (num_longs == 0)
00373         return;
00374     if (FIXNUM_P(val)) {
00375         long i;
00376         long tmp = FIX2LONG(val);
00377         buf[0] = (unsigned long)tmp;
00378         tmp = tmp < 0 ? ~0L : 0;
00379         for (i = 1; i < num_longs; i++)
00380             buf[i] = (unsigned long)tmp;
00381         return;
00382     }
00383     else {
00384         long len = RBIGNUM_LEN(val);
00385         BDIGIT *ds = BDIGITS(val), *dend = ds + len;
00386         long i, j;
00387         for (i = 0; i < num_longs && ds < dend; i++) {
00388             unsigned long l = 0;
00389             for (j = 0; j < DIGSPERLONG && ds < dend; j++, ds++) {
00390                 l |= ((unsigned long)*ds << (j * BITSPERDIG));
00391             }
00392             buf[i] = l;
00393         }
00394         for (; i < num_longs; i++)
00395             buf[i] = 0;
00396         if (RBIGNUM_NEGATIVE_P(val)) {
00397             for (i = 0; i < num_longs; i++) {
00398                 buf[i] = ~buf[i];
00399             }
00400             for (i = 0; i < num_longs; i++) {
00401                 buf[i]++;
00402                 if (buf[i] != 0)
00403                     return;
00404             }
00405         }
00406     }
00407 }
00408 
00409 /* See rb_big_pack comment for endianness and sign of buf. */
00410 VALUE
00411 rb_big_unpack(unsigned long *buf, long num_longs)
00412 {
00413     while (2 <= num_longs) {
00414         if (buf[num_longs-1] == 0 && (long)buf[num_longs-2] >= 0)
00415             num_longs--;
00416         else if (buf[num_longs-1] == ~0UL && (long)buf[num_longs-2] < 0)
00417             num_longs--;
00418         else
00419             break;
00420     }
00421     if (num_longs == 0)
00422         return INT2FIX(0);
00423     else if (num_longs == 1)
00424         return LONG2NUM((long)buf[0]);
00425     else {
00426         VALUE big;
00427         BDIGIT *ds;
00428         long len = num_longs * DIGSPERLONG;
00429         long i;
00430         big = bignew(len, 1);
00431         ds = BDIGITS(big);
00432         for (i = 0; i < num_longs; i++) {
00433             unsigned long d = buf[i];
00434 #if SIZEOF_LONG == SIZEOF_BDIGITS
00435             *ds++ = d;
00436 #else
00437             int j;
00438             for (j = 0; j < DIGSPERLONG; j++) {
00439                 *ds++ = BIGLO(d);
00440                 d = BIGDN(d);
00441             }
00442 #endif
00443         }
00444         if ((long)buf[num_longs-1] < 0) {
00445             get2comp(big);
00446             RBIGNUM_SET_SIGN(big, 0);
00447         }
00448         return bignorm(big);
00449     }
00450 }
00451 
00452 #define QUAD_SIZE 8
00453 
00454 #if SIZEOF_LONG_LONG == QUAD_SIZE && SIZEOF_BDIGITS*2 == SIZEOF_LONG_LONG
00455 
00456 void
00457 rb_quad_pack(char *buf, VALUE val)
00458 {
00459     LONG_LONG q;
00460 
00461     val = rb_to_int(val);
00462     if (FIXNUM_P(val)) {
00463         q = FIX2LONG(val);
00464     }
00465     else {
00466         long len = RBIGNUM_LEN(val);
00467         BDIGIT *ds;
00468 
00469         if (len > SIZEOF_LONG_LONG/SIZEOF_BDIGITS) {
00470             len = SIZEOF_LONG_LONG/SIZEOF_BDIGITS;
00471         }
00472         ds = BDIGITS(val);
00473         q = 0;
00474         while (len--) {
00475             q = BIGUP(q);
00476             q += ds[len];
00477         }
00478         if (!RBIGNUM_SIGN(val)) q = -q;
00479     }
00480     memcpy(buf, (char*)&q, SIZEOF_LONG_LONG);
00481 }
00482 
00483 VALUE
00484 rb_quad_unpack(const char *buf, int sign)
00485 {
00486     unsigned LONG_LONG q;
00487     long neg = 0;
00488     long i;
00489     BDIGIT *digits;
00490     VALUE big;
00491 
00492     memcpy(&q, buf, SIZEOF_LONG_LONG);
00493     if (sign) {
00494         if (FIXABLE((LONG_LONG)q)) return LONG2FIX((LONG_LONG)q);
00495         if ((LONG_LONG)q < 0) {
00496             q = -(LONG_LONG)q;
00497             neg = 1;
00498         }
00499     }
00500     else {
00501         if (POSFIXABLE(q)) return LONG2FIX(q);
00502     }
00503 
00504     i = 0;
00505     big = bignew(DIGSPERLL, 1);
00506     digits = BDIGITS(big);
00507     while (i < DIGSPERLL) {
00508         digits[i++] = BIGLO(q);
00509         q = BIGDN(q);
00510     }
00511 
00512     i = DIGSPERLL;
00513     while (i-- && !digits[i]) ;
00514     RBIGNUM_SET_LEN(big, i+1);
00515 
00516     if (neg) {
00517         RBIGNUM_SET_SIGN(big, 0);
00518     }
00519     return bignorm(big);
00520 }
00521 
00522 #else
00523 
00524 static int
00525 quad_buf_complement(char *buf, size_t len)
00526 {
00527     size_t i;
00528     for (i = 0; i < len; i++)
00529         buf[i] = ~buf[i];
00530     for (i = 0; i < len; i++) {
00531         buf[i]++;
00532         if (buf[i] != 0)
00533             return 0;
00534     }
00535     return 1;
00536 }
00537 
00538 void
00539 rb_quad_pack(char *buf, VALUE val)
00540 {
00541     long len;
00542 
00543     memset(buf, 0, QUAD_SIZE);
00544     val = rb_to_int(val);
00545     if (FIXNUM_P(val)) {
00546         val = rb_int2big(FIX2LONG(val));
00547     }
00548     len = RBIGNUM_LEN(val) * SIZEOF_BDIGITS;
00549     if (len > QUAD_SIZE) {
00550         len = QUAD_SIZE;
00551     }
00552     memcpy(buf, (char*)BDIGITS(val), len);
00553     if (RBIGNUM_NEGATIVE_P(val)) {
00554         quad_buf_complement(buf, QUAD_SIZE);
00555     }
00556 }
00557 
00558 #define BNEG(b) (RSHIFT(((BDIGIT*)(b))[QUAD_SIZE/SIZEOF_BDIGITS-1],BITSPERDIG-1) != 0)
00559 
00560 VALUE
00561 rb_quad_unpack(const char *buf, int sign)
00562 {
00563     VALUE big = bignew(QUAD_SIZE/SIZEOF_BDIGITS, 1);
00564 
00565     memcpy((char*)BDIGITS(big), buf, QUAD_SIZE);
00566     if (sign && BNEG(buf)) {
00567         char *tmp = (char*)BDIGITS(big);
00568 
00569         RBIGNUM_SET_SIGN(big, 0);
00570         quad_buf_complement(tmp, QUAD_SIZE);
00571     }
00572 
00573     return bignorm(big);
00574 }
00575 
00576 #endif
00577 
00578 VALUE
00579 rb_cstr_to_inum(const char *str, int base, int badcheck)
00580 {
00581     const char *s = str;
00582     char *end;
00583     char sign = 1, nondigit = 0;
00584     int c;
00585     BDIGIT_DBL num;
00586     long len, blen = 1;
00587     long i;
00588     VALUE z;
00589     BDIGIT *zds;
00590 
00591 #undef ISDIGIT
00592 #define ISDIGIT(c) ('0' <= (c) && (c) <= '9')
00593 #define conv_digit(c) \
00594     (!ISASCII(c) ? -1 : \
00595      ISDIGIT(c) ? ((c) - '0') : \
00596      ISLOWER(c) ? ((c) - 'a' + 10) : \
00597      ISUPPER(c) ? ((c) - 'A' + 10) : \
00598      -1)
00599 
00600     if (!str) {
00601         if (badcheck) goto bad;
00602         return INT2FIX(0);
00603     }
00604     while (ISSPACE(*str)) str++;
00605 
00606     if (str[0] == '+') {
00607         str++;
00608     }
00609     else if (str[0] == '-') {
00610         str++;
00611         sign = 0;
00612     }
00613     if (str[0] == '+' || str[0] == '-') {
00614         if (badcheck) goto bad;
00615         return INT2FIX(0);
00616     }
00617     if (base <= 0) {
00618         if (str[0] == '0') {
00619             switch (str[1]) {
00620               case 'x': case 'X':
00621                 base = 16;
00622                 break;
00623               case 'b': case 'B':
00624                 base = 2;
00625                 break;
00626               case 'o': case 'O':
00627                 base = 8;
00628                 break;
00629               case 'd': case 'D':
00630                 base = 10;
00631                 break;
00632               default:
00633                 base = 8;
00634             }
00635         }
00636         else if (base < -1) {
00637             base = -base;
00638         }
00639         else {
00640             base = 10;
00641         }
00642     }
00643     switch (base) {
00644       case 2:
00645         len = 1;
00646         if (str[0] == '0' && (str[1] == 'b'||str[1] == 'B')) {
00647             str += 2;
00648         }
00649         break;
00650       case 3:
00651         len = 2;
00652         break;
00653       case 8:
00654         if (str[0] == '0' && (str[1] == 'o'||str[1] == 'O')) {
00655             str += 2;
00656         }
00657       case 4: case 5: case 6: case 7:
00658         len = 3;
00659         break;
00660       case 10:
00661         if (str[0] == '0' && (str[1] == 'd'||str[1] == 'D')) {
00662             str += 2;
00663         }
00664       case 9: case 11: case 12: case 13: case 14: case 15:
00665         len = 4;
00666         break;
00667       case 16:
00668         len = 4;
00669         if (str[0] == '0' && (str[1] == 'x'||str[1] == 'X')) {
00670             str += 2;
00671         }
00672         break;
00673       default:
00674         if (base < 2 || 36 < base) {
00675             rb_raise(rb_eArgError, "invalid radix %d", base);
00676         }
00677         if (base <= 32) {
00678             len = 5;
00679         }
00680         else {
00681             len = 6;
00682         }
00683         break;
00684     }
00685     if (*str == '0') {          /* squeeze preceding 0s */
00686         int us = 0;
00687         while ((c = *++str) == '0' || c == '_') {
00688             if (c == '_') {
00689                 if (++us >= 2)
00690                     break;
00691             } else
00692                 us = 0;
00693         }
00694         if (!(c = *str) || ISSPACE(c)) --str;
00695     }
00696     c = *str;
00697     c = conv_digit(c);
00698     if (c < 0 || c >= base) {
00699         if (badcheck) goto bad;
00700         return INT2FIX(0);
00701     }
00702     len *= strlen(str)*sizeof(char);
00703 
00704     if ((size_t)len <= (sizeof(long)*CHAR_BIT)) {
00705         unsigned long val = STRTOUL(str, &end, base);
00706 
00707         if (str < end && *end == '_') goto bigparse;
00708         if (badcheck) {
00709             if (end == str) goto bad; /* no number */
00710             while (*end && ISSPACE(*end)) end++;
00711             if (*end) goto bad;       /* trailing garbage */
00712         }
00713 
00714         if (POSFIXABLE(val)) {
00715             if (sign) return LONG2FIX(val);
00716             else {
00717                 long result = -(long)val;
00718                 return LONG2FIX(result);
00719             }
00720         }
00721         else {
00722             VALUE big = rb_uint2big(val);
00723             RBIGNUM_SET_SIGN(big, sign);
00724             return bignorm(big);
00725         }
00726     }
00727   bigparse:
00728     len = (len/BITSPERDIG)+1;
00729     if (badcheck && *str == '_') goto bad;
00730 
00731     z = bignew(len, sign);
00732     zds = BDIGITS(z);
00733     for (i=len;i--;) zds[i]=0;
00734     while ((c = *str++) != 0) {
00735         if (c == '_') {
00736             if (nondigit) {
00737                 if (badcheck) goto bad;
00738                 break;
00739             }
00740             nondigit = (char) c;
00741             continue;
00742         }
00743         else if ((c = conv_digit(c)) < 0) {
00744             break;
00745         }
00746         if (c >= base) break;
00747         nondigit = 0;
00748         i = 0;
00749         num = c;
00750         for (;;) {
00751             while (i<blen) {
00752                 num += (BDIGIT_DBL)zds[i]*base;
00753                 zds[i++] = BIGLO(num);
00754                 num = BIGDN(num);
00755             }
00756             if (num) {
00757                 blen++;
00758                 continue;
00759             }
00760             break;
00761         }
00762     }
00763     if (badcheck) {
00764         str--;
00765         if (s+1 < str && str[-1] == '_') goto bad;
00766         while (*str && ISSPACE(*str)) str++;
00767         if (*str) {
00768           bad:
00769             rb_invalid_str(s, "Integer()");
00770         }
00771     }
00772 
00773     return bignorm(z);
00774 }
00775 
00776 VALUE
00777 rb_str_to_inum(VALUE str, int base, int badcheck)
00778 {
00779     char *s;
00780     long len;
00781     VALUE v = 0;
00782     VALUE ret;
00783 
00784     StringValue(str);
00785     rb_must_asciicompat(str);
00786     if (badcheck) {
00787         s = StringValueCStr(str);
00788     }
00789     else {
00790         s = RSTRING_PTR(str);
00791     }
00792     if (s) {
00793         len = RSTRING_LEN(str);
00794         if (s[len]) {           /* no sentinel somehow */
00795             char *p = ALLOCV(v, len+1);
00796 
00797             MEMCPY(p, s, char, len);
00798             p[len] = '\0';
00799             s = p;
00800         }
00801     }
00802     ret = rb_cstr_to_inum(s, base, badcheck);
00803     if (v)
00804         ALLOCV_END(v);
00805     return ret;
00806 }
00807 
00808 #if HAVE_LONG_LONG
00809 
00810 static VALUE
00811 rb_ull2big(unsigned LONG_LONG n)
00812 {
00813     BDIGIT_DBL num = n;
00814     long i = 0;
00815     BDIGIT *digits;
00816     VALUE big;
00817 
00818     big = bignew(DIGSPERLL, 1);
00819     digits = BDIGITS(big);
00820     while (i < DIGSPERLL) {
00821         digits[i++] = BIGLO(num);
00822         num = BIGDN(num);
00823     }
00824 
00825     i = DIGSPERLL;
00826     while (i-- && !digits[i]) ;
00827     RBIGNUM_SET_LEN(big, i+1);
00828     return big;
00829 }
00830 
00831 static VALUE
00832 rb_ll2big(LONG_LONG n)
00833 {
00834     long neg = 0;
00835     VALUE big;
00836 
00837     if (n < 0) {
00838         n = -n;
00839         neg = 1;
00840     }
00841     big = rb_ull2big(n);
00842     if (neg) {
00843         RBIGNUM_SET_SIGN(big, 0);
00844     }
00845     return big;
00846 }
00847 
00848 VALUE
00849 rb_ull2inum(unsigned LONG_LONG n)
00850 {
00851     if (POSFIXABLE(n)) return LONG2FIX(n);
00852     return rb_ull2big(n);
00853 }
00854 
00855 VALUE
00856 rb_ll2inum(LONG_LONG n)
00857 {
00858     if (FIXABLE(n)) return LONG2FIX(n);
00859     return rb_ll2big(n);
00860 }
00861 
00862 #endif  /* HAVE_LONG_LONG */
00863 
00864 VALUE
00865 rb_cstr2inum(const char *str, int base)
00866 {
00867     return rb_cstr_to_inum(str, base, base==0);
00868 }
00869 
00870 VALUE
00871 rb_str2inum(VALUE str, int base)
00872 {
00873     return rb_str_to_inum(str, base, base==0);
00874 }
00875 
00876 const char ruby_digitmap[] = "0123456789abcdefghijklmnopqrstuvwxyz";
00877 
00878 static VALUE bigsqr(VALUE x);
00879 static void bigdivmod(VALUE x, VALUE y, volatile VALUE *divp, volatile VALUE *modp);
00880 
00881 #define POW2_P(x) (((x)&((x)-1))==0)
00882 
00883 static inline int
00884 ones(register unsigned long x)
00885 {
00886 #if SIZEOF_LONG == 8
00887 # define MASK_55 0x5555555555555555UL
00888 # define MASK_33 0x3333333333333333UL
00889 # define MASK_0f 0x0f0f0f0f0f0f0f0fUL
00890 #else
00891 # define MASK_55 0x55555555UL
00892 # define MASK_33 0x33333333UL
00893 # define MASK_0f 0x0f0f0f0fUL
00894 #endif
00895     x -= (x >> 1) & MASK_55;
00896     x = ((x >> 2) & MASK_33) + (x & MASK_33);
00897     x = ((x >> 4) + x) & MASK_0f;
00898     x += (x >> 8);
00899     x += (x >> 16);
00900 #if SIZEOF_LONG == 8
00901     x += (x >> 32);
00902 #endif
00903     return (int)(x & 0x7f);
00904 #undef MASK_0f
00905 #undef MASK_33
00906 #undef MASK_55
00907 }
00908 
00909 static inline unsigned long
00910 next_pow2(register unsigned long x)
00911 {
00912     x |= x >> 1;
00913     x |= x >> 2;
00914     x |= x >> 4;
00915     x |= x >> 8;
00916     x |= x >> 16;
00917 #if SIZEOF_LONG == 8
00918     x |= x >> 32;
00919 #endif
00920     return x + 1;
00921 }
00922 
00923 static inline int
00924 floor_log2(register unsigned long x)
00925 {
00926     x |= x >> 1;
00927     x |= x >> 2;
00928     x |= x >> 4;
00929     x |= x >> 8;
00930     x |= x >> 16;
00931 #if SIZEOF_LONG == 8
00932     x |= x >> 32;
00933 #endif
00934     return (int)ones(x) - 1;
00935 }
00936 
00937 static inline int
00938 ceil_log2(register unsigned long x)
00939 {
00940     return floor_log2(x) + !POW2_P(x);
00941 }
00942 
00943 #define LOG2_KARATSUBA_DIGITS 7
00944 #define KARATSUBA_DIGITS (1L<<LOG2_KARATSUBA_DIGITS)
00945 #define MAX_BIG2STR_TABLE_ENTRIES 64
00946 
00947 static VALUE big2str_power_cache[35][MAX_BIG2STR_TABLE_ENTRIES];
00948 
00949 static void
00950 power_cache_init(void)
00951 {
00952     int i, j;
00953     for (i = 0; i < 35; ++i) {
00954         for (j = 0; j < MAX_BIG2STR_TABLE_ENTRIES; ++j) {
00955             big2str_power_cache[i][j] = Qnil;
00956         }
00957     }
00958 }
00959 
00960 static inline VALUE
00961 power_cache_get_power0(int base, int i)
00962 {
00963     if (NIL_P(big2str_power_cache[base - 2][i])) {
00964         big2str_power_cache[base - 2][i] =
00965             i == 0 ? rb_big_pow(rb_int2big(base), INT2FIX(KARATSUBA_DIGITS))
00966                    : bigsqr(power_cache_get_power0(base, i - 1));
00967         rb_gc_register_mark_object(big2str_power_cache[base - 2][i]);
00968     }
00969     return big2str_power_cache[base - 2][i];
00970 }
00971 
00972 static VALUE
00973 power_cache_get_power(int base, long n1, long* m1)
00974 {
00975     int i, m;
00976     long j;
00977     VALUE t;
00978 
00979     if (n1 <= KARATSUBA_DIGITS)
00980         rb_bug("n1 > KARATSUBA_DIGITS");
00981 
00982     m = ceil_log2(n1);
00983     if (m1) *m1 = 1 << m;
00984     i = m - LOG2_KARATSUBA_DIGITS;
00985     if (i >= MAX_BIG2STR_TABLE_ENTRIES)
00986         i = MAX_BIG2STR_TABLE_ENTRIES - 1;
00987     t = power_cache_get_power0(base, i);
00988 
00989     j = KARATSUBA_DIGITS*(1 << i);
00990     while (n1 > j) {
00991         t = bigsqr(t);
00992         j *= 2;
00993     }
00994     return t;
00995 }
00996 
00997 /* big2str_muraken_find_n1
00998  *
00999  * Let a natural number x is given by:
01000  * x = 2^0 * x_0 + 2^1 * x_1 + ... + 2^(B*n_0 - 1) * x_{B*n_0 - 1},
01001  * where B is BITSPERDIG (i.e. BDIGITS*CHAR_BIT) and n_0 is
01002  * RBIGNUM_LEN(x).
01003  *
01004  * Now, we assume n_1 = min_n \{ n | 2^(B*n_0/2) <= b_1^(n_1) \}, so
01005  * it is realized that 2^(B*n_0) <= {b_1}^{2*n_1}, where b_1 is a
01006  * given radix number. And then, we have n_1 <= (B*n_0) /
01007  * (2*log_2(b_1)), therefore n_1 is given by ceil((B*n_0) /
01008  * (2*log_2(b_1))).
01009  */
01010 static long
01011 big2str_find_n1(VALUE x, int base)
01012 {
01013     static const double log_2[] = {
01014         1.0,              1.58496250072116, 2.0,
01015         2.32192809488736, 2.58496250072116, 2.8073549220576,
01016         3.0,              3.16992500144231, 3.32192809488736,
01017         3.4594316186373,  3.58496250072116, 3.70043971814109,
01018         3.8073549220576,  3.90689059560852, 4.0,
01019         4.08746284125034, 4.16992500144231, 4.24792751344359,
01020         4.32192809488736, 4.39231742277876, 4.4594316186373,
01021         4.52356195605701, 4.58496250072116, 4.64385618977472,
01022         4.70043971814109, 4.75488750216347, 4.8073549220576,
01023         4.85798099512757, 4.90689059560852, 4.95419631038688,
01024         5.0,              5.04439411935845, 5.08746284125034,
01025         5.12928301694497, 5.16992500144231
01026     };
01027     long bits;
01028 
01029     if (base < 2 || 36 < base)
01030         rb_bug("invalid radix %d", base);
01031 
01032     if (FIXNUM_P(x)) {
01033         bits = (SIZEOF_LONG*CHAR_BIT - 1)/2 + 1;
01034     }
01035     else if (BIGZEROP(x)) {
01036         return 0;
01037     }
01038     else if (RBIGNUM_LEN(x) >= LONG_MAX/BITSPERDIG) {
01039         rb_raise(rb_eRangeError, "bignum too big to convert into `string'");
01040     }
01041     else {
01042         bits = BITSPERDIG*RBIGNUM_LEN(x);
01043     }
01044 
01045     /* @shyouhei note: vvvvvvvvvvvvv this cast is suspicious.  But I believe it is OK, because if that cast loses data, this x value is too big, and should have raised RangeError. */
01046     return (long)ceil(((double)bits)/log_2[base - 2]);
01047 }
01048 
01049 static long
01050 big2str_orig(VALUE x, int base, char* ptr, long len, long hbase, int trim)
01051 {
01052     long i = RBIGNUM_LEN(x), j = len;
01053     BDIGIT* ds = BDIGITS(x);
01054 
01055     while (i && j > 0) {
01056         long k = i;
01057         BDIGIT_DBL num = 0;
01058 
01059         while (k--) {               /* x / hbase */
01060             num = BIGUP(num) + ds[k];
01061             ds[k] = (BDIGIT)(num / hbase);
01062             num %= hbase;
01063         }
01064         if (trim && ds[i-1] == 0) i--;
01065         k = SIZEOF_BDIGITS;
01066         while (k--) {
01067             ptr[--j] = ruby_digitmap[num % base];
01068             num /= base;
01069             if (j <= 0) break;
01070             if (trim && i == 0 && num == 0) break;
01071         }
01072     }
01073     if (trim) {
01074         while (j < len && ptr[j] == '0') j++;
01075         MEMMOVE(ptr, ptr + j, char, len - j);
01076         len -= j;
01077     }
01078     return len;
01079 }
01080 
01081 static long
01082 big2str_karatsuba(VALUE x, int base, char* ptr,
01083                   long n1, long len, long hbase, int trim)
01084 {
01085     long lh, ll, m1;
01086     VALUE b, q, r;
01087 
01088     if (BIGZEROP(x)) {
01089         if (trim) return 0;
01090         else {
01091             memset(ptr, '0', len);
01092             return len;
01093         }
01094     }
01095 
01096     if (n1 <= KARATSUBA_DIGITS) {
01097         return big2str_orig(x, base, ptr, len, hbase, trim);
01098     }
01099 
01100     b = power_cache_get_power(base, n1, &m1);
01101     bigdivmod(x, b, &q, &r);
01102     lh = big2str_karatsuba(q, base, ptr, (len - m1)/2,
01103                            len - m1, hbase, trim);
01104     rb_big_resize(q, 0);
01105     ll = big2str_karatsuba(r, base, ptr + lh, m1/2,
01106                            m1, hbase, !lh && trim);
01107     rb_big_resize(r, 0);
01108 
01109     return lh + ll;
01110 }
01111 
01112 VALUE
01113 rb_big2str0(VALUE x, int base, int trim)
01114 {
01115     int off;
01116     VALUE ss, xx;
01117     long n1, n2, len, hbase;
01118     char* ptr;
01119 
01120     if (FIXNUM_P(x)) {
01121         return rb_fix2str(x, base);
01122     }
01123     if (BIGZEROP(x)) {
01124         return rb_usascii_str_new2("0");
01125     }
01126 
01127     if (base < 2 || 36 < base)
01128         rb_raise(rb_eArgError, "invalid radix %d", base);
01129 
01130     n2 = big2str_find_n1(x, base);
01131     n1 = (n2 + 1) / 2;
01132     ss = rb_usascii_str_new(0, n2 + 1); /* plus one for sign */
01133     ptr = RSTRING_PTR(ss);
01134     ptr[0] = RBIGNUM_SIGN(x) ? '+' : '-';
01135 
01136     hbase = base*base;
01137 #if SIZEOF_BDIGITS > 2
01138     hbase *= hbase;
01139 #endif
01140     off = !(trim && RBIGNUM_SIGN(x)); /* erase plus sign if trim */
01141     xx = rb_big_clone(x);
01142     RBIGNUM_SET_SIGN(xx, 1);
01143     if (n1 <= KARATSUBA_DIGITS) {
01144         len = off + big2str_orig(xx, base, ptr + off, n2, hbase, trim);
01145     }
01146     else {
01147         len = off + big2str_karatsuba(xx, base, ptr + off, n1,
01148                                       n2, hbase, trim);
01149     }
01150     rb_big_resize(xx, 0);
01151 
01152     ptr[len] = '\0';
01153     rb_str_resize(ss, len);
01154 
01155     return ss;
01156 }
01157 
01158 VALUE
01159 rb_big2str(VALUE x, int base)
01160 {
01161     return rb_big2str0(x, base, 1);
01162 }
01163 
01164 /*
01165  *  call-seq:
01166  *     big.to_s(base=10)   ->  string
01167  *
01168  *  Returns a string containing the representation of <i>big</i> radix
01169  *  <i>base</i> (2 through 36).
01170  *
01171  *     12345654321.to_s         #=> "12345654321"
01172  *     12345654321.to_s(2)      #=> "1011011111110110111011110000110001"
01173  *     12345654321.to_s(8)      #=> "133766736061"
01174  *     12345654321.to_s(16)     #=> "2dfdbbc31"
01175  *     78546939656932.to_s(36)  #=> "rubyrules"
01176  */
01177 
01178 static VALUE
01179 rb_big_to_s(int argc, VALUE *argv, VALUE x)
01180 {
01181     int base;
01182 
01183     if (argc == 0) base = 10;
01184     else {
01185         VALUE b;
01186 
01187         rb_scan_args(argc, argv, "01", &b);
01188         base = NUM2INT(b);
01189     }
01190     return rb_big2str(x, base);
01191 }
01192 
01193 static VALUE
01194 big2ulong(VALUE x, const char *type, int check)
01195 {
01196     long len = RBIGNUM_LEN(x);
01197     BDIGIT_DBL num;
01198     BDIGIT *ds;
01199 
01200     if (len > DIGSPERLONG) {
01201         if (check)
01202             rb_raise(rb_eRangeError, "bignum too big to convert into `%s'", type);
01203         len = DIGSPERLONG;
01204     }
01205     ds = BDIGITS(x);
01206     num = 0;
01207     while (len--) {
01208         num = BIGUP(num);
01209         num += ds[len];
01210     }
01211     return (VALUE)num;
01212 }
01213 
01214 VALUE
01215 rb_big2ulong_pack(VALUE x)
01216 {
01217     VALUE num = big2ulong(x, "unsigned long", FALSE);
01218     if (!RBIGNUM_SIGN(x)) {
01219         return (VALUE)(-(SIGNED_VALUE)num);
01220     }
01221     return num;
01222 }
01223 
01224 VALUE
01225 rb_big2ulong(VALUE x)
01226 {
01227     VALUE num = big2ulong(x, "unsigned long", TRUE);
01228 
01229     if (RBIGNUM_POSITIVE_P(x)) {
01230         return num;
01231     }
01232     else {
01233         if (num <= LONG_MAX)
01234             return -(long)num;
01235         if (num == 1+(unsigned long)(-(LONG_MIN+1)))
01236             return LONG_MIN;
01237         rb_raise(rb_eRangeError, "bignum out of range of unsigned long");
01238     }
01239     return num;
01240 }
01241 
01242 SIGNED_VALUE
01243 rb_big2long(VALUE x)
01244 {
01245     VALUE num = big2ulong(x, "long", TRUE);
01246 
01247     if (RBIGNUM_POSITIVE_P(x)) {
01248         if (LONG_MAX < num)
01249             rb_raise(rb_eRangeError, "bignum too big to convert into `long'");
01250         return num;
01251     }
01252     else {
01253         if (num <= LONG_MAX)
01254             return -(long)num;
01255         if (num == 1+(unsigned long)(-(LONG_MIN+1)))
01256             return LONG_MIN;
01257         rb_raise(rb_eRangeError, "bignum too big to convert into `long'");
01258     }
01259 }
01260 
01261 #if HAVE_LONG_LONG
01262 
01263 static unsigned LONG_LONG
01264 big2ull(VALUE x, const char *type)
01265 {
01266     long len = RBIGNUM_LEN(x);
01267     BDIGIT_DBL num;
01268     BDIGIT *ds;
01269 
01270     if (len > SIZEOF_LONG_LONG/SIZEOF_BDIGITS)
01271         rb_raise(rb_eRangeError, "bignum too big to convert into `%s'", type);
01272     ds = BDIGITS(x);
01273     num = 0;
01274     while (len--) {
01275         num = BIGUP(num);
01276         num += ds[len];
01277     }
01278     return num;
01279 }
01280 
01281 unsigned LONG_LONG
01282 rb_big2ull(VALUE x)
01283 {
01284     unsigned LONG_LONG num = big2ull(x, "unsigned long long");
01285 
01286     if (RBIGNUM_POSITIVE_P(x)) {
01287         return num;
01288     }
01289     else {
01290         if (num <= LLONG_MAX)
01291             return -(LONG_LONG)num;
01292         if (num == 1+(unsigned LONG_LONG)(-(LLONG_MIN+1)))
01293             return LLONG_MIN;
01294         rb_raise(rb_eRangeError, "bignum out of range of unsigned long long");
01295     }
01296     return num;
01297 }
01298 
01299 LONG_LONG
01300 rb_big2ll(VALUE x)
01301 {
01302     unsigned LONG_LONG num = big2ull(x, "long long");
01303 
01304     if (RBIGNUM_POSITIVE_P(x)) {
01305         if (LLONG_MAX < num)
01306             rb_raise(rb_eRangeError, "bignum too big to convert into `long long'");
01307         return num;
01308     }
01309     else {
01310         if (num <= LLONG_MAX)
01311             return -(LONG_LONG)num;
01312         if (num == 1+(unsigned LONG_LONG)(-(LLONG_MIN+1)))
01313             return LLONG_MIN;
01314         rb_raise(rb_eRangeError, "bignum too big to convert into `long long'");
01315     }
01316 }
01317 
01318 #endif  /* HAVE_LONG_LONG */
01319 
01320 static VALUE
01321 dbl2big(double d)
01322 {
01323     long i = 0;
01324     BDIGIT c;
01325     BDIGIT *digits;
01326     VALUE z;
01327     double u = (d < 0)?-d:d;
01328 
01329     if (isinf(d)) {
01330         rb_raise(rb_eFloatDomainError, d < 0 ? "-Infinity" : "Infinity");
01331     }
01332     if (isnan(d)) {
01333         rb_raise(rb_eFloatDomainError, "NaN");
01334     }
01335 
01336     while (!POSFIXABLE(u) || 0 != (long)u) {
01337         u /= (double)(BIGRAD);
01338         i++;
01339     }
01340     z = bignew(i, d>=0);
01341     digits = BDIGITS(z);
01342     while (i--) {
01343         u *= BIGRAD;
01344         c = (BDIGIT)u;
01345         u -= c;
01346         digits[i] = c;
01347     }
01348 
01349     return z;
01350 }
01351 
01352 VALUE
01353 rb_dbl2big(double d)
01354 {
01355     return bignorm(dbl2big(d));
01356 }
01357 
01358 static int
01359 nlz(BDIGIT x)
01360 {
01361     BDIGIT y;
01362     int n = BITSPERDIG;
01363 #if BITSPERDIG > 64
01364     y = x >> 64; if (y) {n -= 64; x = y;}
01365 #endif
01366 #if BITSPERDIG > 32
01367     y = x >> 32; if (y) {n -= 32; x = y;}
01368 #endif
01369 #if BITSPERDIG > 16
01370     y = x >> 16; if (y) {n -= 16; x = y;}
01371 #endif
01372     y = x >>  8; if (y) {n -=  8; x = y;}
01373     y = x >>  4; if (y) {n -=  4; x = y;}
01374     y = x >>  2; if (y) {n -=  2; x = y;}
01375     y = x >>  1; if (y) {return n - 2;}
01376     return n - x;
01377 }
01378 
01379 static double
01380 big2dbl(VALUE x)
01381 {
01382     double d = 0.0;
01383     long i = (bigtrunc(x), RBIGNUM_LEN(x)), lo = 0, bits;
01384     BDIGIT *ds = BDIGITS(x), dl;
01385 
01386     if (i) {
01387         bits = i * BITSPERDIG - nlz(ds[i-1]);
01388         if (bits > DBL_MANT_DIG+DBL_MAX_EXP) {
01389             d = HUGE_VAL;
01390         }
01391         else {
01392             if (bits > DBL_MANT_DIG+1)
01393                 lo = (bits -= DBL_MANT_DIG+1) / BITSPERDIG;
01394             else
01395                 bits = 0;
01396             while (--i > lo) {
01397                 d = ds[i] + BIGRAD*d;
01398             }
01399             dl = ds[i];
01400             if (bits && (dl & (1UL << (bits %= BITSPERDIG)))) {
01401                 int carry = dl & ~(~(BDIGIT)0 << bits);
01402                 if (!carry) {
01403                     while (i-- > 0) {
01404                         if ((carry = ds[i]) != 0) break;
01405                     }
01406                 }
01407                 if (carry) {
01408                     dl &= (BDIGIT)~0 << bits;
01409                     dl += (BDIGIT)1 << bits;
01410                     if (!dl) d += 1;
01411                 }
01412             }
01413             d = dl + BIGRAD*d;
01414             if (lo) {
01415                 if (lo > INT_MAX / BITSPERDIG)
01416                     d = HUGE_VAL;
01417                 else if (lo < INT_MIN / BITSPERDIG)
01418                     d = 0.0;
01419                 else
01420                     d = ldexp(d, (int)(lo * BITSPERDIG));
01421             }
01422         }
01423     }
01424     if (!RBIGNUM_SIGN(x)) d = -d;
01425     return d;
01426 }
01427 
01428 double
01429 rb_big2dbl(VALUE x)
01430 {
01431     double d = big2dbl(x);
01432 
01433     if (isinf(d)) {
01434         rb_warning("Bignum out of Float range");
01435         if (d < 0.0)
01436             d = -HUGE_VAL;
01437         else
01438             d = HUGE_VAL;
01439     }
01440     return d;
01441 }
01442 
01443 /*
01444  *  call-seq:
01445  *     big.to_f -> float
01446  *
01447  *  Converts <i>big</i> to a <code>Float</code>. If <i>big</i> doesn't
01448  *  fit in a <code>Float</code>, the result is infinity.
01449  *
01450  */
01451 
01452 static VALUE
01453 rb_big_to_f(VALUE x)
01454 {
01455     return DBL2NUM(rb_big2dbl(x));
01456 }
01457 
01458 VALUE
01459 rb_integer_float_cmp(VALUE x, VALUE y)
01460 {
01461     double yd = RFLOAT_VALUE(y);
01462     double yi, yf;
01463     VALUE rel;
01464 
01465     if (isnan(yd))
01466         return Qnil;
01467     if (isinf(yd)) {
01468         if (yd > 0.0) return INT2FIX(-1);
01469         else return INT2FIX(1);
01470     }
01471     yf = modf(yd, &yi);
01472     if (FIXNUM_P(x)) {
01473 #if SIZEOF_LONG * CHAR_BIT < DBL_MANT_DIG /* assume FLT_RADIX == 2 */
01474         double xd = (double)FIX2LONG(x);
01475         if (xd < yd)
01476             return INT2FIX(-1);
01477         if (xd > yd)
01478             return INT2FIX(1);
01479         return INT2FIX(0);
01480 #else
01481         long xl, yl;
01482         if (yi < FIXNUM_MIN)
01483             return INT2FIX(1);
01484         if (FIXNUM_MAX+1 <= yi)
01485             return INT2FIX(-1);
01486         xl = FIX2LONG(x);
01487         yl = (long)yi;
01488         if (xl < yl)
01489             return INT2FIX(-1);
01490         if (xl > yl)
01491             return INT2FIX(1);
01492         if (yf < 0.0)
01493             return INT2FIX(1);
01494         if (0.0 < yf)
01495             return INT2FIX(-1);
01496         return INT2FIX(0);
01497 #endif
01498     }
01499     y = rb_dbl2big(yi);
01500     rel = rb_big_cmp(x, y);
01501     if (yf == 0.0 || rel != INT2FIX(0))
01502         return rel;
01503     if (yf < 0.0)
01504         return INT2FIX(1);
01505     return INT2FIX(-1);
01506 }
01507 
01508 VALUE
01509 rb_integer_float_eq(VALUE x, VALUE y)
01510 {
01511     double yd = RFLOAT_VALUE(y);
01512     double yi, yf;
01513 
01514     if (isnan(yd) || isinf(yd))
01515         return Qfalse;
01516     yf = modf(yd, &yi);
01517     if (yf != 0)
01518         return Qfalse;
01519     if (FIXNUM_P(x)) {
01520 #if SIZEOF_LONG * CHAR_BIT < DBL_MANT_DIG /* assume FLT_RADIX == 2 */
01521         double xd = (double)FIX2LONG(x);
01522         if (xd != yd)
01523             return Qfalse;
01524         return Qtrue;
01525 #else
01526         long xl, yl;
01527         if (yi < LONG_MIN || LONG_MAX < yi)
01528             return Qfalse;
01529         xl = FIX2LONG(x);
01530         yl = (long)yi;
01531         if (xl != yl)
01532             return Qfalse;
01533         return Qtrue;
01534 #endif
01535     }
01536     y = rb_dbl2big(yi);
01537     return rb_big_eq(x, y);
01538 }
01539 
01540 /*
01541  *  call-seq:
01542  *     big <=> numeric   -> -1, 0, +1 or nil
01543  *
01544  *  Comparison---Returns -1, 0, or +1 depending on whether +big+ is
01545  *  less than, equal to, or greater than +numeric+. This is the
01546  *  basis for the tests in Comparable.
01547  *
01548  *  +nil+ is returned if the two values are incomparable.
01549  *
01550  */
01551 
01552 VALUE
01553 rb_big_cmp(VALUE x, VALUE y)
01554 {
01555     long xlen = RBIGNUM_LEN(x);
01556     BDIGIT *xds, *yds;
01557 
01558     switch (TYPE(y)) {
01559       case T_FIXNUM:
01560         y = rb_int2big(FIX2LONG(y));
01561         break;
01562 
01563       case T_BIGNUM:
01564         break;
01565 
01566       case T_FLOAT:
01567         return rb_integer_float_cmp(x, y);
01568 
01569       default:
01570         return rb_num_coerce_cmp(x, y, rb_intern("<=>"));
01571     }
01572 
01573     if (RBIGNUM_SIGN(x) > RBIGNUM_SIGN(y)) return INT2FIX(1);
01574     if (RBIGNUM_SIGN(x) < RBIGNUM_SIGN(y)) return INT2FIX(-1);
01575     if (xlen < RBIGNUM_LEN(y))
01576         return (RBIGNUM_SIGN(x)) ? INT2FIX(-1) : INT2FIX(1);
01577     if (xlen > RBIGNUM_LEN(y))
01578         return (RBIGNUM_SIGN(x)) ? INT2FIX(1) : INT2FIX(-1);
01579 
01580     xds = BDIGITS(x);
01581     yds = BDIGITS(y);
01582 
01583     while (xlen-- && (xds[xlen]==yds[xlen]));
01584     if (-1 == xlen) return INT2FIX(0);
01585     return (xds[xlen] > yds[xlen]) ?
01586         (RBIGNUM_SIGN(x) ? INT2FIX(1) : INT2FIX(-1)) :
01587             (RBIGNUM_SIGN(x) ? INT2FIX(-1) : INT2FIX(1));
01588 }
01589 
01590 enum big_op_t {
01591     big_op_gt,
01592     big_op_ge,
01593     big_op_lt,
01594     big_op_le
01595 };
01596 
01597 static VALUE
01598 big_op(VALUE x, VALUE y, enum big_op_t op)
01599 {
01600     VALUE rel;
01601     int n;
01602 
01603     switch (TYPE(y)) {
01604       case T_FIXNUM:
01605       case T_BIGNUM:
01606         rel = rb_big_cmp(x, y);
01607         break;
01608 
01609       case T_FLOAT:
01610         rel = rb_integer_float_cmp(x, y);
01611         break;
01612 
01613       default:
01614         {
01615             ID id = 0;
01616             switch (op) {
01617                 case big_op_gt: id = '>'; break;
01618                 case big_op_ge: id = rb_intern(">="); break;
01619                 case big_op_lt: id = '<'; break;
01620                 case big_op_le: id = rb_intern("<="); break;
01621             }
01622             return rb_num_coerce_relop(x, y, id);
01623         }
01624     }
01625 
01626     if (NIL_P(rel)) return Qfalse;
01627     n = FIX2INT(rel);
01628 
01629     switch (op) {
01630         case big_op_gt: return n >  0 ? Qtrue : Qfalse;
01631         case big_op_ge: return n >= 0 ? Qtrue : Qfalse;
01632         case big_op_lt: return n <  0 ? Qtrue : Qfalse;
01633         case big_op_le: return n <= 0 ? Qtrue : Qfalse;
01634     }
01635     return Qundef;
01636 }
01637 
01638 /*
01639  * call-seq:
01640  *   big > real  ->  true or false
01641  *
01642  * Returns <code>true</code> if the value of <code>big</code> is
01643  * greater than that of <code>real</code>.
01644  */
01645 
01646 static VALUE
01647 big_gt(VALUE x, VALUE y)
01648 {
01649     return big_op(x, y, big_op_gt);
01650 }
01651 
01652 /*
01653  * call-seq:
01654  *   big >= real  ->  true or false
01655  *
01656  * Returns <code>true</code> if the value of <code>big</code> is
01657  * greater than or equal to that of <code>real</code>.
01658  */
01659 
01660 static VALUE
01661 big_ge(VALUE x, VALUE y)
01662 {
01663     return big_op(x, y, big_op_ge);
01664 }
01665 
01666 /*
01667  * call-seq:
01668  *   big < real  ->  true or false
01669  *
01670  * Returns <code>true</code> if the value of <code>big</code> is
01671  * less than that of <code>real</code>.
01672  */
01673 
01674 static VALUE
01675 big_lt(VALUE x, VALUE y)
01676 {
01677     return big_op(x, y, big_op_lt);
01678 }
01679 
01680 /*
01681  * call-seq:
01682  *   big <= real  ->  true or false
01683  *
01684  * Returns <code>true</code> if the value of <code>big</code> is
01685  * less than or equal to that of <code>real</code>.
01686  */
01687 
01688 static VALUE
01689 big_le(VALUE x, VALUE y)
01690 {
01691     return big_op(x, y, big_op_le);
01692 }
01693 
01694 /*
01695  *  call-seq:
01696  *     big == obj  -> true or false
01697  *
01698  *  Returns <code>true</code> only if <i>obj</i> has the same value
01699  *  as <i>big</i>. Contrast this with <code>Bignum#eql?</code>, which
01700  *  requires <i>obj</i> to be a <code>Bignum</code>.
01701  *
01702  *     68719476736 == 68719476736.0   #=> true
01703  */
01704 
01705 VALUE
01706 rb_big_eq(VALUE x, VALUE y)
01707 {
01708     switch (TYPE(y)) {
01709       case T_FIXNUM:
01710         y = rb_int2big(FIX2LONG(y));
01711         break;
01712       case T_BIGNUM:
01713         break;
01714       case T_FLOAT:
01715         return rb_integer_float_eq(x, y);
01716       default:
01717         return rb_equal(y, x);
01718     }
01719     if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y)) return Qfalse;
01720     if (RBIGNUM_LEN(x) != RBIGNUM_LEN(y)) return Qfalse;
01721     if (MEMCMP(BDIGITS(x),BDIGITS(y),BDIGIT,RBIGNUM_LEN(y)) != 0) return Qfalse;
01722     return Qtrue;
01723 }
01724 
01725 /*
01726  *  call-seq:
01727  *     big.eql?(obj)   -> true or false
01728  *
01729  *  Returns <code>true</code> only if <i>obj</i> is a
01730  *  <code>Bignum</code> with the same value as <i>big</i>. Contrast this
01731  *  with <code>Bignum#==</code>, which performs type conversions.
01732  *
01733  *     68719476736.eql?(68719476736.0)   #=> false
01734  */
01735 
01736 VALUE
01737 rb_big_eql(VALUE x, VALUE y)
01738 {
01739     if (!RB_TYPE_P(y, T_BIGNUM)) return Qfalse;
01740     if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y)) return Qfalse;
01741     if (RBIGNUM_LEN(x) != RBIGNUM_LEN(y)) return Qfalse;
01742     if (MEMCMP(BDIGITS(x),BDIGITS(y),BDIGIT,RBIGNUM_LEN(y)) != 0) return Qfalse;
01743     return Qtrue;
01744 }
01745 
01746 /*
01747  * call-seq:
01748  *    -big   ->  integer
01749  *
01750  * Unary minus (returns an integer whose value is 0-big)
01751  */
01752 
01753 VALUE
01754 rb_big_uminus(VALUE x)
01755 {
01756     VALUE z = rb_big_clone(x);
01757 
01758     RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(x));
01759 
01760     return bignorm(z);
01761 }
01762 
01763 /*
01764  * call-seq:
01765  *     ~big  ->  integer
01766  *
01767  * Inverts the bits in big. As Bignums are conceptually infinite
01768  * length, the result acts as if it had an infinite number of one
01769  * bits to the left. In hex representations, this is displayed
01770  * as two periods to the left of the digits.
01771  *
01772  *   sprintf("%X", ~0x1122334455)    #=> "..FEEDDCCBBAA"
01773  */
01774 
01775 static VALUE
01776 rb_big_neg(VALUE x)
01777 {
01778     VALUE z = rb_big_clone(x);
01779     BDIGIT *ds;
01780     long i;
01781 
01782     if (!RBIGNUM_SIGN(x)) get2comp(z);
01783     ds = BDIGITS(z);
01784     i = RBIGNUM_LEN(x);
01785     if (!i) return INT2FIX(~(SIGNED_VALUE)0);
01786     while (i--) {
01787         ds[i] = ~ds[i];
01788     }
01789     RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(z));
01790     if (RBIGNUM_SIGN(x)) get2comp(z);
01791 
01792     return bignorm(z);
01793 }
01794 
01795 static void
01796 bigsub_core(BDIGIT *xds, long xn, BDIGIT *yds, long yn, BDIGIT *zds, long zn)
01797 {
01798     BDIGIT_DBL_SIGNED num;
01799     long i;
01800 
01801     for (i = 0, num = 0; i < yn; i++) {
01802         num += (BDIGIT_DBL_SIGNED)xds[i] - yds[i];
01803         zds[i] = BIGLO(num);
01804         num = BIGDN(num);
01805     }
01806     while (num && i < xn) {
01807         num += xds[i];
01808         zds[i++] = BIGLO(num);
01809         num = BIGDN(num);
01810     }
01811     while (i < xn) {
01812         zds[i] = xds[i];
01813         i++;
01814     }
01815     assert(i <= zn);
01816     while (i < zn) {
01817         zds[i++] = 0;
01818     }
01819 }
01820 
01821 static VALUE
01822 bigsub(VALUE x, VALUE y)
01823 {
01824     VALUE z = 0;
01825     long i = RBIGNUM_LEN(x);
01826     BDIGIT *xds, *yds;
01827 
01828     /* if x is smaller than y, swap */
01829     if (RBIGNUM_LEN(x) < RBIGNUM_LEN(y)) {
01830         z = x; x = y; y = z;    /* swap x y */
01831     }
01832     else if (RBIGNUM_LEN(x) == RBIGNUM_LEN(y)) {
01833         xds = BDIGITS(x);
01834         yds = BDIGITS(y);
01835         while (i > 0) {
01836             i--;
01837             if (xds[i] > yds[i]) {
01838                 break;
01839             }
01840             if (xds[i] < yds[i]) {
01841                 z = x; x = y; y = z;    /* swap x y */
01842                 break;
01843             }
01844         }
01845     }
01846 
01847     z = bignew(RBIGNUM_LEN(x), z==0);
01848     bigsub_core(BDIGITS(x), RBIGNUM_LEN(x),
01849                 BDIGITS(y), RBIGNUM_LEN(y),
01850                 BDIGITS(z), RBIGNUM_LEN(z));
01851 
01852     return z;
01853 }
01854 
01855 static VALUE bigadd_int(VALUE x, long y);
01856 
01857 static VALUE
01858 bigsub_int(VALUE x, long y0)
01859 {
01860     VALUE z;
01861     BDIGIT *xds, *zds;
01862     long xn;
01863     BDIGIT_DBL_SIGNED num;
01864     long i, y;
01865 
01866     y = y0;
01867     xds = BDIGITS(x);
01868     xn = RBIGNUM_LEN(x);
01869 
01870     z = bignew(xn, RBIGNUM_SIGN(x));
01871     zds = BDIGITS(z);
01872 
01873 #if SIZEOF_BDIGITS == SIZEOF_LONG
01874     num = (BDIGIT_DBL_SIGNED)xds[0] - y;
01875     if (xn == 1 && num < 0) {
01876         RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(x));
01877         zds[0] = (BDIGIT)-num;
01878         RB_GC_GUARD(x);
01879         return bignorm(z);
01880     }
01881     zds[0] = BIGLO(num);
01882     num = BIGDN(num);
01883     i = 1;
01884 #else
01885     num = 0;
01886     for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) {
01887         num += (BDIGIT_DBL_SIGNED)xds[i] - BIGLO(y);
01888         zds[i] = BIGLO(num);
01889         num = BIGDN(num);
01890         y = BIGDN(y);
01891     }
01892 #endif
01893     while (num && i < xn) {
01894         num += xds[i];
01895         zds[i++] = BIGLO(num);
01896         num = BIGDN(num);
01897     }
01898     while (i < xn) {
01899         zds[i] = xds[i];
01900         i++;
01901     }
01902     if (num < 0) {
01903         z = bigsub(x, rb_int2big(y0));
01904     }
01905     RB_GC_GUARD(x);
01906     return bignorm(z);
01907 }
01908 
01909 static VALUE
01910 bigadd_int(VALUE x, long y)
01911 {
01912     VALUE z;
01913     BDIGIT *xds, *zds;
01914     long xn, zn;
01915     BDIGIT_DBL num;
01916     long i;
01917 
01918     xds = BDIGITS(x);
01919     xn = RBIGNUM_LEN(x);
01920 
01921     if (xn < 2) {
01922         zn = 3;
01923     }
01924     else {
01925         zn = xn + 1;
01926     }
01927     z = bignew(zn, RBIGNUM_SIGN(x));
01928     zds = BDIGITS(z);
01929 
01930 #if SIZEOF_BDIGITS == SIZEOF_LONG
01931     num = (BDIGIT_DBL)xds[0] + y;
01932     zds[0] = BIGLO(num);
01933     num = BIGDN(num);
01934     i = 1;
01935 #else
01936     num = 0;
01937     for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) {
01938         num += (BDIGIT_DBL)xds[i] + BIGLO(y);
01939         zds[i] = BIGLO(num);
01940         num = BIGDN(num);
01941         y = BIGDN(y);
01942     }
01943 #endif
01944     while (num && i < xn) {
01945         num += xds[i];
01946         zds[i++] = BIGLO(num);
01947         num = BIGDN(num);
01948     }
01949     if (num) zds[i++] = (BDIGIT)num;
01950     else while (i < xn) {
01951         zds[i] = xds[i];
01952         i++;
01953     }
01954     assert(i <= zn);
01955     while (i < zn) {
01956         zds[i++] = 0;
01957     }
01958     RB_GC_GUARD(x);
01959     return bignorm(z);
01960 }
01961 
01962 static void
01963 bigadd_core(BDIGIT *xds, long xn, BDIGIT *yds, long yn, BDIGIT *zds, long zn)
01964 {
01965     BDIGIT_DBL num = 0;
01966     long i;
01967 
01968     if (xn > yn) {
01969         BDIGIT *tds;
01970         tds = xds; xds = yds; yds = tds;
01971         i = xn; xn = yn; yn = i;
01972     }
01973 
01974     i = 0;
01975     while (i < xn) {
01976         num += (BDIGIT_DBL)xds[i] + yds[i];
01977         zds[i++] = BIGLO(num);
01978         num = BIGDN(num);
01979     }
01980     while (num && i < yn) {
01981         num += yds[i];
01982         zds[i++] = BIGLO(num);
01983         num = BIGDN(num);
01984     }
01985     while (i < yn) {
01986         zds[i] = yds[i];
01987         i++;
01988     }
01989     if (num) zds[i++] = (BDIGIT)num;
01990     assert(i <= zn);
01991     while (i < zn) {
01992         zds[i++] = 0;
01993     }
01994 }
01995 
01996 static VALUE
01997 bigadd(VALUE x, VALUE y, int sign)
01998 {
01999     VALUE z;
02000     long len;
02001 
02002     sign = (sign == RBIGNUM_SIGN(y));
02003     if (RBIGNUM_SIGN(x) != sign) {
02004         if (sign) return bigsub(y, x);
02005         return bigsub(x, y);
02006     }
02007 
02008     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
02009         len = RBIGNUM_LEN(x) + 1;
02010     }
02011     else {
02012         len = RBIGNUM_LEN(y) + 1;
02013     }
02014     z = bignew(len, sign);
02015 
02016     bigadd_core(BDIGITS(x), RBIGNUM_LEN(x),
02017                 BDIGITS(y), RBIGNUM_LEN(y),
02018                 BDIGITS(z), RBIGNUM_LEN(z));
02019 
02020     return z;
02021 }
02022 
02023 /*
02024  *  call-seq:
02025  *     big + other  -> Numeric
02026  *
02027  *  Adds big and other, returning the result.
02028  */
02029 
02030 VALUE
02031 rb_big_plus(VALUE x, VALUE y)
02032 {
02033     long n;
02034 
02035     switch (TYPE(y)) {
02036       case T_FIXNUM:
02037         n = FIX2LONG(y);
02038         if ((n > 0) != RBIGNUM_SIGN(x)) {
02039             if (n < 0) {
02040                 n = -n;
02041             }
02042             return bigsub_int(x, n);
02043         }
02044         if (n < 0) {
02045             n = -n;
02046         }
02047         return bigadd_int(x, n);
02048 
02049       case T_BIGNUM:
02050         return bignorm(bigadd(x, y, 1));
02051 
02052       case T_FLOAT:
02053         return DBL2NUM(rb_big2dbl(x) + RFLOAT_VALUE(y));
02054 
02055       default:
02056         return rb_num_coerce_bin(x, y, '+');
02057     }
02058 }
02059 
02060 /*
02061  *  call-seq:
02062  *     big - other  -> Numeric
02063  *
02064  *  Subtracts other from big, returning the result.
02065  */
02066 
02067 VALUE
02068 rb_big_minus(VALUE x, VALUE y)
02069 {
02070     long n;
02071 
02072     switch (TYPE(y)) {
02073       case T_FIXNUM:
02074         n = FIX2LONG(y);
02075         if ((n > 0) != RBIGNUM_SIGN(x)) {
02076             if (n < 0) {
02077                 n = -n;
02078             }
02079             return bigadd_int(x, n);
02080         }
02081         if (n < 0) {
02082             n = -n;
02083         }
02084         return bigsub_int(x, n);
02085 
02086       case T_BIGNUM:
02087         return bignorm(bigadd(x, y, 0));
02088 
02089       case T_FLOAT:
02090         return DBL2NUM(rb_big2dbl(x) - RFLOAT_VALUE(y));
02091 
02092       default:
02093         return rb_num_coerce_bin(x, y, '-');
02094     }
02095 }
02096 
02097 static long
02098 big_real_len(VALUE x)
02099 {
02100     long i = RBIGNUM_LEN(x);
02101     BDIGIT *xds = BDIGITS(x);
02102     while (--i && !xds[i]);
02103     return i + 1;
02104 }
02105 
02106 static VALUE
02107 bigmul1_single(VALUE x, VALUE y)
02108 {
02109     BDIGIT_DBL n;
02110     VALUE z = bignew(2, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
02111     BDIGIT *xds, *yds, *zds;
02112 
02113     xds = BDIGITS(x);
02114     yds = BDIGITS(y);
02115     zds = BDIGITS(z);
02116 
02117     n = (BDIGIT_DBL)xds[0] * yds[0];
02118     zds[0] = BIGLO(n);
02119     zds[1] = (BDIGIT)BIGDN(n);
02120 
02121     return z;
02122 }
02123 
02124 static VALUE
02125 bigmul1_normal(VALUE x, VALUE y)
02126 {
02127     long xl = RBIGNUM_LEN(x), yl = RBIGNUM_LEN(y), i, j = xl + yl + 1;
02128     BDIGIT_DBL n = 0;
02129     VALUE z = bignew(j, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
02130     BDIGIT *xds, *yds, *zds;
02131 
02132     xds = BDIGITS(x);
02133     yds = BDIGITS(y);
02134     zds = BDIGITS(z);
02135     while (j--) zds[j] = 0;
02136     for (i = 0; i < xl; i++) {
02137         BDIGIT_DBL dd;
02138         dd = xds[i];
02139         if (dd == 0) continue;
02140         n = 0;
02141         for (j = 0; j < yl; j++) {
02142             BDIGIT_DBL ee = n + (BDIGIT_DBL)dd * yds[j];
02143             n = zds[i + j] + ee;
02144             if (ee) zds[i + j] = BIGLO(n);
02145             n = BIGDN(n);
02146         }
02147         if (n) {
02148             zds[i + j] = (BDIGIT)n;
02149         }
02150     }
02151     rb_thread_check_ints();
02152     return z;
02153 }
02154 
02155 static VALUE bigmul0(VALUE x, VALUE y);
02156 
02157 /* balancing multiplication by slicing larger argument */
02158 static VALUE
02159 bigmul1_balance(VALUE x, VALUE y)
02160 {
02161     VALUE z, t1, t2;
02162     long i, xn, yn, r, n;
02163     BDIGIT *yds, *zds, *t1ds;
02164 
02165     xn = RBIGNUM_LEN(x);
02166     yn = RBIGNUM_LEN(y);
02167     assert(2 * xn <= yn || 3 * xn <= 2*(yn+2));
02168 
02169     z = bignew(xn + yn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
02170     t1 = bignew(xn, 1);
02171 
02172     yds = BDIGITS(y);
02173     zds = BDIGITS(z);
02174     t1ds = BDIGITS(t1);
02175 
02176     for (i = 0; i < xn + yn; i++) zds[i] = 0;
02177 
02178     n = 0;
02179     while (yn > 0) {
02180         r = xn > yn ? yn : xn;
02181         MEMCPY(t1ds, yds + n, BDIGIT, r);
02182         RBIGNUM_SET_LEN(t1, r);
02183         t2 = bigmul0(x, t1);
02184         bigadd_core(zds + n, RBIGNUM_LEN(z) - n,
02185                     BDIGITS(t2), big_real_len(t2),
02186                     zds + n, RBIGNUM_LEN(z) - n);
02187         yn -= r;
02188         n += r;
02189     }
02190 
02191     return z;
02192 }
02193 
02194 /* split a bignum into high and low bignums */
02195 static void
02196 big_split(VALUE v, long n, volatile VALUE *ph, volatile VALUE *pl)
02197 {
02198     long hn = 0, ln = RBIGNUM_LEN(v);
02199     VALUE h, l;
02200     BDIGIT *vds = BDIGITS(v);
02201 
02202     if (ln > n) {
02203         hn = ln - n;
02204         ln = n;
02205     }
02206 
02207     if (!hn) {
02208         h = rb_uint2big(0);
02209     }
02210     else {
02211         while (--hn && !vds[hn + ln]);
02212         h = bignew(hn += 2, 1);
02213         MEMCPY(BDIGITS(h), vds + ln, BDIGIT, hn - 1);
02214         BDIGITS(h)[hn - 1] = 0; /* margin for carry */
02215     }
02216 
02217     while (--ln && !vds[ln]);
02218     l = bignew(ln += 2, 1);
02219     MEMCPY(BDIGITS(l), vds, BDIGIT, ln - 1);
02220     BDIGITS(l)[ln - 1] = 0; /* margin for carry */
02221 
02222     *pl = l;
02223     *ph = h;
02224 }
02225 
02226 /* multiplication by karatsuba method */
02227 static VALUE
02228 bigmul1_karatsuba(VALUE x, VALUE y)
02229 {
02230     long i, n, xn, yn, t1n, t2n;
02231     VALUE xh, xl, yh, yl, z, t1, t2, t3;
02232     BDIGIT *zds;
02233 
02234     xn = RBIGNUM_LEN(x);
02235     yn = RBIGNUM_LEN(y);
02236     n = yn / 2;
02237     big_split(x, n, &xh, &xl);
02238     if (x == y) {
02239         yh = xh; yl = xl;
02240     }
02241     else big_split(y, n, &yh, &yl);
02242 
02243     /* x = xh * b + xl
02244      * y = yh * b + yl
02245      *
02246      * Karatsuba method:
02247      *   x * y = z2 * b^2 + z1 * b + z0
02248      *   where
02249      *     z2 = xh * yh
02250      *     z0 = xl * yl
02251      *     z1 = (xh + xl) * (yh + yl) - z2 - z0
02252      *
02253      *  ref: http://en.wikipedia.org/wiki/Karatsuba_algorithm
02254      */
02255 
02256     /* allocate a result bignum */
02257     z = bignew(xn + yn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
02258     zds = BDIGITS(z);
02259 
02260     /* t1 <- xh * yh */
02261     t1 = bigmul0(xh, yh);
02262     t1n = big_real_len(t1);
02263 
02264     /* copy t1 into high bytes of the result (z2) */
02265     MEMCPY(zds + 2 * n, BDIGITS(t1), BDIGIT, t1n);
02266     for (i = 2 * n + t1n; i < xn + yn; i++) zds[i] = 0;
02267 
02268     if (!BIGZEROP(xl) && !BIGZEROP(yl)) {
02269         /* t2 <- xl * yl */
02270         t2 = bigmul0(xl, yl);
02271         t2n = big_real_len(t2);
02272 
02273         /* copy t2 into low bytes of the result (z0) */
02274         MEMCPY(zds, BDIGITS(t2), BDIGIT, t2n);
02275         for (i = t2n; i < 2 * n; i++) zds[i] = 0;
02276     }
02277     else {
02278         t2 = Qundef;
02279         t2n = 0;
02280 
02281         /* copy 0 into low bytes of the result (z0) */
02282         for (i = 0; i < 2 * n; i++) zds[i] = 0;
02283     }
02284 
02285     /* xh <- xh + xl */
02286     if (RBIGNUM_LEN(xl) > RBIGNUM_LEN(xh)) {
02287         t3 = xl; xl = xh; xh = t3;
02288     }
02289     /* xh has a margin for carry */
02290     bigadd_core(BDIGITS(xh), RBIGNUM_LEN(xh),
02291                 BDIGITS(xl), RBIGNUM_LEN(xl),
02292                 BDIGITS(xh), RBIGNUM_LEN(xh));
02293 
02294     /* yh <- yh + yl */
02295     if (x != y) {
02296         if (RBIGNUM_LEN(yl) > RBIGNUM_LEN(yh)) {
02297             t3 = yl; yl = yh; yh = t3;
02298         }
02299         /* yh has a margin for carry */
02300         bigadd_core(BDIGITS(yh), RBIGNUM_LEN(yh),
02301                     BDIGITS(yl), RBIGNUM_LEN(yl),
02302                     BDIGITS(yh), RBIGNUM_LEN(yh));
02303     }
02304     else yh = xh;
02305 
02306     /* t3 <- xh * yh */
02307     t3 = bigmul0(xh, yh);
02308 
02309     i = xn + yn - n;
02310     /* subtract t1 from t3 */
02311     bigsub_core(BDIGITS(t3), big_real_len(t3), BDIGITS(t1), t1n, BDIGITS(t3), big_real_len(t3));
02312 
02313     /* subtract t2 from t3; t3 is now the middle term of the product */
02314     if (t2 != Qundef) bigsub_core(BDIGITS(t3), big_real_len(t3), BDIGITS(t2), t2n, BDIGITS(t3), big_real_len(t3));
02315 
02316     /* add t3 to middle bytes of the result (z1) */
02317     bigadd_core(zds + n, i, BDIGITS(t3), big_real_len(t3), zds + n, i);
02318 
02319     return z;
02320 }
02321 
02322 static void
02323 biglsh_bang(BDIGIT *xds, long xn, unsigned long shift)
02324 {
02325     long const s1 = shift/BITSPERDIG;
02326     int const s2 = (int)(shift%BITSPERDIG);
02327     int const s3 = BITSPERDIG-s2;
02328     BDIGIT* zds;
02329     BDIGIT num;
02330     long i;
02331     if (s1 >= xn) {
02332         MEMZERO(xds, BDIGIT, xn);
02333         return;
02334     }
02335     zds = xds + xn - 1;
02336     xn -= s1 + 1;
02337     num = xds[xn]<<s2;
02338     do {
02339         *zds-- = num | xds[--xn]>>s3;
02340         num = xds[xn]<<s2;
02341     }
02342     while (xn > 0);
02343     *zds = num;
02344     for (i = s1; i > 0; --i)
02345         *zds-- = 0;
02346 }
02347 
02348 static void
02349 bigrsh_bang(BDIGIT* xds, long xn, unsigned long shift)
02350 {
02351     long s1 = shift/BITSPERDIG;
02352     int s2 = (int)(shift%BITSPERDIG);
02353     int s3 = BITSPERDIG - s2;
02354     int i;
02355     BDIGIT num;
02356     BDIGIT* zds;
02357     if (s1 >= xn) {
02358         MEMZERO(xds, BDIGIT, xn);
02359         return;
02360     }
02361 
02362     i = 0;
02363     zds = xds + s1;
02364     num = *zds++>>s2;
02365     do {
02366         xds[i++] = (BDIGIT)(*zds<<s3) | num;
02367         num = *zds++>>s2;
02368     }
02369     while (i < xn - s1 - 1);
02370     xds[i] = num;
02371     MEMZERO(xds + xn - s1, BDIGIT, s1);
02372 }
02373 
02374 static void
02375 big_split3(VALUE v, long n, volatile VALUE* p0, volatile VALUE* p1, volatile VALUE* p2)
02376 {
02377     VALUE v0, v12, v1, v2;
02378 
02379     big_split(v, n, &v12, &v0);
02380     big_split(v12, n, &v2, &v1);
02381 
02382     *p0 = bigtrunc(v0);
02383     *p1 = bigtrunc(v1);
02384     *p2 = bigtrunc(v2);
02385 }
02386 
02387 static VALUE big_lshift(VALUE, unsigned long);
02388 static VALUE big_rshift(VALUE, unsigned long);
02389 static VALUE bigdivrem(VALUE, VALUE, volatile VALUE*, volatile VALUE*);
02390 
02391 static VALUE
02392 bigmul1_toom3(VALUE x, VALUE y)
02393 {
02394     long n, xn, yn, zn;
02395     VALUE x0, x1, x2, y0, y1, y2;
02396     VALUE u0, u1, u2, u3, u4, v1, v2, v3;
02397     VALUE z0, z1, z2, z3, z4, z, t;
02398     BDIGIT* zds;
02399 
02400     xn = RBIGNUM_LEN(x);
02401     yn = RBIGNUM_LEN(y);
02402     assert(xn <= yn);  /* assume y >= x */
02403 
02404     n = (yn + 2) / 3;
02405     big_split3(x, n, &x0, &x1, &x2);
02406     if (x == y) {
02407         y0 = x0; y1 = x1; y2 = x2;
02408     }
02409     else big_split3(y, n, &y0, &y1, &y2);
02410 
02411     /*
02412      * ref. http://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication
02413      *
02414      * x(b) = x0 * b^0 + x1 * b^1 + x2 * b^2
02415      * y(b) = y0 * b^0 + y1 * b^1 + y2 * b^2
02416      *
02417      * z(b) = x(b) * y(b)
02418      * z(b) = z0 * b^0 + z1 * b^1 + z2 * b^2 + z3 * b^3 + z4 * b^4
02419      * where:
02420      *   z0 = x0 * y0
02421      *   z1 = x0 * y1 + x1 * y0
02422      *   z2 = x0 * y2 + x1 * y1 + x2 * y0
02423      *   z3 = x1 * y2 + x2 * y1
02424      *   z4 = x2 * y2
02425      *
02426      * Toom3 method (a.k.a. Toom-Cook method):
02427      * (Step1) calculating 5 points z(b0), z(b1), z(b2), z(b3), z(b4),
02428      * where:
02429      *   b0 = 0, b1 = 1, b2 = -1, b3 = -2, b4 = inf,
02430      *   z(0)   = x(0)   * y(0)   = x0 * y0
02431      *   z(1)   = x(1)   * y(1)   = (x0 + x1 + x2) * (y0 + y1 + y2)
02432      *   z(-1)  = x(-1)  * y(-1)  = (x0 - x1 + x2) * (y0 - y1 + y2)
02433      *   z(-2)  = x(-2)  * y(-2)  = (x0 - 2 * (x1 - 2 * x2)) * (y0 - 2 * (y1 - 2 * y2))
02434      *   z(inf) = x(inf) * y(inf) = x2 * y2
02435      *
02436      * (Step2) interpolating z0, z1, z2, z3, z4, and z5.
02437      *
02438      * (Step3) Substituting base value into b of the polynomial z(b),
02439      */
02440 
02441     /*
02442      * [Step1] calculating 5 points z(b0), z(b1), z(b2), z(b3), z(b4)
02443      */
02444 
02445     /* u1 <- x0 + x2 */
02446     u1 = bigtrunc(bigadd(x0, x2, 1));
02447 
02448     /* x(-1) : u2 <- u1 - x1 = x0 - x1 + x2 */
02449     u2 = bigtrunc(bigsub(u1, x1));
02450 
02451     /* x(1) : u1 <- u1 + x1 = x0 + x1 + x2 */
02452     u1 = bigtrunc(bigadd(u1, x1, 1));
02453 
02454     /* x(-2) : u3 <- 2 * (u2 + x2) - x0 = x0 - 2 * (x1 - 2 * x2) */
02455     u3 = bigadd(u2, x2, 1);
02456     if (BDIGITS(u3)[RBIGNUM_LEN(u3)-1] & BIGRAD_HALF) {
02457         rb_big_resize(u3, RBIGNUM_LEN(u3) + 1);
02458         BDIGITS(u3)[RBIGNUM_LEN(u3)-1] = 0;
02459     }
02460     biglsh_bang(BDIGITS(u3), RBIGNUM_LEN(u3), 1);
02461     u3 = bigtrunc(bigadd(bigtrunc(u3), x0, 0));
02462 
02463     if (x == y) {
02464         v1 = u1; v2 = u2; v3 = u3;
02465     }
02466     else {
02467         /* v1 <- y0 + y2 */
02468         v1 = bigtrunc(bigadd(y0, y2, 1));
02469 
02470         /* y(-1) : v2 <- v1 - y1 = y0 - y1 + y2 */
02471         v2 = bigtrunc(bigsub(v1, y1));
02472 
02473         /* y(1) : v1 <- v1 + y1 = y0 + y1 + y2 */
02474         v1 = bigtrunc(bigadd(v1, y1, 1));
02475 
02476         /* y(-2) : v3 <- 2 * (v2 + y2) - y0 = y0 - 2 * (y1 - 2 * y2) */
02477         v3 = bigadd(v2, y2, 1);
02478         if (BDIGITS(v3)[RBIGNUM_LEN(v3)-1] & BIGRAD_HALF) {
02479             rb_big_resize(v3, RBIGNUM_LEN(v3) + 1);
02480             BDIGITS(v3)[RBIGNUM_LEN(v3)-1] = 0;
02481         }
02482         biglsh_bang(BDIGITS(v3), RBIGNUM_LEN(v3), 1);
02483         v3 = bigtrunc(bigadd(bigtrunc(v3), y0, 0));
02484     }
02485 
02486     /* z(0) : u0 <- x0 * y0 */
02487     u0 = bigtrunc(bigmul0(x0, y0));
02488 
02489     /* z(1) : u1 <- u1 * v1 */
02490     u1 = bigtrunc(bigmul0(u1, v1));
02491 
02492     /* z(-1) : u2 <- u2 * v2 */
02493     u2 = bigtrunc(bigmul0(u2, v2));
02494 
02495     /* z(-2) : u3 <- u3 * v3 */
02496     u3 = bigtrunc(bigmul0(u3, v3));
02497 
02498     /* z(inf) : u4 <- x2 * y2 */
02499     u4 = bigtrunc(bigmul0(x2, y2));
02500 
02501     /* for GC */
02502     v1 = v2 = v3 = Qnil;
02503 
02504     /*
02505      * [Step2] interpolating z0, z1, z2, z3, z4, and z5.
02506      */
02507 
02508     /* z0 <- z(0) == u0 */
02509     z0 = u0;
02510 
02511     /* z4 <- z(inf) == u4 */
02512     z4 = u4;
02513 
02514     /* z3 <- (z(-2) - z(1)) / 3 == (u3 - u1) / 3 */
02515     z3 = bigadd(u3, u1, 0);
02516     bigdivrem(z3, big_three, &z3, NULL); /* TODO: optimize */
02517     bigtrunc(z3);
02518 
02519     /* z1 <- (z(1) - z(-1)) / 2 == (u1 - u2) / 2 */
02520     z1 = bigtrunc(bigadd(u1, u2, 0));
02521     bigrsh_bang(BDIGITS(z1), RBIGNUM_LEN(z1), 1);
02522 
02523     /* z2 <- z(-1) - z(0) == u2 - u0 */
02524     z2 = bigtrunc(bigadd(u2, u0, 0));
02525 
02526     /* z3 <- (z2 - z3) / 2 + 2 * z(inf) == (z2 - z3) / 2 + 2 * u4 */
02527     z3 = bigtrunc(bigadd(z2, z3, 0));
02528     bigrsh_bang(BDIGITS(z3), RBIGNUM_LEN(z3), 1);
02529     t = big_lshift(u4, 1); /* TODO: combining with next addition */
02530     z3 = bigtrunc(bigadd(z3, t, 1));
02531 
02532     /* z2 <- z2 + z1 - z(inf) == z2 + z1 - u4 */
02533     z2 = bigtrunc(bigadd(z2, z1, 1));
02534     z2 = bigtrunc(bigadd(z2, u4, 0));
02535 
02536     /* z1 <- z1 - z3 */
02537     z1 = bigtrunc(bigadd(z1, z3, 0));
02538 
02539     /*
02540      * [Step3] Substituting base value into b of the polynomial z(b),
02541      */
02542 
02543     zn = 6*n + 1;
02544     z = bignew(zn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
02545     zds = BDIGITS(z);
02546     MEMCPY(zds, BDIGITS(z0), BDIGIT, RBIGNUM_LEN(z0));
02547     MEMZERO(zds + RBIGNUM_LEN(z0), BDIGIT, zn - RBIGNUM_LEN(z0));
02548     bigadd_core(zds +   n, zn -   n, BDIGITS(z1), big_real_len(z1), zds +   n, zn -   n);
02549     bigadd_core(zds + 2*n, zn - 2*n, BDIGITS(z2), big_real_len(z2), zds + 2*n, zn - 2*n);
02550     bigadd_core(zds + 3*n, zn - 3*n, BDIGITS(z3), big_real_len(z3), zds + 3*n, zn - 3*n);
02551     bigadd_core(zds + 4*n, zn - 4*n, BDIGITS(z4), big_real_len(z4), zds + 4*n, zn - 4*n);
02552     z = bignorm(z);
02553 
02554     return bignorm(z);
02555 }
02556 
02557 /* efficient squaring (2 times faster than normal multiplication)
02558  * ref: Handbook of Applied Cryptography, Algorithm 14.16
02559  *      http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf
02560  */
02561 static VALUE
02562 bigsqr_fast(VALUE x)
02563 {
02564     long len = RBIGNUM_LEN(x), i, j;
02565     VALUE z = bignew(2 * len + 1, 1);
02566     BDIGIT *xds = BDIGITS(x), *zds = BDIGITS(z);
02567     BDIGIT_DBL c, v, w;
02568 
02569     for (i = 2 * len + 1; i--; ) zds[i] = 0;
02570     for (i = 0; i < len; i++) {
02571         v = (BDIGIT_DBL)xds[i];
02572         if (!v) continue;
02573         c = (BDIGIT_DBL)zds[i + i] + v * v;
02574         zds[i + i] = BIGLO(c);
02575         c = BIGDN(c);
02576         v *= 2;
02577         for (j = i + 1; j < len; j++) {
02578             w = (BDIGIT_DBL)xds[j];
02579             c += (BDIGIT_DBL)zds[i + j] + BIGLO(v) * w;
02580             zds[i + j] = BIGLO(c);
02581             c = BIGDN(c);
02582             if (BIGDN(v)) c += w;
02583         }
02584         if (c) {
02585             c += (BDIGIT_DBL)zds[i + len];
02586             zds[i + len] = BIGLO(c);
02587             c = BIGDN(c);
02588         }
02589         if (c) zds[i + len + 1] += (BDIGIT)c;
02590     }
02591     return z;
02592 }
02593 
02594 #define KARATSUBA_MUL_DIGITS 70
02595 #define TOOM3_MUL_DIGITS 150
02596 
02597 
02598 /* determine whether a bignum is sparse or not by random sampling */
02599 static inline VALUE
02600 big_sparse_p(VALUE x)
02601 {
02602     long c = 0, n = RBIGNUM_LEN(x);
02603 
02604     if (          BDIGITS(x)[rb_genrand_ulong_limited(n / 2) + n / 4]) c++;
02605     if (c <= 1 && BDIGITS(x)[rb_genrand_ulong_limited(n / 2) + n / 4]) c++;
02606     if (c <= 1 && BDIGITS(x)[rb_genrand_ulong_limited(n / 2) + n / 4]) c++;
02607 
02608     return (c <= 1) ? Qtrue : Qfalse;
02609 }
02610 
02611 static VALUE
02612 bigmul0(VALUE x, VALUE y)
02613 {
02614     long xn, yn;
02615 
02616     xn = RBIGNUM_LEN(x);
02617     yn = RBIGNUM_LEN(y);
02618 
02619     /* make sure that y is longer than x */
02620     if (xn > yn) {
02621         VALUE t;
02622         long tn;
02623         t = x; x = y; y = t;
02624         tn = xn; xn = yn; yn = tn;
02625     }
02626     assert(xn <= yn);
02627 
02628     /* normal multiplication when x is small */
02629     if (xn < KARATSUBA_MUL_DIGITS) {
02630       normal:
02631         if (x == y) return bigsqr_fast(x);
02632         if (xn == 1 && yn == 1) return bigmul1_single(x, y);
02633         return bigmul1_normal(x, y);
02634     }
02635 
02636     /* normal multiplication when x or y is a sparse bignum */
02637     if (big_sparse_p(x)) goto normal;
02638     if (big_sparse_p(y)) return bigmul1_normal(y, x);
02639 
02640     /* balance multiplication by slicing y when x is much smaller than y */
02641     if (2 * xn <= yn) return bigmul1_balance(x, y);
02642 
02643     if (xn < TOOM3_MUL_DIGITS) {
02644         /* multiplication by karatsuba method */
02645         return bigmul1_karatsuba(x, y);
02646     }
02647     else if (3*xn <= 2*(yn + 2))
02648         return bigmul1_balance(x, y);
02649     return bigmul1_toom3(x, y);
02650 }
02651 
02652 /*
02653  *  call-seq:
02654  *     big * other  -> Numeric
02655  *
02656  *  Multiplies big and other, returning the result.
02657  */
02658 
02659 VALUE
02660 rb_big_mul(VALUE x, VALUE y)
02661 {
02662     switch (TYPE(y)) {
02663       case T_FIXNUM:
02664         y = rb_int2big(FIX2LONG(y));
02665         break;
02666 
02667       case T_BIGNUM:
02668         break;
02669 
02670       case T_FLOAT:
02671         return DBL2NUM(rb_big2dbl(x) * RFLOAT_VALUE(y));
02672 
02673       default:
02674         return rb_num_coerce_bin(x, y, '*');
02675     }
02676 
02677     return bignorm(bigmul0(x, y));
02678 }
02679 
02680 struct big_div_struct {
02681     long nx, ny, j, nyzero;
02682     BDIGIT *yds, *zds;
02683     volatile VALUE stop;
02684 };
02685 
02686 static void *
02687 bigdivrem1(void *ptr)
02688 {
02689     struct big_div_struct *bds = (struct big_div_struct*)ptr;
02690     long ny = bds->ny;
02691     long i, j;
02692     BDIGIT *yds = bds->yds, *zds = bds->zds;
02693     BDIGIT_DBL t2;
02694     BDIGIT_DBL_SIGNED num;
02695     BDIGIT q;
02696 
02697     j = bds->j;
02698     do {
02699         if (bds->stop) {
02700             bds->j = j;
02701             return 0;
02702         }
02703         if (zds[j] ==  yds[ny-1]) q = (BDIGIT)BIGRAD-1;
02704         else q = (BDIGIT)((BIGUP(zds[j]) + zds[j-1])/yds[ny-1]);
02705         if (q) {
02706            i = bds->nyzero; num = 0; t2 = 0;
02707             do {                        /* multiply and subtract */
02708                 BDIGIT_DBL ee;
02709                 t2 += (BDIGIT_DBL)yds[i] * q;
02710                 ee = num - BIGLO(t2);
02711                 num = (BDIGIT_DBL)zds[j - ny + i] + ee;
02712                 if (ee) zds[j - ny + i] = BIGLO(num);
02713                 num = BIGDN(num);
02714                 t2 = BIGDN(t2);
02715             } while (++i < ny);
02716             num += zds[j - ny + i] - t2;/* borrow from high digit; don't update */
02717             while (num) {               /* "add back" required */
02718                 i = 0; num = 0; q--;
02719                 do {
02720                     BDIGIT_DBL ee = num + yds[i];
02721                     num = (BDIGIT_DBL)zds[j - ny + i] + ee;
02722                     if (ee) zds[j - ny + i] = BIGLO(num);
02723                     num = BIGDN(num);
02724                 } while (++i < ny);
02725                 num--;
02726             }
02727         }
02728         zds[j] = q;
02729     } while (--j >= ny);
02730     return 0;
02731 }
02732 
02733 static void
02734 rb_big_stop(void *ptr)
02735 {
02736     struct big_div_struct *bds = ptr;
02737     bds->stop = Qtrue;
02738 }
02739 
02740 static VALUE
02741 bigdivrem(VALUE x, VALUE y, volatile VALUE *divp, volatile VALUE *modp)
02742 {
02743     struct big_div_struct bds;
02744     long nx = RBIGNUM_LEN(x), ny = RBIGNUM_LEN(y);
02745     long i, j;
02746     VALUE z, yy, zz;
02747     BDIGIT *xds, *yds, *zds, *tds;
02748     BDIGIT_DBL t2;
02749     BDIGIT dd, q;
02750 
02751     if (BIGZEROP(y)) rb_num_zerodiv();
02752     xds = BDIGITS(x);
02753     yds = BDIGITS(y);
02754     if (nx < ny || (nx == ny && xds[nx - 1] < yds[ny - 1])) {
02755         if (divp) *divp = rb_int2big(0);
02756         if (modp) *modp = x;
02757         return Qnil;
02758     }
02759     if (ny == 1) {
02760         dd = yds[0];
02761         z = rb_big_clone(x);
02762         zds = BDIGITS(z);
02763         t2 = 0; i = nx;
02764         while (i--) {
02765             t2 = BIGUP(t2) + zds[i];
02766             zds[i] = (BDIGIT)(t2 / dd);
02767             t2 %= dd;
02768         }
02769         RBIGNUM_SET_SIGN(z, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
02770         if (modp) {
02771             *modp = rb_uint2big((VALUE)t2);
02772             RBIGNUM_SET_SIGN(*modp, RBIGNUM_SIGN(x));
02773         }
02774         if (divp) *divp = z;
02775         return Qnil;
02776     }
02777 
02778     z = bignew(nx==ny?nx+2:nx+1, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
02779     zds = BDIGITS(z);
02780     if (nx==ny) zds[nx+1] = 0;
02781     while (!yds[ny-1]) ny--;
02782 
02783     dd = 0;
02784     q = yds[ny-1];
02785     while ((q & (BDIGIT)(1UL<<(BITSPERDIG-1))) == 0) {
02786         q <<= 1UL;
02787         dd++;
02788     }
02789     if (dd) {
02790         yy = rb_big_clone(y);
02791         tds = BDIGITS(yy);
02792         j = 0;
02793         t2 = 0;
02794         while (j<ny) {
02795             t2 += (BDIGIT_DBL)yds[j]<<dd;
02796             tds[j++] = BIGLO(t2);
02797             t2 = BIGDN(t2);
02798         }
02799         yds = tds;
02800         RB_GC_GUARD(y) = yy;
02801         j = 0;
02802         t2 = 0;
02803         while (j<nx) {
02804             t2 += (BDIGIT_DBL)xds[j]<<dd;
02805             zds[j++] = BIGLO(t2);
02806             t2 = BIGDN(t2);
02807         }
02808         zds[j] = (BDIGIT)t2;
02809     }
02810     else {
02811         zds[nx] = 0;
02812         j = nx;
02813         while (j--) zds[j] = xds[j];
02814     }
02815 
02816     bds.nx = nx;
02817     bds.ny = ny;
02818     bds.zds = zds;
02819     bds.yds = yds;
02820     bds.stop = Qfalse;
02821     bds.j = nx==ny?nx+1:nx;
02822     for (bds.nyzero = 0; !yds[bds.nyzero]; bds.nyzero++);
02823     if (nx > 10000 || ny > 10000) {
02824       retry:
02825         bds.stop = Qfalse;
02826         rb_thread_call_without_gvl(bigdivrem1, &bds, rb_big_stop, &bds);
02827 
02828         if (bds.stop == Qtrue) {
02829             /* execute trap handler, but exception was not raised. */
02830             goto retry;
02831         }
02832     }
02833     else {
02834         bigdivrem1(&bds);
02835     }
02836 
02837     if (divp) {                 /* move quotient down in z */
02838         *divp = zz = rb_big_clone(z);
02839         zds = BDIGITS(zz);
02840         j = (nx==ny ? nx+2 : nx+1) - ny;
02841         for (i = 0;i < j;i++) zds[i] = zds[i+ny];
02842         if (!zds[i-1]) i--;
02843         RBIGNUM_SET_LEN(zz, i);
02844     }
02845     if (modp) {                 /* normalize remainder */
02846         *modp = zz = rb_big_clone(z);
02847         zds = BDIGITS(zz);
02848         while (ny > 1 && !zds[ny-1]) --ny;
02849         if (dd) {
02850             t2 = 0; i = ny;
02851             while (i--) {
02852                 t2 = (t2 | zds[i]) >> dd;
02853                 q = zds[i];
02854                 zds[i] = BIGLO(t2);
02855                 t2 = BIGUP(q);
02856             }
02857         }
02858         if (!zds[ny-1]) ny--;
02859         RBIGNUM_SET_LEN(zz, ny);
02860         RBIGNUM_SET_SIGN(zz, RBIGNUM_SIGN(x));
02861     }
02862     return z;
02863 }
02864 
02865 static void
02866 bigdivmod(VALUE x, VALUE y, volatile VALUE *divp, volatile VALUE *modp)
02867 {
02868     VALUE mod;
02869 
02870     bigdivrem(x, y, divp, &mod);
02871     if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y) && !BIGZEROP(mod)) {
02872         if (divp) *divp = bigadd(*divp, rb_int2big(1), 0);
02873         if (modp) *modp = bigadd(mod, y, 1);
02874     }
02875     else if (modp) {
02876         *modp = mod;
02877     }
02878 }
02879 
02880 
02881 static VALUE
02882 rb_big_divide(VALUE x, VALUE y, ID op)
02883 {
02884     VALUE z;
02885 
02886     switch (TYPE(y)) {
02887       case T_FIXNUM:
02888         y = rb_int2big(FIX2LONG(y));
02889         break;
02890 
02891       case T_BIGNUM:
02892         break;
02893 
02894       case T_FLOAT:
02895         {
02896             if (op == '/') {
02897                 return DBL2NUM(rb_big2dbl(x) / RFLOAT_VALUE(y));
02898             }
02899             else {
02900                 double dy = RFLOAT_VALUE(y);
02901                 if (dy == 0.0) rb_num_zerodiv();
02902                 return rb_dbl2big(rb_big2dbl(x) / dy);
02903             }
02904         }
02905 
02906       default:
02907         return rb_num_coerce_bin(x, y, op);
02908     }
02909     bigdivmod(x, y, &z, 0);
02910 
02911     return bignorm(z);
02912 }
02913 
02914 /*
02915  *  call-seq:
02916  *     big / other     -> Numeric
02917  *
02918  * Performs division: the class of the resulting object depends on
02919  * the class of <code>numeric</code> and on the magnitude of the
02920  * result.
02921  */
02922 
02923 VALUE
02924 rb_big_div(VALUE x, VALUE y)
02925 {
02926     return rb_big_divide(x, y, '/');
02927 }
02928 
02929 /*
02930  *  call-seq:
02931  *     big.div(other)  -> integer
02932  *
02933  * Performs integer division: returns integer value.
02934  */
02935 
02936 VALUE
02937 rb_big_idiv(VALUE x, VALUE y)
02938 {
02939     return rb_big_divide(x, y, rb_intern("div"));
02940 }
02941 
02942 /*
02943  *  call-seq:
02944  *     big % other         -> Numeric
02945  *     big.modulo(other)   -> Numeric
02946  *
02947  *  Returns big modulo other. See Numeric.divmod for more
02948  *  information.
02949  */
02950 
02951 VALUE
02952 rb_big_modulo(VALUE x, VALUE y)
02953 {
02954     VALUE z;
02955 
02956     switch (TYPE(y)) {
02957       case T_FIXNUM:
02958         y = rb_int2big(FIX2LONG(y));
02959         break;
02960 
02961       case T_BIGNUM:
02962         break;
02963 
02964       default:
02965         return rb_num_coerce_bin(x, y, '%');
02966     }
02967     bigdivmod(x, y, 0, &z);
02968 
02969     return bignorm(z);
02970 }
02971 
02972 /*
02973  *  call-seq:
02974  *     big.remainder(numeric)    -> number
02975  *
02976  *  Returns the remainder after dividing <i>big</i> by <i>numeric</i>.
02977  *
02978  *     -1234567890987654321.remainder(13731)      #=> -6966
02979  *     -1234567890987654321.remainder(13731.24)   #=> -9906.22531493148
02980  */
02981 static VALUE
02982 rb_big_remainder(VALUE x, VALUE y)
02983 {
02984     VALUE z;
02985 
02986     switch (TYPE(y)) {
02987       case T_FIXNUM:
02988         y = rb_int2big(FIX2LONG(y));
02989         break;
02990 
02991       case T_BIGNUM:
02992         break;
02993 
02994       default:
02995         return rb_num_coerce_bin(x, y, rb_intern("remainder"));
02996     }
02997     bigdivrem(x, y, 0, &z);
02998 
02999     return bignorm(z);
03000 }
03001 
03002 /*
03003  *  call-seq:
03004  *     big.divmod(numeric)   -> array
03005  *
03006  *  See <code>Numeric#divmod</code>.
03007  *
03008  */
03009 VALUE
03010 rb_big_divmod(VALUE x, VALUE y)
03011 {
03012     VALUE div, mod;
03013 
03014     switch (TYPE(y)) {
03015       case T_FIXNUM:
03016         y = rb_int2big(FIX2LONG(y));
03017         break;
03018 
03019       case T_BIGNUM:
03020         break;
03021 
03022       default:
03023         return rb_num_coerce_bin(x, y, rb_intern("divmod"));
03024     }
03025     bigdivmod(x, y, &div, &mod);
03026 
03027     return rb_assoc_new(bignorm(div), bignorm(mod));
03028 }
03029 
03030 static int
03031 bdigbitsize(BDIGIT x)
03032 {
03033     int size = 1;
03034     int nb = BITSPERDIG / 2;
03035     BDIGIT bits = (~0 << nb);
03036 
03037     if (!x) return 0;
03038     while (x > 1) {
03039         if (x & bits) {
03040             size += nb;
03041             x >>= nb;
03042         }
03043         x &= ~bits;
03044         nb /= 2;
03045         bits >>= nb;
03046     }
03047 
03048     return size;
03049 }
03050 
03051 static VALUE big_lshift(VALUE, unsigned long);
03052 static VALUE big_rshift(VALUE, unsigned long);
03053 
03054 static VALUE
03055 big_shift(VALUE x, long n)
03056 {
03057     if (n < 0)
03058         return big_lshift(x, (unsigned long)-n);
03059     else if (n > 0)
03060         return big_rshift(x, (unsigned long)n);
03061     return x;
03062 }
03063 
03064 static VALUE
03065 big_fdiv(VALUE x, VALUE y)
03066 {
03067 #define DBL_BIGDIG ((DBL_MANT_DIG + BITSPERDIG) / BITSPERDIG)
03068     VALUE z;
03069     long l, ex, ey;
03070     int i;
03071 
03072     bigtrunc(x);
03073     l = RBIGNUM_LEN(x) - 1;
03074     ex = l * BITSPERDIG;
03075     ex += bdigbitsize(BDIGITS(x)[l]);
03076     ex -= 2 * DBL_BIGDIG * BITSPERDIG;
03077     if (ex) x = big_shift(x, ex);
03078 
03079     switch (TYPE(y)) {
03080       case T_FIXNUM:
03081         y = rb_int2big(FIX2LONG(y));
03082       case T_BIGNUM:
03083         bigtrunc(y);
03084         l = RBIGNUM_LEN(y) - 1;
03085         ey = l * BITSPERDIG;
03086         ey += bdigbitsize(BDIGITS(y)[l]);
03087         ey -= DBL_BIGDIG * BITSPERDIG;
03088         if (ey) y = big_shift(y, ey);
03089         break;
03090       case T_FLOAT:
03091         y = dbl2big(ldexp(frexp(RFLOAT_VALUE(y), &i), DBL_MANT_DIG));
03092         ey = i - DBL_MANT_DIG;
03093         break;
03094       default:
03095         rb_bug("big_fdiv");
03096     }
03097     bigdivrem(x, y, &z, 0);
03098     l = ex - ey;
03099 #if SIZEOF_LONG > SIZEOF_INT
03100     {
03101         /* Visual C++ can't be here */
03102         if (l > INT_MAX) return DBL2NUM(INFINITY);
03103         if (l < INT_MIN) return DBL2NUM(0.0);
03104     }
03105 #endif
03106     return DBL2NUM(ldexp(big2dbl(z), (int)l));
03107 }
03108 
03109 /*
03110  *  call-seq:
03111   *     big.fdiv(numeric) -> float
03112  *
03113  *  Returns the floating point result of dividing <i>big</i> by
03114  *  <i>numeric</i>.
03115  *
03116  *     -1234567890987654321.fdiv(13731)      #=> -89910996357705.5
03117  *     -1234567890987654321.fdiv(13731.24)   #=> -89909424858035.7
03118  *
03119  */
03120 
03121 
03122 VALUE
03123 rb_big_fdiv(VALUE x, VALUE y)
03124 {
03125     double dx, dy;
03126 
03127     dx = big2dbl(x);
03128     switch (TYPE(y)) {
03129       case T_FIXNUM:
03130         dy = (double)FIX2LONG(y);
03131         if (isinf(dx))
03132             return big_fdiv(x, y);
03133         break;
03134 
03135       case T_BIGNUM:
03136         dy = rb_big2dbl(y);
03137         if (isinf(dx) || isinf(dy))
03138             return big_fdiv(x, y);
03139         break;
03140 
03141       case T_FLOAT:
03142         dy = RFLOAT_VALUE(y);
03143         if (isnan(dy))
03144             return y;
03145         if (isinf(dx))
03146             return big_fdiv(x, y);
03147         break;
03148 
03149       default:
03150         return rb_num_coerce_bin(x, y, rb_intern("fdiv"));
03151     }
03152     return DBL2NUM(dx / dy);
03153 }
03154 
03155 static VALUE
03156 bigsqr(VALUE x)
03157 {
03158     return bigtrunc(bigmul0(x, x));
03159 }
03160 
03161 /*
03162  *  call-seq:
03163  *     big ** exponent   -> numeric
03164  *
03165  *  Raises _big_ to the _exponent_ power (which may be an integer, float,
03166  *  or anything that will coerce to a number). The result may be
03167  *  a Fixnum, Bignum, or Float
03168  *
03169  *    123456789 ** 2      #=> 15241578750190521
03170  *    123456789 ** 1.2    #=> 5126464716.09932
03171  *    123456789 ** -2     #=> 6.5610001194102e-17
03172  */
03173 
03174 VALUE
03175 rb_big_pow(VALUE x, VALUE y)
03176 {
03177     double d;
03178     SIGNED_VALUE yy;
03179 
03180     if (y == INT2FIX(0)) return INT2FIX(1);
03181     switch (TYPE(y)) {
03182       case T_FLOAT:
03183         d = RFLOAT_VALUE(y);
03184         if ((!RBIGNUM_SIGN(x) && !BIGZEROP(x)) && d != round(d))
03185             return rb_funcall(rb_complex_raw1(x), rb_intern("**"), 1, y);
03186         break;
03187 
03188       case T_BIGNUM:
03189         rb_warn("in a**b, b may be too big");
03190         d = rb_big2dbl(y);
03191         break;
03192 
03193       case T_FIXNUM:
03194         yy = FIX2LONG(y);
03195 
03196         if (yy < 0)
03197             return rb_funcall(rb_rational_raw1(x), rb_intern("**"), 1, y);
03198         else {
03199             VALUE z = 0;
03200             SIGNED_VALUE mask;
03201             const long xlen = RBIGNUM_LEN(x) - 1;
03202             const long xbits = ffs(RBIGNUM_DIGITS(x)[xlen]) + SIZEOF_BDIGITS*BITSPERDIG*xlen;
03203             const long BIGLEN_LIMIT = BITSPERDIG*1024*1024;
03204 
03205             if ((xbits > BIGLEN_LIMIT) || (xbits * yy > BIGLEN_LIMIT)) {
03206                 rb_warn("in a**b, b may be too big");
03207                 d = (double)yy;
03208                 break;
03209             }
03210             for (mask = FIXNUM_MAX + 1; mask; mask >>= 1) {
03211                 if (z) z = bigsqr(z);
03212                 if (yy & mask) {
03213                     z = z ? bigtrunc(bigmul0(z, x)) : x;
03214                 }
03215             }
03216             return bignorm(z);
03217         }
03218         /* NOTREACHED */
03219         break;
03220 
03221       default:
03222         return rb_num_coerce_bin(x, y, rb_intern("**"));
03223     }
03224     return DBL2NUM(pow(rb_big2dbl(x), d));
03225 }
03226 
03227 static VALUE
03228 bigand_int(VALUE x, long y)
03229 {
03230     VALUE z;
03231     BDIGIT *xds, *zds;
03232     long xn, zn;
03233     long i;
03234     char sign;
03235 
03236     if (y == 0) return INT2FIX(0);
03237     sign = (y > 0);
03238     xds = BDIGITS(x);
03239     zn = xn = RBIGNUM_LEN(x);
03240 #if SIZEOF_BDIGITS == SIZEOF_LONG
03241     if (sign) {
03242         y &= xds[0];
03243         return LONG2NUM(y);
03244     }
03245 #endif
03246 
03247     z = bignew(zn, RBIGNUM_SIGN(x) || sign);
03248     zds = BDIGITS(z);
03249 
03250 #if SIZEOF_BDIGITS == SIZEOF_LONG
03251     i = 1;
03252     zds[0] = xds[0] & y;
03253 #else
03254     {
03255         BDIGIT_DBL num = y;
03256 
03257         for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) {
03258             zds[i] = xds[i] & BIGLO(num);
03259             num = BIGDN(num);
03260         }
03261     }
03262 #endif
03263     while (i < xn) {
03264         zds[i] = sign?0:xds[i];
03265         i++;
03266     }
03267     if (!RBIGNUM_SIGN(z)) get2comp(z);
03268     return bignorm(z);
03269 }
03270 
03271 /*
03272  * call-seq:
03273  *     big & numeric   ->  integer
03274  *
03275  * Performs bitwise +and+ between _big_ and _numeric_.
03276  */
03277 
03278 VALUE
03279 rb_big_and(VALUE xx, VALUE yy)
03280 {
03281     volatile VALUE x, y, z;
03282     BDIGIT *ds1, *ds2, *zds;
03283     long i, l1, l2;
03284     char sign;
03285 
03286     if (!FIXNUM_P(yy) && !RB_TYPE_P(yy, T_BIGNUM)) {
03287         return rb_num_coerce_bit(xx, yy, '&');
03288     }
03289 
03290     x = xx;
03291     y = yy;
03292 
03293     if (!RBIGNUM_SIGN(x)) {
03294         x = rb_big_clone(x);
03295         get2comp(x);
03296     }
03297     if (FIXNUM_P(y)) {
03298         return bigand_int(x, FIX2LONG(y));
03299     }
03300     if (!RBIGNUM_SIGN(y)) {
03301         y = rb_big_clone(y);
03302         get2comp(y);
03303     }
03304     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
03305         l1 = RBIGNUM_LEN(y);
03306         l2 = RBIGNUM_LEN(x);
03307         ds1 = BDIGITS(y);
03308         ds2 = BDIGITS(x);
03309         sign = RBIGNUM_SIGN(y);
03310     }
03311     else {
03312         l1 = RBIGNUM_LEN(x);
03313         l2 = RBIGNUM_LEN(y);
03314         ds1 = BDIGITS(x);
03315         ds2 = BDIGITS(y);
03316         sign = RBIGNUM_SIGN(x);
03317     }
03318     z = bignew(l2, RBIGNUM_SIGN(x) || RBIGNUM_SIGN(y));
03319     zds = BDIGITS(z);
03320 
03321     for (i=0; i<l1; i++) {
03322         zds[i] = ds1[i] & ds2[i];
03323     }
03324     for (; i<l2; i++) {
03325         zds[i] = sign?0:ds2[i];
03326     }
03327     if (!RBIGNUM_SIGN(z)) get2comp(z);
03328     return bignorm(z);
03329 }
03330 
03331 static VALUE
03332 bigor_int(VALUE x, long y)
03333 {
03334     VALUE z;
03335     BDIGIT *xds, *zds;
03336     long xn, zn;
03337     long i;
03338     char sign;
03339 
03340     sign = (y >= 0);
03341     xds = BDIGITS(x);
03342     zn = xn = RBIGNUM_LEN(x);
03343     z = bignew(zn, RBIGNUM_SIGN(x) && sign);
03344     zds = BDIGITS(z);
03345 
03346 #if SIZEOF_BDIGITS == SIZEOF_LONG
03347     i = 1;
03348     zds[0] = xds[0] | y;
03349 #else
03350     {
03351         BDIGIT_DBL num = y;
03352 
03353         for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) {
03354             zds[i] = xds[i] | BIGLO(num);
03355             num = BIGDN(num);
03356         }
03357     }
03358 #endif
03359     while (i < xn) {
03360         zds[i] = sign?xds[i]:(BDIGIT)(BIGRAD-1);
03361         i++;
03362     }
03363     if (!RBIGNUM_SIGN(z)) get2comp(z);
03364     return bignorm(z);
03365 }
03366 
03367 /*
03368  * call-seq:
03369  *     big | numeric   ->  integer
03370  *
03371  * Performs bitwise +or+ between _big_ and _numeric_.
03372  */
03373 
03374 VALUE
03375 rb_big_or(VALUE xx, VALUE yy)
03376 {
03377     volatile VALUE x, y, z;
03378     BDIGIT *ds1, *ds2, *zds;
03379     long i, l1, l2;
03380     char sign;
03381 
03382     if (!FIXNUM_P(yy) && !RB_TYPE_P(yy, T_BIGNUM)) {
03383         return rb_num_coerce_bit(xx, yy, '|');
03384     }
03385 
03386     x = xx;
03387     y = yy;
03388 
03389     if (!RBIGNUM_SIGN(x)) {
03390         x = rb_big_clone(x);
03391         get2comp(x);
03392     }
03393     if (FIXNUM_P(y)) {
03394         return bigor_int(x, FIX2LONG(y));
03395     }
03396     if (!RBIGNUM_SIGN(y)) {
03397         y = rb_big_clone(y);
03398         get2comp(y);
03399     }
03400     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
03401         l1 = RBIGNUM_LEN(y);
03402         l2 = RBIGNUM_LEN(x);
03403         ds1 = BDIGITS(y);
03404         ds2 = BDIGITS(x);
03405         sign = RBIGNUM_SIGN(y);
03406     }
03407     else {
03408         l1 = RBIGNUM_LEN(x);
03409         l2 = RBIGNUM_LEN(y);
03410         ds1 = BDIGITS(x);
03411         ds2 = BDIGITS(y);
03412         sign = RBIGNUM_SIGN(x);
03413     }
03414     z = bignew(l2, RBIGNUM_SIGN(x) && RBIGNUM_SIGN(y));
03415     zds = BDIGITS(z);
03416 
03417     for (i=0; i<l1; i++) {
03418         zds[i] = ds1[i] | ds2[i];
03419     }
03420     for (; i<l2; i++) {
03421         zds[i] = sign?ds2[i]:(BDIGIT)(BIGRAD-1);
03422     }
03423     if (!RBIGNUM_SIGN(z)) get2comp(z);
03424     return bignorm(z);
03425 }
03426 
03427 static VALUE
03428 bigxor_int(VALUE x, long y)
03429 {
03430     VALUE z;
03431     BDIGIT *xds, *zds;
03432     long xn, zn;
03433     long i;
03434     char sign;
03435 
03436     sign = (y >= 0) ? 1 : 0;
03437     xds = BDIGITS(x);
03438     zn = xn = RBIGNUM_LEN(x);
03439     z = bignew(zn, !(RBIGNUM_SIGN(x) ^ sign));
03440     zds = BDIGITS(z);
03441 
03442 #if SIZEOF_BDIGITS == SIZEOF_LONG
03443     i = 1;
03444     zds[0] = xds[0] ^ y;
03445 #else
03446     {
03447         BDIGIT_DBL num = y;
03448 
03449         for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) {
03450             zds[i] = xds[i] ^ BIGLO(num);
03451             num = BIGDN(num);
03452         }
03453     }
03454 #endif
03455     while (i < xn) {
03456         zds[i] = sign?xds[i]:~xds[i];
03457         i++;
03458     }
03459     if (!RBIGNUM_SIGN(z)) get2comp(z);
03460     return bignorm(z);
03461 }
03462 /*
03463  * call-seq:
03464  *     big ^ numeric   ->  integer
03465  *
03466  * Performs bitwise +exclusive or+ between _big_ and _numeric_.
03467  */
03468 
03469 VALUE
03470 rb_big_xor(VALUE xx, VALUE yy)
03471 {
03472     volatile VALUE x, y;
03473     VALUE z;
03474     BDIGIT *ds1, *ds2, *zds;
03475     long i, l1, l2;
03476     char sign;
03477 
03478     if (!FIXNUM_P(yy) && !RB_TYPE_P(yy, T_BIGNUM)) {
03479         return rb_num_coerce_bit(xx, yy, '^');
03480     }
03481 
03482     x = xx;
03483     y = yy;
03484 
03485     if (!RBIGNUM_SIGN(x)) {
03486         x = rb_big_clone(x);
03487         get2comp(x);
03488     }
03489     if (FIXNUM_P(y)) {
03490         return bigxor_int(x, FIX2LONG(y));
03491     }
03492     if (!RBIGNUM_SIGN(y)) {
03493         y = rb_big_clone(y);
03494         get2comp(y);
03495     }
03496     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
03497         l1 = RBIGNUM_LEN(y);
03498         l2 = RBIGNUM_LEN(x);
03499         ds1 = BDIGITS(y);
03500         ds2 = BDIGITS(x);
03501         sign = RBIGNUM_SIGN(y);
03502     }
03503     else {
03504         l1 = RBIGNUM_LEN(x);
03505         l2 = RBIGNUM_LEN(y);
03506         ds1 = BDIGITS(x);
03507         ds2 = BDIGITS(y);
03508         sign = RBIGNUM_SIGN(x);
03509     }
03510     RBIGNUM_SET_SIGN(x, RBIGNUM_SIGN(x)?1:0);
03511     RBIGNUM_SET_SIGN(y, RBIGNUM_SIGN(y)?1:0);
03512     z = bignew(l2, !(RBIGNUM_SIGN(x) ^ RBIGNUM_SIGN(y)));
03513     zds = BDIGITS(z);
03514 
03515     for (i=0; i<l1; i++) {
03516         zds[i] = ds1[i] ^ ds2[i];
03517     }
03518     for (; i<l2; i++) {
03519         zds[i] = sign?ds2[i]:~ds2[i];
03520     }
03521     if (!RBIGNUM_SIGN(z)) get2comp(z);
03522 
03523     return bignorm(z);
03524 }
03525 
03526 static VALUE
03527 check_shiftdown(VALUE y, VALUE x)
03528 {
03529     if (!RBIGNUM_LEN(x)) return INT2FIX(0);
03530     if (RBIGNUM_LEN(y) > SIZEOF_LONG / SIZEOF_BDIGITS) {
03531         return RBIGNUM_SIGN(x) ? INT2FIX(0) : INT2FIX(-1);
03532     }
03533     return Qnil;
03534 }
03535 
03536 /*
03537  * call-seq:
03538  *     big << numeric   ->  integer
03539  *
03540  * Shifts big left _numeric_ positions (right if _numeric_ is negative).
03541  */
03542 
03543 VALUE
03544 rb_big_lshift(VALUE x, VALUE y)
03545 {
03546     long shift;
03547     int neg = 0;
03548 
03549     for (;;) {
03550         if (FIXNUM_P(y)) {
03551             shift = FIX2LONG(y);
03552             if (shift < 0) {
03553                 neg = 1;
03554                 shift = -shift;
03555             }
03556             break;
03557         }
03558         else if (RB_TYPE_P(y, T_BIGNUM)) {
03559             if (!RBIGNUM_SIGN(y)) {
03560                 VALUE t = check_shiftdown(y, x);
03561                 if (!NIL_P(t)) return t;
03562                 neg = 1;
03563             }
03564             shift = big2ulong(y, "long", TRUE);
03565             break;
03566         }
03567         y = rb_to_int(y);
03568     }
03569 
03570     x = neg ? big_rshift(x, shift) : big_lshift(x, shift);
03571     return bignorm(x);
03572 }
03573 
03574 static VALUE
03575 big_lshift(VALUE x, unsigned long shift)
03576 {
03577     BDIGIT *xds, *zds;
03578     long s1 = shift/BITSPERDIG;
03579     int s2 = (int)(shift%BITSPERDIG);
03580     VALUE z;
03581     BDIGIT_DBL num = 0;
03582     long len, i;
03583 
03584     len = RBIGNUM_LEN(x);
03585     z = bignew(len+s1+1, RBIGNUM_SIGN(x));
03586     zds = BDIGITS(z);
03587     for (i=0; i<s1; i++) {
03588         *zds++ = 0;
03589     }
03590     xds = BDIGITS(x);
03591     for (i=0; i<len; i++) {
03592         num = num | (BDIGIT_DBL)*xds++<<s2;
03593         *zds++ = BIGLO(num);
03594         num = BIGDN(num);
03595     }
03596     *zds = BIGLO(num);
03597     return z;
03598 }
03599 
03600 /*
03601  * call-seq:
03602  *     big >> numeric   ->  integer
03603  *
03604  * Shifts big right _numeric_ positions (left if _numeric_ is negative).
03605  */
03606 
03607 VALUE
03608 rb_big_rshift(VALUE x, VALUE y)
03609 {
03610     long shift;
03611     int neg = 0;
03612 
03613     for (;;) {
03614         if (FIXNUM_P(y)) {
03615             shift = FIX2LONG(y);
03616             if (shift < 0) {
03617                 neg = 1;
03618                 shift = -shift;
03619             }
03620             break;
03621         }
03622         else if (RB_TYPE_P(y, T_BIGNUM)) {
03623             if (RBIGNUM_SIGN(y)) {
03624                 VALUE t = check_shiftdown(y, x);
03625                 if (!NIL_P(t)) return t;
03626             }
03627             else {
03628                 neg = 1;
03629             }
03630             shift = big2ulong(y, "long", TRUE);
03631             break;
03632         }
03633         y = rb_to_int(y);
03634     }
03635 
03636     x = neg ? big_lshift(x, shift) : big_rshift(x, shift);
03637     return bignorm(x);
03638 }
03639 
03640 static VALUE
03641 big_rshift(VALUE x, unsigned long shift)
03642 {
03643     BDIGIT *xds, *zds;
03644     long s1 = shift/BITSPERDIG;
03645     int s2 = (int)(shift%BITSPERDIG);
03646     VALUE z;
03647     BDIGIT_DBL num = 0;
03648     long i, j;
03649     volatile VALUE save_x;
03650 
03651     if (s1 > RBIGNUM_LEN(x)) {
03652         if (RBIGNUM_SIGN(x))
03653             return INT2FIX(0);
03654         else
03655             return INT2FIX(-1);
03656     }
03657     if (!RBIGNUM_SIGN(x)) {
03658         x = rb_big_clone(x);
03659         get2comp(x);
03660     }
03661     save_x = x;
03662     xds = BDIGITS(x);
03663     i = RBIGNUM_LEN(x); j = i - s1;
03664     if (j == 0) {
03665         if (RBIGNUM_SIGN(x)) return INT2FIX(0);
03666         else return INT2FIX(-1);
03667     }
03668     z = bignew(j, RBIGNUM_SIGN(x));
03669     if (!RBIGNUM_SIGN(x)) {
03670         num = ((BDIGIT_DBL)~0) << BITSPERDIG;
03671     }
03672     zds = BDIGITS(z);
03673     while (i--, j--) {
03674         num = (num | xds[i]) >> s2;
03675         zds[j] = BIGLO(num);
03676         num = BIGUP(xds[i]);
03677     }
03678     if (!RBIGNUM_SIGN(x)) {
03679         get2comp(z);
03680     }
03681     RB_GC_GUARD(save_x);
03682     return z;
03683 }
03684 
03685 /*
03686  *  call-seq:
03687  *     big[n] -> 0, 1
03688  *
03689  *  Bit Reference---Returns the <em>n</em>th bit in the (assumed) binary
03690  *  representation of <i>big</i>, where <i>big</i>[0] is the least
03691  *  significant bit.
03692  *
03693  *     a = 9**15
03694  *     50.downto(0) do |n|
03695  *       print a[n]
03696  *     end
03697  *
03698  *  <em>produces:</em>
03699  *
03700  *     000101110110100000111000011110010100111100010111001
03701  *
03702  */
03703 
03704 static VALUE
03705 rb_big_aref(VALUE x, VALUE y)
03706 {
03707     BDIGIT *xds;
03708     BDIGIT_DBL num;
03709     VALUE shift;
03710     long i, s1, s2;
03711 
03712     if (RB_TYPE_P(y, T_BIGNUM)) {
03713         if (!RBIGNUM_SIGN(y))
03714             return INT2FIX(0);
03715         bigtrunc(y);
03716         if (RBIGNUM_LEN(y) > DIGSPERLONG) {
03717           out_of_range:
03718             return RBIGNUM_SIGN(x) ? INT2FIX(0) : INT2FIX(1);
03719         }
03720         shift = big2ulong(y, "long", FALSE);
03721     }
03722     else {
03723         i = NUM2LONG(y);
03724         if (i < 0) return INT2FIX(0);
03725         shift = (VALUE)i;
03726     }
03727     s1 = shift/BITSPERDIG;
03728     s2 = shift%BITSPERDIG;
03729 
03730     if (s1 >= RBIGNUM_LEN(x)) goto out_of_range;
03731     if (!RBIGNUM_SIGN(x)) {
03732         xds = BDIGITS(x);
03733         i = 0; num = 1;
03734         while (num += ~xds[i], ++i <= s1) {
03735             num = BIGDN(num);
03736         }
03737     }
03738     else {
03739         num = BDIGITS(x)[s1];
03740     }
03741     if (num & ((BDIGIT_DBL)1<<s2))
03742         return INT2FIX(1);
03743     return INT2FIX(0);
03744 }
03745 
03746 /*
03747  * call-seq:
03748  *   big.hash   -> fixnum
03749  *
03750  * Compute a hash based on the value of _big_.
03751  */
03752 
03753 static VALUE
03754 rb_big_hash(VALUE x)
03755 {
03756     st_index_t hash;
03757 
03758     hash = rb_memhash(BDIGITS(x), sizeof(BDIGIT)*RBIGNUM_LEN(x)) ^ RBIGNUM_SIGN(x);
03759     return INT2FIX(hash);
03760 }
03761 
03762 /*
03763  * MISSING: documentation
03764  */
03765 
03766 static VALUE
03767 rb_big_coerce(VALUE x, VALUE y)
03768 {
03769     if (FIXNUM_P(y)) {
03770         y = rb_int2big(FIX2LONG(y));
03771     }
03772     else if (!RB_TYPE_P(y, T_BIGNUM)) {
03773         rb_raise(rb_eTypeError, "can't coerce %s to Bignum",
03774                  rb_obj_classname(y));
03775     }
03776     return rb_assoc_new(y, x);
03777 }
03778 
03779 /*
03780  *  call-seq:
03781  *     big.abs -> aBignum
03782  *
03783  *  Returns the absolute value of <i>big</i>.
03784  *
03785  *     -1234567890987654321.abs   #=> 1234567890987654321
03786  */
03787 
03788 static VALUE
03789 rb_big_abs(VALUE x)
03790 {
03791     if (!RBIGNUM_SIGN(x)) {
03792         x = rb_big_clone(x);
03793         RBIGNUM_SET_SIGN(x, 1);
03794     }
03795     return x;
03796 }
03797 
03798 /*
03799  *  call-seq:
03800  *     big.size -> integer
03801  *
03802  *  Returns the number of bytes in the machine representation of
03803  *  <i>big</i>.
03804  *
03805  *     (256**10 - 1).size   #=> 12
03806  *     (256**20 - 1).size   #=> 20
03807  *     (256**40 - 1).size   #=> 40
03808  */
03809 
03810 static VALUE
03811 rb_big_size(VALUE big)
03812 {
03813     return LONG2FIX(RBIGNUM_LEN(big)*SIZEOF_BDIGITS);
03814 }
03815 
03816 /*
03817  *  call-seq:
03818  *     big.odd? -> true or false
03819  *
03820  *  Returns <code>true</code> if <i>big</i> is an odd number.
03821  */
03822 
03823 static VALUE
03824 rb_big_odd_p(VALUE num)
03825 {
03826     if (BDIGITS(num)[0] & 1) {
03827         return Qtrue;
03828     }
03829     return Qfalse;
03830 }
03831 
03832 /*
03833  *  call-seq:
03834  *     big.even? -> true or false
03835  *
03836  *  Returns <code>true</code> if <i>big</i> is an even number.
03837  */
03838 
03839 static VALUE
03840 rb_big_even_p(VALUE num)
03841 {
03842     if (BDIGITS(num)[0] & 1) {
03843         return Qfalse;
03844     }
03845     return Qtrue;
03846 }
03847 
03848 /*
03849  *  Bignum objects hold integers outside the range of
03850  *  Fixnum. Bignum objects are created
03851  *  automatically when integer calculations would otherwise overflow a
03852  *  Fixnum. When a calculation involving
03853  *  Bignum objects returns a result that will fit in a
03854  *  Fixnum, the result is automatically converted.
03855  *
03856  *  For the purposes of the bitwise operations and <code>[]</code>, a
03857  *  Bignum is treated as if it were an infinite-length
03858  *  bitstring with 2's complement representation.
03859  *
03860  *  While Fixnum values are immediate, Bignum
03861  *  objects are not---assignment and parameter passing work with
03862  *  references to objects, not the objects themselves.
03863  *
03864  */
03865 
03866 void
03867 Init_Bignum(void)
03868 {
03869     rb_cBignum = rb_define_class("Bignum", rb_cInteger);
03870 
03871     rb_define_method(rb_cBignum, "to_s", rb_big_to_s, -1);
03872     rb_define_alias(rb_cBignum, "inspect", "to_s");
03873     rb_define_method(rb_cBignum, "coerce", rb_big_coerce, 1);
03874     rb_define_method(rb_cBignum, "-@", rb_big_uminus, 0);
03875     rb_define_method(rb_cBignum, "+", rb_big_plus, 1);
03876     rb_define_method(rb_cBignum, "-", rb_big_minus, 1);
03877     rb_define_method(rb_cBignum, "*", rb_big_mul, 1);
03878     rb_define_method(rb_cBignum, "/", rb_big_div, 1);
03879     rb_define_method(rb_cBignum, "%", rb_big_modulo, 1);
03880     rb_define_method(rb_cBignum, "div", rb_big_idiv, 1);
03881     rb_define_method(rb_cBignum, "divmod", rb_big_divmod, 1);
03882     rb_define_method(rb_cBignum, "modulo", rb_big_modulo, 1);
03883     rb_define_method(rb_cBignum, "remainder", rb_big_remainder, 1);
03884     rb_define_method(rb_cBignum, "fdiv", rb_big_fdiv, 1);
03885     rb_define_method(rb_cBignum, "**", rb_big_pow, 1);
03886     rb_define_method(rb_cBignum, "&", rb_big_and, 1);
03887     rb_define_method(rb_cBignum, "|", rb_big_or, 1);
03888     rb_define_method(rb_cBignum, "^", rb_big_xor, 1);
03889     rb_define_method(rb_cBignum, "~", rb_big_neg, 0);
03890     rb_define_method(rb_cBignum, "<<", rb_big_lshift, 1);
03891     rb_define_method(rb_cBignum, ">>", rb_big_rshift, 1);
03892     rb_define_method(rb_cBignum, "[]", rb_big_aref, 1);
03893 
03894     rb_define_method(rb_cBignum, "<=>", rb_big_cmp, 1);
03895     rb_define_method(rb_cBignum, "==", rb_big_eq, 1);
03896     rb_define_method(rb_cBignum, ">", big_gt, 1);
03897     rb_define_method(rb_cBignum, ">=", big_ge, 1);
03898     rb_define_method(rb_cBignum, "<", big_lt, 1);
03899     rb_define_method(rb_cBignum, "<=", big_le, 1);
03900     rb_define_method(rb_cBignum, "===", rb_big_eq, 1);
03901     rb_define_method(rb_cBignum, "eql?", rb_big_eql, 1);
03902     rb_define_method(rb_cBignum, "hash", rb_big_hash, 0);
03903     rb_define_method(rb_cBignum, "to_f", rb_big_to_f, 0);
03904     rb_define_method(rb_cBignum, "abs", rb_big_abs, 0);
03905     rb_define_method(rb_cBignum, "magnitude", rb_big_abs, 0);
03906     rb_define_method(rb_cBignum, "size", rb_big_size, 0);
03907     rb_define_method(rb_cBignum, "odd?", rb_big_odd_p, 0);
03908     rb_define_method(rb_cBignum, "even?", rb_big_even_p, 0);
03909 
03910     power_cache_init();
03911 
03912     big_three = rb_uint2big(3);
03913     rb_gc_register_mark_object(big_three);
03914 }
03915