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