bignum.c

   1 /**********************************************************************
   2
   3   bignum.c -
   4
   5   $Author$
   6   created at: Fri Jun 10 00:48:55 JST 1994
   7
   8   Copyright (C) 1993-2007 Yukihiro Matsumoto
   9
  10 **********************************************************************/
  11
  12 #include "ruby/ruby.h"
  13
  14 #include <math.h>
  15 #include <float.h>
  16 #include <ctype.h>
  17 #ifdef HAVE_IEEEFP_H
  18 #include <ieeefp.h>
  19 #endif
  20
  21 VALUE rb_cBignum;
  22
  23 #if defined __MINGW32__
  24 #define USHORT _USHORT
  25 #endif
  26
  27 #define BDIGITS(x) (RBIGNUM_DIGITS(x))
  28 #define BITSPERDIG (SIZEOF_BDIGITS*CHAR_BIT)
  29 #define BIGRAD ((BDIGIT_DBL)1 << BITSPERDIG)
  30 #define DIGSPERLONG ((unsigned int)(SIZEOF_LONG/SIZEOF_BDIGITS))
  31 #if HAVE_LONG_LONG
  32 # define DIGSPERLL ((unsigned int)(SIZEOF_LONG_LONG/SIZEOF_BDIGITS))
  33 #endif
  34 #define BIGUP(x) ((BDIGIT_DBL)(x) << BITSPERDIG)
  35 #define BIGDN(x) RSHIFT(x,BITSPERDIG)
  36 #define BIGLO(x) ((BDIGIT)((x) & (BIGRAD-1)))
  37 #define BDIGMAX ((BDIGIT)-1)
  38
  39 #define BIGZEROP(x) (RBIGNUM_LEN(x) == 0 || \
  40                      (BDIGITS(x)[0] == 0 && \
  41                       (RBIGNUM_LEN(x) == 1 || bigzero_p(x))))
  42
  43 static int
  44 bigzero_p(VALUE x)
  45 {
  46     long i;
  47     for (i = RBIGNUM_LEN(x) - 1; 0 <= i; i--) {
  48         if (BDIGITS(x)[i]) return 0;
  49     }
  50     return 1;
  51 }
  52
  53 int
  54 rb_cmpint(VALUE val, VALUE a, VALUE b)
  55 {
  56     if (NIL_P(val)) {
  57         rb_cmperr(a, b);
  58     }
  59     if (FIXNUM_P(val)) return FIX2INT(val);
  60     if (TYPE(val) == T_BIGNUM) {
  61         if (BIGZEROP(val)) return 0;
  62         if (RBIGNUM_SIGN(val)) return 1;
  63         return -1;
  64     }
  65     if (RTEST(rb_funcall(val, '>', 1, INT2FIX(0)))) return 1;
  66     if (RTEST(rb_funcall(val, '<', 1, INT2FIX(0)))) return -1;
  67     return 0;
  68 }
  69
  70 #define RBIGNUM_SET_LEN(b,l) \
  71   ((RBASIC(b)->flags & RBIGNUM_EMBED_FLAG) ? \
  72    (RBASIC(b)->flags = (RBASIC(b)->flags & ~RBIGNUM_EMBED_LEN_MASK) | \
  73       ((l) << RBIGNUM_EMBED_LEN_SHIFT)) : \
  74    (RBIGNUM(b)->as.heap.len = (l)))
  75
  76 static void
  77 rb_big_realloc(VALUE big, long len)
  78 {
  79     BDIGIT *ds;
  80     if (RBASIC(big)->flags & RBIGNUM_EMBED_FLAG) {
  81         if (RBIGNUM_EMBED_LEN_MAX < len) {
  82             ds = ALLOC_N(BDIGIT, len);
  83             MEMCPY(ds, RBIGNUM(big)->as.ary, BDIGIT, RBIGNUM_EMBED_LEN_MAX);
  84             RBIGNUM(big)->as.heap.len = RBIGNUM_LEN(big);
  85             RBIGNUM(big)->as.heap.digits = ds;
  86             RBASIC(big)->flags &= ~RBIGNUM_EMBED_FLAG;
  87         }
  88     }
  89     else {
  90         if (len <= RBIGNUM_EMBED_LEN_MAX) {
  91             ds = RBIGNUM(big)->as.heap.digits;
  92             RBASIC(big)->flags |= RBIGNUM_EMBED_FLAG;
  93             RBIGNUM_SET_LEN(big, len);
  94             if (ds) {
  95                 MEMCPY(RBIGNUM(big)->as.ary, ds, BDIGIT, len);
  96                 xfree(ds);
  97             }
  98         }
  99         else {
 100             if (RBIGNUM_LEN(big) == 0) {
 101                 RBIGNUM(big)->as.heap.digits = ALLOC_N(BDIGIT, len);
 102             }
 103             else {
 104                 REALLOC_N(RBIGNUM(big)->as.heap.digits, BDIGIT, len);
 105             }
 106         }
 107     }
 108 }
 109
 110 void
 111 rb_big_resize(VALUE big, long len)
 112 {
 113     rb_big_realloc(big, len);
 114     RBIGNUM_SET_LEN(big, len);
 115 }
 116
 117 static VALUE
 118 bignew_1(VALUE klass, long len, int sign)
 119 {
 120     NEWOBJ(big, struct RBignum);
 121     OBJSETUP(big, klass, T_BIGNUM);
 122     RBIGNUM_SET_SIGN(big, sign?1:0);
 123     if (len <= RBIGNUM_EMBED_LEN_MAX) {
 124         RBASIC(big)->flags |= RBIGNUM_EMBED_FLAG;
 125         RBIGNUM_SET_LEN(big, len);
 126     }
 127     else {
 128         rb_big_resize((VALUE)big, len);
 129     }
 130
 131     return (VALUE)big;
 132 }
 133
 134 #define bignew(len,sign) bignew_1(rb_cBignum,len,sign)
 135
 136 VALUE
 137 rb_big_clone(VALUE x)
 138 {
 139     VALUE z = bignew_1(CLASS_OF(x), RBIGNUM_LEN(x), RBIGNUM_SIGN(x));
 140
 141     MEMCPY(BDIGITS(z), BDIGITS(x), BDIGIT, RBIGNUM_LEN(x));
 142     return z;
 143 }
 144
 145 /* modify a bignum by 2's complement */
 146 static void
 147 get2comp(VALUE x)
 148 {
 149     long i = RBIGNUM_LEN(x);
 150     BDIGIT *ds = BDIGITS(x);
 151     BDIGIT_DBL num;
 152
 153     if (!i) return;
 154     while (i--) ds[i] = ~ds[i];
 155     i = 0; num = 1;
 156     do {
 157         num += ds[i];
 158         ds[i++] = BIGLO(num);
 159         num = BIGDN(num);
 160     } while (i < RBIGNUM_LEN(x));
 161     if (num != 0) {
 162         rb_big_resize(x, RBIGNUM_LEN(x)+1);
 163         ds = BDIGITS(x);
 164         ds[RBIGNUM_LEN(x)-1] = 1;
 165     }
 166 }
 167
 168 void
 169 rb_big_2comp(VALUE x)                   /* get 2's complement */
 170 {
 171     get2comp(x);
 172 }
 173
 174 static VALUE
 175 bigtrunc(VALUE x)
 176 {
 177     long len = RBIGNUM_LEN(x);
 178     BDIGIT *ds = BDIGITS(x);
 179
 180     if (len == 0) return x;
 181     while (--len && !ds[len]);
 182     rb_big_resize(x, len+1);
 183     return x;
 184 }
 185
 186 static VALUE
 187 bigfixize(VALUE x)
 188 {
 189     long len = RBIGNUM_LEN(x);
 190     BDIGIT *ds = BDIGITS(x);
 191
 192     if (len*SIZEOF_BDIGITS <= sizeof(long)) {
 193         long num = 0;
 194         while (len--) {
 195             num = BIGUP(num) + ds[len];
 196         }
 197         if (num >= 0) {
 198             if (RBIGNUM_SIGN(x)) {
 199                 if (POSFIXABLE(num)) return LONG2FIX(num);
 200             }
 201             else {
 202                 if (NEGFIXABLE(-(long)num)) return LONG2FIX(-(long)num);
 203             }
 204         }
 205     }
 206     return x;
 207 }
 208
 209 static VALUE
 210 bignorm(VALUE x)
 211 {
 212     if (!FIXNUM_P(x) && TYPE(x) == T_BIGNUM) {
 213         x = bigfixize(bigtrunc(x));
 214     }
 215     return x;
 216 }
 217
 218 VALUE
 219 rb_big_norm(VALUE x)
 220 {
 221     return bignorm(x);
 222 }
 223
 224 VALUE
 225 rb_uint2big(VALUE n)
 226 {
 227     BDIGIT_DBL num = n;
 228     long i = 0;
 229     BDIGIT *digits;
 230     VALUE big;
 231
 232     big = bignew(DIGSPERLONG, 1);
 233     digits = BDIGITS(big);
 234     while (i < DIGSPERLONG) {
 235         digits[i++] = BIGLO(num);
 236         num = BIGDN(num);
 237     }
 238
 239     i = DIGSPERLONG;
 240     while (--i && !digits[i]) ;
 241     RBIGNUM_SET_LEN(big, i+1);
 242     return big;
 243 }
 244
 245 VALUE
 246 rb_int2big(SIGNED_VALUE n)
 247 {
 248     long neg = 0;
 249     VALUE big;
 250
 251     if (n < 0) {
 252         n = -n;
 253         neg = 1;
 254     }
 255     big = rb_uint2big(n);
 256     if (neg) {
 257         RBIGNUM_SET_SIGN(big, 0);
 258     }
 259     return big;
 260 }
 261
 262 VALUE
 263 rb_uint2inum(VALUE n)
 264 {
 265     if (POSFIXABLE(n)) return LONG2FIX(n);
 266     return rb_uint2big(n);
 267 }
 268
 269 VALUE
 270 rb_int2inum(SIGNED_VALUE n)
 271 {
 272     if (FIXABLE(n)) return LONG2FIX(n);
 273     return rb_int2big(n);
 274 }
 275
 276 #ifdef HAVE_LONG_LONG
 277
 278 void
 279 rb_quad_pack(char *buf, VALUE val)
 280 {
 281     LONG_LONG q;
 282
 283     val = rb_to_int(val);
 284     if (FIXNUM_P(val)) {
 285         q = FIX2LONG(val);
 286     }
 287     else {
 288         long len = RBIGNUM_LEN(val);
 289         BDIGIT *ds;
 290
 291         if (len > SIZEOF_LONG_LONG/SIZEOF_BDIGITS) {
 292             len = SIZEOF_LONG_LONG/SIZEOF_BDIGITS;
 293         }
 294         ds = BDIGITS(val);
 295         q = 0;
 296         while (len--) {
 297             q = BIGUP(q);
 298             q += ds[len];
 299         }
 300         if (!RBIGNUM_SIGN(val)) q = -q;
 301     }
 302     memcpy(buf, (char*)&q, SIZEOF_LONG_LONG);
 303 }
 304
 305 VALUE
 306 rb_quad_unpack(const char *buf, int sign)
 307 {
 308     unsigned LONG_LONG q;
 309     long neg = 0;
 310     long i;
 311     BDIGIT *digits;
 312     VALUE big;
 313
 314     memcpy(&q, buf, SIZEOF_LONG_LONG);
 315     if (sign) {
 316         if (FIXABLE((LONG_LONG)q)) return LONG2FIX((LONG_LONG)q);
 317         if ((LONG_LONG)q < 0) {
 318             q = -(LONG_LONG)q;
 319             neg = 1;
 320         }
 321     }
 322     else {
 323         if (POSFIXABLE(q)) return LONG2FIX(q);
 324     }
 325
 326     i = 0;
 327     big = bignew(DIGSPERLL, 1);
 328     digits = BDIGITS(big);
 329     while (i < DIGSPERLL) {
 330         digits[i++] = BIGLO(q);
 331         q = BIGDN(q);
 332     }
 333
 334     i = DIGSPERLL;
 335     while (i-- && !digits[i]) ;
 336     RBIGNUM_SET_LEN(big, i+1);
 337
 338     if (neg) {
 339         RBIGNUM_SET_SIGN(big, 0);
 340     }
 341     return bignorm(big);
 342 }
 343
 344 #else
 345
 346 #define QUAD_SIZE 8
 347
 348 void
 349 rb_quad_pack(char *buf, VALUE val)
 350 {
 351     long len;
 352
 353     memset(buf, 0, QUAD_SIZE);
 354     val = rb_to_int(val);
 355     if (FIXNUM_P(val)) {
 356         val = rb_int2big(FIX2LONG(val));
 357     }
 358     len = RBIGNUM_LEN(val) * SIZEOF_BDIGITS;
 359     if (len > QUAD_SIZE) {
 360         rb_raise(rb_eRangeError, "bignum too big to convert into `quad int'");
 361     }
 362     memcpy(buf, (char*)BDIGITS(val), len);
 363     if (!RBIGNUM_SIGN(val)) {
 364         len = QUAD_SIZE;
 365         while (len--) {
 366             *buf = ~*buf;
 367             buf++;
 368         }
 369     }
 370 }
 371
 372 #define BNEG(b) (RSHIFT(((BDIGIT*)b)[QUAD_SIZE/SIZEOF_BDIGITS-1],BITSPERDIG-1) != 0)
 373
 374 VALUE
 375 rb_quad_unpack(const char *buf, int sign)
 376 {
 377     VALUE big = bignew(QUAD_SIZE/SIZEOF_BDIGITS, 1);
 378
 379     memcpy((char*)BDIGITS(big), buf, QUAD_SIZE);
 380     if (sign && BNEG(buf)) {
 381         long len = QUAD_SIZE;
 382         char *tmp = (char*)BDIGITS(big);
 383
 384         RBIGNUM_SET_SIGN(big, 0);
 385         while (len--) {
 386             *tmp = ~*tmp;
 387             tmp++;
 388         }
 389     }
 390
 391     return bignorm(big);
 392 }
 393
 394 #endif
 395
 396 VALUE
 397 rb_cstr_to_inum(const char *str, int base, int badcheck)
 398 {
 399     const char *s = str;
 400     char *end;
 401     char sign = 1, nondigit = 0;
 402     int c;
 403     BDIGIT_DBL num;
 404     long len, blen = 1;
 405     long i;
 406     VALUE z;
 407     BDIGIT *zds;
 408
 409 #define conv_digit(c) \
 410     (!ISASCII(c) ? -1 : \
 411      ISDIGIT(c) ? ((c) - '0') : \
 412      ISLOWER(c) ? ((c) - 'a' + 10) : \
 413      ISUPPER(c) ? ((c) - 'A' + 10) : \
 414      -1)
 415
 416     if (!str) {
 417         if (badcheck) goto bad;
 418         return INT2FIX(0);
 419     }
 420     while (ISSPACE(*str)) str++;
 421
 422     if (str[0] == '+') {
 423         str++;
 424     }
 425     else if (str[0] == '-') {
 426         str++;
 427         sign = 0;
 428     }
 429     if (str[0] == '+' || str[0] == '-') {
 430         if (badcheck) goto bad;
 431         return INT2FIX(0);
 432     }
 433     if (base <= 0) {
 434         if (str[0] == '0') {
 435             switch (str[1]) {
 436               case 'x': case 'X':
 437                 base = 16;
 438                 break;
 439               case 'b': case 'B':
 440                 base = 2;
 441                 break;
 442               case 'o': case 'O':
 443                 base = 8;
 444                 break;
 445               case 'd': case 'D':
 446                 base = 10;
 447                 break;
 448               default:
 449                 base = 8;
 450             }
 451         }
 452         else if (base < -1) {
 453             base = -base;
 454         }
 455         else {
 456             base = 10;
 457         }
 458     }
 459     switch (base) {
 460       case 2:
 461         len = 1;
 462         if (str[0] == '0' && (str[1] == 'b'||str[1] == 'B')) {
 463             str += 2;
 464         }
 465         break;
 466       case 3:
 467         len = 2;
 468         break;
 469       case 8:
 470         if (str[0] == '0' && (str[1] == 'o'||str[1] == 'O')) {
 471             str += 2;
 472         }
 473       case 4: case 5: case 6: case 7:
 474         len = 3;
 475         break;
 476       case 10:
 477         if (str[0] == '0' && (str[1] == 'd'||str[1] == 'D')) {
 478             str += 2;
 479         }
 480       case 9: case 11: case 12: case 13: case 14: case 15:
 481         len = 4;
 482         break;
 483       case 16:
 484         len = 4;
 485         if (str[0] == '0' && (str[1] == 'x'||str[1] == 'X')) {
 486             str += 2;
 487         }
 488         break;
 489       default:
 490         if (base < 2 || 36 < base) {
 491             rb_raise(rb_eArgError, "invalid radix %d", base);
 492         }
 493         if (base <= 32) {
 494             len = 5;
 495         }
 496         else {
 497             len = 6;
 498         }
 499         break;
 500     }
 501     if (*str == '0') {          /* squeeze preceding 0s */
 502         int us = 0;
 503         while ((c = *++str) == '0' || c == '_') {
 504             if (c == '_') {
 505                 if (++us >= 2)
 506                     break;
 507             } else
 508                 us = 0;
 509         }
 510         if (!(c = *str) || ISSPACE(c)) --str;
 511     }
 512     c = *str;
 513     c = conv_digit(c);
 514     if (c < 0 || c >= base) {
 515         if (badcheck) goto bad;
 516         return INT2FIX(0);
 517     }
 518     len *= strlen(str)*sizeof(char);
 519
 520     if (len <= (sizeof(long)*CHAR_BIT)) {
 521         unsigned long val = STRTOUL(str, &end, base);
 522
 523         if (str < end && *end == '_') goto bigparse;
 524         if (badcheck) {
 525             if (end == str) goto bad; /* no number */
 526             while (*end && ISSPACE(*end)) end++;
 527             if (*end) goto bad;       /* trailing garbage */
 528         }
 529
 530         if (POSFIXABLE(val)) {
 531             if (sign) return LONG2FIX(val);
 532             else {
 533                 long result = -(long)val;
 534                 return LONG2FIX(result);
 535             }
 536         }
 537         else {
 538             VALUE big = rb_uint2big(val);
 539             RBIGNUM_SET_SIGN(big, sign);
 540             return bignorm(big);
 541         }
 542     }
 543   bigparse:
 544     len = (len/BITSPERDIG)+1;
 545     if (badcheck && *str == '_') goto bad;
 546
 547     z = bignew(len, sign);
 548     zds = BDIGITS(z);
 549     for (i=len;i--;) zds[i]=0;
 550     while ((c = *str++) != 0) {
 551         if (c == '_') {
 552             if (nondigit) {
 553                 if (badcheck) goto bad;
 554                 break;
 555             }
 556             nondigit = c;
 557             continue;
 558         }
 559         else if ((c = conv_digit(c)) < 0) {
 560             break;
 561         }
 562         if (c >= base) break;
 563         nondigit = 0;
 564         i = 0;
 565         num = c;
 566         for (;;) {
 567             while (i<blen) {
 568                 num += (BDIGIT_DBL)zds[i]*base;
 569                 zds[i++] = BIGLO(num);
 570                 num = BIGDN(num);
 571             }
 572             if (num) {
 573                 blen++;
 574                 continue;
 575             }
 576             break;
 577         }
 578     }
 579     if (badcheck) {
 580         str--;
 581         if (s+1 < str && str[-1] == '_') goto bad;
 582         while (*str && ISSPACE(*str)) str++;
 583         if (*str) {
 584           bad:
 585             rb_invalid_str(s, "Integer");
 586         }
 587     }
 588
 589     return bignorm(z);
 590 }
 591
 592 VALUE
 593 rb_str_to_inum(VALUE str, int base, int badcheck)
 594 {
 595     char *s;
 596     long len;
 597
 598     StringValue(str);
 599     if (badcheck) {
 600         s = StringValueCStr(str);
 601     }
 602     else {
 603         s = RSTRING_PTR(str);
 604     }
 605     if (s) {
 606         len = RSTRING_LEN(str);
 607         if (s[len]) {           /* no sentinel somehow */
 608             char *p = ALLOCA_N(char, len+1);
 609
 610             MEMCPY(p, s, char, len);
 611             p[len] = '\0';
 612             s = p;
 613         }
 614     }
 615     return rb_cstr_to_inum(s, base, badcheck);
 616 }
 617
 618 #if HAVE_LONG_LONG
 619
 620 static VALUE
 621 rb_ull2big(unsigned LONG_LONG n)
 622 {
 623     BDIGIT_DBL num = n;
 624     long i = 0;
 625     BDIGIT *digits;
 626     VALUE big;
 627
 628     big = bignew(DIGSPERLL, 1);
 629     digits = BDIGITS(big);
 630     while (i < DIGSPERLL) {
 631         digits[i++] = BIGLO(num);
 632         num = BIGDN(num);
 633     }
 634
 635     i = DIGSPERLL;
 636     while (i-- && !digits[i]) ;
 637     RBIGNUM_SET_LEN(big, i+1);
 638     return big;
 639 }
 640
 641 static VALUE
 642 rb_ll2big(LONG_LONG n)
 643 {
 644     long neg = 0;
 645     VALUE big;
 646
 647     if (n < 0) {
 648         n = -n;
 649         neg = 1;
 650     }
 651     big = rb_ull2big(n);
 652     if (neg) {
 653         RBIGNUM_SET_SIGN(big, 0);
 654     }
 655     return big;
 656 }
 657
 658 VALUE
 659 rb_ull2inum(unsigned LONG_LONG n)
 660 {
 661     if (POSFIXABLE(n)) return LONG2FIX(n);
 662     return rb_ull2big(n);
 663 }
 664
 665 VALUE
 666 rb_ll2inum(LONG_LONG n)
 667 {
 668     if (FIXABLE(n)) return LONG2FIX(n);
 669     return rb_ll2big(n);
 670 }
 671
 672 #endif  /* HAVE_LONG_LONG */
 673
 674 VALUE
 675 rb_cstr2inum(const char *str, int base)
 676 {
 677     return rb_cstr_to_inum(str, base, base==0);
 678 }
 679
 680 VALUE
 681 rb_str2inum(VALUE str, int base)
 682 {
 683     return rb_str_to_inum(str, base, base==0);
 684 }
 685
 686 const char ruby_digitmap[] = "0123456789abcdefghijklmnopqrstuvwxyz";
 687
 688 static VALUE bigsqr(VALUE x);
 689 static void bigdivmod(VALUE x, VALUE y, VALUE *divp, VALUE *modp);
 690
 691 #define POW2_P(x) (((x)&((x)-1))==0)
 692
 693 static inline int
 694 ones(register unsigned long x)
 695 {
 696 #if SIZEOF_LONG == 8
 697 # define MASK_55 0x5555555555555555UL
 698 # define MASK_33 0x3333333333333333UL
 699 # define MASK_0f 0x0f0f0f0f0f0f0f0fUL
 700 #else
 701 # define MASK_55 0x55555555UL
 702 # define MASK_33 0x33333333UL
 703 # define MASK_0f 0x0f0f0f0fUL
 704 #endif
 705     x -= (x >> 1) & MASK_55;
 706     x = ((x >> 2) & MASK_33) + (x & MASK_33);
 707     x = ((x >> 4) + x) & MASK_0f;
 708     x += (x >> 8);
 709     x += (x >> 16);
 710 #if SIZEOF_LONG == 8
 711     x += (x >> 32);
 712 #endif
 713     return (int)(x & 0x7f);
 714 #undef MASK_0f
 715 #undef MASK_33
 716 #undef MASK_55
 717 }
 718
 719 static inline unsigned long
 720 next_pow2(register unsigned long x)
 721 {
 722     x |= x >> 1;
 723     x |= x >> 2;
 724     x |= x >> 4;
 725     x |= x >> 8;
 726     x |= x >> 16;
 727 #if SIZEOF_LONG == 8
 728     x |= x >> 32;
 729 #endif
 730     return x + 1;
 731 }
 732
 733 static inline int
 734 floor_log2(register unsigned long x)
 735 {
 736     x |= x >> 1;
 737     x |= x >> 2;
 738     x |= x >> 4;
 739     x |= x >> 8;
 740     x |= x >> 16;
 741 #if SIZEOF_LONG == 8
 742     x |= x >> 32;
 743 #endif
 744     return (int)ones(x) - 1;
 745 }
 746
 747 static inline int
 748 ceil_log2(register unsigned long x)
 749 {
 750     return floor_log2(x) + !POW2_P(x);
 751 }
 752
 753 #define LOG2_KARATSUBA_DIGITS 7
 754 #define KARATSUBA_DIGITS (1L<<LOG2_KARATSUBA_DIGITS)
 755 #define MAX_BIG2STR_TABLE_ENTRIES 64
 756
 757 static VALUE big2str_power_cache[35][MAX_BIG2STR_TABLE_ENTRIES];
 758
 759 static void
 760 power_cache_init(void)
 761 {
 762     int i, j;
 763     for (i = 0; i < 35; ++i) {
 764         for (j = 0; j < MAX_BIG2STR_TABLE_ENTRIES; ++j) {
 765             big2str_power_cache[i][j] = Qnil;
 766         }
 767     }
 768 }
 769
 770 static inline VALUE
 771 power_cache_get_power0(int base, int i)
 772 {
 773     if (NIL_P(big2str_power_cache[base - 2][i])) {
 774         big2str_power_cache[base - 2][i] =
 775             i == 0 ? rb_big_pow(rb_int2big(base), INT2FIX(KARATSUBA_DIGITS))
 776                    : bigsqr(power_cache_get_power0(base, i - 1));
 777         rb_global_variable(&big2str_power_cache[base - 2][i]);
 778     }
 779     return big2str_power_cache[base - 2][i];
 780 }
 781
 782 static VALUE
 783 power_cache_get_power(int base, long n1, long* m1)
 784 {
 785     long i, j, m;
 786     VALUE t;
 787
 788     if (n1 <= KARATSUBA_DIGITS)
 789         rb_bug("n1 > KARATSUBA_DIGITS");
 790
 791     m = ceil_log2(n1);
 792     if (m1) *m1 = 1 << m;
 793     i = m - LOG2_KARATSUBA_DIGITS;
 794     if (i >= MAX_BIG2STR_TABLE_ENTRIES)
 795         i = MAX_BIG2STR_TABLE_ENTRIES - 1;
 796     t = power_cache_get_power0(base, i);
 797
 798     j = KARATSUBA_DIGITS*(1 << i);
 799     while (n1 > j) {
 800         t = bigsqr(t);
 801         j *= 2;
 802     }
 803     return t;
 804 }
 805
 806 /* big2str_muraken_find_n1
 807  *
 808  * Let a natural number x is given by:
 809  * x = 2^0 * x_0 + 2^1 * x_1 + ... + 2^(B*n_0 - 1) * x_{B*n_0 - 1},
 810  * where B is BITSPERDIG (i.e. BDIGITS*CHAR_BIT) and n_0 is
 811  * RBIGNUM_LEN(x).
 812  *
 813  * Now, we assume n_1 = min_n \{ n | 2^(B*n_0/2) <= b_1^(n_1) \}, so
 814  * it is realized that 2^(B*n_0) <= {b_1}^{2*n_1}, where b_1 is a
 815  * given radix number. And then, we have n_1 <= (B*n_0) /
 816  * (2*log_2(b_1)), therefore n_1 is given by ceil((B*n_0) /
 817  * (2*log_2(b_1))).
 818  */
 819 static long
 820 big2str_find_n1(VALUE x, int base)
 821 {
 822     static const double log_2[] = {
 823         1.0,              1.58496250072116, 2.0,
 824         2.32192809488736, 2.58496250072116, 2.8073549220576,
 825         3.0,              3.16992500144231, 3.32192809488736,
 826         3.4594316186373,  3.58496250072116, 3.70043971814109,
 827         3.8073549220576,  3.90689059560852, 4.0,
 828         4.08746284125034, 4.16992500144231, 4.24792751344359,
 829         4.32192809488736, 4.39231742277876, 4.4594316186373,
 830         4.52356195605701, 4.58496250072116, 4.64385618977472,
 831         4.70043971814109, 4.75488750216347, 4.8073549220576,
 832         4.85798099512757, 4.90689059560852, 4.95419631038688,
 833         5.0,              5.04439411935845, 5.08746284125034,
 834         5.12928301694497, 5.16992500144231
 835     };
 836     long bits;
 837
 838     if (base < 2 || 36 < base)
 839         rb_bug("invalid radix %d", base);
 840
 841     if (FIXNUM_P(x)) {
 842         bits = (SIZEOF_LONG*CHAR_BIT - 1)/2 + 1;
 843     }
 844     else if (BIGZEROP(x)) {
 845         return 0;
 846     }
 847     else if (RBIGNUM_LEN(x) >= LONG_MAX/BITSPERDIG) {
 848         rb_raise(rb_eRangeError, "bignum too big to convert into `string'");
 849     }
 850     else {
 851         bits = BITSPERDIG*RBIGNUM_LEN(x);
 852     }
 853
 854     return (long)ceil(bits/log_2[base - 2]);
 855 }
 856
 857 static long
 858 big2str_orig(VALUE x, int base, char* ptr, long len, long hbase, int trim)
 859 {
 860     long i = RBIGNUM_LEN(x), j = len;
 861     BDIGIT* ds = BDIGITS(x);
 862
 863     while (i && j > 0) {
 864         long k = i;
 865         BDIGIT_DBL num = 0;
 866
 867         while (k--) {               /* x / hbase */
 868             num = BIGUP(num) + ds[k];
 869             ds[k] = (BDIGIT)(num / hbase);
 870             num %= hbase;
 871         }
 872         if (trim && ds[i-1] == 0) i--;
 873         k = SIZEOF_BDIGITS;
 874         while (k--) {
 875             ptr[--j] = ruby_digitmap[num % base];
 876             num /= base;
 877             if (j <= 0) break;
 878             if (trim && i == 0 && num == 0) break;
 879         }
 880     }
 881     if (trim) {
 882         while (j < len && ptr[j] == '0') j++;
 883         MEMMOVE(ptr, ptr + j, char, len - j);
 884         len -= j;
 885     }
 886     return len;
 887 }
 888
 889 static long
 890 big2str_karatsuba(VALUE x, int base, char* ptr,
 891                   long n1, long len, long hbase, int trim)
 892 {
 893     long lh, ll, m1;
 894     VALUE b, q, r;
 895
 896     if (FIXNUM_P(x)) {
 897         VALUE str = rb_fix2str(x, base);
 898         char* str_ptr = RSTRING_PTR(str);
 899         long str_len = RSTRING_LEN(str);
 900         if (trim) {
 901             if (FIX2INT(x) == 0) return 0;
 902             MEMCPY(ptr, str_ptr, char, str_len);
 903             return str_len;
 904         }
 905         else {
 906             memset(ptr, '0', len - str_len);
 907             MEMCPY(ptr + len - str_len, str_ptr, char, str_len);
 908             return len;
 909         }
 910     }
 911     if (BIGZEROP(x)) {
 912         if (trim) return 0;
 913         else {
 914             memset(ptr, '0', len);
 915             return len;
 916         }
 917     }
 918
 919     if (n1 <= KARATSUBA_DIGITS) {
 920         return big2str_orig(x, base, ptr, len, hbase, trim);
 921     }
 922
 923     b = power_cache_get_power(base, n1, &m1);
 924     bigdivmod(x, b, &q, &r);
 925     lh = big2str_karatsuba(q, base, ptr,      (len - m1)/2,
 926                            len - m1, hbase, trim);
 927     ll = big2str_karatsuba(r, base, ptr + lh, m1/2,
 928                            m1,       hbase, !lh && trim);
 929
 930     return lh + ll;
 931 }
 932
 933 VALUE
 934 rb_big2str0(VALUE x, int base, int trim)
 935 {
 936     int off;
 937     VALUE ss, xx;
 938     long n1, n2, len, hbase;
 939     char* ptr;
 940
 941     if (FIXNUM_P(x)) {
 942         return rb_fix2str(x, base);
 943     }
 944     if (BIGZEROP(x)) {
 945         return rb_usascii_str_new2("0");
 946     }
 947
 948     if (base < 2 || 36 < base)
 949         rb_raise(rb_eArgError, "invalid radix %d", base);
 950
 951     n2 = big2str_find_n1(x, base);
 952     n1 = (n2 + 1) / 2;
 953     ss = rb_usascii_str_new(0, n2 + 1); /* plus one for sign */
 954     ptr = RSTRING_PTR(ss);
 955     ptr[0] = RBIGNUM_SIGN(x) ? '+' : '-';
 956
 957     hbase = base*base;
 958 #if SIZEOF_BDIGITS > 2
 959     hbase *= hbase;
 960 #endif
 961     off = !(trim && RBIGNUM_SIGN(x)); /* erase plus sign if trim */
 962     xx = rb_big_clone(x);
 963     RBIGNUM_SET_SIGN(xx, 1);
 964     if (n1 <= KARATSUBA_DIGITS) {
 965         len = off + big2str_orig(xx, base, ptr + off, n2, hbase, trim);
 966     }
 967     else {
 968         len = off + big2str_karatsuba(xx, base, ptr + off, n1,
 969                                       n2, hbase, trim);
 970     }
 971
 972     ptr[len] = '\0';
 973     rb_str_resize(ss, len);
 974
 975     return ss;
 976 }
 977
 978 VALUE
 979 rb_big2str(VALUE x, int base)
 980 {
 981     return rb_big2str0(x, base, 1);
 982 }
 983
 984 /*
 985  *  call-seq:
 986  *     big.to_s(base=10)   =>  string
 987  *
 988  *  Returns a string containing the representation of <i>big</i> radix
 989  *  <i>base</i> (2 through 36).
 990  *
 991  *     12345654321.to_s         #=> "12345654321"
 992  *     12345654321.to_s(2)      #=> "1011011111110110111011110000110001"
 993  *     12345654321.to_s(8)      #=> "133766736061"
 994  *     12345654321.to_s(16)     #=> "2dfdbbc31"
 995  *     78546939656932.to_s(36)  #=> "rubyrules"
 996  */
 997
 998 static VALUE
 999 rb_big_to_s(int argc, VALUE *argv, VALUE x)
1000 {
1001     int base;
1002
1003     if (argc == 0) base = 10;
1004     else {
1005         VALUE b;
1006
1007         rb_scan_args(argc, argv, "01", &b);
1008         base = NUM2INT(b);
1009     }
1010     return rb_big2str(x, base);
1011 }
1012
1013 static VALUE
1014 big2ulong(VALUE x, const char *type, int check)
1015 {
1016     long len = RBIGNUM_LEN(x);
1017     BDIGIT_DBL num;
1018     BDIGIT *ds;
1019
1020     if (len > DIGSPERLONG) {
1021         if (check)
1022             rb_raise(rb_eRangeError, "bignum too big to convert into `%s'", type);
1023         len = DIGSPERLONG;
1024     }
1025     ds = BDIGITS(x);
1026     num = 0;
1027     while (len--) {
1028         num = BIGUP(num);
1029         num += ds[len];
1030     }
1031     return num;
1032 }
1033
1034 VALUE
1035 rb_big2ulong_pack(VALUE x)
1036 {
1037     VALUE num = big2ulong(x, "unsigned long", Qfalse);
1038     if (!RBIGNUM_SIGN(x)) {
1039         return -num;
1040     }
1041     return num;
1042 }
1043
1044 VALUE
1045 rb_big2ulong(VALUE x)
1046 {
1047     VALUE num = big2ulong(x, "unsigned long", Qtrue);
1048
1049     if (!RBIGNUM_SIGN(x)) {
1050         if ((SIGNED_VALUE)num < 0) {
1051             rb_raise(rb_eRangeError, "bignum out of range of unsigned long");
1052         }
1053         return -num;
1054     }
1055     return num;
1056 }
1057
1058 SIGNED_VALUE
1059 rb_big2long(VALUE x)
1060 {
1061     VALUE num = big2ulong(x, "long", Qtrue);
1062
1063     if ((SIGNED_VALUE)num < 0 &&
1064         (RBIGNUM_SIGN(x) || (SIGNED_VALUE)num != LONG_MIN)) {
1065         rb_raise(rb_eRangeError, "bignum too big to convert into `long'");
1066     }
1067     if (!RBIGNUM_SIGN(x)) return -(SIGNED_VALUE)num;
1068     return num;
1069 }
1070
1071 #if HAVE_LONG_LONG
1072
1073 static unsigned LONG_LONG
1074 big2ull(VALUE x, const char *type)
1075 {
1076     long len = RBIGNUM_LEN(x);
1077     BDIGIT_DBL num;
1078     BDIGIT *ds;
1079
1080     if (len > SIZEOF_LONG_LONG/SIZEOF_BDIGITS)
1081         rb_raise(rb_eRangeError, "bignum too big to convert into `%s'", type);
1082     ds = BDIGITS(x);
1083     num = 0;
1084     while (len--) {
1085         num = BIGUP(num);
1086         num += ds[len];
1087     }
1088     return num;
1089 }
1090
1091 unsigned LONG_LONG
1092 rb_big2ull(VALUE x)
1093 {
1094     unsigned LONG_LONG num = big2ull(x, "unsigned long long");
1095
1096     if (!RBIGNUM_SIGN(x)) return -num;
1097     return num;
1098 }
1099
1100 LONG_LONG
1101 rb_big2ll(VALUE x)
1102 {
1103     unsigned LONG_LONG num = big2ull(x, "long long");
1104
1105     if ((LONG_LONG)num < 0 && (RBIGNUM_SIGN(x)
1106                                || (LONG_LONG)num != LLONG_MIN)) {
1107         rb_raise(rb_eRangeError, "bignum too big to convert into `long long'");
1108     }
1109     if (!RBIGNUM_SIGN(x)) return -(LONG_LONG)num;
1110     return num;
1111 }
1112
1113 #endif  /* HAVE_LONG_LONG */
1114
1115 static VALUE
1116 dbl2big(double d)
1117 {
1118     long i = 0;
1119     BDIGIT c;
1120     BDIGIT *digits;
1121     VALUE z;
1122     double u = (d < 0)?-d:d;
1123
1124     if (isinf(d)) {
1125         rb_raise(rb_eFloatDomainError, d < 0 ? "-Infinity" : "Infinity");
1126     }
1127     if (isnan(d)) {
1128         rb_raise(rb_eFloatDomainError, "NaN");
1129     }
1130
1131     while (!POSFIXABLE(u) || 0 != (long)u) {
1132         u /= (double)(BIGRAD);
1133         i++;
1134     }
1135     z = bignew(i, d>=0);
1136     digits = BDIGITS(z);
1137     while (i--) {
1138         u *= BIGRAD;
1139         c = (BDIGIT)u;
1140         u -= c;
1141         digits[i] = c;
1142     }
1143
1144     return z;
1145 }
1146
1147 VALUE
1148 rb_dbl2big(double d)
1149 {
1150     return bignorm(dbl2big(d));
1151 }
1152
1153 static int
1154 nlz(BDIGIT x)
1155 {
1156     BDIGIT y;
1157     int n = BITSPERDIG;
1158 #if BITSPERDIG > 64
1159     y = x >> 64; if (y) {n -= 64; x = y;}
1160 #endif
1161 #if BITSPERDIG > 32
1162     y = x >> 32; if (y) {n -= 32; x = y;}
1163 #endif
1164 #if BITSPERDIG > 16
1165     y = x >> 16; if (y) {n -= 16; x = y;}
1166 #endif
1167     y = x >>  8; if (y) {n -=  8; x = y;}
1168     y = x >>  4; if (y) {n -=  4; x = y;}
1169     y = x >>  2; if (y) {n -=  2; x = y;}
1170     y = x >>  1; if (y) {return n - 2;}
1171     return n - x;
1172 }
1173
1174 static double
1175 big2dbl(VALUE x)
1176 {
1177     double d = 0.0;
1178     long i = RBIGNUM_LEN(x), lo = 0, bits;
1179     BDIGIT *ds = BDIGITS(x), dl;
1180
1181     if (i) {
1182         bits = i * BITSPERDIG - nlz(ds[i-1]);
1183         if (bits > DBL_MANT_DIG+DBL_MAX_EXP) {
1184             d = HUGE_VAL;
1185         }
1186         else {
1187             if (bits > DBL_MANT_DIG+1)
1188                 lo = (bits -= DBL_MANT_DIG+1) / BITSPERDIG;
1189             else
1190                 bits = 0;
1191             while (--i > lo) {
1192                 d = ds[i] + BIGRAD*d;
1193             }
1194             dl = ds[i];
1195             if (bits && (dl & (1UL << (bits %= BITSPERDIG)))) {
1196                 int carry = dl & ~(~0UL << bits);
1197                 if (!carry) {
1198                     while (i-- > 0) {
1199                         if ((carry = ds[i]) != 0) break;
1200                     }
1201                 }
1202                 if (carry) {
1203                     dl &= ~0UL << bits;
1204                     dl += 1UL << bits;
1205                     if (!dl) d += 1;
1206                 }
1207             }
1208             d = dl + BIGRAD*d;
1209             if (lo) d = ldexp(d, lo * BITSPERDIG);
1210         }
1211     }
1212     if (!RBIGNUM_SIGN(x)) d = -d;
1213     return d;
1214 }
1215
1216 double
1217 rb_big2dbl(VALUE x)
1218 {
1219     double d = big2dbl(x);
1220
1221     if (isinf(d)) {
1222         rb_warning("Bignum out of Float range");
1223         d = HUGE_VAL;
1224     }
1225     return d;
1226 }
1227
1228 /*
1229  *  call-seq:
1230  *     big.to_f -> float
1231  *
1232  *  Converts <i>big</i> to a <code>Float</code>. If <i>big</i> doesn't
1233  *  fit in a <code>Float</code>, the result is infinity.
1234  *
1235  */
1236
1237 static VALUE
1238 rb_big_to_f(VALUE x)
1239 {
1240     return DOUBLE2NUM(rb_big2dbl(x));
1241 }
1242
1243 /*
1244  *  call-seq:
1245  *     big <=> numeric   => -1, 0, +1
1246  *
1247  *  Comparison---Returns -1, 0, or +1 depending on whether <i>big</i> is
1248  *  less than, equal to, or greater than <i>numeric</i>. This is the
1249  *  basis for the tests in <code>Comparable</code>.
1250  *
1251  */
1252
1253 VALUE
1254 rb_big_cmp(VALUE x, VALUE y)
1255 {
1256     long xlen = RBIGNUM_LEN(x);
1257
1258     switch (TYPE(y)) {
1259       case T_FIXNUM:
1260         y = rb_int2big(FIX2LONG(y));
1261         break;
1262
1263       case T_BIGNUM:
1264         break;
1265
1266       case T_FLOAT:
1267         return rb_dbl_cmp(rb_big2dbl(x), RFLOAT_VALUE(y));
1268
1269       default:
1270         return rb_num_coerce_cmp(x, y, rb_intern("<=>"));
1271     }
1272
1273     if (RBIGNUM_SIGN(x) > RBIGNUM_SIGN(y)) return INT2FIX(1);
1274     if (RBIGNUM_SIGN(x) < RBIGNUM_SIGN(y)) return INT2FIX(-1);
1275     if (xlen < RBIGNUM_LEN(y))
1276         return (RBIGNUM_SIGN(x)) ? INT2FIX(-1) : INT2FIX(1);
1277     if (xlen > RBIGNUM_LEN(y))
1278         return (RBIGNUM_SIGN(x)) ? INT2FIX(1) : INT2FIX(-1);
1279
1280     while(xlen-- && (BDIGITS(x)[xlen]==BDIGITS(y)[xlen]));
1281     if (-1 == xlen) return INT2FIX(0);
1282     return (BDIGITS(x)[xlen] > BDIGITS(y)[xlen]) ?
1283         (RBIGNUM_SIGN(x) ? INT2FIX(1) : INT2FIX(-1)) :
1284             (RBIGNUM_SIGN(x) ? INT2FIX(-1) : INT2FIX(1));
1285 }
1286
1287 /*
1288  *  call-seq:
1289  *     big == obj  => true or false
1290  *
1291  *  Returns <code>true</code> only if <i>obj</i> has the same value
1292  *  as <i>big</i>. Contrast this with <code>Bignum#eql?</code>, which
1293  *  requires <i>obj</i> to be a <code>Bignum</code>.
1294  *
1295  *     68719476736 == 68719476736.0   #=> true
1296  */
1297
1298 VALUE
1299 rb_big_eq(VALUE x, VALUE y)
1300 {
1301     switch (TYPE(y)) {
1302       case T_FIXNUM:
1303         y = rb_int2big(FIX2LONG(y));
1304         break;
1305       case T_BIGNUM:
1306         break;
1307       case T_FLOAT:
1308         {
1309             volatile double a, b;
1310
1311             a = RFLOAT_VALUE(y);
1312             if (isnan(a)) return Qfalse;
1313             b = rb_big2dbl(x);
1314             return (a == b)?Qtrue:Qfalse;
1315         }
1316       default:
1317         return rb_equal(y, x);
1318     }
1319     if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y)) return Qfalse;
1320     if (RBIGNUM_LEN(x) != RBIGNUM_LEN(y)) return Qfalse;
1321     if (MEMCMP(BDIGITS(x),BDIGITS(y),BDIGIT,RBIGNUM_LEN(y)) != 0) return Qfalse;
1322     return Qtrue;
1323 }
1324
1325 /*
1326  *  call-seq:
1327  *     big.eql?(obj)   => true or false
1328  *
1329  *  Returns <code>true</code> only if <i>obj</i> is a
1330  *  <code>Bignum</code> with the same value as <i>big</i>. Contrast this
1331  *  with <code>Bignum#==</code>, which performs type conversions.
1332  *
1333  *     68719476736.eql?(68719476736.0)   #=> false
1334  */
1335
1336 static VALUE
1337 rb_big_eql(VALUE x, VALUE y)
1338 {
1339     if (TYPE(y) != T_BIGNUM) return Qfalse;
1340     if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y)) return Qfalse;
1341     if (RBIGNUM_LEN(x) != RBIGNUM_LEN(y)) return Qfalse;
1342     if (MEMCMP(BDIGITS(x),BDIGITS(y),BDIGIT,RBIGNUM_LEN(y)) != 0) return Qfalse;
1343     return Qtrue;
1344 }
1345
1346 /*
1347  * call-seq:
1348  *    -big   =>  other_big
1349  *
1350  * Unary minus (returns a new Bignum whose value is 0-big)
1351  */
1352
1353 static VALUE
1354 rb_big_uminus(VALUE x)
1355 {
1356     VALUE z = rb_big_clone(x);
1357
1358     RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(x));
1359
1360     return bignorm(z);
1361 }
1362
1363 /*
1364  * call-seq:
1365  *     ~big  =>  integer
1366  *
1367  * Inverts the bits in big. As Bignums are conceptually infinite
1368  * length, the result acts as if it had an infinite number of one
1369  * bits to the left. In hex representations, this is displayed
1370  * as two periods to the left of the digits.
1371  *
1372  *   sprintf("%X", ~0x1122334455)    #=> "..FEEDDCCBBAA"
1373  */
1374
1375 static VALUE
1376 rb_big_neg(VALUE x)
1377 {
1378     VALUE z = rb_big_clone(x);
1379     BDIGIT *ds;
1380     long i;
1381
1382     if (!RBIGNUM_SIGN(x)) get2comp(z);
1383     ds = BDIGITS(z);
1384     i = RBIGNUM_LEN(x);
1385     if (!i) return INT2FIX(~(SIGNED_VALUE)0);
1386     while (i--) {
1387         ds[i] = ~ds[i];
1388     }
1389     RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(z));
1390     if (RBIGNUM_SIGN(x)) get2comp(z);
1391
1392     return bignorm(z);
1393 }
1394
1395 static VALUE
1396 bigsub(VALUE x, VALUE y)
1397 {
1398     VALUE z = 0;
1399     BDIGIT *zds;
1400     BDIGIT_DBL_SIGNED num;
1401     long i = RBIGNUM_LEN(x);
1402
1403     /* if x is larger than y, swap */
1404     if (RBIGNUM_LEN(x) < RBIGNUM_LEN(y)) {
1405         z = x; x = y; y = z;    /* swap x y */
1406     }
1407     else if (RBIGNUM_LEN(x) == RBIGNUM_LEN(y)) {
1408         while (i > 0) {
1409             i--;
1410             if (BDIGITS(x)[i] > BDIGITS(y)[i]) {
1411                 break;
1412             }
1413             if (BDIGITS(x)[i] < BDIGITS(y)[i]) {
1414                 z = x; x = y; y = z;    /* swap x y */
1415                 break;
1416             }
1417         }
1418     }
1419
1420     z = bignew(RBIGNUM_LEN(x), z==0);
1421     zds = BDIGITS(z);
1422
1423     for (i = 0, num = 0; i < RBIGNUM_LEN(y); i++) {
1424         num += (BDIGIT_DBL_SIGNED)BDIGITS(x)[i] - BDIGITS(y)[i];
1425         zds[i] = BIGLO(num);
1426         num = BIGDN(num);
1427     }
1428     while (num && i < RBIGNUM_LEN(x)) {
1429         num += BDIGITS(x)[i];
1430         zds[i++] = BIGLO(num);
1431         num = BIGDN(num);
1432     }
1433     while (i < RBIGNUM_LEN(x)) {
1434         zds[i] = BDIGITS(x)[i];
1435         i++;
1436     }
1437
1438     return z;
1439 }
1440
1441 static VALUE
1442 bigadd(VALUE x, VALUE y, int sign)
1443 {
1444     VALUE z;
1445     BDIGIT_DBL num;
1446     long i, len;
1447
1448     sign = (sign == RBIGNUM_SIGN(y));
1449     if (RBIGNUM_SIGN(x) != sign) {
1450         if (sign) return bigsub(y, x);
1451         return bigsub(x, y);
1452     }
1453
1454     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
1455         len = RBIGNUM_LEN(x) + 1;
1456         z = x; x = y; y = z;
1457     }
1458     else {
1459         len = RBIGNUM_LEN(y) + 1;
1460     }
1461     z = bignew(len, sign);
1462
1463     len = RBIGNUM_LEN(x);
1464     for (i = 0, num = 0; i < len; i++) {
1465         num += (BDIGIT_DBL)BDIGITS(x)[i] + BDIGITS(y)[i];
1466         BDIGITS(z)[i] = BIGLO(num);
1467         num = BIGDN(num);
1468     }
1469     len = RBIGNUM_LEN(y);
1470     while (num && i < len) {
1471         num += BDIGITS(y)[i];
1472         BDIGITS(z)[i++] = BIGLO(num);
1473         num = BIGDN(num);
1474     }
1475     while (i < len) {
1476         BDIGITS(z)[i] = BDIGITS(y)[i];
1477         i++;
1478     }
1479     BDIGITS(z)[i] = (BDIGIT)num;
1480
1481     return z;
1482 }
1483
1484 /*
1485  *  call-seq:
1486  *     big + other  => Numeric
1487  *
1488  *  Adds big and other, returning the result.
1489  */
1490
1491 VALUE
1492 rb_big_plus(VALUE x, VALUE y)
1493 {
1494     switch (TYPE(y)) {
1495       case T_FIXNUM:
1496         y = rb_int2big(FIX2LONG(y));
1497         /* fall through */
1498       case T_BIGNUM:
1499         return bignorm(bigadd(x, y, 1));
1500
1501       case T_FLOAT:
1502         return DOUBLE2NUM(rb_big2dbl(x) + RFLOAT_VALUE(y));
1503
1504       default:
1505         return rb_num_coerce_bin(x, y, '+');
1506     }
1507 }
1508
1509 /*
1510  *  call-seq:
1511  *     big - other  => Numeric
1512  *
1513  *  Subtracts other from big, returning the result.
1514  */
1515
1516 VALUE
1517 rb_big_minus(VALUE x, VALUE y)
1518 {
1519     switch (TYPE(y)) {
1520       case T_FIXNUM:
1521         y = rb_int2big(FIX2LONG(y));
1522         /* fall through */
1523       case T_BIGNUM:
1524         return bignorm(bigadd(x, y, 0));
1525
1526       case T_FLOAT:
1527         return DOUBLE2NUM(rb_big2dbl(x) - RFLOAT_VALUE(y));
1528
1529       default:
1530         return rb_num_coerce_bin(x, y, '-');
1531     }
1532 }
1533
1534 static void
1535 rb_big_stop(void *ptr)
1536 {
1537     VALUE *stop = (VALUE*)ptr;
1538     *stop = Qtrue;
1539 }
1540
1541 struct big_mul_struct {
1542     VALUE x, y, z, stop;
1543 };
1544
1545 static VALUE
1546 bigmul1(void *ptr)
1547 {
1548     struct big_mul_struct *bms = (struct big_mul_struct*)ptr;
1549     long i, j;
1550     BDIGIT_DBL n = 0;
1551     VALUE x = bms->x, y = bms->y, z = bms->z;
1552     BDIGIT *zds;
1553
1554     j = RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1;
1555     zds = BDIGITS(z);
1556     while (j--) zds[j] = 0;
1557     for (i = 0; i < RBIGNUM_LEN(x); i++) {
1558         BDIGIT_DBL dd;
1559         if (bms->stop) return Qnil;
1560         dd = BDIGITS(x)[i];
1561         if (dd == 0) continue;
1562         n = 0;
1563         for (j = 0; j < RBIGNUM_LEN(y); j++) {
1564             BDIGIT_DBL ee = n + (BDIGIT_DBL)dd * BDIGITS(y)[j];
1565             n = zds[i + j] + ee;
1566             if (ee) zds[i + j] = BIGLO(n);
1567             n = BIGDN(n);
1568         }
1569         if (n) {
1570             zds[i + j] = n;
1571         }
1572     }
1573     return z;
1574 }
1575
1576 static VALUE
1577 rb_big_mul0(VALUE x, VALUE y)
1578 {
1579     struct big_mul_struct bms;
1580     volatile VALUE z;
1581
1582     switch (TYPE(y)) {
1583       case T_FIXNUM:
1584         y = rb_int2big(FIX2LONG(y));
1585         break;
1586
1587       case T_BIGNUM:
1588         break;
1589
1590       case T_FLOAT:
1591         return DOUBLE2NUM(rb_big2dbl(x) * RFLOAT_VALUE(y));
1592
1593       default:
1594         return rb_num_coerce_bin(x, y, '*');
1595     }
1596
1597     bms.x = x;
1598     bms.y = y;
1599     bms.z = bignew(RBIGNUM_LEN(x) + RBIGNUM_LEN(y) + 1, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
1600     bms.stop = Qfalse;
1601
1602     if (RBIGNUM_LEN(x) + RBIGNUM_LEN(y) > 10000) {
1603         z = rb_thread_blocking_region(bigmul1, &bms, rb_big_stop, &bms.stop);
1604     }
1605     else {
1606         z = bigmul1(&bms);
1607     }
1608
1609     return z;
1610 }
1611
1612 /*
1613  *  call-seq:
1614  *     big * other  => Numeric
1615  *
1616  *  Multiplies big and other, returning the result.
1617  */
1618
1619 VALUE
1620 rb_big_mul(VALUE x, VALUE y)
1621 {
1622     return bignorm(rb_big_mul0(x, y));
1623 }
1624
1625 struct big_div_struct {
1626     long nx, ny;
1627     BDIGIT *yds, *zds;
1628     VALUE stop;
1629 };
1630
1631 static VALUE
1632 bigdivrem1(void *ptr)
1633 {
1634     struct big_div_struct *bds = (struct big_div_struct*)ptr;
1635     long nx = bds->nx, ny = bds->ny;
1636     long i, j;
1637     BDIGIT *yds = bds->yds, *zds = bds->zds;
1638     BDIGIT_DBL t2;
1639     BDIGIT_DBL_SIGNED num;
1640     BDIGIT q;
1641
1642     j = nx==ny?nx+1:nx;
1643     do {
1644         if (bds->stop) return Qnil;
1645         if (zds[j] ==  yds[ny-1]) q = BIGRAD-1;
1646         else q = (BDIGIT)((BIGUP(zds[j]) + zds[j-1])/yds[ny-1]);
1647         if (q) {
1648             i = 0; num = 0; t2 = 0;
1649             do {                        /* multiply and subtract */
1650                 BDIGIT_DBL ee;
1651                 t2 += (BDIGIT_DBL)yds[i] * q;
1652                 ee = num - BIGLO(t2);
1653                 num = (BDIGIT_DBL)zds[j - ny + i] + ee;
1654                 if (ee) zds[j - ny + i] = BIGLO(num);
1655                 num = BIGDN(num);
1656                 t2 = BIGDN(t2);
1657             } while (++i < ny);
1658             num += zds[j - ny + i] - t2;/* borrow from high digit; don't update */
1659             while (num) {               /* "add back" required */
1660                 i = 0; num = 0; q--;
1661                 do {
1662                     BDIGIT_DBL ee = num + yds[i];
1663                     num = (BDIGIT_DBL)zds[j - ny + i] + ee;
1664                     if (ee) zds[j - ny + i] = BIGLO(num);
1665                     num = BIGDN(num);
1666                 } while (++i < ny);
1667                 num--;
1668             }
1669         }
1670         zds[j] = q;
1671     } while (--j >= ny);
1672     return Qnil;
1673 }
1674
1675 static VALUE
1676 bigdivrem(VALUE x, VALUE y, VALUE *divp, VALUE *modp)
1677 {
1678     struct big_div_struct bds;
1679     long nx = RBIGNUM_LEN(x), ny = RBIGNUM_LEN(y);
1680     long i, j;
1681     volatile VALUE yy, z;
1682     BDIGIT *xds, *yds, *zds, *tds;
1683     BDIGIT_DBL t2;
1684     BDIGIT dd, q;
1685
1686     if (BIGZEROP(y)) rb_num_zerodiv();
1687     yds = BDIGITS(y);
1688     if (nx < ny || (nx == ny && BDIGITS(x)[nx - 1] < BDIGITS(y)[ny - 1])) {
1689         if (divp) *divp = rb_int2big(0);
1690         if (modp) *modp = x;
1691         return Qnil;
1692     }
1693     xds = BDIGITS(x);
1694     if (ny == 1) {
1695         dd = yds[0];
1696         z = rb_big_clone(x);
1697         zds = BDIGITS(z);
1698         t2 = 0; i = nx;
1699         while (i--) {
1700             t2 = BIGUP(t2) + zds[i];
1701             zds[i] = (BDIGIT)(t2 / dd);
1702             t2 %= dd;
1703         }
1704         RBIGNUM_SET_SIGN(z, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
1705         if (modp) {
1706             *modp = rb_uint2big((VALUE)t2);
1707             RBIGNUM_SET_SIGN(*modp, RBIGNUM_SIGN(x));
1708         }
1709         if (divp) *divp = z;
1710         return Qnil;
1711     }
1712     z = bignew(nx==ny?nx+2:nx+1, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y));
1713     zds = BDIGITS(z);
1714     if (nx==ny) zds[nx+1] = 0;
1715     while (!yds[ny-1]) ny--;
1716
1717     dd = 0;
1718     q = yds[ny-1];
1719     while ((q & (1UL<<(BITSPERDIG-1))) == 0) {
1720         q <<= 1UL;
1721         dd++;
1722     }
1723     if (dd) {
1724         yy = rb_big_clone(y);
1725         tds = BDIGITS(yy);
1726         j = 0;
1727         t2 = 0;
1728         while (j<ny) {
1729             t2 += (BDIGIT_DBL)yds[j]<<dd;
1730             tds[j++] = BIGLO(t2);
1731             t2 = BIGDN(t2);
1732         }
1733         yds = tds;
1734         j = 0;
1735         t2 = 0;
1736         while (j<nx) {
1737             t2 += (BDIGIT_DBL)xds[j]<<dd;
1738             zds[j++] = BIGLO(t2);
1739             t2 = BIGDN(t2);
1740         }
1741         zds[j] = (BDIGIT)t2;
1742     }
1743     else {
1744         zds[nx] = 0;
1745         j = nx;
1746         while (j--) zds[j] = xds[j];
1747     }
1748
1749     bds.nx = nx;
1750     bds.ny = ny;
1751     bds.zds = zds;
1752     bds.yds = yds;
1753     bds.stop = Qfalse;
1754     if (RBIGNUM_LEN(x) > 10000 || RBIGNUM_LEN(y) > 10000) {
1755         rb_thread_blocking_region(bigdivrem1, &bds, rb_big_stop, &bds.stop);
1756     }
1757     else {
1758         bigdivrem1(&bds);
1759     }
1760
1761     if (divp) {                 /* move quotient down in z */
1762         *divp = rb_big_clone(z);
1763         zds = BDIGITS(*divp);
1764         j = (nx==ny ? nx+2 : nx+1) - ny;
1765         for (i = 0;i < j;i++) zds[i] = zds[i+ny];
1766         RBIGNUM_SET_LEN(*divp, i);
1767     }
1768     if (modp) {                 /* normalize remainder */
1769         *modp = rb_big_clone(z);
1770         zds = BDIGITS(*modp);
1771         while (--ny && !zds[ny]); ++ny;
1772         if (dd) {
1773             t2 = 0; i = ny;
1774             while(i--) {
1775                 t2 = (t2 | zds[i]) >> dd;
1776                 q = zds[i];
1777                 zds[i] = BIGLO(t2);
1778                 t2 = BIGUP(q);
1779             }
1780         }
1781         RBIGNUM_SET_LEN(*modp, ny);
1782         RBIGNUM_SET_SIGN(*modp, RBIGNUM_SIGN(x));
1783     }
1784     return z;
1785 }
1786
1787 static void
1788 bigdivmod(VALUE x, VALUE y, VALUE *divp, VALUE *modp)
1789 {
1790     VALUE mod;
1791
1792     bigdivrem(x, y, divp, &mod);
1793     if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y) && !BIGZEROP(mod)) {
1794         if (divp) *divp = bigadd(*divp, rb_int2big(1), 0);
1795         if (modp) *modp = bigadd(mod, y, 1);
1796     }
1797     else {
1798         if (divp) *divp = *divp;
1799         if (modp) *modp = mod;
1800     }
1801 }
1802
1803
1804 static VALUE
1805 rb_big_divide(VALUE x, VALUE y, ID op)
1806 {
1807     VALUE z;
1808
1809     switch (TYPE(y)) {
1810       case T_FIXNUM:
1811         y = rb_int2big(FIX2LONG(y));
1812         break;
1813
1814       case T_BIGNUM:
1815         break;
1816
1817       case T_FLOAT:
1818         {
1819             double div = rb_big2dbl(x) / RFLOAT_VALUE(y);
1820             if (op == '/') {
1821                 return DOUBLE2NUM(div);
1822             }
1823             else {
1824                 return rb_dbl2big(div);
1825             }
1826         }
1827
1828       default:
1829         return rb_num_coerce_bin(x, y, op);
1830     }
1831     bigdivmod(x, y, &z, 0);
1832
1833     return bignorm(z);
1834 }
1835
1836 /*
1837  *  call-seq:
1838  *     big / other     => Numeric
1839  *
1840  *  Divides big by other, returning the result.
1841  */
1842
1843 VALUE
1844 rb_big_div(VALUE x, VALUE y)
1845 {
1846   return rb_big_divide(x, y, '/');
1847 }
1848
1849 VALUE
1850 rb_big_idiv(VALUE x, VALUE y)
1851 {
1852   return rb_big_divide(x, y, rb_intern("div"));
1853 }
1854
1855 /*
1856  *  call-seq:
1857  *     big % other         => Numeric
1858  *     big.modulo(other)   => Numeric
1859  *
1860  *  Returns big modulo other. See Numeric.divmod for more
1861  *  information.
1862  */
1863
1864 VALUE
1865 rb_big_modulo(VALUE x, VALUE y)
1866 {
1867     VALUE z;
1868
1869     switch (TYPE(y)) {
1870       case T_FIXNUM:
1871         y = rb_int2big(FIX2LONG(y));
1872         break;
1873
1874       case T_BIGNUM:
1875         break;
1876
1877       default:
1878         return rb_num_coerce_bin(x, y, '%');
1879     }
1880     bigdivmod(x, y, 0, &z);
1881
1882     return bignorm(z);
1883 }
1884
1885 /*
1886  *  call-seq:
1887  *     big.remainder(numeric)    => number
1888  *
1889  *  Returns the remainder after dividing <i>big</i> by <i>numeric</i>.
1890  *
1891  *     -1234567890987654321.remainder(13731)      #=> -6966
1892  *     -1234567890987654321.remainder(13731.24)   #=> -9906.22531493148
1893  */
1894 static VALUE
1895 rb_big_remainder(VALUE x, VALUE y)
1896 {
1897     VALUE z;
1898
1899     switch (TYPE(y)) {
1900       case T_FIXNUM:
1901         y = rb_int2big(FIX2LONG(y));
1902         break;
1903
1904       case T_BIGNUM:
1905         break;
1906
1907       default:
1908         return rb_num_coerce_bin(x, y, rb_intern("remainder"));
1909     }
1910     bigdivrem(x, y, 0, &z);
1911
1912     return bignorm(z);
1913 }
1914
1915 /*
1916  *  call-seq:
1917  *     big.divmod(numeric)   => array
1918  *
1919  *  See <code>Numeric#divmod</code>.
1920  *
1921  */
1922 VALUE
1923 rb_big_divmod(VALUE x, VALUE y)
1924 {
1925     VALUE div, mod;
1926
1927     switch (TYPE(y)) {
1928       case T_FIXNUM:
1929         y = rb_int2big(FIX2LONG(y));
1930         break;
1931
1932       case T_BIGNUM:
1933         break;
1934
1935       default:
1936         return rb_num_coerce_bin(x, y, rb_intern("divmod"));
1937     }
1938     bigdivmod(x, y, &div, &mod);
1939
1940     return rb_assoc_new(bignorm(div), bignorm(mod));
1941 }
1942
1943 static int
1944 bdigbitsize(BDIGIT x)
1945 {
1946     int size = 1;
1947     int nb = BITSPERDIG / 2;
1948     BDIGIT bits = (~0 << nb);
1949
1950     if (!x) return 0;
1951     while (x > 1) {
1952         if (x & bits) {
1953             size += nb;
1954             x >>= nb;
1955         }
1956         x &= ~bits;
1957         nb /= 2;
1958         bits >>= nb;
1959     }
1960
1961     return size;
1962 }
1963
1964 static VALUE big_lshift(VALUE, unsigned long);
1965 static VALUE big_rshift(VALUE, unsigned long);
1966
1967 static VALUE big_shift(VALUE x, int n)
1968 {
1969     if (n < 0)
1970         return big_lshift(x, (unsigned int)-n);
1971     else if (n > 0)
1972         return big_rshift(x, (unsigned int)n);
1973     return x;
1974 }
1975
1976 /*
1977  *  call-seq:
1978   *     big.fdiv(numeric) -> float
1979  *
1980  *  Returns the floating point result of dividing <i>big</i> by
1981  *  <i>numeric</i>.
1982  *
1983  *     -1234567890987654321.fdiv(13731)      #=> -89910996357705.5
1984  *     -1234567890987654321.fdiv(13731.24)   #=> -89909424858035.7
1985  *
1986  */
1987
1988 static VALUE
1989 rb_big_fdiv(VALUE x, VALUE y)
1990 {
1991     double dx = big2dbl(x);
1992     double dy;
1993
1994     if (isinf(dx)) {
1995 #define DBL_BIGDIG ((DBL_MANT_DIG + BITSPERDIG) / BITSPERDIG)
1996         VALUE z;
1997         int ex, ey;
1998
1999         ex = (RBIGNUM_LEN(bigtrunc(x)) - 1) * BITSPERDIG;
2000         ex += bdigbitsize(BDIGITS(x)[RBIGNUM_LEN(x) - 1]);
2001         ex -= 2 * DBL_BIGDIG * BITSPERDIG;
2002         if (ex) x = big_shift(x, ex);
2003
2004         switch (TYPE(y)) {
2005           case T_FIXNUM:
2006             y = rb_int2big(FIX2LONG(y));
2007           case T_BIGNUM: {
2008             ey = (RBIGNUM_LEN(bigtrunc(y)) - 1) * BITSPERDIG;
2009             ey += bdigbitsize(BDIGITS(y)[RBIGNUM_LEN(y) - 1]);
2010             ey -= DBL_BIGDIG * BITSPERDIG;
2011             if (ey) y = big_shift(y, ey);
2012           bignum:
2013             bigdivrem(x, y, &z, 0);
2014             return DOUBLE2NUM(ldexp(big2dbl(z), ex - ey));
2015           }
2016           case T_FLOAT:
2017             if (isnan(RFLOAT_VALUE(y))) return y;
2018             y = dbl2big(ldexp(frexp(RFLOAT_VALUE(y), &ey), DBL_MANT_DIG));
2019             ey -= DBL_MANT_DIG;
2020             goto bignum;
2021         }
2022     }
2023     switch (TYPE(y)) {
2024       case T_FIXNUM:
2025         dy = (double)FIX2LONG(y);
2026         break;
2027
2028       case T_BIGNUM:
2029         dy = rb_big2dbl(y);
2030         break;
2031
2032       case T_FLOAT:
2033         dy = RFLOAT_VALUE(y);
2034         break;
2035
2036       default:
2037         return rb_num_coerce_bin(x, y, rb_intern("fdiv"));
2038     }
2039     return DOUBLE2NUM(dx / dy);
2040 }
2041
2042 static VALUE
2043 bigsqr(VALUE x)
2044 {
2045     long len = RBIGNUM_LEN(x), k = len / 2, i;
2046     VALUE a, b, a2, z;
2047     BDIGIT_DBL num;
2048
2049     if (len < 4000 / BITSPERDIG) {
2050         return bigtrunc(rb_big_mul0(x, x));
2051     }
2052
2053     a = bignew(len - k, 1);
2054     MEMCPY(BDIGITS(a), BDIGITS(x) + k, BDIGIT, len - k);
2055     b = bignew(k, 1);
2056     MEMCPY(BDIGITS(b), BDIGITS(x), BDIGIT, k);
2057
2058     a2 = bigtrunc(bigsqr(a));
2059     z = bigsqr(b);
2060     rb_big_realloc(z, (len = 2 * k + RBIGNUM_LEN(a2)) + 1);
2061     while (RBIGNUM_LEN(z) < 2 * k) {
2062         BDIGITS(z)[RBIGNUM_LEN(z)] = 0;
2063         RBIGNUM_SET_LEN(z, RBIGNUM_LEN(z)+1);
2064     }
2065     MEMCPY(BDIGITS(z) + 2 * k, BDIGITS(a2), BDIGIT, RBIGNUM_LEN(a2));
2066     RBIGNUM_SET_LEN(z, len);
2067     a2 = bigtrunc(rb_big_mul0(a, b));
2068     len = RBIGNUM_LEN(a2);
2069     for (i = 0, num = 0; i < len; i++) {
2070         num += (BDIGIT_DBL)BDIGITS(z)[i + k] + ((BDIGIT_DBL)BDIGITS(a2)[i] << 1);
2071         BDIGITS(z)[i + k] = BIGLO(num);
2072         num = BIGDN(num);
2073     }
2074     if (num) {
2075         len = RBIGNUM_LEN(z);
2076         for (i += k; i < len && num; ++i) {
2077             num += (BDIGIT_DBL)BDIGITS(z)[i];
2078             BDIGITS(z)[i] = BIGLO(num);
2079             num = BIGDN(num);
2080         }
2081         if (num) {
2082             BDIGITS(z)[RBIGNUM_LEN(z)] = BIGLO(num);
2083             RBIGNUM_SET_LEN(z, RBIGNUM_LEN(z)+1);
2084         }
2085     }
2086     return bigtrunc(z);
2087 }
2088
2089 /*
2090  *  call-seq:
2091  *     big ** exponent   => numeric
2092  *
2093  *  Raises _big_ to the _exponent_ power (which may be an integer, float,
2094  *  or anything that will coerce to a number). The result may be
2095  *  a Fixnum, Bignum, or Float
2096  *
2097  *    123456789 ** 2      #=> 15241578750190521
2098  *    123456789 ** 1.2    #=> 5126464716.09932
2099  *    123456789 ** -2     #=> 6.5610001194102e-17
2100  */
2101
2102 VALUE
2103 rb_big_pow(VALUE x, VALUE y)
2104 {
2105     double d;
2106     SIGNED_VALUE yy;
2107
2108     if (y == INT2FIX(0)) return INT2FIX(1);
2109     switch (TYPE(y)) {
2110       case T_FLOAT:
2111         d = RFLOAT_VALUE(y);
2112         break;
2113
2114       case T_BIGNUM:
2115         if (rb_funcall(y, '<', 1, INT2FIX(0)))
2116           return rb_funcall(rb_rational_raw1(x), rb_intern("**"), 1, y);
2117
2118         rb_warn("in a**b, b may be too big");
2119         d = rb_big2dbl(y);
2120         break;
2121
2122       case T_FIXNUM:
2123         yy = FIX2LONG(y);
2124
2125         if (yy < 0)
2126           return rb_funcall(rb_rational_raw1(x), rb_intern("**"), 1, y);
2127         else {
2128             VALUE z = 0;
2129             SIGNED_VALUE mask;
2130             const long BIGLEN_LIMIT = 1024*1024 / SIZEOF_BDIGITS;
2131
2132             if ((RBIGNUM_LEN(x) > BIGLEN_LIMIT) ||
2133                 (RBIGNUM_LEN(x) > BIGLEN_LIMIT / yy)) {
2134                 rb_warn("in a**b, b may be too big");
2135                 d = (double)yy;
2136                 break;
2137             }
2138             for (mask = FIXNUM_MAX + 1; mask; mask >>= 1) {
2139                 if (z) z = bigtrunc(bigsqr(z));
2140                 if (yy & mask) {
2141                     z = z ? bigtrunc(rb_big_mul0(z, x)) : x;
2142                 }
2143             }
2144             return bignorm(z);
2145         }
2146         /* NOTREACHED */
2147         break;
2148
2149       default:
2150         return rb_num_coerce_bin(x, y, rb_intern("**"));
2151     }
2152     return DOUBLE2NUM(pow(rb_big2dbl(x), d));
2153 }
2154
2155 static VALUE
2156 bit_coerce(VALUE x)
2157 {
2158     while (!FIXNUM_P(x) && TYPE(x) != T_BIGNUM) {
2159         if (TYPE(x) == T_FLOAT) {
2160             rb_raise(rb_eTypeError, "can't convert Float into Integer");
2161         }
2162         x = rb_to_int(x);
2163     }
2164     return x;
2165 }
2166
2167 /*
2168  * call-seq:
2169  *     big & numeric   =>  integer
2170  *
2171  * Performs bitwise +and+ between _big_ and _numeric_.
2172  */
2173
2174 VALUE
2175 rb_big_and(VALUE xx, VALUE yy)
2176 {
2177     volatile VALUE x, y, z;
2178     BDIGIT *ds1, *ds2, *zds;
2179     long i, l1, l2;
2180     char sign;
2181
2182     x = xx;
2183     y = bit_coerce(yy);
2184     if (FIXNUM_P(y)) {
2185         y = rb_int2big(FIX2LONG(y));
2186     }
2187     if (!RBIGNUM_SIGN(y)) {
2188         y = rb_big_clone(y);
2189         get2comp(y);
2190     }
2191     if (!RBIGNUM_SIGN(x)) {
2192         x = rb_big_clone(x);
2193         get2comp(x);
2194     }
2195     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
2196         l1 = RBIGNUM_LEN(y);
2197         l2 = RBIGNUM_LEN(x);
2198         ds1 = BDIGITS(y);
2199         ds2 = BDIGITS(x);
2200         sign = RBIGNUM_SIGN(y);
2201     }
2202     else {
2203         l1 = RBIGNUM_LEN(x);
2204         l2 = RBIGNUM_LEN(y);
2205         ds1 = BDIGITS(x);
2206         ds2 = BDIGITS(y);
2207         sign = RBIGNUM_SIGN(x);
2208     }
2209     z = bignew(l2, RBIGNUM_SIGN(x) || RBIGNUM_SIGN(y));
2210     zds = BDIGITS(z);
2211
2212     for (i=0; i<l1; i++) {
2213         zds[i] = ds1[i] & ds2[i];
2214     }
2215     for (; i<l2; i++) {
2216         zds[i] = sign?0:ds2[i];
2217     }
2218     if (!RBIGNUM_SIGN(z)) get2comp(z);
2219     return bignorm(z);
2220 }
2221
2222 /*
2223  * call-seq:
2224  *     big | numeric   =>  integer
2225  *
2226  * Performs bitwise +or+ between _big_ and _numeric_.
2227  */
2228
2229 VALUE
2230 rb_big_or(VALUE xx, VALUE yy)
2231 {
2232     volatile VALUE x, y, z;
2233     BDIGIT *ds1, *ds2, *zds;
2234     long i, l1, l2;
2235     char sign;
2236
2237     x = xx;
2238     y = bit_coerce(yy);
2239     if (FIXNUM_P(y)) {
2240         y = rb_int2big(FIX2LONG(y));
2241     }
2242
2243     if (!RBIGNUM_SIGN(y)) {
2244         y = rb_big_clone(y);
2245         get2comp(y);
2246     }
2247     if (!RBIGNUM_SIGN(x)) {
2248         x = rb_big_clone(x);
2249         get2comp(x);
2250     }
2251     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
2252         l1 = RBIGNUM_LEN(y);
2253         l2 = RBIGNUM_LEN(x);
2254         ds1 = BDIGITS(y);
2255         ds2 = BDIGITS(x);
2256         sign = RBIGNUM_SIGN(y);
2257     }
2258     else {
2259         l1 = RBIGNUM_LEN(x);
2260         l2 = RBIGNUM_LEN(y);
2261         ds1 = BDIGITS(x);
2262         ds2 = BDIGITS(y);
2263         sign = RBIGNUM_SIGN(x);
2264     }
2265     z = bignew(l2, RBIGNUM_SIGN(x) && RBIGNUM_SIGN(y));
2266     zds = BDIGITS(z);
2267
2268     for (i=0; i<l1; i++) {
2269         zds[i] = ds1[i] | ds2[i];
2270     }
2271     for (; i<l2; i++) {
2272         zds[i] = sign?ds2[i]:(BIGRAD-1);
2273     }
2274     if (!RBIGNUM_SIGN(z)) get2comp(z);
2275
2276     return bignorm(z);
2277 }
2278
2279 /*
2280  * call-seq:
2281  *     big ^ numeric   =>  integer
2282  *
2283  * Performs bitwise +exclusive or+ between _big_ and _numeric_.
2284  */
2285
2286 VALUE
2287 rb_big_xor(VALUE xx, VALUE yy)
2288 {
2289     volatile VALUE x, y;
2290     VALUE z;
2291     BDIGIT *ds1, *ds2, *zds;
2292     long i, l1, l2;
2293     char sign;
2294
2295     x = xx;
2296     y = bit_coerce(yy);
2297     if (FIXNUM_P(y)) {
2298         y = rb_int2big(FIX2LONG(y));
2299     }
2300
2301     if (!RBIGNUM_SIGN(y)) {
2302         y = rb_big_clone(y);
2303         get2comp(y);
2304     }
2305     if (!RBIGNUM_SIGN(x)) {
2306         x = rb_big_clone(x);
2307         get2comp(x);
2308     }
2309     if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) {
2310         l1 = RBIGNUM_LEN(y);
2311         l2 = RBIGNUM_LEN(x);
2312         ds1 = BDIGITS(y);
2313         ds2 = BDIGITS(x);
2314         sign = RBIGNUM_SIGN(y);
2315     }
2316     else {
2317         l1 = RBIGNUM_LEN(x);
2318         l2 = RBIGNUM_LEN(y);
2319         ds1 = BDIGITS(x);
2320         ds2 = BDIGITS(y);
2321         sign = RBIGNUM_SIGN(x);
2322     }
2323     RBIGNUM_SET_SIGN(x, RBIGNUM_SIGN(x)?1:0);
2324     RBIGNUM_SET_SIGN(y, RBIGNUM_SIGN(y)?1:0);
2325     z = bignew(l2, !(RBIGNUM_SIGN(x) ^ RBIGNUM_SIGN(y)));
2326     zds = BDIGITS(z);
2327
2328     for (i=0; i<l1; i++) {
2329         zds[i] = ds1[i] ^ ds2[i];
2330     }
2331     for (; i<l2; i++) {
2332         zds[i] = sign?ds2[i]:~ds2[i];
2333     }
2334     if (!RBIGNUM_SIGN(z)) get2comp(z);
2335
2336     return bignorm(z);
2337 }
2338
2339 static VALUE
2340 check_shiftdown(VALUE y, VALUE x)
2341 {
2342     if (!RBIGNUM_LEN(x)) return INT2FIX(0);
2343     if (RBIGNUM_LEN(y) > SIZEOF_LONG / SIZEOF_BDIGITS) {
2344         return RBIGNUM_SIGN(x) ? INT2FIX(0) : INT2FIX(-1);
2345     }
2346     return Qnil;
2347 }
2348
2349 /*
2350  * call-seq:
2351  *     big << numeric   =>  integer
2352  *
2353  * Shifts big left _numeric_ positions (right if _numeric_ is negative).
2354  */
2355
2356 VALUE
2357 rb_big_lshift(VALUE x, VALUE y)
2358 {
2359     long shift;
2360     int neg = 0;
2361
2362     for (;;) {
2363         if (FIXNUM_P(y)) {
2364             shift = FIX2LONG(y);
2365             if (shift < 0) {
2366                 neg = 1;
2367                 shift = -shift;
2368             }
2369             break;
2370         }
2371         else if (TYPE(y) == T_BIGNUM) {
2372             if (!RBIGNUM_SIGN(y)) {
2373                 VALUE t = check_shiftdown(y, x);
2374                 if (!NIL_P(t)) return t;
2375                 neg = 1;
2376             }
2377             shift = big2ulong(y, "long", Qtrue);
2378             break;
2379         }
2380         y = rb_to_int(y);
2381     }
2382
2383     if (neg) return big_rshift(x, shift);
2384     return big_lshift(x, shift);
2385 }
2386
2387 static VALUE
2388 big_lshift(VALUE x, unsigned long shift)
2389 {
2390     BDIGIT *xds, *zds;
2391     long s1 = shift/BITSPERDIG;
2392     int s2 = shift%BITSPERDIG;
2393     VALUE z;
2394     BDIGIT_DBL num = 0;
2395     long len, i;
2396
2397     len = RBIGNUM_LEN(x);
2398     z = bignew(len+s1+1, RBIGNUM_SIGN(x));
2399     zds = BDIGITS(z);
2400     for (i=0; i<s1; i++) {
2401         *zds++ = 0;
2402     }
2403     xds = BDIGITS(x);
2404     for (i=0; i<len; i++) {
2405         num = num | (BDIGIT_DBL)*xds++<<s2;
2406         *zds++ = BIGLO(num);
2407         num = BIGDN(num);
2408     }
2409     *zds = BIGLO(num);
2410     return bignorm(z);
2411 }
2412
2413 /*
2414  * call-seq:
2415  *     big >> numeric   =>  integer
2416  *
2417  * Shifts big right _numeric_ positions (left if _numeric_ is negative).
2418  */
2419
2420 VALUE
2421 rb_big_rshift(VALUE x, VALUE y)
2422 {
2423     long shift;
2424     int neg = 0;
2425
2426     for (;;) {
2427         if (FIXNUM_P(y)) {
2428             shift = FIX2LONG(y);
2429             if (shift < 0) {
2430                 neg = 1;
2431                 shift = -shift;
2432             }
2433             break;
2434         }
2435         else if (TYPE(y) == T_BIGNUM) {
2436             if (RBIGNUM_SIGN(y)) {
2437                 VALUE t = check_shiftdown(y, x);
2438                 if (!NIL_P(t)) return t;
2439             }
2440             else {
2441                 neg = 1;
2442             }
2443             shift = big2ulong(y, "long", Qtrue);
2444             break;
2445         }
2446         y = rb_to_int(y);
2447     }
2448
2449     if (neg) return big_lshift(x, shift);
2450     return big_rshift(x, shift);
2451 }
2452
2453 static VALUE
2454 big_rshift(VALUE x, unsigned long shift)
2455 {
2456     BDIGIT *xds, *zds;
2457     long s1 = shift/BITSPERDIG;
2458     int s2 = shift%BITSPERDIG;
2459     VALUE z;
2460     BDIGIT_DBL num = 0;
2461     long i, j;
2462     volatile VALUE save_x;
2463
2464     if (s1 > RBIGNUM_LEN(x)) {
2465         if (RBIGNUM_SIGN(x))
2466             return INT2FIX(0);
2467         else
2468             return INT2FIX(-1);
2469     }
2470     if (!RBIGNUM_SIGN(x)) {
2471         save_x = x = rb_big_clone(x);
2472         get2comp(x);
2473     }
2474     xds = BDIGITS(x);
2475     i = RBIGNUM_LEN(x); j = i - s1;
2476     if (j == 0) {
2477         if (RBIGNUM_SIGN(x)) return INT2FIX(0);
2478         else return INT2FIX(-1);
2479     }
2480     z = bignew(j, RBIGNUM_SIGN(x));
2481     if (!RBIGNUM_SIGN(x)) {
2482         num = ((BDIGIT_DBL)~0) << BITSPERDIG;
2483     }
2484     zds = BDIGITS(z);
2485     while (i--, j--) {
2486         num = (num | xds[i]) >> s2;
2487         zds[j] = BIGLO(num);
2488         num = BIGUP(xds[i]);
2489     }
2490     if (!RBIGNUM_SIGN(x)) {
2491         get2comp(z);
2492     }
2493     return bignorm(z);
2494 }
2495
2496 /*
2497  *  call-seq:
2498  *     big[n] -> 0, 1
2499  *
2500  *  Bit Reference---Returns the <em>n</em>th bit in the (assumed) binary
2501  *  representation of <i>big</i>, where <i>big</i>[0] is the least
2502  *  significant bit.
2503  *
2504  *     a = 9**15
2505  *     50.downto(0) do |n|
2506  *       print a[n]
2507  *     end
2508  *
2509  *  <em>produces:</em>
2510  *
2511  *     000101110110100000111000011110010100111100010111001
2512  *
2513  */
2514
2515 static VALUE
2516 rb_big_aref(VALUE x, VALUE y)
2517 {
2518     BDIGIT *xds;
2519     BDIGIT_DBL num;
2520     VALUE shift;
2521     long i, s1, s2;
2522
2523     if (TYPE(y) == T_BIGNUM) {
2524         if (!RBIGNUM_SIGN(y))
2525             return INT2FIX(0);
2526         if (RBIGNUM_LEN(bigtrunc(y)) > DIGSPERLONG) {
2527           out_of_range:
2528             return RBIGNUM_SIGN(x) ? INT2FIX(0) : INT2FIX(1);
2529         }
2530         shift = big2ulong(y, "long", Qfalse);
2531     }
2532     else {
2533         i = NUM2LONG(y);
2534         if (i < 0) return INT2FIX(0);
2535         shift = (VALUE)i;
2536     }
2537     s1 = shift/BITSPERDIG;
2538     s2 = shift%BITSPERDIG;
2539
2540     if (s1 >= RBIGNUM_LEN(x)) goto out_of_range;
2541     if (!RBIGNUM_SIGN(x)) {
2542         xds = BDIGITS(x);
2543         i = 0; num = 1;
2544         while (num += ~xds[i], ++i <= s1) {
2545             num = BIGDN(num);
2546         }
2547     }
2548     else {
2549         num = BDIGITS(x)[s1];
2550     }
2551     if (num & ((BDIGIT_DBL)1<<s2))
2552         return INT2FIX(1);
2553     return INT2FIX(0);
2554 }
2555
2556 /*
2557  * call-seq:
2558  *   big.hash   => fixnum
2559  *
2560  * Compute a hash based on the value of _big_.
2561  */
2562
2563 static VALUE
2564 rb_big_hash(VALUE x)
2565 {
2566     int hash;
2567
2568     hash = rb_memhash(BDIGITS(x), sizeof(BDIGIT)*RBIGNUM_LEN(x)) ^ RBIGNUM_SIGN(x);
2569     return INT2FIX(hash);
2570 }
2571
2572 /*
2573  * MISSING: documentation
2574  */
2575
2576 static VALUE
2577 rb_big_coerce(VALUE x, VALUE y)
2578 {
2579     if (FIXNUM_P(y)) {
2580         return rb_assoc_new(rb_int2big(FIX2LONG(y)), x);
2581     }
2582     else if (TYPE(y) == T_BIGNUM) {
2583        return rb_assoc_new(y, x);
2584     }
2585     else {
2586         rb_raise(rb_eTypeError, "can't coerce %s to Bignum",
2587                  rb_obj_classname(y));
2588     }
2589     /* not reached */
2590     return Qnil;
2591 }
2592
2593 /*
2594  *  call-seq:
2595  *     big.abs -> aBignum
2596  *
2597  *  Returns the absolute value of <i>big</i>.
2598  *
2599  *     -1234567890987654321.abs   #=> 1234567890987654321
2600  */
2601
2602 static VALUE
2603 rb_big_abs(VALUE x)
2604 {
2605     if (!RBIGNUM_SIGN(x)) {
2606         x = rb_big_clone(x);
2607         RBIGNUM_SET_SIGN(x, 1);
2608     }
2609     return x;
2610 }
2611
2612 /*
2613  *  call-seq:
2614  *     big.size -> integer
2615  *
2616  *  Returns the number of bytes in the machine representation of
2617  *  <i>big</i>.
2618  *
2619  *     (256**10 - 1).size   #=> 12
2620  *     (256**20 - 1).size   #=> 20
2621  *     (256**40 - 1).size   #=> 40
2622  */
2623
2624 static VALUE
2625 rb_big_size(VALUE big)
2626 {
2627     return LONG2FIX(RBIGNUM_LEN(big)*SIZEOF_BDIGITS);
2628 }
2629
2630 /*
2631  *  call-seq:
2632  *     big.odd? -> true or false
2633  *
2634  *  Returns <code>true</code> if <i>big</i> is an odd number.
2635  */
2636
2637 static VALUE
2638 rb_big_odd_p(VALUE num)
2639 {
2640     if (BDIGITS(num)[0] & 1) {
2641         return Qtrue;
2642     }
2643     return Qfalse;
2644 }
2645
2646 /*
2647  *  call-seq:
2648  *     big.even? -> true or false
2649  *
2650  *  Returns <code>true</code> if <i>big</i> is an even number.
2651  */
2652
2653 static VALUE
2654 rb_big_even_p(VALUE num)
2655 {
2656     if (BDIGITS(num)[0] & 1) {
2657         return Qfalse;
2658     }
2659     return Qtrue;
2660 }
2661
2662 /*
2663  *  Bignum objects hold integers outside the range of
2664  *  Fixnum. Bignum objects are created
2665  *  automatically when integer calculations would otherwise overflow a
2666  *  Fixnum. When a calculation involving
2667  *  Bignum objects returns a result that will fit in a
2668  *  Fixnum, the result is automatically converted.
2669  *
2670  *  For the purposes of the bitwise operations and <code>[]</code>, a
2671  *  Bignum is treated as if it were an infinite-length
2672  *  bitstring with 2's complement representation.
2673  *
2674  *  While Fixnum values are immediate, Bignum
2675  *  objects are not---assignment and parameter passing work with
2676  *  references to objects, not the objects themselves.
2677  *
2678  */
2679
2680 void
2681 Init_Bignum(void)
2682 {
2683     rb_cBignum = rb_define_class("Bignum", rb_cInteger);
2684
2685     rb_define_method(rb_cBignum, "to_s", rb_big_to_s, -1);
2686     rb_define_method(rb_cBignum, "coerce", rb_big_coerce, 1);
2687     rb_define_method(rb_cBignum, "-@", rb_big_uminus, 0);
2688     rb_define_method(rb_cBignum, "+", rb_big_plus, 1);
2689     rb_define_method(rb_cBignum, "-", rb_big_minus, 1);
2690     rb_define_method(rb_cBignum, "*", rb_big_mul, 1);
2691     rb_define_method(rb_cBignum, "/", rb_big_div, 1);
2692     rb_define_method(rb_cBignum, "%", rb_big_modulo, 1);
2693     rb_define_method(rb_cBignum, "div", rb_big_idiv, 1);
2694     rb_define_method(rb_cBignum, "divmod", rb_big_divmod, 1);
2695     rb_define_method(rb_cBignum, "modulo", rb_big_modulo, 1);
2696     rb_define_method(rb_cBignum, "remainder", rb_big_remainder, 1);
2697     rb_define_method(rb_cBignum, "fdiv", rb_big_fdiv, 1);
2698     rb_define_method(rb_cBignum, "**", rb_big_pow, 1);
2699     rb_define_method(rb_cBignum, "&", rb_big_and, 1);
2700     rb_define_method(rb_cBignum, "|", rb_big_or, 1);
2701     rb_define_method(rb_cBignum, "^", rb_big_xor, 1);
2702     rb_define_method(rb_cBignum, "~", rb_big_neg, 0);
2703     rb_define_method(rb_cBignum, "<<", rb_big_lshift, 1);
2704     rb_define_method(rb_cBignum, ">>", rb_big_rshift, 1);
2705     rb_define_method(rb_cBignum, "[]", rb_big_aref, 1);
2706
2707     rb_define_method(rb_cBignum, "<=>", rb_big_cmp, 1);
2708     rb_define_method(rb_cBignum, "==", rb_big_eq, 1);
2709     rb_define_method(rb_cBignum, "eql?", rb_big_eql, 1);
2710     rb_define_method(rb_cBignum, "hash", rb_big_hash, 0);
2711     rb_define_method(rb_cBignum, "to_f", rb_big_to_f, 0);
2712     rb_define_method(rb_cBignum, "abs", rb_big_abs, 0);
2713     rb_define_method(rb_cBignum, "size", rb_big_size, 0);
2714     rb_define_method(rb_cBignum, "odd?", rb_big_odd_p, 0);
2715     rb_define_method(rb_cBignum, "even?", rb_big_even_p, 0);
2716
2717     power_cache_init();
2718 }