1/********************************************************************** 2 3 bignum.c - 4 5 $Author: nagachika $ 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#include "ruby/thread.h" 14#include "ruby/util.h" 15#include "internal.h" 16 17#ifdef HAVE_STRINGS_H 18#include <strings.h> 19#endif 20#include <math.h> 21#include <float.h> 22#include <ctype.h> 23#ifdef HAVE_IEEEFP_H 24#include <ieeefp.h> 25#endif 26#include <assert.h> 27 28VALUE rb_cBignum; 29 30static VALUE big_three = Qnil; 31 32#if defined __MINGW32__ 33#define USHORT _USHORT 34#endif 35 36#define BDIGITS(x) (RBIGNUM_DIGITS(x)) 37#define BITSPERDIG (SIZEOF_BDIGITS*CHAR_BIT) 38#define BIGRAD ((BDIGIT_DBL)1 << BITSPERDIG) 39#define BIGRAD_HALF ((BDIGIT)(BIGRAD >> 1)) 40#define DIGSPERLONG (SIZEOF_LONG/SIZEOF_BDIGITS) 41#if HAVE_LONG_LONG 42# define DIGSPERLL (SIZEOF_LONG_LONG/SIZEOF_BDIGITS) 43#endif 44#define BIGUP(x) ((BDIGIT_DBL)(x) << BITSPERDIG) 45#define BIGDN(x) RSHIFT((x),BITSPERDIG) 46#define BIGLO(x) ((BDIGIT)((x) & (BIGRAD-1))) 47#define BDIGMAX ((BDIGIT)-1) 48 49#define BIGZEROP(x) (RBIGNUM_LEN(x) == 0 || \ 50 (BDIGITS(x)[0] == 0 && \ 51 (RBIGNUM_LEN(x) == 1 || bigzero_p(x)))) 52 53#define BIGNUM_DEBUG 0 54#if BIGNUM_DEBUG 55#define ON_DEBUG(x) do { x; } while (0) 56static void 57dump_bignum(VALUE x) 58{ 59 long i; 60 printf("%c0x0", RBIGNUM_SIGN(x) ? '+' : '-'); 61 for (i = RBIGNUM_LEN(x); i--; ) { 62 printf("_%08"PRIxBDIGIT, BDIGITS(x)[i]); 63 } 64 printf(", len=%lu", RBIGNUM_LEN(x)); 65 puts(""); 66} 67 68static VALUE 69rb_big_dump(VALUE x) 70{ 71 dump_bignum(x); 72 return x; 73} 74#else 75#define ON_DEBUG(x) 76#endif 77 78static int 79bigzero_p(VALUE x) 80{ 81 long i; 82 BDIGIT *ds = BDIGITS(x); 83 84 for (i = RBIGNUM_LEN(x) - 1; 0 <= i; i--) { 85 if (ds[i]) return 0; 86 } 87 return 1; 88} 89 90int 91rb_bigzero_p(VALUE x) 92{ 93 return BIGZEROP(x); 94} 95 96int 97rb_cmpint(VALUE val, VALUE a, VALUE b) 98{ 99 if (NIL_P(val)) { 100 rb_cmperr(a, b); 101 } 102 if (FIXNUM_P(val)) { 103 long l = FIX2LONG(val); 104 if (l > 0) return 1; 105 if (l < 0) return -1; 106 return 0; 107 } 108 if (RB_TYPE_P(val, T_BIGNUM)) { 109 if (BIGZEROP(val)) return 0; 110 if (RBIGNUM_SIGN(val)) return 1; 111 return -1; 112 } 113 if (RTEST(rb_funcall(val, '>', 1, INT2FIX(0)))) return 1; 114 if (RTEST(rb_funcall(val, '<', 1, INT2FIX(0)))) return -1; 115 return 0; 116} 117 118#define RBIGNUM_SET_LEN(b,l) \ 119 ((RBASIC(b)->flags & RBIGNUM_EMBED_FLAG) ? \ 120 (void)(RBASIC(b)->flags = \ 121 (RBASIC(b)->flags & ~RBIGNUM_EMBED_LEN_MASK) | \ 122 ((l) << RBIGNUM_EMBED_LEN_SHIFT)) : \ 123 (void)(RBIGNUM(b)->as.heap.len = (l))) 124 125static void 126rb_big_realloc(VALUE big, long len) 127{ 128 BDIGIT *ds; 129 if (RBASIC(big)->flags & RBIGNUM_EMBED_FLAG) { 130 if (RBIGNUM_EMBED_LEN_MAX < len) { 131 ds = ALLOC_N(BDIGIT, len); 132 MEMCPY(ds, RBIGNUM(big)->as.ary, BDIGIT, RBIGNUM_EMBED_LEN_MAX); 133 RBIGNUM(big)->as.heap.len = RBIGNUM_LEN(big); 134 RBIGNUM(big)->as.heap.digits = ds; 135 RBASIC(big)->flags &= ~RBIGNUM_EMBED_FLAG; 136 } 137 } 138 else { 139 if (len <= RBIGNUM_EMBED_LEN_MAX) { 140 ds = RBIGNUM(big)->as.heap.digits; 141 RBASIC(big)->flags |= RBIGNUM_EMBED_FLAG; 142 RBIGNUM_SET_LEN(big, len); 143 if (ds) { 144 MEMCPY(RBIGNUM(big)->as.ary, ds, BDIGIT, len); 145 xfree(ds); 146 } 147 } 148 else { 149 if (RBIGNUM_LEN(big) == 0) { 150 RBIGNUM(big)->as.heap.digits = ALLOC_N(BDIGIT, len); 151 } 152 else { 153 REALLOC_N(RBIGNUM(big)->as.heap.digits, BDIGIT, len); 154 } 155 } 156 } 157} 158 159void 160rb_big_resize(VALUE big, long len) 161{ 162 rb_big_realloc(big, len); 163 RBIGNUM_SET_LEN(big, len); 164} 165 166static VALUE 167bignew_1(VALUE klass, long len, int sign) 168{ 169 NEWOBJ_OF(big, struct RBignum, klass, T_BIGNUM); 170 RBIGNUM_SET_SIGN(big, sign?1:0); 171 if (len <= RBIGNUM_EMBED_LEN_MAX) { 172 RBASIC(big)->flags |= RBIGNUM_EMBED_FLAG; 173 RBIGNUM_SET_LEN(big, len); 174 } 175 else { 176 RBIGNUM(big)->as.heap.digits = ALLOC_N(BDIGIT, len); 177 RBIGNUM(big)->as.heap.len = len; 178 } 179 OBJ_FREEZE(big); 180 return (VALUE)big; 181} 182 183#define bignew(len,sign) bignew_1(rb_cBignum,(len),(sign)) 184 185VALUE 186rb_big_new(long len, int sign) 187{ 188 return bignew(len, sign != 0); 189} 190 191VALUE 192rb_big_clone(VALUE x) 193{ 194 long len = RBIGNUM_LEN(x); 195 VALUE z = bignew_1(CLASS_OF(x), len, RBIGNUM_SIGN(x)); 196 197 MEMCPY(BDIGITS(z), BDIGITS(x), BDIGIT, len); 198 return z; 199} 200 201/* modify a bignum by 2's complement */ 202static void 203get2comp(VALUE x) 204{ 205 long i = RBIGNUM_LEN(x); 206 BDIGIT *ds = BDIGITS(x); 207 BDIGIT_DBL num; 208 209 if (!i) return; 210 while (i--) ds[i] = ~ds[i]; 211 i = 0; num = 1; 212 do { 213 num += ds[i]; 214 ds[i++] = BIGLO(num); 215 num = BIGDN(num); 216 } while (i < RBIGNUM_LEN(x)); 217 if (num != 0) { 218 rb_big_resize(x, RBIGNUM_LEN(x)+1); 219 ds = BDIGITS(x); 220 ds[RBIGNUM_LEN(x)-1] = 1; 221 } 222} 223 224void 225rb_big_2comp(VALUE x) /* get 2's complement */ 226{ 227 get2comp(x); 228} 229 230static inline VALUE 231bigtrunc(VALUE x) 232{ 233 long len = RBIGNUM_LEN(x); 234 BDIGIT *ds = BDIGITS(x); 235 236 if (len == 0) return x; 237 while (--len && !ds[len]); 238 if (RBIGNUM_LEN(x) > len+1) { 239 rb_big_resize(x, len+1); 240 } 241 return x; 242} 243 244static inline VALUE 245bigfixize(VALUE x) 246{ 247 long len = RBIGNUM_LEN(x); 248 BDIGIT *ds = BDIGITS(x); 249 250 if (len == 0) return INT2FIX(0); 251 if ((size_t)(len*SIZEOF_BDIGITS) <= sizeof(long)) { 252 long num = 0; 253#if 2*SIZEOF_BDIGITS > SIZEOF_LONG 254 num = (long)ds[0]; 255#else 256 while (len--) { 257 num = (long)(BIGUP(num) + ds[len]); 258 } 259#endif 260 if (num >= 0) { 261 if (RBIGNUM_SIGN(x)) { 262 if (POSFIXABLE(num)) return LONG2FIX(num); 263 } 264 else { 265 if (NEGFIXABLE(-num)) return LONG2FIX(-num); 266 } 267 } 268 } 269 return x; 270} 271 272static VALUE 273bignorm(VALUE x) 274{ 275 if (RB_TYPE_P(x, T_BIGNUM)) { 276 x = bigfixize(bigtrunc(x)); 277 } 278 return x; 279} 280 281VALUE 282rb_big_norm(VALUE x) 283{ 284 return bignorm(x); 285} 286 287VALUE 288rb_uint2big(VALUE n) 289{ 290 BDIGIT_DBL num = n; 291 long i = 0; 292 BDIGIT *digits; 293 VALUE big; 294 295 big = bignew(DIGSPERLONG, 1); 296 digits = BDIGITS(big); 297 while (i < DIGSPERLONG) { 298 digits[i++] = BIGLO(num); 299 num = BIGDN(num); 300 } 301 302 i = DIGSPERLONG; 303 while (--i && !digits[i]) ; 304 RBIGNUM_SET_LEN(big, i+1); 305 return big; 306} 307 308VALUE 309rb_int2big(SIGNED_VALUE n) 310{ 311 long neg = 0; 312 VALUE u; 313 VALUE big; 314 315 if (n < 0) { 316 u = 1 + (VALUE)(-(n + 1)); /* u = -n avoiding overflow */ 317 neg = 1; 318 } 319 else { 320 u = n; 321 } 322 big = rb_uint2big(u); 323 if (neg) { 324 RBIGNUM_SET_SIGN(big, 0); 325 } 326 return big; 327} 328 329VALUE 330rb_uint2inum(VALUE n) 331{ 332 if (POSFIXABLE(n)) return LONG2FIX(n); 333 return rb_uint2big(n); 334} 335 336VALUE 337rb_int2inum(SIGNED_VALUE n) 338{ 339 if (FIXABLE(n)) return LONG2FIX(n); 340 return rb_int2big(n); 341} 342 343#if SIZEOF_LONG % SIZEOF_BDIGITS != 0 344# error unexpected SIZEOF_LONG : SIZEOF_BDIGITS ratio 345#endif 346 347/* 348 * buf is an array of long integers. 349 * buf is ordered from least significant word to most significant word. 350 * buf[0] is the least significant word and 351 * buf[num_longs-1] is the most significant word. 352 * This means words in buf is little endian. 353 * However each word in buf is native endian. 354 * (buf[i]&1) is the least significant bit and 355 * (buf[i]&(1<<(SIZEOF_LONG*CHAR_BIT-1))) is the most significant bit 356 * for each 0 <= i < num_longs. 357 * So buf is little endian at whole on a little endian machine. 358 * But buf is mixed endian on a big endian machine. 359 * 360 * The buf represents negative integers as two's complement. 361 * So, the most significant bit of the most significant word, 362 * (buf[num_longs-1]>>(SIZEOF_LONG*CHAR_BIT-1)), 363 * is the sign bit: 1 means negative and 0 means zero or positive. 364 * 365 * If given size of buf (num_longs) is not enough to represent val, 366 * higier words (including a sign bit) are ignored. 367 */ 368void 369rb_big_pack(VALUE val, unsigned long *buf, long num_longs) 370{ 371 val = rb_to_int(val); 372 if (num_longs == 0) 373 return; 374 if (FIXNUM_P(val)) { 375 long i; 376 long tmp = FIX2LONG(val); 377 buf[0] = (unsigned long)tmp; 378 tmp = tmp < 0 ? ~0L : 0; 379 for (i = 1; i < num_longs; i++) 380 buf[i] = (unsigned long)tmp; 381 return; 382 } 383 else { 384 long len = RBIGNUM_LEN(val); 385 BDIGIT *ds = BDIGITS(val), *dend = ds + len; 386 long i, j; 387 for (i = 0; i < num_longs && ds < dend; i++) { 388 unsigned long l = 0; 389 for (j = 0; j < DIGSPERLONG && ds < dend; j++, ds++) { 390 l |= ((unsigned long)*ds << (j * BITSPERDIG)); 391 } 392 buf[i] = l; 393 } 394 for (; i < num_longs; i++) 395 buf[i] = 0; 396 if (RBIGNUM_NEGATIVE_P(val)) { 397 for (i = 0; i < num_longs; i++) { 398 buf[i] = ~buf[i]; 399 } 400 for (i = 0; i < num_longs; i++) { 401 buf[i]++; 402 if (buf[i] != 0) 403 return; 404 } 405 } 406 } 407} 408 409/* See rb_big_pack comment for endianness and sign of buf. */ 410VALUE 411rb_big_unpack(unsigned long *buf, long num_longs) 412{ 413 while (2 <= num_longs) { 414 if (buf[num_longs-1] == 0 && (long)buf[num_longs-2] >= 0) 415 num_longs--; 416 else if (buf[num_longs-1] == ~0UL && (long)buf[num_longs-2] < 0) 417 num_longs--; 418 else 419 break; 420 } 421 if (num_longs == 0) 422 return INT2FIX(0); 423 else if (num_longs == 1) 424 return LONG2NUM((long)buf[0]); 425 else { 426 VALUE big; 427 BDIGIT *ds; 428 long len = num_longs * DIGSPERLONG; 429 long i; 430 big = bignew(len, 1); 431 ds = BDIGITS(big); 432 for (i = 0; i < num_longs; i++) { 433 unsigned long d = buf[i]; 434#if SIZEOF_LONG == SIZEOF_BDIGITS 435 *ds++ = d; 436#else 437 int j; 438 for (j = 0; j < DIGSPERLONG; j++) { 439 *ds++ = BIGLO(d); 440 d = BIGDN(d); 441 } 442#endif 443 } 444 if ((long)buf[num_longs-1] < 0) { 445 get2comp(big); 446 RBIGNUM_SET_SIGN(big, 0); 447 } 448 return bignorm(big); 449 } 450} 451 452#define QUAD_SIZE 8 453 454#if SIZEOF_LONG_LONG == QUAD_SIZE && SIZEOF_BDIGITS*2 == SIZEOF_LONG_LONG 455 456void 457rb_quad_pack(char *buf, VALUE val) 458{ 459 LONG_LONG q; 460 461 val = rb_to_int(val); 462 if (FIXNUM_P(val)) { 463 q = FIX2LONG(val); 464 } 465 else { 466 long len = RBIGNUM_LEN(val); 467 BDIGIT *ds; 468 469 if (len > SIZEOF_LONG_LONG/SIZEOF_BDIGITS) { 470 len = SIZEOF_LONG_LONG/SIZEOF_BDIGITS; 471 } 472 ds = BDIGITS(val); 473 q = 0; 474 while (len--) { 475 q = BIGUP(q); 476 q += ds[len]; 477 } 478 if (!RBIGNUM_SIGN(val)) q = -q; 479 } 480 memcpy(buf, (char*)&q, SIZEOF_LONG_LONG); 481} 482 483VALUE 484rb_quad_unpack(const char *buf, int sign) 485{ 486 unsigned LONG_LONG q; 487 long neg = 0; 488 long i; 489 BDIGIT *digits; 490 VALUE big; 491 492 memcpy(&q, buf, SIZEOF_LONG_LONG); 493 if (sign) { 494 if (FIXABLE((LONG_LONG)q)) return LONG2FIX((LONG_LONG)q); 495 if ((LONG_LONG)q < 0) { 496 q = -(LONG_LONG)q; 497 neg = 1; 498 } 499 } 500 else { 501 if (POSFIXABLE(q)) return LONG2FIX(q); 502 } 503 504 i = 0; 505 big = bignew(DIGSPERLL, 1); 506 digits = BDIGITS(big); 507 while (i < DIGSPERLL) { 508 digits[i++] = BIGLO(q); 509 q = BIGDN(q); 510 } 511 512 i = DIGSPERLL; 513 while (i-- && !digits[i]) ; 514 RBIGNUM_SET_LEN(big, i+1); 515 516 if (neg) { 517 RBIGNUM_SET_SIGN(big, 0); 518 } 519 return bignorm(big); 520} 521 522#else 523 524static int 525quad_buf_complement(char *buf, size_t len) 526{ 527 size_t i; 528 for (i = 0; i < len; i++) 529 buf[i] = ~buf[i]; 530 for (i = 0; i < len; i++) { 531 buf[i]++; 532 if (buf[i] != 0) 533 return 0; 534 } 535 return 1; 536} 537 538void 539rb_quad_pack(char *buf, VALUE val) 540{ 541 long len; 542 543 memset(buf, 0, QUAD_SIZE); 544 val = rb_to_int(val); 545 if (FIXNUM_P(val)) { 546 val = rb_int2big(FIX2LONG(val)); 547 } 548 len = RBIGNUM_LEN(val) * SIZEOF_BDIGITS; 549 if (len > QUAD_SIZE) { 550 len = QUAD_SIZE; 551 } 552 memcpy(buf, (char*)BDIGITS(val), len); 553 if (RBIGNUM_NEGATIVE_P(val)) { 554 quad_buf_complement(buf, QUAD_SIZE); 555 } 556} 557 558#define BNEG(b) (RSHIFT(((BDIGIT*)(b))[QUAD_SIZE/SIZEOF_BDIGITS-1],BITSPERDIG-1) != 0) 559 560VALUE 561rb_quad_unpack(const char *buf, int sign) 562{ 563 VALUE big = bignew(QUAD_SIZE/SIZEOF_BDIGITS, 1); 564 565 memcpy((char*)BDIGITS(big), buf, QUAD_SIZE); 566 if (sign && BNEG(buf)) { 567 char *tmp = (char*)BDIGITS(big); 568 569 RBIGNUM_SET_SIGN(big, 0); 570 quad_buf_complement(tmp, QUAD_SIZE); 571 } 572 573 return bignorm(big); 574} 575 576#endif 577 578VALUE 579rb_cstr_to_inum(const char *str, int base, int badcheck) 580{ 581 const char *s = str; 582 char *end; 583 char sign = 1, nondigit = 0; 584 int c; 585 BDIGIT_DBL num; 586 long len, blen = 1; 587 long i; 588 VALUE z; 589 BDIGIT *zds; 590 591#undef ISDIGIT 592#define ISDIGIT(c) ('0' <= (c) && (c) <= '9') 593#define conv_digit(c) \ 594 (!ISASCII(c) ? -1 : \ 595 ISDIGIT(c) ? ((c) - '0') : \ 596 ISLOWER(c) ? ((c) - 'a' + 10) : \ 597 ISUPPER(c) ? ((c) - 'A' + 10) : \ 598 -1) 599 600 if (!str) { 601 if (badcheck) goto bad; 602 return INT2FIX(0); 603 } 604 while (ISSPACE(*str)) str++; 605 606 if (str[0] == '+') { 607 str++; 608 } 609 else if (str[0] == '-') { 610 str++; 611 sign = 0; 612 } 613 if (str[0] == '+' || str[0] == '-') { 614 if (badcheck) goto bad; 615 return INT2FIX(0); 616 } 617 if (base <= 0) { 618 if (str[0] == '0') { 619 switch (str[1]) { 620 case 'x': case 'X': 621 base = 16; 622 break; 623 case 'b': case 'B': 624 base = 2; 625 break; 626 case 'o': case 'O': 627 base = 8; 628 break; 629 case 'd': case 'D': 630 base = 10; 631 break; 632 default: 633 base = 8; 634 } 635 } 636 else if (base < -1) { 637 base = -base; 638 } 639 else { 640 base = 10; 641 } 642 } 643 switch (base) { 644 case 2: 645 len = 1; 646 if (str[0] == '0' && (str[1] == 'b'||str[1] == 'B')) { 647 str += 2; 648 } 649 break; 650 case 3: 651 len = 2; 652 break; 653 case 8: 654 if (str[0] == '0' && (str[1] == 'o'||str[1] == 'O')) { 655 str += 2; 656 } 657 case 4: case 5: case 6: case 7: 658 len = 3; 659 break; 660 case 10: 661 if (str[0] == '0' && (str[1] == 'd'||str[1] == 'D')) { 662 str += 2; 663 } 664 case 9: case 11: case 12: case 13: case 14: case 15: 665 len = 4; 666 break; 667 case 16: 668 len = 4; 669 if (str[0] == '0' && (str[1] == 'x'||str[1] == 'X')) { 670 str += 2; 671 } 672 break; 673 default: 674 if (base < 2 || 36 < base) { 675 rb_raise(rb_eArgError, "invalid radix %d", base); 676 } 677 if (base <= 32) { 678 len = 5; 679 } 680 else { 681 len = 6; 682 } 683 break; 684 } 685 if (*str == '0') { /* squeeze preceding 0s */ 686 int us = 0; 687 while ((c = *++str) == '0' || c == '_') { 688 if (c == '_') { 689 if (++us >= 2) 690 break; 691 } else 692 us = 0; 693 } 694 if (!(c = *str) || ISSPACE(c)) --str; 695 } 696 c = *str; 697 c = conv_digit(c); 698 if (c < 0 || c >= base) { 699 if (badcheck) goto bad; 700 return INT2FIX(0); 701 } 702 len *= strlen(str)*sizeof(char); 703 704 if ((size_t)len <= (sizeof(long)*CHAR_BIT)) { 705 unsigned long val = STRTOUL(str, &end, base); 706 707 if (str < end && *end == '_') goto bigparse; 708 if (badcheck) { 709 if (end == str) goto bad; /* no number */ 710 while (*end && ISSPACE(*end)) end++; 711 if (*end) goto bad; /* trailing garbage */ 712 } 713 714 if (POSFIXABLE(val)) { 715 if (sign) return LONG2FIX(val); 716 else { 717 long result = -(long)val; 718 return LONG2FIX(result); 719 } 720 } 721 else { 722 VALUE big = rb_uint2big(val); 723 RBIGNUM_SET_SIGN(big, sign); 724 return bignorm(big); 725 } 726 } 727 bigparse: 728 len = (len/BITSPERDIG)+1; 729 if (badcheck && *str == '_') goto bad; 730 731 z = bignew(len, sign); 732 zds = BDIGITS(z); 733 for (i=len;i--;) zds[i]=0; 734 while ((c = *str++) != 0) { 735 if (c == '_') { 736 if (nondigit) { 737 if (badcheck) goto bad; 738 break; 739 } 740 nondigit = (char) c; 741 continue; 742 } 743 else if ((c = conv_digit(c)) < 0) { 744 break; 745 } 746 if (c >= base) break; 747 nondigit = 0; 748 i = 0; 749 num = c; 750 for (;;) { 751 while (i<blen) { 752 num += (BDIGIT_DBL)zds[i]*base; 753 zds[i++] = BIGLO(num); 754 num = BIGDN(num); 755 } 756 if (num) { 757 blen++; 758 continue; 759 } 760 break; 761 } 762 } 763 if (badcheck) { 764 str--; 765 if (s+1 < str && str[-1] == '_') goto bad; 766 while (*str && ISSPACE(*str)) str++; 767 if (*str) { 768 bad: 769 rb_invalid_str(s, "Integer()"); 770 } 771 } 772 773 return bignorm(z); 774} 775 776VALUE 777rb_str_to_inum(VALUE str, int base, int badcheck) 778{ 779 char *s; 780 long len; 781 VALUE v = 0; 782 VALUE ret; 783 784 StringValue(str); 785 rb_must_asciicompat(str); 786 if (badcheck) { 787 s = StringValueCStr(str); 788 } 789 else { 790 s = RSTRING_PTR(str); 791 } 792 if (s) { 793 len = RSTRING_LEN(str); 794 if (s[len]) { /* no sentinel somehow */ 795 char *p = ALLOCV(v, len+1); 796 797 MEMCPY(p, s, char, len); 798 p[len] = '\0'; 799 s = p; 800 } 801 } 802 ret = rb_cstr_to_inum(s, base, badcheck); 803 if (v) 804 ALLOCV_END(v); 805 return ret; 806} 807 808#if HAVE_LONG_LONG 809 810static VALUE 811rb_ull2big(unsigned LONG_LONG n) 812{ 813 BDIGIT_DBL num = n; 814 long i = 0; 815 BDIGIT *digits; 816 VALUE big; 817 818 big = bignew(DIGSPERLL, 1); 819 digits = BDIGITS(big); 820 while (i < DIGSPERLL) { 821 digits[i++] = BIGLO(num); 822 num = BIGDN(num); 823 } 824 825 i = DIGSPERLL; 826 while (i-- && !digits[i]) ; 827 RBIGNUM_SET_LEN(big, i+1); 828 return big; 829} 830 831static VALUE 832rb_ll2big(LONG_LONG n) 833{ 834 long neg = 0; 835 VALUE big; 836 837 if (n < 0) { 838 n = -n; 839 neg = 1; 840 } 841 big = rb_ull2big(n); 842 if (neg) { 843 RBIGNUM_SET_SIGN(big, 0); 844 } 845 return big; 846} 847 848VALUE 849rb_ull2inum(unsigned LONG_LONG n) 850{ 851 if (POSFIXABLE(n)) return LONG2FIX(n); 852 return rb_ull2big(n); 853} 854 855VALUE 856rb_ll2inum(LONG_LONG n) 857{ 858 if (FIXABLE(n)) return LONG2FIX(n); 859 return rb_ll2big(n); 860} 861 862#endif /* HAVE_LONG_LONG */ 863 864VALUE 865rb_cstr2inum(const char *str, int base) 866{ 867 return rb_cstr_to_inum(str, base, base==0); 868} 869 870VALUE 871rb_str2inum(VALUE str, int base) 872{ 873 return rb_str_to_inum(str, base, base==0); 874} 875 876const char ruby_digitmap[] = "0123456789abcdefghijklmnopqrstuvwxyz"; 877 878static VALUE bigsqr(VALUE x); 879static void bigdivmod(VALUE x, VALUE y, volatile VALUE *divp, volatile VALUE *modp); 880 881#define POW2_P(x) (((x)&((x)-1))==0) 882 883static inline int 884ones(register unsigned long x) 885{ 886#if SIZEOF_LONG == 8 887# define MASK_55 0x5555555555555555UL 888# define MASK_33 0x3333333333333333UL 889# define MASK_0f 0x0f0f0f0f0f0f0f0fUL 890#else 891# define MASK_55 0x55555555UL 892# define MASK_33 0x33333333UL 893# define MASK_0f 0x0f0f0f0fUL 894#endif 895 x -= (x >> 1) & MASK_55; 896 x = ((x >> 2) & MASK_33) + (x & MASK_33); 897 x = ((x >> 4) + x) & MASK_0f; 898 x += (x >> 8); 899 x += (x >> 16); 900#if SIZEOF_LONG == 8 901 x += (x >> 32); 902#endif 903 return (int)(x & 0x7f); 904#undef MASK_0f 905#undef MASK_33 906#undef MASK_55 907} 908 909static inline unsigned long 910next_pow2(register unsigned long x) 911{ 912 x |= x >> 1; 913 x |= x >> 2; 914 x |= x >> 4; 915 x |= x >> 8; 916 x |= x >> 16; 917#if SIZEOF_LONG == 8 918 x |= x >> 32; 919#endif 920 return x + 1; 921} 922 923static inline int 924floor_log2(register unsigned long x) 925{ 926 x |= x >> 1; 927 x |= x >> 2; 928 x |= x >> 4; 929 x |= x >> 8; 930 x |= x >> 16; 931#if SIZEOF_LONG == 8 932 x |= x >> 32; 933#endif 934 return (int)ones(x) - 1; 935} 936 937static inline int 938ceil_log2(register unsigned long x) 939{ 940 return floor_log2(x) + !POW2_P(x); 941} 942 943#define LOG2_KARATSUBA_DIGITS 7 944#define KARATSUBA_DIGITS (1L<<LOG2_KARATSUBA_DIGITS) 945#define MAX_BIG2STR_TABLE_ENTRIES 64 946 947static VALUE big2str_power_cache[35][MAX_BIG2STR_TABLE_ENTRIES]; 948 949static void 950power_cache_init(void) 951{ 952 int i, j; 953 for (i = 0; i < 35; ++i) { 954 for (j = 0; j < MAX_BIG2STR_TABLE_ENTRIES; ++j) { 955 big2str_power_cache[i][j] = Qnil; 956 } 957 } 958} 959 960static inline VALUE 961power_cache_get_power0(int base, int i) 962{ 963 if (NIL_P(big2str_power_cache[base - 2][i])) { 964 big2str_power_cache[base - 2][i] = 965 i == 0 ? rb_big_pow(rb_int2big(base), INT2FIX(KARATSUBA_DIGITS)) 966 : bigsqr(power_cache_get_power0(base, i - 1)); 967 rb_gc_register_mark_object(big2str_power_cache[base - 2][i]); 968 } 969 return big2str_power_cache[base - 2][i]; 970} 971 972static VALUE 973power_cache_get_power(int base, long n1, long* m1) 974{ 975 int i, m; 976 long j; 977 VALUE t; 978 979 if (n1 <= KARATSUBA_DIGITS) 980 rb_bug("n1 > KARATSUBA_DIGITS"); 981 982 m = ceil_log2(n1); 983 if (m1) *m1 = 1 << m; 984 i = m - LOG2_KARATSUBA_DIGITS; 985 if (i >= MAX_BIG2STR_TABLE_ENTRIES) 986 i = MAX_BIG2STR_TABLE_ENTRIES - 1; 987 t = power_cache_get_power0(base, i); 988 989 j = KARATSUBA_DIGITS*(1 << i); 990 while (n1 > j) { 991 t = bigsqr(t); 992 j *= 2; 993 } 994 return t; 995} 996 997/* big2str_muraken_find_n1 998 * 999 * Let a natural number x is given by: 1000 * x = 2^0 * x_0 + 2^1 * x_1 + ... + 2^(B*n_0 - 1) * x_{B*n_0 - 1}, 1001 * where B is BITSPERDIG (i.e. BDIGITS*CHAR_BIT) and n_0 is 1002 * RBIGNUM_LEN(x). 1003 * 1004 * Now, we assume n_1 = min_n \{ n | 2^(B*n_0/2) <= b_1^(n_1) \}, so 1005 * it is realized that 2^(B*n_0) <= {b_1}^{2*n_1}, where b_1 is a 1006 * given radix number. And then, we have n_1 <= (B*n_0) / 1007 * (2*log_2(b_1)), therefore n_1 is given by ceil((B*n_0) / 1008 * (2*log_2(b_1))). 1009 */ 1010static long 1011big2str_find_n1(VALUE x, int base) 1012{ 1013 static const double log_2[] = { 1014 1.0, 1.58496250072116, 2.0, 1015 2.32192809488736, 2.58496250072116, 2.8073549220576, 1016 3.0, 3.16992500144231, 3.32192809488736, 1017 3.4594316186373, 3.58496250072116, 3.70043971814109, 1018 3.8073549220576, 3.90689059560852, 4.0, 1019 4.08746284125034, 4.16992500144231, 4.24792751344359, 1020 4.32192809488736, 4.39231742277876, 4.4594316186373, 1021 4.52356195605701, 4.58496250072116, 4.64385618977472, 1022 4.70043971814109, 4.75488750216347, 4.8073549220576, 1023 4.85798099512757, 4.90689059560852, 4.95419631038688, 1024 5.0, 5.04439411935845, 5.08746284125034, 1025 5.12928301694497, 5.16992500144231 1026 }; 1027 long bits; 1028 1029 if (base < 2 || 36 < base) 1030 rb_bug("invalid radix %d", base); 1031 1032 if (FIXNUM_P(x)) { 1033 bits = (SIZEOF_LONG*CHAR_BIT - 1)/2 + 1; 1034 } 1035 else if (BIGZEROP(x)) { 1036 return 0; 1037 } 1038 else if (RBIGNUM_LEN(x) >= LONG_MAX/BITSPERDIG) { 1039 rb_raise(rb_eRangeError, "bignum too big to convert into `string'"); 1040 } 1041 else { 1042 bits = BITSPERDIG*RBIGNUM_LEN(x); 1043 } 1044 1045 /* @shyouhei note: vvvvvvvvvvvvv this cast is suspicious. But I believe it is OK, because if that cast loses data, this x value is too big, and should have raised RangeError. */ 1046 return (long)ceil(((double)bits)/log_2[base - 2]); 1047} 1048 1049static long 1050big2str_orig(VALUE x, int base, char* ptr, long len, long hbase, int trim) 1051{ 1052 long i = RBIGNUM_LEN(x), j = len; 1053 BDIGIT* ds = BDIGITS(x); 1054 1055 while (i && j > 0) { 1056 long k = i; 1057 BDIGIT_DBL num = 0; 1058 1059 while (k--) { /* x / hbase */ 1060 num = BIGUP(num) + ds[k]; 1061 ds[k] = (BDIGIT)(num / hbase); 1062 num %= hbase; 1063 } 1064 if (trim && ds[i-1] == 0) i--; 1065 k = SIZEOF_BDIGITS; 1066 while (k--) { 1067 ptr[--j] = ruby_digitmap[num % base]; 1068 num /= base; 1069 if (j <= 0) break; 1070 if (trim && i == 0 && num == 0) break; 1071 } 1072 } 1073 if (trim) { 1074 while (j < len && ptr[j] == '0') j++; 1075 MEMMOVE(ptr, ptr + j, char, len - j); 1076 len -= j; 1077 } 1078 return len; 1079} 1080 1081static long 1082big2str_karatsuba(VALUE x, int base, char* ptr, 1083 long n1, long len, long hbase, int trim) 1084{ 1085 long lh, ll, m1; 1086 VALUE b, q, r; 1087 1088 if (BIGZEROP(x)) { 1089 if (trim) return 0; 1090 else { 1091 memset(ptr, '0', len); 1092 return len; 1093 } 1094 } 1095 1096 if (n1 <= KARATSUBA_DIGITS) { 1097 return big2str_orig(x, base, ptr, len, hbase, trim); 1098 } 1099 1100 b = power_cache_get_power(base, n1, &m1); 1101 bigdivmod(x, b, &q, &r); 1102 lh = big2str_karatsuba(q, base, ptr, (len - m1)/2, 1103 len - m1, hbase, trim); 1104 rb_big_resize(q, 0); 1105 ll = big2str_karatsuba(r, base, ptr + lh, m1/2, 1106 m1, hbase, !lh && trim); 1107 rb_big_resize(r, 0); 1108 1109 return lh + ll; 1110} 1111 1112VALUE 1113rb_big2str0(VALUE x, int base, int trim) 1114{ 1115 int off; 1116 VALUE ss, xx; 1117 long n1, n2, len, hbase; 1118 char* ptr; 1119 1120 if (FIXNUM_P(x)) { 1121 return rb_fix2str(x, base); 1122 } 1123 if (BIGZEROP(x)) { 1124 return rb_usascii_str_new2("0"); 1125 } 1126 1127 if (base < 2 || 36 < base) 1128 rb_raise(rb_eArgError, "invalid radix %d", base); 1129 1130 n2 = big2str_find_n1(x, base); 1131 n1 = (n2 + 1) / 2; 1132 ss = rb_usascii_str_new(0, n2 + 1); /* plus one for sign */ 1133 ptr = RSTRING_PTR(ss); 1134 ptr[0] = RBIGNUM_SIGN(x) ? '+' : '-'; 1135 1136 hbase = base*base; 1137#if SIZEOF_BDIGITS > 2 1138 hbase *= hbase; 1139#endif 1140 off = !(trim && RBIGNUM_SIGN(x)); /* erase plus sign if trim */ 1141 xx = rb_big_clone(x); 1142 RBIGNUM_SET_SIGN(xx, 1); 1143 if (n1 <= KARATSUBA_DIGITS) { 1144 len = off + big2str_orig(xx, base, ptr + off, n2, hbase, trim); 1145 } 1146 else { 1147 len = off + big2str_karatsuba(xx, base, ptr + off, n1, 1148 n2, hbase, trim); 1149 } 1150 rb_big_resize(xx, 0); 1151 1152 ptr[len] = '\0'; 1153 rb_str_resize(ss, len); 1154 1155 return ss; 1156} 1157 1158VALUE 1159rb_big2str(VALUE x, int base) 1160{ 1161 return rb_big2str0(x, base, 1); 1162} 1163 1164/* 1165 * call-seq: 1166 * big.to_s(base=10) -> string 1167 * 1168 * Returns a string containing the representation of <i>big</i> radix 1169 * <i>base</i> (2 through 36). 1170 * 1171 * 12345654321.to_s #=> "12345654321" 1172 * 12345654321.to_s(2) #=> "1011011111110110111011110000110001" 1173 * 12345654321.to_s(8) #=> "133766736061" 1174 * 12345654321.to_s(16) #=> "2dfdbbc31" 1175 * 78546939656932.to_s(36) #=> "rubyrules" 1176 */ 1177 1178static VALUE 1179rb_big_to_s(int argc, VALUE *argv, VALUE x) 1180{ 1181 int base; 1182 1183 if (argc == 0) base = 10; 1184 else { 1185 VALUE b; 1186 1187 rb_scan_args(argc, argv, "01", &b); 1188 base = NUM2INT(b); 1189 } 1190 return rb_big2str(x, base); 1191} 1192 1193static VALUE 1194big2ulong(VALUE x, const char *type, int check) 1195{ 1196 long len = RBIGNUM_LEN(x); 1197 BDIGIT_DBL num; 1198 BDIGIT *ds; 1199 1200 if (len > DIGSPERLONG) { 1201 if (check) 1202 rb_raise(rb_eRangeError, "bignum too big to convert into `%s'", type); 1203 len = DIGSPERLONG; 1204 } 1205 ds = BDIGITS(x); 1206 num = 0; 1207 while (len--) { 1208 num = BIGUP(num); 1209 num += ds[len]; 1210 } 1211 return (VALUE)num; 1212} 1213 1214VALUE 1215rb_big2ulong_pack(VALUE x) 1216{ 1217 VALUE num = big2ulong(x, "unsigned long", FALSE); 1218 if (!RBIGNUM_SIGN(x)) { 1219 return (VALUE)(-(SIGNED_VALUE)num); 1220 } 1221 return num; 1222} 1223 1224VALUE 1225rb_big2ulong(VALUE x) 1226{ 1227 VALUE num = big2ulong(x, "unsigned long", TRUE); 1228 1229 if (RBIGNUM_POSITIVE_P(x)) { 1230 return num; 1231 } 1232 else { 1233 if (num <= LONG_MAX) 1234 return -(long)num; 1235 if (num == 1+(unsigned long)(-(LONG_MIN+1))) 1236 return LONG_MIN; 1237 rb_raise(rb_eRangeError, "bignum out of range of unsigned long"); 1238 } 1239 return num; 1240} 1241 1242SIGNED_VALUE 1243rb_big2long(VALUE x) 1244{ 1245 VALUE num = big2ulong(x, "long", TRUE); 1246 1247 if (RBIGNUM_POSITIVE_P(x)) { 1248 if (LONG_MAX < num) 1249 rb_raise(rb_eRangeError, "bignum too big to convert into `long'"); 1250 return num; 1251 } 1252 else { 1253 if (num <= LONG_MAX) 1254 return -(long)num; 1255 if (num == 1+(unsigned long)(-(LONG_MIN+1))) 1256 return LONG_MIN; 1257 rb_raise(rb_eRangeError, "bignum too big to convert into `long'"); 1258 } 1259} 1260 1261#if HAVE_LONG_LONG 1262 1263static unsigned LONG_LONG 1264big2ull(VALUE x, const char *type) 1265{ 1266 long len = RBIGNUM_LEN(x); 1267 BDIGIT_DBL num; 1268 BDIGIT *ds; 1269 1270 if (len > SIZEOF_LONG_LONG/SIZEOF_BDIGITS) 1271 rb_raise(rb_eRangeError, "bignum too big to convert into `%s'", type); 1272 ds = BDIGITS(x); 1273 num = 0; 1274 while (len--) { 1275 num = BIGUP(num); 1276 num += ds[len]; 1277 } 1278 return num; 1279} 1280 1281unsigned LONG_LONG 1282rb_big2ull(VALUE x) 1283{ 1284 unsigned LONG_LONG num = big2ull(x, "unsigned long long"); 1285 1286 if (RBIGNUM_POSITIVE_P(x)) { 1287 return num; 1288 } 1289 else { 1290 if (num <= LLONG_MAX) 1291 return -(LONG_LONG)num; 1292 if (num == 1+(unsigned LONG_LONG)(-(LLONG_MIN+1))) 1293 return LLONG_MIN; 1294 rb_raise(rb_eRangeError, "bignum out of range of unsigned long long"); 1295 } 1296 return num; 1297} 1298 1299LONG_LONG 1300rb_big2ll(VALUE x) 1301{ 1302 unsigned LONG_LONG num = big2ull(x, "long long"); 1303 1304 if (RBIGNUM_POSITIVE_P(x)) { 1305 if (LLONG_MAX < num) 1306 rb_raise(rb_eRangeError, "bignum too big to convert into `long long'"); 1307 return num; 1308 } 1309 else { 1310 if (num <= LLONG_MAX) 1311 return -(LONG_LONG)num; 1312 if (num == 1+(unsigned LONG_LONG)(-(LLONG_MIN+1))) 1313 return LLONG_MIN; 1314 rb_raise(rb_eRangeError, "bignum too big to convert into `long long'"); 1315 } 1316} 1317 1318#endif /* HAVE_LONG_LONG */ 1319 1320static VALUE 1321dbl2big(double d) 1322{ 1323 long i = 0; 1324 BDIGIT c; 1325 BDIGIT *digits; 1326 VALUE z; 1327 double u = (d < 0)?-d:d; 1328 1329 if (isinf(d)) { 1330 rb_raise(rb_eFloatDomainError, d < 0 ? "-Infinity" : "Infinity"); 1331 } 1332 if (isnan(d)) { 1333 rb_raise(rb_eFloatDomainError, "NaN"); 1334 } 1335 1336 while (!POSFIXABLE(u) || 0 != (long)u) { 1337 u /= (double)(BIGRAD); 1338 i++; 1339 } 1340 z = bignew(i, d>=0); 1341 digits = BDIGITS(z); 1342 while (i--) { 1343 u *= BIGRAD; 1344 c = (BDIGIT)u; 1345 u -= c; 1346 digits[i] = c; 1347 } 1348 1349 return z; 1350} 1351 1352VALUE 1353rb_dbl2big(double d) 1354{ 1355 return bignorm(dbl2big(d)); 1356} 1357 1358static int 1359nlz(BDIGIT x) 1360{ 1361 BDIGIT y; 1362 int n = BITSPERDIG; 1363#if BITSPERDIG > 64 1364 y = x >> 64; if (y) {n -= 64; x = y;} 1365#endif 1366#if BITSPERDIG > 32 1367 y = x >> 32; if (y) {n -= 32; x = y;} 1368#endif 1369#if BITSPERDIG > 16 1370 y = x >> 16; if (y) {n -= 16; x = y;} 1371#endif 1372 y = x >> 8; if (y) {n -= 8; x = y;} 1373 y = x >> 4; if (y) {n -= 4; x = y;} 1374 y = x >> 2; if (y) {n -= 2; x = y;} 1375 y = x >> 1; if (y) {return n - 2;} 1376 return n - x; 1377} 1378 1379static double 1380big2dbl(VALUE x) 1381{ 1382 double d = 0.0; 1383 long i = (bigtrunc(x), RBIGNUM_LEN(x)), lo = 0, bits; 1384 BDIGIT *ds = BDIGITS(x), dl; 1385 1386 if (i) { 1387 bits = i * BITSPERDIG - nlz(ds[i-1]); 1388 if (bits > DBL_MANT_DIG+DBL_MAX_EXP) { 1389 d = HUGE_VAL; 1390 } 1391 else { 1392 if (bits > DBL_MANT_DIG+1) 1393 lo = (bits -= DBL_MANT_DIG+1) / BITSPERDIG; 1394 else 1395 bits = 0; 1396 while (--i > lo) { 1397 d = ds[i] + BIGRAD*d; 1398 } 1399 dl = ds[i]; 1400 if (bits && (dl & (1UL << (bits %= BITSPERDIG)))) { 1401 int carry = dl & ~(~(BDIGIT)0 << bits); 1402 if (!carry) { 1403 while (i-- > 0) { 1404 if ((carry = ds[i]) != 0) break; 1405 } 1406 } 1407 if (carry) { 1408 dl &= (BDIGIT)~0 << bits; 1409 dl += (BDIGIT)1 << bits; 1410 if (!dl) d += 1; 1411 } 1412 } 1413 d = dl + BIGRAD*d; 1414 if (lo) { 1415 if (lo > INT_MAX / BITSPERDIG) 1416 d = HUGE_VAL; 1417 else if (lo < INT_MIN / BITSPERDIG) 1418 d = 0.0; 1419 else 1420 d = ldexp(d, (int)(lo * BITSPERDIG)); 1421 } 1422 } 1423 } 1424 if (!RBIGNUM_SIGN(x)) d = -d; 1425 return d; 1426} 1427 1428double 1429rb_big2dbl(VALUE x) 1430{ 1431 double d = big2dbl(x); 1432 1433 if (isinf(d)) { 1434 rb_warning("Bignum out of Float range"); 1435 if (d < 0.0) 1436 d = -HUGE_VAL; 1437 else 1438 d = HUGE_VAL; 1439 } 1440 return d; 1441} 1442 1443/* 1444 * call-seq: 1445 * big.to_f -> float 1446 * 1447 * Converts <i>big</i> to a <code>Float</code>. If <i>big</i> doesn't 1448 * fit in a <code>Float</code>, the result is infinity. 1449 * 1450 */ 1451 1452static VALUE 1453rb_big_to_f(VALUE x) 1454{ 1455 return DBL2NUM(rb_big2dbl(x)); 1456} 1457 1458VALUE 1459rb_integer_float_cmp(VALUE x, VALUE y) 1460{ 1461 double yd = RFLOAT_VALUE(y); 1462 double yi, yf; 1463 VALUE rel; 1464 1465 if (isnan(yd)) 1466 return Qnil; 1467 if (isinf(yd)) { 1468 if (yd > 0.0) return INT2FIX(-1); 1469 else return INT2FIX(1); 1470 } 1471 yf = modf(yd, &yi); 1472 if (FIXNUM_P(x)) { 1473#if SIZEOF_LONG * CHAR_BIT < DBL_MANT_DIG /* assume FLT_RADIX == 2 */ 1474 double xd = (double)FIX2LONG(x); 1475 if (xd < yd) 1476 return INT2FIX(-1); 1477 if (xd > yd) 1478 return INT2FIX(1); 1479 return INT2FIX(0); 1480#else 1481 long xl, yl; 1482 if (yi < FIXNUM_MIN) 1483 return INT2FIX(1); 1484 if (FIXNUM_MAX+1 <= yi) 1485 return INT2FIX(-1); 1486 xl = FIX2LONG(x); 1487 yl = (long)yi; 1488 if (xl < yl) 1489 return INT2FIX(-1); 1490 if (xl > yl) 1491 return INT2FIX(1); 1492 if (yf < 0.0) 1493 return INT2FIX(1); 1494 if (0.0 < yf) 1495 return INT2FIX(-1); 1496 return INT2FIX(0); 1497#endif 1498 } 1499 y = rb_dbl2big(yi); 1500 rel = rb_big_cmp(x, y); 1501 if (yf == 0.0 || rel != INT2FIX(0)) 1502 return rel; 1503 if (yf < 0.0) 1504 return INT2FIX(1); 1505 return INT2FIX(-1); 1506} 1507 1508VALUE 1509rb_integer_float_eq(VALUE x, VALUE y) 1510{ 1511 double yd = RFLOAT_VALUE(y); 1512 double yi, yf; 1513 1514 if (isnan(yd) || isinf(yd)) 1515 return Qfalse; 1516 yf = modf(yd, &yi); 1517 if (yf != 0) 1518 return Qfalse; 1519 if (FIXNUM_P(x)) { 1520#if SIZEOF_LONG * CHAR_BIT < DBL_MANT_DIG /* assume FLT_RADIX == 2 */ 1521 double xd = (double)FIX2LONG(x); 1522 if (xd != yd) 1523 return Qfalse; 1524 return Qtrue; 1525#else 1526 long xl, yl; 1527 if (yi < LONG_MIN || LONG_MAX < yi) 1528 return Qfalse; 1529 xl = FIX2LONG(x); 1530 yl = (long)yi; 1531 if (xl != yl) 1532 return Qfalse; 1533 return Qtrue; 1534#endif 1535 } 1536 y = rb_dbl2big(yi); 1537 return rb_big_eq(x, y); 1538} 1539 1540/* 1541 * call-seq: 1542 * big <=> numeric -> -1, 0, +1 or nil 1543 * 1544 * Comparison---Returns -1, 0, or +1 depending on whether +big+ is 1545 * less than, equal to, or greater than +numeric+. This is the 1546 * basis for the tests in Comparable. 1547 * 1548 * +nil+ is returned if the two values are incomparable. 1549 * 1550 */ 1551 1552VALUE 1553rb_big_cmp(VALUE x, VALUE y) 1554{ 1555 long xlen = RBIGNUM_LEN(x); 1556 BDIGIT *xds, *yds; 1557 1558 switch (TYPE(y)) { 1559 case T_FIXNUM: 1560 y = rb_int2big(FIX2LONG(y)); 1561 break; 1562 1563 case T_BIGNUM: 1564 break; 1565 1566 case T_FLOAT: 1567 return rb_integer_float_cmp(x, y); 1568 1569 default: 1570 return rb_num_coerce_cmp(x, y, rb_intern("<=>")); 1571 } 1572 1573 if (RBIGNUM_SIGN(x) > RBIGNUM_SIGN(y)) return INT2FIX(1); 1574 if (RBIGNUM_SIGN(x) < RBIGNUM_SIGN(y)) return INT2FIX(-1); 1575 if (xlen < RBIGNUM_LEN(y)) 1576 return (RBIGNUM_SIGN(x)) ? INT2FIX(-1) : INT2FIX(1); 1577 if (xlen > RBIGNUM_LEN(y)) 1578 return (RBIGNUM_SIGN(x)) ? INT2FIX(1) : INT2FIX(-1); 1579 1580 xds = BDIGITS(x); 1581 yds = BDIGITS(y); 1582 1583 while (xlen-- && (xds[xlen]==yds[xlen])); 1584 if (-1 == xlen) return INT2FIX(0); 1585 return (xds[xlen] > yds[xlen]) ? 1586 (RBIGNUM_SIGN(x) ? INT2FIX(1) : INT2FIX(-1)) : 1587 (RBIGNUM_SIGN(x) ? INT2FIX(-1) : INT2FIX(1)); 1588} 1589 1590enum big_op_t { 1591 big_op_gt, 1592 big_op_ge, 1593 big_op_lt, 1594 big_op_le 1595}; 1596 1597static VALUE 1598big_op(VALUE x, VALUE y, enum big_op_t op) 1599{ 1600 VALUE rel; 1601 int n; 1602 1603 switch (TYPE(y)) { 1604 case T_FIXNUM: 1605 case T_BIGNUM: 1606 rel = rb_big_cmp(x, y); 1607 break; 1608 1609 case T_FLOAT: 1610 rel = rb_integer_float_cmp(x, y); 1611 break; 1612 1613 default: 1614 { 1615 ID id = 0; 1616 switch (op) { 1617 case big_op_gt: id = '>'; break; 1618 case big_op_ge: id = rb_intern(">="); break; 1619 case big_op_lt: id = '<'; break; 1620 case big_op_le: id = rb_intern("<="); break; 1621 } 1622 return rb_num_coerce_relop(x, y, id); 1623 } 1624 } 1625 1626 if (NIL_P(rel)) return Qfalse; 1627 n = FIX2INT(rel); 1628 1629 switch (op) { 1630 case big_op_gt: return n > 0 ? Qtrue : Qfalse; 1631 case big_op_ge: return n >= 0 ? Qtrue : Qfalse; 1632 case big_op_lt: return n < 0 ? Qtrue : Qfalse; 1633 case big_op_le: return n <= 0 ? Qtrue : Qfalse; 1634 } 1635 return Qundef; 1636} 1637 1638/* 1639 * call-seq: 1640 * big > real -> true or false 1641 * 1642 * Returns <code>true</code> if the value of <code>big</code> is 1643 * greater than that of <code>real</code>. 1644 */ 1645 1646static VALUE 1647big_gt(VALUE x, VALUE y) 1648{ 1649 return big_op(x, y, big_op_gt); 1650} 1651 1652/* 1653 * call-seq: 1654 * big >= real -> true or false 1655 * 1656 * Returns <code>true</code> if the value of <code>big</code> is 1657 * greater than or equal to that of <code>real</code>. 1658 */ 1659 1660static VALUE 1661big_ge(VALUE x, VALUE y) 1662{ 1663 return big_op(x, y, big_op_ge); 1664} 1665 1666/* 1667 * call-seq: 1668 * big < real -> true or false 1669 * 1670 * Returns <code>true</code> if the value of <code>big</code> is 1671 * less than that of <code>real</code>. 1672 */ 1673 1674static VALUE 1675big_lt(VALUE x, VALUE y) 1676{ 1677 return big_op(x, y, big_op_lt); 1678} 1679 1680/* 1681 * call-seq: 1682 * big <= real -> true or false 1683 * 1684 * Returns <code>true</code> if the value of <code>big</code> is 1685 * less than or equal to that of <code>real</code>. 1686 */ 1687 1688static VALUE 1689big_le(VALUE x, VALUE y) 1690{ 1691 return big_op(x, y, big_op_le); 1692} 1693 1694/* 1695 * call-seq: 1696 * big == obj -> true or false 1697 * 1698 * Returns <code>true</code> only if <i>obj</i> has the same value 1699 * as <i>big</i>. Contrast this with <code>Bignum#eql?</code>, which 1700 * requires <i>obj</i> to be a <code>Bignum</code>. 1701 * 1702 * 68719476736 == 68719476736.0 #=> true 1703 */ 1704 1705VALUE 1706rb_big_eq(VALUE x, VALUE y) 1707{ 1708 switch (TYPE(y)) { 1709 case T_FIXNUM: 1710 y = rb_int2big(FIX2LONG(y)); 1711 break; 1712 case T_BIGNUM: 1713 break; 1714 case T_FLOAT: 1715 return rb_integer_float_eq(x, y); 1716 default: 1717 return rb_equal(y, x); 1718 } 1719 if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y)) return Qfalse; 1720 if (RBIGNUM_LEN(x) != RBIGNUM_LEN(y)) return Qfalse; 1721 if (MEMCMP(BDIGITS(x),BDIGITS(y),BDIGIT,RBIGNUM_LEN(y)) != 0) return Qfalse; 1722 return Qtrue; 1723} 1724 1725/* 1726 * call-seq: 1727 * big.eql?(obj) -> true or false 1728 * 1729 * Returns <code>true</code> only if <i>obj</i> is a 1730 * <code>Bignum</code> with the same value as <i>big</i>. Contrast this 1731 * with <code>Bignum#==</code>, which performs type conversions. 1732 * 1733 * 68719476736.eql?(68719476736.0) #=> false 1734 */ 1735 1736VALUE 1737rb_big_eql(VALUE x, VALUE y) 1738{ 1739 if (!RB_TYPE_P(y, T_BIGNUM)) return Qfalse; 1740 if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y)) return Qfalse; 1741 if (RBIGNUM_LEN(x) != RBIGNUM_LEN(y)) return Qfalse; 1742 if (MEMCMP(BDIGITS(x),BDIGITS(y),BDIGIT,RBIGNUM_LEN(y)) != 0) return Qfalse; 1743 return Qtrue; 1744} 1745 1746/* 1747 * call-seq: 1748 * -big -> integer 1749 * 1750 * Unary minus (returns an integer whose value is 0-big) 1751 */ 1752 1753VALUE 1754rb_big_uminus(VALUE x) 1755{ 1756 VALUE z = rb_big_clone(x); 1757 1758 RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(x)); 1759 1760 return bignorm(z); 1761} 1762 1763/* 1764 * call-seq: 1765 * ~big -> integer 1766 * 1767 * Inverts the bits in big. As Bignums are conceptually infinite 1768 * length, the result acts as if it had an infinite number of one 1769 * bits to the left. In hex representations, this is displayed 1770 * as two periods to the left of the digits. 1771 * 1772 * sprintf("%X", ~0x1122334455) #=> "..FEEDDCCBBAA" 1773 */ 1774 1775static VALUE 1776rb_big_neg(VALUE x) 1777{ 1778 VALUE z = rb_big_clone(x); 1779 BDIGIT *ds; 1780 long i; 1781 1782 if (!RBIGNUM_SIGN(x)) get2comp(z); 1783 ds = BDIGITS(z); 1784 i = RBIGNUM_LEN(x); 1785 if (!i) return INT2FIX(~(SIGNED_VALUE)0); 1786 while (i--) { 1787 ds[i] = ~ds[i]; 1788 } 1789 RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(z)); 1790 if (RBIGNUM_SIGN(x)) get2comp(z); 1791 1792 return bignorm(z); 1793} 1794 1795static void 1796bigsub_core(BDIGIT *xds, long xn, BDIGIT *yds, long yn, BDIGIT *zds, long zn) 1797{ 1798 BDIGIT_DBL_SIGNED num; 1799 long i; 1800 1801 for (i = 0, num = 0; i < yn; i++) { 1802 num += (BDIGIT_DBL_SIGNED)xds[i] - yds[i]; 1803 zds[i] = BIGLO(num); 1804 num = BIGDN(num); 1805 } 1806 while (num && i < xn) { 1807 num += xds[i]; 1808 zds[i++] = BIGLO(num); 1809 num = BIGDN(num); 1810 } 1811 while (i < xn) { 1812 zds[i] = xds[i]; 1813 i++; 1814 } 1815 assert(i <= zn); 1816 while (i < zn) { 1817 zds[i++] = 0; 1818 } 1819} 1820 1821static VALUE 1822bigsub(VALUE x, VALUE y) 1823{ 1824 VALUE z = 0; 1825 long i = RBIGNUM_LEN(x); 1826 BDIGIT *xds, *yds; 1827 1828 /* if x is smaller than y, swap */ 1829 if (RBIGNUM_LEN(x) < RBIGNUM_LEN(y)) { 1830 z = x; x = y; y = z; /* swap x y */ 1831 } 1832 else if (RBIGNUM_LEN(x) == RBIGNUM_LEN(y)) { 1833 xds = BDIGITS(x); 1834 yds = BDIGITS(y); 1835 while (i > 0) { 1836 i--; 1837 if (xds[i] > yds[i]) { 1838 break; 1839 } 1840 if (xds[i] < yds[i]) { 1841 z = x; x = y; y = z; /* swap x y */ 1842 break; 1843 } 1844 } 1845 } 1846 1847 z = bignew(RBIGNUM_LEN(x), z==0); 1848 bigsub_core(BDIGITS(x), RBIGNUM_LEN(x), 1849 BDIGITS(y), RBIGNUM_LEN(y), 1850 BDIGITS(z), RBIGNUM_LEN(z)); 1851 1852 return z; 1853} 1854 1855static VALUE bigadd_int(VALUE x, long y); 1856 1857static VALUE 1858bigsub_int(VALUE x, long y0) 1859{ 1860 VALUE z; 1861 BDIGIT *xds, *zds; 1862 long xn; 1863 BDIGIT_DBL_SIGNED num; 1864 long i, y; 1865 1866 y = y0; 1867 xds = BDIGITS(x); 1868 xn = RBIGNUM_LEN(x); 1869 1870 z = bignew(xn, RBIGNUM_SIGN(x)); 1871 zds = BDIGITS(z); 1872 1873#if SIZEOF_BDIGITS == SIZEOF_LONG 1874 num = (BDIGIT_DBL_SIGNED)xds[0] - y; 1875 if (xn == 1 && num < 0) { 1876 RBIGNUM_SET_SIGN(z, !RBIGNUM_SIGN(x)); 1877 zds[0] = (BDIGIT)-num; 1878 RB_GC_GUARD(x); 1879 return bignorm(z); 1880 } 1881 zds[0] = BIGLO(num); 1882 num = BIGDN(num); 1883 i = 1; 1884#else 1885 num = 0; 1886 for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) { 1887 num += (BDIGIT_DBL_SIGNED)xds[i] - BIGLO(y); 1888 zds[i] = BIGLO(num); 1889 num = BIGDN(num); 1890 y = BIGDN(y); 1891 } 1892#endif 1893 while (num && i < xn) { 1894 num += xds[i]; 1895 zds[i++] = BIGLO(num); 1896 num = BIGDN(num); 1897 } 1898 while (i < xn) { 1899 zds[i] = xds[i]; 1900 i++; 1901 } 1902 if (num < 0) { 1903 z = bigsub(x, rb_int2big(y0)); 1904 } 1905 RB_GC_GUARD(x); 1906 return bignorm(z); 1907} 1908 1909static VALUE 1910bigadd_int(VALUE x, long y) 1911{ 1912 VALUE z; 1913 BDIGIT *xds, *zds; 1914 long xn, zn; 1915 BDIGIT_DBL num; 1916 long i; 1917 1918 xds = BDIGITS(x); 1919 xn = RBIGNUM_LEN(x); 1920 1921 if (xn < 2) { 1922 zn = 3; 1923 } 1924 else { 1925 zn = xn + 1; 1926 } 1927 z = bignew(zn, RBIGNUM_SIGN(x)); 1928 zds = BDIGITS(z); 1929 1930#if SIZEOF_BDIGITS == SIZEOF_LONG 1931 num = (BDIGIT_DBL)xds[0] + y; 1932 zds[0] = BIGLO(num); 1933 num = BIGDN(num); 1934 i = 1; 1935#else 1936 num = 0; 1937 for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) { 1938 num += (BDIGIT_DBL)xds[i] + BIGLO(y); 1939 zds[i] = BIGLO(num); 1940 num = BIGDN(num); 1941 y = BIGDN(y); 1942 } 1943#endif 1944 while (num && i < xn) { 1945 num += xds[i]; 1946 zds[i++] = BIGLO(num); 1947 num = BIGDN(num); 1948 } 1949 if (num) zds[i++] = (BDIGIT)num; 1950 else while (i < xn) { 1951 zds[i] = xds[i]; 1952 i++; 1953 } 1954 assert(i <= zn); 1955 while (i < zn) { 1956 zds[i++] = 0; 1957 } 1958 RB_GC_GUARD(x); 1959 return bignorm(z); 1960} 1961 1962static void 1963bigadd_core(BDIGIT *xds, long xn, BDIGIT *yds, long yn, BDIGIT *zds, long zn) 1964{ 1965 BDIGIT_DBL num = 0; 1966 long i; 1967 1968 if (xn > yn) { 1969 BDIGIT *tds; 1970 tds = xds; xds = yds; yds = tds; 1971 i = xn; xn = yn; yn = i; 1972 } 1973 1974 i = 0; 1975 while (i < xn) { 1976 num += (BDIGIT_DBL)xds[i] + yds[i]; 1977 zds[i++] = BIGLO(num); 1978 num = BIGDN(num); 1979 } 1980 while (num && i < yn) { 1981 num += yds[i]; 1982 zds[i++] = BIGLO(num); 1983 num = BIGDN(num); 1984 } 1985 while (i < yn) { 1986 zds[i] = yds[i]; 1987 i++; 1988 } 1989 if (num) zds[i++] = (BDIGIT)num; 1990 assert(i <= zn); 1991 while (i < zn) { 1992 zds[i++] = 0; 1993 } 1994} 1995 1996static VALUE 1997bigadd(VALUE x, VALUE y, int sign) 1998{ 1999 VALUE z; 2000 long len; 2001 2002 sign = (sign == RBIGNUM_SIGN(y)); 2003 if (RBIGNUM_SIGN(x) != sign) { 2004 if (sign) return bigsub(y, x); 2005 return bigsub(x, y); 2006 } 2007 2008 if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) { 2009 len = RBIGNUM_LEN(x) + 1; 2010 } 2011 else { 2012 len = RBIGNUM_LEN(y) + 1; 2013 } 2014 z = bignew(len, sign); 2015 2016 bigadd_core(BDIGITS(x), RBIGNUM_LEN(x), 2017 BDIGITS(y), RBIGNUM_LEN(y), 2018 BDIGITS(z), RBIGNUM_LEN(z)); 2019 2020 return z; 2021} 2022 2023/* 2024 * call-seq: 2025 * big + other -> Numeric 2026 * 2027 * Adds big and other, returning the result. 2028 */ 2029 2030VALUE 2031rb_big_plus(VALUE x, VALUE y) 2032{ 2033 long n; 2034 2035 switch (TYPE(y)) { 2036 case T_FIXNUM: 2037 n = FIX2LONG(y); 2038 if ((n > 0) != RBIGNUM_SIGN(x)) { 2039 if (n < 0) { 2040 n = -n; 2041 } 2042 return bigsub_int(x, n); 2043 } 2044 if (n < 0) { 2045 n = -n; 2046 } 2047 return bigadd_int(x, n); 2048 2049 case T_BIGNUM: 2050 return bignorm(bigadd(x, y, 1)); 2051 2052 case T_FLOAT: 2053 return DBL2NUM(rb_big2dbl(x) + RFLOAT_VALUE(y)); 2054 2055 default: 2056 return rb_num_coerce_bin(x, y, '+'); 2057 } 2058} 2059 2060/* 2061 * call-seq: 2062 * big - other -> Numeric 2063 * 2064 * Subtracts other from big, returning the result. 2065 */ 2066 2067VALUE 2068rb_big_minus(VALUE x, VALUE y) 2069{ 2070 long n; 2071 2072 switch (TYPE(y)) { 2073 case T_FIXNUM: 2074 n = FIX2LONG(y); 2075 if ((n > 0) != RBIGNUM_SIGN(x)) { 2076 if (n < 0) { 2077 n = -n; 2078 } 2079 return bigadd_int(x, n); 2080 } 2081 if (n < 0) { 2082 n = -n; 2083 } 2084 return bigsub_int(x, n); 2085 2086 case T_BIGNUM: 2087 return bignorm(bigadd(x, y, 0)); 2088 2089 case T_FLOAT: 2090 return DBL2NUM(rb_big2dbl(x) - RFLOAT_VALUE(y)); 2091 2092 default: 2093 return rb_num_coerce_bin(x, y, '-'); 2094 } 2095} 2096 2097static long 2098big_real_len(VALUE x) 2099{ 2100 long i = RBIGNUM_LEN(x); 2101 BDIGIT *xds = BDIGITS(x); 2102 while (--i && !xds[i]); 2103 return i + 1; 2104} 2105 2106static VALUE 2107bigmul1_single(VALUE x, VALUE y) 2108{ 2109 BDIGIT_DBL n; 2110 VALUE z = bignew(2, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y)); 2111 BDIGIT *xds, *yds, *zds; 2112 2113 xds = BDIGITS(x); 2114 yds = BDIGITS(y); 2115 zds = BDIGITS(z); 2116 2117 n = (BDIGIT_DBL)xds[0] * yds[0]; 2118 zds[0] = BIGLO(n); 2119 zds[1] = (BDIGIT)BIGDN(n); 2120 2121 return z; 2122} 2123 2124static VALUE 2125bigmul1_normal(VALUE x, VALUE y) 2126{ 2127 long xl = RBIGNUM_LEN(x), yl = RBIGNUM_LEN(y), i, j = xl + yl + 1; 2128 BDIGIT_DBL n = 0; 2129 VALUE z = bignew(j, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y)); 2130 BDIGIT *xds, *yds, *zds; 2131 2132 xds = BDIGITS(x); 2133 yds = BDIGITS(y); 2134 zds = BDIGITS(z); 2135 while (j--) zds[j] = 0; 2136 for (i = 0; i < xl; i++) { 2137 BDIGIT_DBL dd; 2138 dd = xds[i]; 2139 if (dd == 0) continue; 2140 n = 0; 2141 for (j = 0; j < yl; j++) { 2142 BDIGIT_DBL ee = n + (BDIGIT_DBL)dd * yds[j]; 2143 n = zds[i + j] + ee; 2144 if (ee) zds[i + j] = BIGLO(n); 2145 n = BIGDN(n); 2146 } 2147 if (n) { 2148 zds[i + j] = (BDIGIT)n; 2149 } 2150 } 2151 rb_thread_check_ints(); 2152 return z; 2153} 2154 2155static VALUE bigmul0(VALUE x, VALUE y); 2156 2157/* balancing multiplication by slicing larger argument */ 2158static VALUE 2159bigmul1_balance(VALUE x, VALUE y) 2160{ 2161 VALUE z, t1, t2; 2162 long i, xn, yn, r, n; 2163 BDIGIT *yds, *zds, *t1ds; 2164 2165 xn = RBIGNUM_LEN(x); 2166 yn = RBIGNUM_LEN(y); 2167 assert(2 * xn <= yn || 3 * xn <= 2*(yn+2)); 2168 2169 z = bignew(xn + yn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y)); 2170 t1 = bignew(xn, 1); 2171 2172 yds = BDIGITS(y); 2173 zds = BDIGITS(z); 2174 t1ds = BDIGITS(t1); 2175 2176 for (i = 0; i < xn + yn; i++) zds[i] = 0; 2177 2178 n = 0; 2179 while (yn > 0) { 2180 r = xn > yn ? yn : xn; 2181 MEMCPY(t1ds, yds + n, BDIGIT, r); 2182 RBIGNUM_SET_LEN(t1, r); 2183 t2 = bigmul0(x, t1); 2184 bigadd_core(zds + n, RBIGNUM_LEN(z) - n, 2185 BDIGITS(t2), big_real_len(t2), 2186 zds + n, RBIGNUM_LEN(z) - n); 2187 yn -= r; 2188 n += r; 2189 } 2190 2191 return z; 2192} 2193 2194/* split a bignum into high and low bignums */ 2195static void 2196big_split(VALUE v, long n, volatile VALUE *ph, volatile VALUE *pl) 2197{ 2198 long hn = 0, ln = RBIGNUM_LEN(v); 2199 VALUE h, l; 2200 BDIGIT *vds = BDIGITS(v); 2201 2202 if (ln > n) { 2203 hn = ln - n; 2204 ln = n; 2205 } 2206 2207 if (!hn) { 2208 h = rb_uint2big(0); 2209 } 2210 else { 2211 while (--hn && !vds[hn + ln]); 2212 h = bignew(hn += 2, 1); 2213 MEMCPY(BDIGITS(h), vds + ln, BDIGIT, hn - 1); 2214 BDIGITS(h)[hn - 1] = 0; /* margin for carry */ 2215 } 2216 2217 while (--ln && !vds[ln]); 2218 l = bignew(ln += 2, 1); 2219 MEMCPY(BDIGITS(l), vds, BDIGIT, ln - 1); 2220 BDIGITS(l)[ln - 1] = 0; /* margin for carry */ 2221 2222 *pl = l; 2223 *ph = h; 2224} 2225 2226/* multiplication by karatsuba method */ 2227static VALUE 2228bigmul1_karatsuba(VALUE x, VALUE y) 2229{ 2230 long i, n, xn, yn, t1n, t2n; 2231 VALUE xh, xl, yh, yl, z, t1, t2, t3; 2232 BDIGIT *zds; 2233 2234 xn = RBIGNUM_LEN(x); 2235 yn = RBIGNUM_LEN(y); 2236 n = yn / 2; 2237 big_split(x, n, &xh, &xl); 2238 if (x == y) { 2239 yh = xh; yl = xl; 2240 } 2241 else big_split(y, n, &yh, &yl); 2242 2243 /* x = xh * b + xl 2244 * y = yh * b + yl 2245 * 2246 * Karatsuba method: 2247 * x * y = z2 * b^2 + z1 * b + z0 2248 * where 2249 * z2 = xh * yh 2250 * z0 = xl * yl 2251 * z1 = (xh + xl) * (yh + yl) - z2 - z0 2252 * 2253 * ref: http://en.wikipedia.org/wiki/Karatsuba_algorithm 2254 */ 2255 2256 /* allocate a result bignum */ 2257 z = bignew(xn + yn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y)); 2258 zds = BDIGITS(z); 2259 2260 /* t1 <- xh * yh */ 2261 t1 = bigmul0(xh, yh); 2262 t1n = big_real_len(t1); 2263 2264 /* copy t1 into high bytes of the result (z2) */ 2265 MEMCPY(zds + 2 * n, BDIGITS(t1), BDIGIT, t1n); 2266 for (i = 2 * n + t1n; i < xn + yn; i++) zds[i] = 0; 2267 2268 if (!BIGZEROP(xl) && !BIGZEROP(yl)) { 2269 /* t2 <- xl * yl */ 2270 t2 = bigmul0(xl, yl); 2271 t2n = big_real_len(t2); 2272 2273 /* copy t2 into low bytes of the result (z0) */ 2274 MEMCPY(zds, BDIGITS(t2), BDIGIT, t2n); 2275 for (i = t2n; i < 2 * n; i++) zds[i] = 0; 2276 } 2277 else { 2278 t2 = Qundef; 2279 t2n = 0; 2280 2281 /* copy 0 into low bytes of the result (z0) */ 2282 for (i = 0; i < 2 * n; i++) zds[i] = 0; 2283 } 2284 2285 /* xh <- xh + xl */ 2286 if (RBIGNUM_LEN(xl) > RBIGNUM_LEN(xh)) { 2287 t3 = xl; xl = xh; xh = t3; 2288 } 2289 /* xh has a margin for carry */ 2290 bigadd_core(BDIGITS(xh), RBIGNUM_LEN(xh), 2291 BDIGITS(xl), RBIGNUM_LEN(xl), 2292 BDIGITS(xh), RBIGNUM_LEN(xh)); 2293 2294 /* yh <- yh + yl */ 2295 if (x != y) { 2296 if (RBIGNUM_LEN(yl) > RBIGNUM_LEN(yh)) { 2297 t3 = yl; yl = yh; yh = t3; 2298 } 2299 /* yh has a margin for carry */ 2300 bigadd_core(BDIGITS(yh), RBIGNUM_LEN(yh), 2301 BDIGITS(yl), RBIGNUM_LEN(yl), 2302 BDIGITS(yh), RBIGNUM_LEN(yh)); 2303 } 2304 else yh = xh; 2305 2306 /* t3 <- xh * yh */ 2307 t3 = bigmul0(xh, yh); 2308 2309 i = xn + yn - n; 2310 /* subtract t1 from t3 */ 2311 bigsub_core(BDIGITS(t3), big_real_len(t3), BDIGITS(t1), t1n, BDIGITS(t3), big_real_len(t3)); 2312 2313 /* subtract t2 from t3; t3 is now the middle term of the product */ 2314 if (t2 != Qundef) bigsub_core(BDIGITS(t3), big_real_len(t3), BDIGITS(t2), t2n, BDIGITS(t3), big_real_len(t3)); 2315 2316 /* add t3 to middle bytes of the result (z1) */ 2317 bigadd_core(zds + n, i, BDIGITS(t3), big_real_len(t3), zds + n, i); 2318 2319 return z; 2320} 2321 2322static void 2323biglsh_bang(BDIGIT *xds, long xn, unsigned long shift) 2324{ 2325 long const s1 = shift/BITSPERDIG; 2326 int const s2 = (int)(shift%BITSPERDIG); 2327 int const s3 = BITSPERDIG-s2; 2328 BDIGIT* zds; 2329 BDIGIT num; 2330 long i; 2331 if (s1 >= xn) { 2332 MEMZERO(xds, BDIGIT, xn); 2333 return; 2334 } 2335 zds = xds + xn - 1; 2336 xn -= s1 + 1; 2337 num = xds[xn]<<s2; 2338 while (0 < xn) { 2339 *zds-- = num | xds[--xn]>>s3; 2340 num = xds[xn]<<s2; 2341 } 2342 assert(xds <= zds); 2343 *zds = num; 2344 for (i = s1; i > 0; --i) 2345 *zds-- = 0; 2346} 2347 2348static void 2349bigrsh_bang(BDIGIT* xds, long xn, unsigned long shift) 2350{ 2351 long s1 = shift/BITSPERDIG; 2352 int s2 = (int)(shift%BITSPERDIG); 2353 int s3 = BITSPERDIG - s2; 2354 int i; 2355 BDIGIT num; 2356 BDIGIT* zds; 2357 if (s1 >= xn) { 2358 MEMZERO(xds, BDIGIT, xn); 2359 return; 2360 } 2361 2362 i = 0; 2363 zds = xds + s1; 2364 num = *zds++>>s2; 2365 while (i < xn - s1 - 1) { 2366 xds[i++] = (BDIGIT)(*zds<<s3) | num; 2367 num = *zds++>>s2; 2368 } 2369 assert(i < xn); 2370 xds[i] = num; 2371 MEMZERO(xds + xn - s1, BDIGIT, s1); 2372} 2373 2374static void 2375big_split3(VALUE v, long n, volatile VALUE* p0, volatile VALUE* p1, volatile VALUE* p2) 2376{ 2377 VALUE v0, v12, v1, v2; 2378 2379 big_split(v, n, &v12, &v0); 2380 big_split(v12, n, &v2, &v1); 2381 2382 *p0 = bigtrunc(v0); 2383 *p1 = bigtrunc(v1); 2384 *p2 = bigtrunc(v2); 2385} 2386 2387static VALUE big_lshift(VALUE, unsigned long); 2388static VALUE big_rshift(VALUE, unsigned long); 2389static VALUE bigdivrem(VALUE, VALUE, volatile VALUE*, volatile VALUE*); 2390 2391static VALUE 2392bigmul1_toom3(VALUE x, VALUE y) 2393{ 2394 long n, xn, yn, zn; 2395 VALUE x0, x1, x2, y0, y1, y2; 2396 VALUE u0, u1, u2, u3, u4, v1, v2, v3; 2397 VALUE z0, z1, z2, z3, z4, z, t; 2398 BDIGIT* zds; 2399 2400 xn = RBIGNUM_LEN(x); 2401 yn = RBIGNUM_LEN(y); 2402 assert(xn <= yn); /* assume y >= x */ 2403 2404 n = (yn + 2) / 3; 2405 big_split3(x, n, &x0, &x1, &x2); 2406 if (x == y) { 2407 y0 = x0; y1 = x1; y2 = x2; 2408 } 2409 else big_split3(y, n, &y0, &y1, &y2); 2410 2411 /* 2412 * ref. http://en.wikipedia.org/wiki/Toom%E2%80%93Cook_multiplication 2413 * 2414 * x(b) = x0 * b^0 + x1 * b^1 + x2 * b^2 2415 * y(b) = y0 * b^0 + y1 * b^1 + y2 * b^2 2416 * 2417 * z(b) = x(b) * y(b) 2418 * z(b) = z0 * b^0 + z1 * b^1 + z2 * b^2 + z3 * b^3 + z4 * b^4 2419 * where: 2420 * z0 = x0 * y0 2421 * z1 = x0 * y1 + x1 * y0 2422 * z2 = x0 * y2 + x1 * y1 + x2 * y0 2423 * z3 = x1 * y2 + x2 * y1 2424 * z4 = x2 * y2 2425 * 2426 * Toom3 method (a.k.a. Toom-Cook method): 2427 * (Step1) calculating 5 points z(b0), z(b1), z(b2), z(b3), z(b4), 2428 * where: 2429 * b0 = 0, b1 = 1, b2 = -1, b3 = -2, b4 = inf, 2430 * z(0) = x(0) * y(0) = x0 * y0 2431 * z(1) = x(1) * y(1) = (x0 + x1 + x2) * (y0 + y1 + y2) 2432 * z(-1) = x(-1) * y(-1) = (x0 - x1 + x2) * (y0 - y1 + y2) 2433 * z(-2) = x(-2) * y(-2) = (x0 - 2 * (x1 - 2 * x2)) * (y0 - 2 * (y1 - 2 * y2)) 2434 * z(inf) = x(inf) * y(inf) = x2 * y2 2435 * 2436 * (Step2) interpolating z0, z1, z2, z3, z4, and z5. 2437 * 2438 * (Step3) Substituting base value into b of the polynomial z(b), 2439 */ 2440 2441 /* 2442 * [Step1] calculating 5 points z(b0), z(b1), z(b2), z(b3), z(b4) 2443 */ 2444 2445 /* u1 <- x0 + x2 */ 2446 u1 = bigtrunc(bigadd(x0, x2, 1)); 2447 2448 /* x(-1) : u2 <- u1 - x1 = x0 - x1 + x2 */ 2449 u2 = bigtrunc(bigsub(u1, x1)); 2450 2451 /* x(1) : u1 <- u1 + x1 = x0 + x1 + x2 */ 2452 u1 = bigtrunc(bigadd(u1, x1, 1)); 2453 2454 /* x(-2) : u3 <- 2 * (u2 + x2) - x0 = x0 - 2 * (x1 - 2 * x2) */ 2455 u3 = bigadd(u2, x2, 1); 2456 if (BDIGITS(u3)[RBIGNUM_LEN(u3)-1] & BIGRAD_HALF) { 2457 rb_big_resize(u3, RBIGNUM_LEN(u3) + 1); 2458 BDIGITS(u3)[RBIGNUM_LEN(u3)-1] = 0; 2459 } 2460 biglsh_bang(BDIGITS(u3), RBIGNUM_LEN(u3), 1); 2461 u3 = bigtrunc(bigadd(bigtrunc(u3), x0, 0)); 2462 2463 if (x == y) { 2464 v1 = u1; v2 = u2; v3 = u3; 2465 } 2466 else { 2467 /* v1 <- y0 + y2 */ 2468 v1 = bigtrunc(bigadd(y0, y2, 1)); 2469 2470 /* y(-1) : v2 <- v1 - y1 = y0 - y1 + y2 */ 2471 v2 = bigtrunc(bigsub(v1, y1)); 2472 2473 /* y(1) : v1 <- v1 + y1 = y0 + y1 + y2 */ 2474 v1 = bigtrunc(bigadd(v1, y1, 1)); 2475 2476 /* y(-2) : v3 <- 2 * (v2 + y2) - y0 = y0 - 2 * (y1 - 2 * y2) */ 2477 v3 = bigadd(v2, y2, 1); 2478 if (BDIGITS(v3)[RBIGNUM_LEN(v3)-1] & BIGRAD_HALF) { 2479 rb_big_resize(v3, RBIGNUM_LEN(v3) + 1); 2480 BDIGITS(v3)[RBIGNUM_LEN(v3)-1] = 0; 2481 } 2482 biglsh_bang(BDIGITS(v3), RBIGNUM_LEN(v3), 1); 2483 v3 = bigtrunc(bigadd(bigtrunc(v3), y0, 0)); 2484 } 2485 2486 /* z(0) : u0 <- x0 * y0 */ 2487 u0 = bigtrunc(bigmul0(x0, y0)); 2488 2489 /* z(1) : u1 <- u1 * v1 */ 2490 u1 = bigtrunc(bigmul0(u1, v1)); 2491 2492 /* z(-1) : u2 <- u2 * v2 */ 2493 u2 = bigtrunc(bigmul0(u2, v2)); 2494 2495 /* z(-2) : u3 <- u3 * v3 */ 2496 u3 = bigtrunc(bigmul0(u3, v3)); 2497 2498 /* z(inf) : u4 <- x2 * y2 */ 2499 u4 = bigtrunc(bigmul0(x2, y2)); 2500 2501 /* for GC */ 2502 v1 = v2 = v3 = Qnil; 2503 2504 /* 2505 * [Step2] interpolating z0, z1, z2, z3, z4, and z5. 2506 */ 2507 2508 /* z0 <- z(0) == u0 */ 2509 z0 = u0; 2510 2511 /* z4 <- z(inf) == u4 */ 2512 z4 = u4; 2513 2514 /* z3 <- (z(-2) - z(1)) / 3 == (u3 - u1) / 3 */ 2515 z3 = bigadd(u3, u1, 0); 2516 bigdivrem(z3, big_three, &z3, NULL); /* TODO: optimize */ 2517 bigtrunc(z3); 2518 2519 /* z1 <- (z(1) - z(-1)) / 2 == (u1 - u2) / 2 */ 2520 z1 = bigtrunc(bigadd(u1, u2, 0)); 2521 bigrsh_bang(BDIGITS(z1), RBIGNUM_LEN(z1), 1); 2522 2523 /* z2 <- z(-1) - z(0) == u2 - u0 */ 2524 z2 = bigtrunc(bigadd(u2, u0, 0)); 2525 2526 /* z3 <- (z2 - z3) / 2 + 2 * z(inf) == (z2 - z3) / 2 + 2 * u4 */ 2527 z3 = bigtrunc(bigadd(z2, z3, 0)); 2528 bigrsh_bang(BDIGITS(z3), RBIGNUM_LEN(z3), 1); 2529 t = big_lshift(u4, 1); /* TODO: combining with next addition */ 2530 z3 = bigtrunc(bigadd(z3, t, 1)); 2531 2532 /* z2 <- z2 + z1 - z(inf) == z2 + z1 - u4 */ 2533 z2 = bigtrunc(bigadd(z2, z1, 1)); 2534 z2 = bigtrunc(bigadd(z2, u4, 0)); 2535 2536 /* z1 <- z1 - z3 */ 2537 z1 = bigtrunc(bigadd(z1, z3, 0)); 2538 2539 /* 2540 * [Step3] Substituting base value into b of the polynomial z(b), 2541 */ 2542 2543 zn = 6*n + 1; 2544 z = bignew(zn, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y)); 2545 zds = BDIGITS(z); 2546 MEMCPY(zds, BDIGITS(z0), BDIGIT, RBIGNUM_LEN(z0)); 2547 MEMZERO(zds + RBIGNUM_LEN(z0), BDIGIT, zn - RBIGNUM_LEN(z0)); 2548 bigadd_core(zds + n, zn - n, BDIGITS(z1), big_real_len(z1), zds + n, zn - n); 2549 bigadd_core(zds + 2*n, zn - 2*n, BDIGITS(z2), big_real_len(z2), zds + 2*n, zn - 2*n); 2550 bigadd_core(zds + 3*n, zn - 3*n, BDIGITS(z3), big_real_len(z3), zds + 3*n, zn - 3*n); 2551 bigadd_core(zds + 4*n, zn - 4*n, BDIGITS(z4), big_real_len(z4), zds + 4*n, zn - 4*n); 2552 z = bignorm(z); 2553 2554 return bignorm(z); 2555} 2556 2557/* efficient squaring (2 times faster than normal multiplication) 2558 * ref: Handbook of Applied Cryptography, Algorithm 14.16 2559 * http://www.cacr.math.uwaterloo.ca/hac/about/chap14.pdf 2560 */ 2561static VALUE 2562bigsqr_fast(VALUE x) 2563{ 2564 long len = RBIGNUM_LEN(x), i, j; 2565 VALUE z = bignew(2 * len + 1, 1); 2566 BDIGIT *xds = BDIGITS(x), *zds = BDIGITS(z); 2567 BDIGIT_DBL c, v, w; 2568 2569 for (i = 2 * len + 1; i--; ) zds[i] = 0; 2570 for (i = 0; i < len; i++) { 2571 v = (BDIGIT_DBL)xds[i]; 2572 if (!v) continue; 2573 c = (BDIGIT_DBL)zds[i + i] + v * v; 2574 zds[i + i] = BIGLO(c); 2575 c = BIGDN(c); 2576 v *= 2; 2577 for (j = i + 1; j < len; j++) { 2578 w = (BDIGIT_DBL)xds[j]; 2579 c += (BDIGIT_DBL)zds[i + j] + BIGLO(v) * w; 2580 zds[i + j] = BIGLO(c); 2581 c = BIGDN(c); 2582 if (BIGDN(v)) c += w; 2583 } 2584 if (c) { 2585 c += (BDIGIT_DBL)zds[i + len]; 2586 zds[i + len] = BIGLO(c); 2587 c = BIGDN(c); 2588 } 2589 if (c) zds[i + len + 1] += (BDIGIT)c; 2590 } 2591 return z; 2592} 2593 2594#define KARATSUBA_MUL_DIGITS 70 2595#define TOOM3_MUL_DIGITS 150 2596 2597 2598/* determine whether a bignum is sparse or not by random sampling */ 2599static inline VALUE 2600big_sparse_p(VALUE x) 2601{ 2602 long c = 0, n = RBIGNUM_LEN(x); 2603 2604 if ( BDIGITS(x)[rb_genrand_ulong_limited(n / 2) + n / 4]) c++; 2605 if (c <= 1 && BDIGITS(x)[rb_genrand_ulong_limited(n / 2) + n / 4]) c++; 2606 if (c <= 1 && BDIGITS(x)[rb_genrand_ulong_limited(n / 2) + n / 4]) c++; 2607 2608 return (c <= 1) ? Qtrue : Qfalse; 2609} 2610 2611static VALUE 2612bigmul0(VALUE x, VALUE y) 2613{ 2614 long xn, yn; 2615 2616 xn = RBIGNUM_LEN(x); 2617 yn = RBIGNUM_LEN(y); 2618 2619 /* make sure that y is longer than x */ 2620 if (xn > yn) { 2621 VALUE t; 2622 long tn; 2623 t = x; x = y; y = t; 2624 tn = xn; xn = yn; yn = tn; 2625 } 2626 assert(xn <= yn); 2627 2628 /* normal multiplication when x is small */ 2629 if (xn < KARATSUBA_MUL_DIGITS) { 2630 normal: 2631 if (x == y) return bigsqr_fast(x); 2632 if (xn == 1 && yn == 1) return bigmul1_single(x, y); 2633 return bigmul1_normal(x, y); 2634 } 2635 2636 /* normal multiplication when x or y is a sparse bignum */ 2637 if (big_sparse_p(x)) goto normal; 2638 if (big_sparse_p(y)) return bigmul1_normal(y, x); 2639 2640 /* balance multiplication by slicing y when x is much smaller than y */ 2641 if (2 * xn <= yn) return bigmul1_balance(x, y); 2642 2643 if (xn < TOOM3_MUL_DIGITS) { 2644 /* multiplication by karatsuba method */ 2645 return bigmul1_karatsuba(x, y); 2646 } 2647 else if (3*xn <= 2*(yn + 2)) 2648 return bigmul1_balance(x, y); 2649 return bigmul1_toom3(x, y); 2650} 2651 2652/* 2653 * call-seq: 2654 * big * other -> Numeric 2655 * 2656 * Multiplies big and other, returning the result. 2657 */ 2658 2659VALUE 2660rb_big_mul(VALUE x, VALUE y) 2661{ 2662 switch (TYPE(y)) { 2663 case T_FIXNUM: 2664 y = rb_int2big(FIX2LONG(y)); 2665 break; 2666 2667 case T_BIGNUM: 2668 break; 2669 2670 case T_FLOAT: 2671 return DBL2NUM(rb_big2dbl(x) * RFLOAT_VALUE(y)); 2672 2673 default: 2674 return rb_num_coerce_bin(x, y, '*'); 2675 } 2676 2677 return bignorm(bigmul0(x, y)); 2678} 2679 2680struct big_div_struct { 2681 long nx, ny, j, nyzero; 2682 BDIGIT *yds, *zds; 2683 volatile VALUE stop; 2684}; 2685 2686static void * 2687bigdivrem1(void *ptr) 2688{ 2689 struct big_div_struct *bds = (struct big_div_struct*)ptr; 2690 long ny = bds->ny; 2691 long i, j; 2692 BDIGIT *yds = bds->yds, *zds = bds->zds; 2693 BDIGIT_DBL t2; 2694 BDIGIT_DBL_SIGNED num; 2695 BDIGIT q; 2696 2697 j = bds->j; 2698 do { 2699 if (bds->stop) { 2700 bds->j = j; 2701 return 0; 2702 } 2703 if (zds[j] == yds[ny-1]) q = (BDIGIT)BIGRAD-1; 2704 else q = (BDIGIT)((BIGUP(zds[j]) + zds[j-1])/yds[ny-1]); 2705 if (q) { 2706 i = bds->nyzero; num = 0; t2 = 0; 2707 do { /* multiply and subtract */ 2708 BDIGIT_DBL ee; 2709 t2 += (BDIGIT_DBL)yds[i] * q; 2710 ee = num - BIGLO(t2); 2711 num = (BDIGIT_DBL)zds[j - ny + i] + ee; 2712 if (ee) zds[j - ny + i] = BIGLO(num); 2713 num = BIGDN(num); 2714 t2 = BIGDN(t2); 2715 } while (++i < ny); 2716 num += zds[j - ny + i] - t2;/* borrow from high digit; don't update */ 2717 while (num) { /* "add back" required */ 2718 i = 0; num = 0; q--; 2719 do { 2720 BDIGIT_DBL ee = num + yds[i]; 2721 num = (BDIGIT_DBL)zds[j - ny + i] + ee; 2722 if (ee) zds[j - ny + i] = BIGLO(num); 2723 num = BIGDN(num); 2724 } while (++i < ny); 2725 num--; 2726 } 2727 } 2728 zds[j] = q; 2729 } while (--j >= ny); 2730 return 0; 2731} 2732 2733static void 2734rb_big_stop(void *ptr) 2735{ 2736 struct big_div_struct *bds = ptr; 2737 bds->stop = Qtrue; 2738} 2739 2740static VALUE 2741bigdivrem(VALUE x, VALUE y, volatile VALUE *divp, volatile VALUE *modp) 2742{ 2743 struct big_div_struct bds; 2744 long nx = RBIGNUM_LEN(x), ny = RBIGNUM_LEN(y); 2745 long i, j; 2746 VALUE z, yy, zz; 2747 BDIGIT *xds, *yds, *zds, *tds; 2748 BDIGIT_DBL t2; 2749 BDIGIT dd, q; 2750 2751 if (BIGZEROP(y)) rb_num_zerodiv(); 2752 xds = BDIGITS(x); 2753 yds = BDIGITS(y); 2754 if (nx < ny || (nx == ny && xds[nx - 1] < yds[ny - 1])) { 2755 if (divp) *divp = rb_int2big(0); 2756 if (modp) *modp = x; 2757 return Qnil; 2758 } 2759 if (ny == 1) { 2760 dd = yds[0]; 2761 z = rb_big_clone(x); 2762 zds = BDIGITS(z); 2763 t2 = 0; i = nx; 2764 while (i--) { 2765 t2 = BIGUP(t2) + zds[i]; 2766 zds[i] = (BDIGIT)(t2 / dd); 2767 t2 %= dd; 2768 } 2769 RBIGNUM_SET_SIGN(z, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y)); 2770 if (modp) { 2771 *modp = rb_uint2big((VALUE)t2); 2772 RBIGNUM_SET_SIGN(*modp, RBIGNUM_SIGN(x)); 2773 } 2774 if (divp) *divp = z; 2775 return Qnil; 2776 } 2777 2778 z = bignew(nx==ny?nx+2:nx+1, RBIGNUM_SIGN(x)==RBIGNUM_SIGN(y)); 2779 zds = BDIGITS(z); 2780 if (nx==ny) zds[nx+1] = 0; 2781 while (!yds[ny-1]) ny--; 2782 2783 dd = 0; 2784 q = yds[ny-1]; 2785 while ((q & (BDIGIT)(1UL<<(BITSPERDIG-1))) == 0) { 2786 q <<= 1UL; 2787 dd++; 2788 } 2789 if (dd) { 2790 yy = rb_big_clone(y); 2791 tds = BDIGITS(yy); 2792 j = 0; 2793 t2 = 0; 2794 while (j<ny) { 2795 t2 += (BDIGIT_DBL)yds[j]<<dd; 2796 tds[j++] = BIGLO(t2); 2797 t2 = BIGDN(t2); 2798 } 2799 yds = tds; 2800 RB_GC_GUARD(y) = yy; 2801 j = 0; 2802 t2 = 0; 2803 while (j<nx) { 2804 t2 += (BDIGIT_DBL)xds[j]<<dd; 2805 zds[j++] = BIGLO(t2); 2806 t2 = BIGDN(t2); 2807 } 2808 zds[j] = (BDIGIT)t2; 2809 } 2810 else { 2811 zds[nx] = 0; 2812 j = nx; 2813 while (j--) zds[j] = xds[j]; 2814 } 2815 2816 bds.nx = nx; 2817 bds.ny = ny; 2818 bds.zds = zds; 2819 bds.yds = yds; 2820 bds.stop = Qfalse; 2821 bds.j = nx==ny?nx+1:nx; 2822 for (bds.nyzero = 0; !yds[bds.nyzero]; bds.nyzero++); 2823 if (nx > 10000 || ny > 10000) { 2824 retry: 2825 bds.stop = Qfalse; 2826 rb_thread_call_without_gvl(bigdivrem1, &bds, rb_big_stop, &bds); 2827 2828 if (bds.stop == Qtrue) { 2829 /* execute trap handler, but exception was not raised. */ 2830 goto retry; 2831 } 2832 } 2833 else { 2834 bigdivrem1(&bds); 2835 } 2836 2837 if (divp) { /* move quotient down in z */ 2838 *divp = zz = rb_big_clone(z); 2839 zds = BDIGITS(zz); 2840 j = (nx==ny ? nx+2 : nx+1) - ny; 2841 for (i = 0;i < j;i++) zds[i] = zds[i+ny]; 2842 if (!zds[i-1]) i--; 2843 RBIGNUM_SET_LEN(zz, i); 2844 } 2845 if (modp) { /* normalize remainder */ 2846 *modp = zz = rb_big_clone(z); 2847 zds = BDIGITS(zz); 2848 while (ny > 1 && !zds[ny-1]) --ny; 2849 if (dd) { 2850 t2 = 0; i = ny; 2851 while (i--) { 2852 t2 = (t2 | zds[i]) >> dd; 2853 q = zds[i]; 2854 zds[i] = BIGLO(t2); 2855 t2 = BIGUP(q); 2856 } 2857 } 2858 if (!zds[ny-1]) ny--; 2859 RBIGNUM_SET_LEN(zz, ny); 2860 RBIGNUM_SET_SIGN(zz, RBIGNUM_SIGN(x)); 2861 } 2862 return z; 2863} 2864 2865static void 2866bigdivmod(VALUE x, VALUE y, volatile VALUE *divp, volatile VALUE *modp) 2867{ 2868 VALUE mod; 2869 2870 bigdivrem(x, y, divp, &mod); 2871 if (RBIGNUM_SIGN(x) != RBIGNUM_SIGN(y) && !BIGZEROP(mod)) { 2872 if (divp) *divp = bigadd(*divp, rb_int2big(1), 0); 2873 if (modp) *modp = bigadd(mod, y, 1); 2874 } 2875 else if (modp) { 2876 *modp = mod; 2877 } 2878} 2879 2880 2881static VALUE 2882rb_big_divide(VALUE x, VALUE y, ID op) 2883{ 2884 VALUE z; 2885 2886 switch (TYPE(y)) { 2887 case T_FIXNUM: 2888 y = rb_int2big(FIX2LONG(y)); 2889 break; 2890 2891 case T_BIGNUM: 2892 break; 2893 2894 case T_FLOAT: 2895 { 2896 if (op == '/') { 2897 return DBL2NUM(rb_big2dbl(x) / RFLOAT_VALUE(y)); 2898 } 2899 else { 2900 double dy = RFLOAT_VALUE(y); 2901 if (dy == 0.0) rb_num_zerodiv(); 2902 return rb_dbl2big(rb_big2dbl(x) / dy); 2903 } 2904 } 2905 2906 default: 2907 return rb_num_coerce_bin(x, y, op); 2908 } 2909 bigdivmod(x, y, &z, 0); 2910 2911 return bignorm(z); 2912} 2913 2914/* 2915 * call-seq: 2916 * big / other -> Numeric 2917 * 2918 * Performs division: the class of the resulting object depends on 2919 * the class of <code>numeric</code> and on the magnitude of the 2920 * result. 2921 */ 2922 2923VALUE 2924rb_big_div(VALUE x, VALUE y) 2925{ 2926 return rb_big_divide(x, y, '/'); 2927} 2928 2929/* 2930 * call-seq: 2931 * big.div(other) -> integer 2932 * 2933 * Performs integer division: returns integer value. 2934 */ 2935 2936VALUE 2937rb_big_idiv(VALUE x, VALUE y) 2938{ 2939 return rb_big_divide(x, y, rb_intern("div")); 2940} 2941 2942/* 2943 * call-seq: 2944 * big % other -> Numeric 2945 * big.modulo(other) -> Numeric 2946 * 2947 * Returns big modulo other. See Numeric.divmod for more 2948 * information. 2949 */ 2950 2951VALUE 2952rb_big_modulo(VALUE x, VALUE y) 2953{ 2954 VALUE z; 2955 2956 switch (TYPE(y)) { 2957 case T_FIXNUM: 2958 y = rb_int2big(FIX2LONG(y)); 2959 break; 2960 2961 case T_BIGNUM: 2962 break; 2963 2964 default: 2965 return rb_num_coerce_bin(x, y, '%'); 2966 } 2967 bigdivmod(x, y, 0, &z); 2968 2969 return bignorm(z); 2970} 2971 2972/* 2973 * call-seq: 2974 * big.remainder(numeric) -> number 2975 * 2976 * Returns the remainder after dividing <i>big</i> by <i>numeric</i>. 2977 * 2978 * -1234567890987654321.remainder(13731) #=> -6966 2979 * -1234567890987654321.remainder(13731.24) #=> -9906.22531493148 2980 */ 2981static VALUE 2982rb_big_remainder(VALUE x, VALUE y) 2983{ 2984 VALUE z; 2985 2986 switch (TYPE(y)) { 2987 case T_FIXNUM: 2988 y = rb_int2big(FIX2LONG(y)); 2989 break; 2990 2991 case T_BIGNUM: 2992 break; 2993 2994 default: 2995 return rb_num_coerce_bin(x, y, rb_intern("remainder")); 2996 } 2997 bigdivrem(x, y, 0, &z); 2998 2999 return bignorm(z); 3000} 3001 3002/* 3003 * call-seq: 3004 * big.divmod(numeric) -> array 3005 * 3006 * See <code>Numeric#divmod</code>. 3007 * 3008 */ 3009VALUE 3010rb_big_divmod(VALUE x, VALUE y) 3011{ 3012 VALUE div, mod; 3013 3014 switch (TYPE(y)) { 3015 case T_FIXNUM: 3016 y = rb_int2big(FIX2LONG(y)); 3017 break; 3018 3019 case T_BIGNUM: 3020 break; 3021 3022 default: 3023 return rb_num_coerce_bin(x, y, rb_intern("divmod")); 3024 } 3025 bigdivmod(x, y, &div, &mod); 3026 3027 return rb_assoc_new(bignorm(div), bignorm(mod)); 3028} 3029 3030static int 3031bdigbitsize(BDIGIT x) 3032{ 3033 int size = 1; 3034 int nb = BITSPERDIG / 2; 3035 BDIGIT bits = (~0 << nb); 3036 3037 if (!x) return 0; 3038 while (x > 1) { 3039 if (x & bits) { 3040 size += nb; 3041 x >>= nb; 3042 } 3043 x &= ~bits; 3044 nb /= 2; 3045 bits >>= nb; 3046 } 3047 3048 return size; 3049} 3050 3051static VALUE big_lshift(VALUE, unsigned long); 3052static VALUE big_rshift(VALUE, unsigned long); 3053 3054static VALUE 3055big_shift(VALUE x, long n) 3056{ 3057 if (n < 0) 3058 return big_lshift(x, (unsigned long)-n); 3059 else if (n > 0) 3060 return big_rshift(x, (unsigned long)n); 3061 return x; 3062} 3063 3064static VALUE 3065big_fdiv(VALUE x, VALUE y) 3066{ 3067#define DBL_BIGDIG ((DBL_MANT_DIG + BITSPERDIG) / BITSPERDIG) 3068 VALUE z; 3069 long l, ex, ey; 3070 int i; 3071 3072 bigtrunc(x); 3073 l = RBIGNUM_LEN(x) - 1; 3074 ex = l * BITSPERDIG; 3075 ex += bdigbitsize(BDIGITS(x)[l]); 3076 ex -= 2 * DBL_BIGDIG * BITSPERDIG; 3077 if (ex) x = big_shift(x, ex); 3078 3079 switch (TYPE(y)) { 3080 case T_FIXNUM: 3081 y = rb_int2big(FIX2LONG(y)); 3082 case T_BIGNUM: 3083 bigtrunc(y); 3084 l = RBIGNUM_LEN(y) - 1; 3085 ey = l * BITSPERDIG; 3086 ey += bdigbitsize(BDIGITS(y)[l]); 3087 ey -= DBL_BIGDIG * BITSPERDIG; 3088 if (ey) y = big_shift(y, ey); 3089 break; 3090 case T_FLOAT: 3091 y = dbl2big(ldexp(frexp(RFLOAT_VALUE(y), &i), DBL_MANT_DIG)); 3092 ey = i - DBL_MANT_DIG; 3093 break; 3094 default: 3095 rb_bug("big_fdiv"); 3096 } 3097 bigdivrem(x, y, &z, 0); 3098 l = ex - ey; 3099#if SIZEOF_LONG > SIZEOF_INT 3100 { 3101 /* Visual C++ can't be here */ 3102 if (l > INT_MAX) return DBL2NUM(INFINITY); 3103 if (l < INT_MIN) return DBL2NUM(0.0); 3104 } 3105#endif 3106 return DBL2NUM(ldexp(big2dbl(z), (int)l)); 3107} 3108 3109/* 3110 * call-seq: 3111 * big.fdiv(numeric) -> float 3112 * 3113 * Returns the floating point result of dividing <i>big</i> by 3114 * <i>numeric</i>. 3115 * 3116 * -1234567890987654321.fdiv(13731) #=> -89910996357705.5 3117 * -1234567890987654321.fdiv(13731.24) #=> -89909424858035.7 3118 * 3119 */ 3120 3121 3122VALUE 3123rb_big_fdiv(VALUE x, VALUE y) 3124{ 3125 double dx, dy; 3126 3127 dx = big2dbl(x); 3128 switch (TYPE(y)) { 3129 case T_FIXNUM: 3130 dy = (double)FIX2LONG(y); 3131 if (isinf(dx)) 3132 return big_fdiv(x, y); 3133 break; 3134 3135 case T_BIGNUM: 3136 dy = rb_big2dbl(y); 3137 if (isinf(dx) || isinf(dy)) 3138 return big_fdiv(x, y); 3139 break; 3140 3141 case T_FLOAT: 3142 dy = RFLOAT_VALUE(y); 3143 if (isnan(dy)) 3144 return y; 3145 if (isinf(dx)) 3146 return big_fdiv(x, y); 3147 break; 3148 3149 default: 3150 return rb_num_coerce_bin(x, y, rb_intern("fdiv")); 3151 } 3152 return DBL2NUM(dx / dy); 3153} 3154 3155static VALUE 3156bigsqr(VALUE x) 3157{ 3158 return bigtrunc(bigmul0(x, x)); 3159} 3160 3161/* 3162 * call-seq: 3163 * big ** exponent -> numeric 3164 * 3165 * Raises _big_ to the _exponent_ power (which may be an integer, float, 3166 * or anything that will coerce to a number). The result may be 3167 * a Fixnum, Bignum, or Float 3168 * 3169 * 123456789 ** 2 #=> 15241578750190521 3170 * 123456789 ** 1.2 #=> 5126464716.09932 3171 * 123456789 ** -2 #=> 6.5610001194102e-17 3172 */ 3173 3174VALUE 3175rb_big_pow(VALUE x, VALUE y) 3176{ 3177 double d; 3178 SIGNED_VALUE yy; 3179 3180 if (y == INT2FIX(0)) return INT2FIX(1); 3181 switch (TYPE(y)) { 3182 case T_FLOAT: 3183 d = RFLOAT_VALUE(y); 3184 if ((!RBIGNUM_SIGN(x) && !BIGZEROP(x)) && d != round(d)) 3185 return rb_funcall(rb_complex_raw1(x), rb_intern("**"), 1, y); 3186 break; 3187 3188 case T_BIGNUM: 3189 rb_warn("in a**b, b may be too big"); 3190 d = rb_big2dbl(y); 3191 break; 3192 3193 case T_FIXNUM: 3194 yy = FIX2LONG(y); 3195 3196 if (yy < 0) 3197 return rb_funcall(rb_rational_raw1(x), rb_intern("**"), 1, y); 3198 else { 3199 VALUE z = 0; 3200 SIGNED_VALUE mask; 3201 const long xlen = RBIGNUM_LEN(x) - 1; 3202 const long xbits = ffs(RBIGNUM_DIGITS(x)[xlen]) + SIZEOF_BDIGITS*BITSPERDIG*xlen; 3203 const long BIGLEN_LIMIT = BITSPERDIG*1024*1024; 3204 3205 if ((xbits > BIGLEN_LIMIT) || (xbits * yy > BIGLEN_LIMIT)) { 3206 rb_warn("in a**b, b may be too big"); 3207 d = (double)yy; 3208 break; 3209 } 3210 for (mask = FIXNUM_MAX + 1; mask; mask >>= 1) { 3211 if (z) z = bigsqr(z); 3212 if (yy & mask) { 3213 z = z ? bigtrunc(bigmul0(z, x)) : x; 3214 } 3215 } 3216 return bignorm(z); 3217 } 3218 /* NOTREACHED */ 3219 break; 3220 3221 default: 3222 return rb_num_coerce_bin(x, y, rb_intern("**")); 3223 } 3224 return DBL2NUM(pow(rb_big2dbl(x), d)); 3225} 3226 3227static VALUE 3228bigand_int(VALUE x, long y) 3229{ 3230 VALUE z; 3231 BDIGIT *xds, *zds; 3232 long xn, zn; 3233 long i; 3234 char sign; 3235 3236 if (y == 0) return INT2FIX(0); 3237 sign = (y > 0); 3238 xds = BDIGITS(x); 3239 zn = xn = RBIGNUM_LEN(x); 3240#if SIZEOF_BDIGITS == SIZEOF_LONG 3241 if (sign) { 3242 y &= xds[0]; 3243 return LONG2NUM(y); 3244 } 3245#endif 3246 3247 z = bignew(zn, RBIGNUM_SIGN(x) || sign); 3248 zds = BDIGITS(z); 3249 3250#if SIZEOF_BDIGITS == SIZEOF_LONG 3251 i = 1; 3252 zds[0] = xds[0] & y; 3253#else 3254 { 3255 BDIGIT_DBL num = y; 3256 3257 for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) { 3258 zds[i] = xds[i] & BIGLO(num); 3259 num = BIGDN(num); 3260 } 3261 } 3262#endif 3263 while (i < xn) { 3264 zds[i] = sign?0:xds[i]; 3265 i++; 3266 } 3267 if (!RBIGNUM_SIGN(z)) get2comp(z); 3268 return bignorm(z); 3269} 3270 3271/* 3272 * call-seq: 3273 * big & numeric -> integer 3274 * 3275 * Performs bitwise +and+ between _big_ and _numeric_. 3276 */ 3277 3278VALUE 3279rb_big_and(VALUE xx, VALUE yy) 3280{ 3281 volatile VALUE x, y, z; 3282 BDIGIT *ds1, *ds2, *zds; 3283 long i, l1, l2; 3284 char sign; 3285 3286 if (!FIXNUM_P(yy) && !RB_TYPE_P(yy, T_BIGNUM)) { 3287 return rb_num_coerce_bit(xx, yy, '&'); 3288 } 3289 3290 x = xx; 3291 y = yy; 3292 3293 if (!RBIGNUM_SIGN(x)) { 3294 x = rb_big_clone(x); 3295 get2comp(x); 3296 } 3297 if (FIXNUM_P(y)) { 3298 return bigand_int(x, FIX2LONG(y)); 3299 } 3300 if (!RBIGNUM_SIGN(y)) { 3301 y = rb_big_clone(y); 3302 get2comp(y); 3303 } 3304 if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) { 3305 l1 = RBIGNUM_LEN(y); 3306 l2 = RBIGNUM_LEN(x); 3307 ds1 = BDIGITS(y); 3308 ds2 = BDIGITS(x); 3309 sign = RBIGNUM_SIGN(y); 3310 } 3311 else { 3312 l1 = RBIGNUM_LEN(x); 3313 l2 = RBIGNUM_LEN(y); 3314 ds1 = BDIGITS(x); 3315 ds2 = BDIGITS(y); 3316 sign = RBIGNUM_SIGN(x); 3317 } 3318 z = bignew(l2, RBIGNUM_SIGN(x) || RBIGNUM_SIGN(y)); 3319 zds = BDIGITS(z); 3320 3321 for (i=0; i<l1; i++) { 3322 zds[i] = ds1[i] & ds2[i]; 3323 } 3324 for (; i<l2; i++) { 3325 zds[i] = sign?0:ds2[i]; 3326 } 3327 if (!RBIGNUM_SIGN(z)) get2comp(z); 3328 return bignorm(z); 3329} 3330 3331static VALUE 3332bigor_int(VALUE x, long y) 3333{ 3334 VALUE z; 3335 BDIGIT *xds, *zds; 3336 long xn, zn; 3337 long i; 3338 char sign; 3339 3340 sign = (y >= 0); 3341 xds = BDIGITS(x); 3342 zn = xn = RBIGNUM_LEN(x); 3343 z = bignew(zn, RBIGNUM_SIGN(x) && sign); 3344 zds = BDIGITS(z); 3345 3346#if SIZEOF_BDIGITS == SIZEOF_LONG 3347 i = 1; 3348 zds[0] = xds[0] | y; 3349#else 3350 { 3351 BDIGIT_DBL num = y; 3352 3353 for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) { 3354 zds[i] = xds[i] | BIGLO(num); 3355 num = BIGDN(num); 3356 } 3357 } 3358#endif 3359 while (i < xn) { 3360 zds[i] = sign?xds[i]:(BDIGIT)(BIGRAD-1); 3361 i++; 3362 } 3363 if (!RBIGNUM_SIGN(z)) get2comp(z); 3364 return bignorm(z); 3365} 3366 3367/* 3368 * call-seq: 3369 * big | numeric -> integer 3370 * 3371 * Performs bitwise +or+ between _big_ and _numeric_. 3372 */ 3373 3374VALUE 3375rb_big_or(VALUE xx, VALUE yy) 3376{ 3377 volatile VALUE x, y, z; 3378 BDIGIT *ds1, *ds2, *zds; 3379 long i, l1, l2; 3380 char sign; 3381 3382 if (!FIXNUM_P(yy) && !RB_TYPE_P(yy, T_BIGNUM)) { 3383 return rb_num_coerce_bit(xx, yy, '|'); 3384 } 3385 3386 x = xx; 3387 y = yy; 3388 3389 if (!RBIGNUM_SIGN(x)) { 3390 x = rb_big_clone(x); 3391 get2comp(x); 3392 } 3393 if (FIXNUM_P(y)) { 3394 return bigor_int(x, FIX2LONG(y)); 3395 } 3396 if (!RBIGNUM_SIGN(y)) { 3397 y = rb_big_clone(y); 3398 get2comp(y); 3399 } 3400 if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) { 3401 l1 = RBIGNUM_LEN(y); 3402 l2 = RBIGNUM_LEN(x); 3403 ds1 = BDIGITS(y); 3404 ds2 = BDIGITS(x); 3405 sign = RBIGNUM_SIGN(y); 3406 } 3407 else { 3408 l1 = RBIGNUM_LEN(x); 3409 l2 = RBIGNUM_LEN(y); 3410 ds1 = BDIGITS(x); 3411 ds2 = BDIGITS(y); 3412 sign = RBIGNUM_SIGN(x); 3413 } 3414 z = bignew(l2, RBIGNUM_SIGN(x) && RBIGNUM_SIGN(y)); 3415 zds = BDIGITS(z); 3416 3417 for (i=0; i<l1; i++) { 3418 zds[i] = ds1[i] | ds2[i]; 3419 } 3420 for (; i<l2; i++) { 3421 zds[i] = sign?ds2[i]:(BDIGIT)(BIGRAD-1); 3422 } 3423 if (!RBIGNUM_SIGN(z)) get2comp(z); 3424 return bignorm(z); 3425} 3426 3427static VALUE 3428bigxor_int(VALUE x, long y) 3429{ 3430 VALUE z; 3431 BDIGIT *xds, *zds; 3432 long xn, zn; 3433 long i; 3434 char sign; 3435 3436 sign = (y >= 0) ? 1 : 0; 3437 xds = BDIGITS(x); 3438 zn = xn = RBIGNUM_LEN(x); 3439 z = bignew(zn, !(RBIGNUM_SIGN(x) ^ sign)); 3440 zds = BDIGITS(z); 3441 3442#if SIZEOF_BDIGITS == SIZEOF_LONG 3443 i = 1; 3444 zds[0] = xds[0] ^ y; 3445#else 3446 { 3447 BDIGIT_DBL num = y; 3448 3449 for (i=0; i<(int)(sizeof(y)/sizeof(BDIGIT)); i++) { 3450 zds[i] = xds[i] ^ BIGLO(num); 3451 num = BIGDN(num); 3452 } 3453 } 3454#endif 3455 while (i < xn) { 3456 zds[i] = sign?xds[i]:~xds[i]; 3457 i++; 3458 } 3459 if (!RBIGNUM_SIGN(z)) get2comp(z); 3460 return bignorm(z); 3461} 3462/* 3463 * call-seq: 3464 * big ^ numeric -> integer 3465 * 3466 * Performs bitwise +exclusive or+ between _big_ and _numeric_. 3467 */ 3468 3469VALUE 3470rb_big_xor(VALUE xx, VALUE yy) 3471{ 3472 volatile VALUE x, y; 3473 VALUE z; 3474 BDIGIT *ds1, *ds2, *zds; 3475 long i, l1, l2; 3476 char sign; 3477 3478 if (!FIXNUM_P(yy) && !RB_TYPE_P(yy, T_BIGNUM)) { 3479 return rb_num_coerce_bit(xx, yy, '^'); 3480 } 3481 3482 x = xx; 3483 y = yy; 3484 3485 if (!RBIGNUM_SIGN(x)) { 3486 x = rb_big_clone(x); 3487 get2comp(x); 3488 } 3489 if (FIXNUM_P(y)) { 3490 return bigxor_int(x, FIX2LONG(y)); 3491 } 3492 if (!RBIGNUM_SIGN(y)) { 3493 y = rb_big_clone(y); 3494 get2comp(y); 3495 } 3496 if (RBIGNUM_LEN(x) > RBIGNUM_LEN(y)) { 3497 l1 = RBIGNUM_LEN(y); 3498 l2 = RBIGNUM_LEN(x); 3499 ds1 = BDIGITS(y); 3500 ds2 = BDIGITS(x); 3501 sign = RBIGNUM_SIGN(y); 3502 } 3503 else { 3504 l1 = RBIGNUM_LEN(x); 3505 l2 = RBIGNUM_LEN(y); 3506 ds1 = BDIGITS(x); 3507 ds2 = BDIGITS(y); 3508 sign = RBIGNUM_SIGN(x); 3509 } 3510 RBIGNUM_SET_SIGN(x, RBIGNUM_SIGN(x)?1:0); 3511 RBIGNUM_SET_SIGN(y, RBIGNUM_SIGN(y)?1:0); 3512 z = bignew(l2, !(RBIGNUM_SIGN(x) ^ RBIGNUM_SIGN(y))); 3513 zds = BDIGITS(z); 3514 3515 for (i=0; i<l1; i++) { 3516 zds[i] = ds1[i] ^ ds2[i]; 3517 } 3518 for (; i<l2; i++) { 3519 zds[i] = sign?ds2[i]:~ds2[i]; 3520 } 3521 if (!RBIGNUM_SIGN(z)) get2comp(z); 3522 3523 return bignorm(z); 3524} 3525 3526static VALUE 3527check_shiftdown(VALUE y, VALUE x) 3528{ 3529 if (!RBIGNUM_LEN(x)) return INT2FIX(0); 3530 if (RBIGNUM_LEN(y) > SIZEOF_LONG / SIZEOF_BDIGITS) { 3531 return RBIGNUM_SIGN(x) ? INT2FIX(0) : INT2FIX(-1); 3532 } 3533 return Qnil; 3534} 3535 3536/* 3537 * call-seq: 3538 * big << numeric -> integer 3539 * 3540 * Shifts big left _numeric_ positions (right if _numeric_ is negative). 3541 */ 3542 3543VALUE 3544rb_big_lshift(VALUE x, VALUE y) 3545{ 3546 long shift; 3547 int neg = 0; 3548 3549 for (;;) { 3550 if (FIXNUM_P(y)) { 3551 shift = FIX2LONG(y); 3552 if (shift < 0) { 3553 neg = 1; 3554 shift = -shift; 3555 } 3556 break; 3557 } 3558 else if (RB_TYPE_P(y, T_BIGNUM)) { 3559 if (!RBIGNUM_SIGN(y)) { 3560 VALUE t = check_shiftdown(y, x); 3561 if (!NIL_P(t)) return t; 3562 neg = 1; 3563 } 3564 shift = big2ulong(y, "long", TRUE); 3565 break; 3566 } 3567 y = rb_to_int(y); 3568 } 3569 3570 x = neg ? big_rshift(x, shift) : big_lshift(x, shift); 3571 return bignorm(x); 3572} 3573 3574static VALUE 3575big_lshift(VALUE x, unsigned long shift) 3576{ 3577 BDIGIT *xds, *zds; 3578 long s1 = shift/BITSPERDIG; 3579 int s2 = (int)(shift%BITSPERDIG); 3580 VALUE z; 3581 BDIGIT_DBL num = 0; 3582 long len, i; 3583 3584 len = RBIGNUM_LEN(x); 3585 z = bignew(len+s1+1, RBIGNUM_SIGN(x)); 3586 zds = BDIGITS(z); 3587 for (i=0; i<s1; i++) { 3588 *zds++ = 0; 3589 } 3590 xds = BDIGITS(x); 3591 for (i=0; i<len; i++) { 3592 num = num | (BDIGIT_DBL)*xds++<<s2; 3593 *zds++ = BIGLO(num); 3594 num = BIGDN(num); 3595 } 3596 *zds = BIGLO(num); 3597 return z; 3598} 3599 3600/* 3601 * call-seq: 3602 * big >> numeric -> integer 3603 * 3604 * Shifts big right _numeric_ positions (left if _numeric_ is negative). 3605 */ 3606 3607VALUE 3608rb_big_rshift(VALUE x, VALUE y) 3609{ 3610 long shift; 3611 int neg = 0; 3612 3613 for (;;) { 3614 if (FIXNUM_P(y)) { 3615 shift = FIX2LONG(y); 3616 if (shift < 0) { 3617 neg = 1; 3618 shift = -shift; 3619 } 3620 break; 3621 } 3622 else if (RB_TYPE_P(y, T_BIGNUM)) { 3623 if (RBIGNUM_SIGN(y)) { 3624 VALUE t = check_shiftdown(y, x); 3625 if (!NIL_P(t)) return t; 3626 } 3627 else { 3628 neg = 1; 3629 } 3630 shift = big2ulong(y, "long", TRUE); 3631 break; 3632 } 3633 y = rb_to_int(y); 3634 } 3635 3636 x = neg ? big_lshift(x, shift) : big_rshift(x, shift); 3637 return bignorm(x); 3638} 3639 3640static VALUE 3641big_rshift(VALUE x, unsigned long shift) 3642{ 3643 BDIGIT *xds, *zds; 3644 long s1 = shift/BITSPERDIG; 3645 int s2 = (int)(shift%BITSPERDIG); 3646 VALUE z; 3647 BDIGIT_DBL num = 0; 3648 long i, j; 3649 volatile VALUE save_x; 3650 3651 if (s1 > RBIGNUM_LEN(x)) { 3652 if (RBIGNUM_SIGN(x)) 3653 return INT2FIX(0); 3654 else 3655 return INT2FIX(-1); 3656 } 3657 if (!RBIGNUM_SIGN(x)) { 3658 x = rb_big_clone(x); 3659 get2comp(x); 3660 } 3661 save_x = x; 3662 xds = BDIGITS(x); 3663 i = RBIGNUM_LEN(x); j = i - s1; 3664 if (j == 0) { 3665 if (RBIGNUM_SIGN(x)) return INT2FIX(0); 3666 else return INT2FIX(-1); 3667 } 3668 z = bignew(j, RBIGNUM_SIGN(x)); 3669 if (!RBIGNUM_SIGN(x)) { 3670 num = ((BDIGIT_DBL)~0) << BITSPERDIG; 3671 } 3672 zds = BDIGITS(z); 3673 while (i--, j--) { 3674 num = (num | xds[i]) >> s2; 3675 zds[j] = BIGLO(num); 3676 num = BIGUP(xds[i]); 3677 } 3678 if (!RBIGNUM_SIGN(x)) { 3679 get2comp(z); 3680 } 3681 RB_GC_GUARD(save_x); 3682 return z; 3683} 3684 3685/* 3686 * call-seq: 3687 * big[n] -> 0, 1 3688 * 3689 * Bit Reference---Returns the <em>n</em>th bit in the (assumed) binary 3690 * representation of <i>big</i>, where <i>big</i>[0] is the least 3691 * significant bit. 3692 * 3693 * a = 9**15 3694 * 50.downto(0) do |n| 3695 * print a[n] 3696 * end 3697 * 3698 * <em>produces:</em> 3699 * 3700 * 000101110110100000111000011110010100111100010111001 3701 * 3702 */ 3703 3704static VALUE 3705rb_big_aref(VALUE x, VALUE y) 3706{ 3707 BDIGIT *xds; 3708 BDIGIT_DBL num; 3709 VALUE shift; 3710 long i, s1, s2; 3711 3712 if (RB_TYPE_P(y, T_BIGNUM)) { 3713 if (!RBIGNUM_SIGN(y)) 3714 return INT2FIX(0); 3715 bigtrunc(y); 3716 if (RBIGNUM_LEN(y) > DIGSPERLONG) { 3717 out_of_range: 3718 return RBIGNUM_SIGN(x) ? INT2FIX(0) : INT2FIX(1); 3719 } 3720 shift = big2ulong(y, "long", FALSE); 3721 } 3722 else { 3723 i = NUM2LONG(y); 3724 if (i < 0) return INT2FIX(0); 3725 shift = (VALUE)i; 3726 } 3727 s1 = shift/BITSPERDIG; 3728 s2 = shift%BITSPERDIG; 3729 3730 if (s1 >= RBIGNUM_LEN(x)) goto out_of_range; 3731 if (!RBIGNUM_SIGN(x)) { 3732 xds = BDIGITS(x); 3733 i = 0; num = 1; 3734 while (num += ~xds[i], ++i <= s1) { 3735 num = BIGDN(num); 3736 } 3737 } 3738 else { 3739 num = BDIGITS(x)[s1]; 3740 } 3741 if (num & ((BDIGIT_DBL)1<<s2)) 3742 return INT2FIX(1); 3743 return INT2FIX(0); 3744} 3745 3746/* 3747 * call-seq: 3748 * big.hash -> fixnum 3749 * 3750 * Compute a hash based on the value of _big_. 3751 */ 3752 3753static VALUE 3754rb_big_hash(VALUE x) 3755{ 3756 st_index_t hash; 3757 3758 hash = rb_memhash(BDIGITS(x), sizeof(BDIGIT)*RBIGNUM_LEN(x)) ^ RBIGNUM_SIGN(x); 3759 return INT2FIX(hash); 3760} 3761 3762/* 3763 * MISSING: documentation 3764 */ 3765 3766static VALUE 3767rb_big_coerce(VALUE x, VALUE y) 3768{ 3769 if (FIXNUM_P(y)) { 3770 y = rb_int2big(FIX2LONG(y)); 3771 } 3772 else if (!RB_TYPE_P(y, T_BIGNUM)) { 3773 rb_raise(rb_eTypeError, "can't coerce %s to Bignum", 3774 rb_obj_classname(y)); 3775 } 3776 return rb_assoc_new(y, x); 3777} 3778 3779/* 3780 * call-seq: 3781 * big.abs -> aBignum 3782 * 3783 * Returns the absolute value of <i>big</i>. 3784 * 3785 * -1234567890987654321.abs #=> 1234567890987654321 3786 */ 3787 3788static VALUE 3789rb_big_abs(VALUE x) 3790{ 3791 if (!RBIGNUM_SIGN(x)) { 3792 x = rb_big_clone(x); 3793 RBIGNUM_SET_SIGN(x, 1); 3794 } 3795 return x; 3796} 3797 3798/* 3799 * call-seq: 3800 * big.size -> integer 3801 * 3802 * Returns the number of bytes in the machine representation of 3803 * <i>big</i>. 3804 * 3805 * (256**10 - 1).size #=> 12 3806 * (256**20 - 1).size #=> 20 3807 * (256**40 - 1).size #=> 40 3808 */ 3809 3810static VALUE 3811rb_big_size(VALUE big) 3812{ 3813 return LONG2FIX(RBIGNUM_LEN(big)*SIZEOF_BDIGITS); 3814} 3815 3816/* 3817 * call-seq: 3818 * big.odd? -> true or false 3819 * 3820 * Returns <code>true</code> if <i>big</i> is an odd number. 3821 */ 3822 3823static VALUE 3824rb_big_odd_p(VALUE num) 3825{ 3826 if (BDIGITS(num)[0] & 1) { 3827 return Qtrue; 3828 } 3829 return Qfalse; 3830} 3831 3832/* 3833 * call-seq: 3834 * big.even? -> true or false 3835 * 3836 * Returns <code>true</code> if <i>big</i> is an even number. 3837 */ 3838 3839static VALUE 3840rb_big_even_p(VALUE num) 3841{ 3842 if (BDIGITS(num)[0] & 1) { 3843 return Qfalse; 3844 } 3845 return Qtrue; 3846} 3847 3848/* 3849 * Bignum objects hold integers outside the range of 3850 * Fixnum. Bignum objects are created 3851 * automatically when integer calculations would otherwise overflow a 3852 * Fixnum. When a calculation involving 3853 * Bignum objects returns a result that will fit in a 3854 * Fixnum, the result is automatically converted. 3855 * 3856 * For the purposes of the bitwise operations and <code>[]</code>, a 3857 * Bignum is treated as if it were an infinite-length 3858 * bitstring with 2's complement representation. 3859 * 3860 * While Fixnum values are immediate, Bignum 3861 * objects are not---assignment and parameter passing work with 3862 * references to objects, not the objects themselves. 3863 * 3864 */ 3865 3866void 3867Init_Bignum(void) 3868{ 3869 rb_cBignum = rb_define_class("Bignum", rb_cInteger); 3870 3871 rb_define_method(rb_cBignum, "to_s", rb_big_to_s, -1); 3872 rb_define_alias(rb_cBignum, "inspect", "to_s"); 3873 rb_define_method(rb_cBignum, "coerce", rb_big_coerce, 1); 3874 rb_define_method(rb_cBignum, "-@", rb_big_uminus, 0); 3875 rb_define_method(rb_cBignum, "+", rb_big_plus, 1); 3876 rb_define_method(rb_cBignum, "-", rb_big_minus, 1); 3877 rb_define_method(rb_cBignum, "*", rb_big_mul, 1); 3878 rb_define_method(rb_cBignum, "/", rb_big_div, 1); 3879 rb_define_method(rb_cBignum, "%", rb_big_modulo, 1); 3880 rb_define_method(rb_cBignum, "div", rb_big_idiv, 1); 3881 rb_define_method(rb_cBignum, "divmod", rb_big_divmod, 1); 3882 rb_define_method(rb_cBignum, "modulo", rb_big_modulo, 1); 3883 rb_define_method(rb_cBignum, "remainder", rb_big_remainder, 1); 3884 rb_define_method(rb_cBignum, "fdiv", rb_big_fdiv, 1); 3885 rb_define_method(rb_cBignum, "**", rb_big_pow, 1); 3886 rb_define_method(rb_cBignum, "&", rb_big_and, 1); 3887 rb_define_method(rb_cBignum, "|", rb_big_or, 1); 3888 rb_define_method(rb_cBignum, "^", rb_big_xor, 1); 3889 rb_define_method(rb_cBignum, "~", rb_big_neg, 0); 3890 rb_define_method(rb_cBignum, "<<", rb_big_lshift, 1); 3891 rb_define_method(rb_cBignum, ">>", rb_big_rshift, 1); 3892 rb_define_method(rb_cBignum, "[]", rb_big_aref, 1); 3893 3894 rb_define_method(rb_cBignum, "<=>", rb_big_cmp, 1); 3895 rb_define_method(rb_cBignum, "==", rb_big_eq, 1); 3896 rb_define_method(rb_cBignum, ">", big_gt, 1); 3897 rb_define_method(rb_cBignum, ">=", big_ge, 1); 3898 rb_define_method(rb_cBignum, "<", big_lt, 1); 3899 rb_define_method(rb_cBignum, "<=", big_le, 1); 3900 rb_define_method(rb_cBignum, "===", rb_big_eq, 1); 3901 rb_define_method(rb_cBignum, "eql?", rb_big_eql, 1); 3902 rb_define_method(rb_cBignum, "hash", rb_big_hash, 0); 3903 rb_define_method(rb_cBignum, "to_f", rb_big_to_f, 0); 3904 rb_define_method(rb_cBignum, "abs", rb_big_abs, 0); 3905 rb_define_method(rb_cBignum, "magnitude", rb_big_abs, 0); 3906 rb_define_method(rb_cBignum, "size", rb_big_size, 0); 3907 rb_define_method(rb_cBignum, "odd?", rb_big_odd_p, 0); 3908 rb_define_method(rb_cBignum, "even?", rb_big_even_p, 0); 3909 3910 power_cache_init(); 3911 3912 big_three = rb_uint2big(3); 3913 rb_gc_register_mark_object(big_three); 3914} 3915