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