vm/bignum.cpp

   1 /*
   2    Copyright (C) 1989-94 Massachusetts Institute of Technology
   3    Portions copyright (C) 2004-2008 Slava Pestov
   4
   5    This material was developed by the Scheme project at the Massachusetts
   6    Institute of Technology, Department of Electrical Engineering and
   7    Computer Science.  Permission to copy and modify this software, to
   8    redistribute either the original software or a modified version, and
   9    to use this software for any purpose is granted, subject to the
  10    following restrictions and understandings.
  11
  12    1. Any copy made of this software must include this copyright notice
  13    in full.
  14
  15    2. Users of this software agree to make their best efforts (a) to
  16    return to the MIT Scheme project any improvements or extensions that
  17    they make, so that these may be included in future releases; and (b)
  18    to inform MIT of noteworthy uses of this software.
  19
  20    3. All materials developed as a consequence of the use of this
  21    software shall duly acknowledge such use, in accordance with the usual
  22    standards of acknowledging credit in academic research.
  23
  24    4. MIT has made no warrantee or representation that the operation of
  25    this software will be error-free, and MIT is under no obligation to
  26    provide any services, by way of maintenance, update, or otherwise.
  27
  28    5. In conjunction with products arising from the use of this material,
  29    there shall be no use of the name of the Massachusetts Institute of
  30    Technology nor of any adaptation thereof in any advertising,
  31    promotional, or sales literature without prior written consent from
  32    MIT in each case. */
  33
  34 /* Changes for Scheme 48:
  35  *  - Converted to ANSI.
  36  *  - Added bitwise operations.
  37  *  - Added s48 to the beginning of all externally visible names.
  38  *  - Cached the bignum representations of -1, 0, and 1.
  39  */
  40
  41 /* Changes for Factor:
  42  *  - Adapt bignumint.h for Factor memory manager
  43  *  - Add more bignum <-> C type conversions
  44  *  - Remove unused functions
  45  *  - Add local variable GC root recording
  46  *  - Remove s48 prefix from function names
  47  *  - Various fixes for Win64
  48  *  - Port to C++
  49  *  - Added bignum_gcd implementation
  50  */
  51
  52 #include "master.hpp"
  53
  54 namespace factor {
  55
  56 /* Exports */
  57
  58 int factor_vm::bignum_equal_p(bignum* x, bignum* y) {
  59   return ((BIGNUM_ZERO_P(x))
  60               ? (BIGNUM_ZERO_P(y))
  61               : ((!(BIGNUM_ZERO_P(y))) &&
  62                  ((BIGNUM_NEGATIVE_P(x)) ? (BIGNUM_NEGATIVE_P(y))
  63                                          : (!(BIGNUM_NEGATIVE_P(y)))) &&
  64                  (bignum_equal_p_unsigned(x, y))));
  65 }
  66
  67 enum bignum_comparison factor_vm::bignum_compare(bignum* x, bignum* y) {
  68   return ((BIGNUM_ZERO_P(x)) ? ((BIGNUM_ZERO_P(y)) ? bignum_comparison_equal
  69                                                    : (BIGNUM_NEGATIVE_P(y))
  70                                     ? bignum_comparison_greater
  71                                     : bignum_comparison_less)
  72                              : (BIGNUM_ZERO_P(y))
  73               ? ((BIGNUM_NEGATIVE_P(x)) ? bignum_comparison_less
  74                                         : bignum_comparison_greater)
  75               : (BIGNUM_NEGATIVE_P(x))
  76               ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_compare_unsigned(y, x))
  77                                         : (bignum_comparison_less))
  78               : ((BIGNUM_NEGATIVE_P(y)) ? (bignum_comparison_greater)
  79                                         : (bignum_compare_unsigned(x, y))));
  80 }
  81
  82 /* Allocates memory */
  83 bignum* factor_vm::bignum_add(bignum* x, bignum* y) {
  84   return (
  85       (BIGNUM_ZERO_P(x)) ? (y) : (BIGNUM_ZERO_P(y))
  86           ? (x)
  87           : ((BIGNUM_NEGATIVE_P(x))
  88                  ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_add_unsigned(x, y, 1))
  89                                            : (bignum_subtract_unsigned(y, x)))
  90                  : ((BIGNUM_NEGATIVE_P(y)) ? (bignum_subtract_unsigned(x, y))
  91                                            : (bignum_add_unsigned(x, y, 0)))));
  92 }
  93
  94 /* Allocates memory */
  95 bignum* factor_vm::bignum_subtract(bignum* x, bignum* y) {
  96   return ((BIGNUM_ZERO_P(x))
  97               ? ((BIGNUM_ZERO_P(y)) ? (y) : (bignum_new_sign(
  98                                                 y, (!(BIGNUM_NEGATIVE_P(y))))))
  99               : ((BIGNUM_ZERO_P(y))
 100                      ? (x)
 101                      : ((BIGNUM_NEGATIVE_P(x))
 102                             ? ((BIGNUM_NEGATIVE_P(y))
 103                                    ? (bignum_subtract_unsigned(y, x))
 104                                    : (bignum_add_unsigned(x, y, 1)))
 105                             : ((BIGNUM_NEGATIVE_P(y))
 106                                    ? (bignum_add_unsigned(x, y, 0))
 107                                    : (bignum_subtract_unsigned(x, y))))));
 108 }
 109
 110 #ifdef _WIN64
 111 bignum *factor_vm::bignum_square(bignum * x)
 112 {
 113     return bignum_multiply(x, x);
 114 }
 115 #else
 116 /* Allocates memory */
 117 bignum *factor_vm::bignum_square(bignum * x)
 118 {
 119     GC_BIGNUM(x);
 120
 121     bignum_length_type length = (BIGNUM_LENGTH (x));
 122     bignum * z = (allot_bignum_zeroed ((length + length), 0));
 123
 124     bignum_digit_type * scan_z = BIGNUM_START_PTR (z);
 125     bignum_digit_type * scan_x = BIGNUM_START_PTR (x);
 126     bignum_digit_type * end_x = scan_x + length;
 127
 128     for (int i = 0; i < length; ++i) {
 129         bignum_twodigit_type carry;
 130         bignum_twodigit_type f = BIGNUM_REF (x, i);
 131         bignum_digit_type *pz = scan_z + (i << 1);
 132         bignum_digit_type *px = scan_x + i + 1;
 133
 134         carry = *pz + f * f;
 135         *pz++ = carry & BIGNUM_DIGIT_MASK;
 136         carry >>= BIGNUM_DIGIT_LENGTH;
 137         BIGNUM_ASSERT (carry <= BIGNUM_DIGIT_MASK);
 138
 139         f <<= 1;
 140         while (px < end_x)
 141         {
 142             carry += *pz + *px++ * f;
 143             *pz++ = carry & BIGNUM_DIGIT_MASK;
 144             carry >>= BIGNUM_DIGIT_LENGTH;
 145             BIGNUM_ASSERT (carry <= (BIGNUM_DIGIT_MASK << 1));
 146         }
 147         if (carry) {
 148             carry += *pz;
 149             *pz++ = carry & BIGNUM_DIGIT_MASK;
 150             carry >>= BIGNUM_DIGIT_LENGTH;
 151         }
 152         if (carry)
 153             *pz += carry & BIGNUM_DIGIT_MASK;
 154         BIGNUM_ASSERT ((carry >> BIGNUM_DIGIT_LENGTH) == 0);
 155     }
 156     return (bignum_trim (z));
 157 }
 158 #endif
 159
 160 /* Allocates memory */
 161 bignum* factor_vm::bignum_multiply(bignum* x, bignum* y) {
 162
 163 #ifndef _WIN64
 164   if (x == y) {
 165     return bignum_square(x);
 166   }
 167 #endif
 168
 169   bignum_length_type x_length = (BIGNUM_LENGTH(x));
 170   bignum_length_type y_length = (BIGNUM_LENGTH(y));
 171   int negative_p = ((BIGNUM_NEGATIVE_P(x)) ? (!(BIGNUM_NEGATIVE_P(y)))
 172                                            : (BIGNUM_NEGATIVE_P(y)));
 173   if (BIGNUM_ZERO_P(x))
 174     return (x);
 175   if (BIGNUM_ZERO_P(y))
 176     return (y);
 177   if (x_length == 1) {
 178     bignum_digit_type digit = (BIGNUM_REF(x, 0));
 179     if (digit == 1)
 180       return (bignum_maybe_new_sign(y, negative_p));
 181     if (digit < BIGNUM_RADIX_ROOT)
 182       return (bignum_multiply_unsigned_small_factor(y, digit, negative_p));
 183   }
 184   if (y_length == 1) {
 185     bignum_digit_type digit = (BIGNUM_REF(y, 0));
 186     if (digit == 1)
 187       return (bignum_maybe_new_sign(x, negative_p));
 188     if (digit < BIGNUM_RADIX_ROOT)
 189       return (bignum_multiply_unsigned_small_factor(x, digit, negative_p));
 190   }
 191   return (bignum_multiply_unsigned(x, y, negative_p));
 192 }
 193
 194 /* Allocates memory */
 195 void factor_vm::bignum_divide(bignum* numerator, bignum* denominator,
 196                               bignum** quotient, bignum** remainder) {
 197   if (BIGNUM_ZERO_P(denominator)) {
 198     divide_by_zero_error();
 199     return;
 200   }
 201   if (BIGNUM_ZERO_P(numerator)) {
 202     (*quotient) = numerator;
 203     (*remainder) = numerator;
 204   } else {
 205     int r_negative_p = (BIGNUM_NEGATIVE_P(numerator));
 206     int q_negative_p =
 207         ((BIGNUM_NEGATIVE_P(denominator)) ? (!r_negative_p) : r_negative_p);
 208     switch (bignum_compare_unsigned(numerator, denominator)) {
 209       case bignum_comparison_equal: {
 210         (*quotient) = (BIGNUM_ONE(q_negative_p));
 211         (*remainder) = (BIGNUM_ZERO());
 212         break;
 213       }
 214       case bignum_comparison_less: {
 215         (*quotient) = (BIGNUM_ZERO());
 216         (*remainder) = numerator;
 217         break;
 218       }
 219       case bignum_comparison_greater: {
 220         if ((BIGNUM_LENGTH(denominator)) == 1) {
 221           bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
 222           if (digit == 1) {
 223             (*quotient) = (bignum_maybe_new_sign(numerator, q_negative_p));
 224             (*remainder) = (BIGNUM_ZERO());
 225             break;
 226           } else if (digit < BIGNUM_RADIX_ROOT) {
 227             bignum_divide_unsigned_small_denominator(numerator, digit, quotient,
 228                                                      remainder, q_negative_p,
 229                                                      r_negative_p);
 230             break;
 231           } else {
 232             bignum_divide_unsigned_medium_denominator(
 233                 numerator, digit, quotient, remainder, q_negative_p,
 234                 r_negative_p);
 235             break;
 236           }
 237         }
 238         bignum_divide_unsigned_large_denominator(
 239             numerator, denominator, quotient, remainder, q_negative_p,
 240             r_negative_p);
 241         break;
 242       }
 243     }
 244   }
 245 }
 246
 247 /* Allocates memory */
 248 bignum* factor_vm::bignum_quotient(bignum* numerator, bignum* denominator) {
 249   if (BIGNUM_ZERO_P(denominator)) {
 250     divide_by_zero_error();
 251     return (BIGNUM_OUT_OF_BAND);
 252   }
 253   if (BIGNUM_ZERO_P(numerator))
 254     return numerator;
 255   {
 256     int q_negative_p =
 257         ((BIGNUM_NEGATIVE_P(denominator)) ? (!(BIGNUM_NEGATIVE_P(numerator)))
 258                                           : (BIGNUM_NEGATIVE_P(numerator)));
 259     switch (bignum_compare_unsigned(numerator, denominator)) {
 260       case bignum_comparison_equal:
 261         return (BIGNUM_ONE(q_negative_p));
 262       case bignum_comparison_less:
 263         return (BIGNUM_ZERO());
 264       case bignum_comparison_greater:
 265       default: /* to appease gcc -Wall */
 266                {
 267         bignum* quotient;
 268         if ((BIGNUM_LENGTH(denominator)) == 1) {
 269           bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
 270           if (digit == 1)
 271             return (bignum_maybe_new_sign(numerator, q_negative_p));
 272           if (digit < BIGNUM_RADIX_ROOT)
 273             bignum_divide_unsigned_small_denominator(
 274                 numerator, digit, (&quotient), ((bignum**)0), q_negative_p, 0);
 275           else
 276             bignum_divide_unsigned_medium_denominator(
 277                 numerator, digit, (&quotient), ((bignum**)0), q_negative_p, 0);
 278         } else
 279           bignum_divide_unsigned_large_denominator(
 280               numerator, denominator, (&quotient), ((bignum**)0), q_negative_p,
 281               0);
 282         return (quotient);
 283       }
 284     }
 285   }
 286 }
 287
 288 /* Allocates memory */
 289 bignum* factor_vm::bignum_remainder(bignum* numerator, bignum* denominator) {
 290   if (BIGNUM_ZERO_P(denominator)) {
 291     divide_by_zero_error();
 292     return (BIGNUM_OUT_OF_BAND);
 293   }
 294   if (BIGNUM_ZERO_P(numerator))
 295     return numerator;
 296   switch (bignum_compare_unsigned(numerator, denominator)) {
 297     case bignum_comparison_equal:
 298       return (BIGNUM_ZERO());
 299     case bignum_comparison_less:
 300       return numerator;
 301     case bignum_comparison_greater:
 302     default: /* to appease gcc -Wall */
 303              {
 304       bignum* remainder;
 305       if ((BIGNUM_LENGTH(denominator)) == 1) {
 306         bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
 307         if (digit == 1)
 308           return (BIGNUM_ZERO());
 309         if (digit < BIGNUM_RADIX_ROOT)
 310           return (bignum_remainder_unsigned_small_denominator(
 311               numerator, digit, (BIGNUM_NEGATIVE_P(numerator))));
 312         bignum_divide_unsigned_medium_denominator(
 313             numerator, digit, ((bignum**)0), (&remainder), 0,
 314             (BIGNUM_NEGATIVE_P(numerator)));
 315       } else
 316         bignum_divide_unsigned_large_denominator(
 317             numerator, denominator, ((bignum**)0), (&remainder), 0,
 318             (BIGNUM_NEGATIVE_P(numerator)));
 319       return (remainder);
 320     }
 321   }
 322 }
 323
 324 /* cell_to_bignum, fixnum_to_bignum, long_long_to_bignum, ulong_long_to_bignum
 325  */
 326 /* Allocates memory */
 327 #define FOO_TO_BIGNUM(name, type, stype, utype)                       \
 328   bignum* factor_vm::name##_to_bignum(type n) {                       \
 329     int negative_p;                                                   \
 330     /* Special cases win when these small constants are cached. */    \
 331     if (n == 0)                                                       \
 332       return (BIGNUM_ZERO());                                         \
 333     if (n == 1)                                                       \
 334       return (BIGNUM_ONE(0));                                         \
 335     if (n < (type) 0 && n == (type) - 1)                              \
 336       return (BIGNUM_ONE(1));                                         \
 337     {                                                                 \
 338       utype accumulator =                                             \
 339           ((negative_p = (n < (type) 0)) ? ((type)(-(stype) n)) : n); \
 340       if (accumulator < BIGNUM_RADIX)                                 \
 341       {                                                               \
 342         bignum* result = allot_bignum(1, negative_p);                 \
 343         bignum_digit_type* scan = (BIGNUM_START_PTR(result));         \
 344         *scan = (accumulator & BIGNUM_DIGIT_MASK);                    \
 345         return (result);                                              \
 346       } else {                                                        \
 347         BIGNUM_ASSERT(BIGNUM_DIGITS_FOR(type) == 2);                  \
 348         bignum* result = allot_bignum(2, negative_p);                 \
 349         bignum_digit_type* scan = (BIGNUM_START_PTR(result));         \
 350         (*scan++) = (accumulator & BIGNUM_DIGIT_MASK);                \
 351         accumulator >>= BIGNUM_DIGIT_LENGTH;                          \
 352         *scan = (accumulator & BIGNUM_DIGIT_MASK);                    \
 353         BIGNUM_ASSERT((accumulator >> BIGNUM_DIGIT_LENGTH) == 0);     \
 354         return (result);                                              \
 355       }                                                               \
 356     }                                                                 \
 357   }
 358
 359 FOO_TO_BIGNUM(cell, cell, fixnum, cell)
 360 FOO_TO_BIGNUM(fixnum, fixnum, fixnum, cell)
 361 FOO_TO_BIGNUM(long_long, int64_t, int64_t, uint64_t)
 362 FOO_TO_BIGNUM(ulong_long, uint64_t, int64_t, uint64_t)
 363
 364 /* cannot allocate memory */
 365 /* bignum_to_cell, fixnum_to_cell, long_long_to_cell, ulong_long_to_cell */
 366 #define BIGNUM_TO_FOO(name, type, stype, utype)                            \
 367   type factor_vm::bignum_to_##name(bignum * bignum) {                      \
 368     if (BIGNUM_ZERO_P(bignum))                                             \
 369       return (0);                                                          \
 370     {                                                                      \
 371       utype accumulator = 0;                                               \
 372       bignum_digit_type* start = (BIGNUM_START_PTR(bignum));               \
 373       bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bignum)));         \
 374       while (start < scan)                                                 \
 375         accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));  \
 376       return ((BIGNUM_NEGATIVE_P(bignum)) ? ((type)(-(stype) accumulator)) \
 377                                           : accumulator);                  \
 378     }                                                                      \
 379   }
 380
 381 BIGNUM_TO_FOO(cell, cell, fixnum, cell)
 382 BIGNUM_TO_FOO(fixnum, fixnum, fixnum, cell)
 383 BIGNUM_TO_FOO(long_long, int64_t, int64_t, uint64_t)
 384 BIGNUM_TO_FOO(ulong_long, uint64_t, int64_t, uint64_t)
 385
 386 #define DTB_WRITE_DIGIT(factor)                \
 387   {                                            \
 388     significand *= (factor);                   \
 389     digit = ((bignum_digit_type) significand); \
 390     (*--scan) = digit;                         \
 391     significand -= ((double)digit);            \
 392   }
 393
 394 #define inf std::numeric_limits<double>::infinity()
 395
 396 /* Allocates memory */
 397 bignum* factor_vm::double_to_bignum(double x) {
 398   if (x == inf || x == -inf || x != x)
 399     return (BIGNUM_ZERO());
 400   int exponent;
 401   double significand = (frexp(x, (&exponent)));
 402   if (exponent <= 0)
 403     return (BIGNUM_ZERO());
 404   if (exponent == 1)
 405     return (BIGNUM_ONE(x < 0));
 406   if (significand < 0)
 407     significand = (-significand);
 408   {
 409     bignum_length_type length = (BIGNUM_BITS_TO_DIGITS(exponent));
 410     bignum* result = (allot_bignum(length, (x < 0)));
 411     bignum_digit_type* start = (BIGNUM_START_PTR(result));
 412     bignum_digit_type* scan = (start + length);
 413     bignum_digit_type digit;
 414     int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
 415     if (odd_bits > 0)
 416       DTB_WRITE_DIGIT((fixnum)1 << odd_bits);
 417     while (start < scan) {
 418       if (significand == 0) {
 419         while (start < scan)
 420           (*--scan) = 0;
 421         break;
 422       }
 423       DTB_WRITE_DIGIT(BIGNUM_RADIX);
 424     }
 425     return (result);
 426   }
 427 }
 428
 429 #undef DTB_WRITE_DIGIT
 430
 431 /* Comparisons */
 432
 433 int factor_vm::bignum_equal_p_unsigned(bignum* x, bignum* y) {
 434   bignum_length_type length = (BIGNUM_LENGTH(x));
 435   if (length != (BIGNUM_LENGTH(y)))
 436     return (0);
 437   else {
 438     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 439     bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 440     bignum_digit_type* end_x = (scan_x + length);
 441     while (scan_x < end_x)
 442       if ((*scan_x++) != (*scan_y++))
 443         return (0);
 444     return (1);
 445   }
 446 }
 447
 448 enum bignum_comparison factor_vm::bignum_compare_unsigned(bignum* x,
 449                                                           bignum* y) {
 450   bignum_length_type x_length = (BIGNUM_LENGTH(x));
 451   bignum_length_type y_length = (BIGNUM_LENGTH(y));
 452   if (x_length < y_length)
 453     return (bignum_comparison_less);
 454   if (x_length > y_length)
 455     return (bignum_comparison_greater);
 456   {
 457     bignum_digit_type* start_x = (BIGNUM_START_PTR(x));
 458     bignum_digit_type* scan_x = (start_x + x_length);
 459     bignum_digit_type* scan_y = ((BIGNUM_START_PTR(y)) + y_length);
 460     while (start_x < scan_x) {
 461       bignum_digit_type digit_x = (*--scan_x);
 462       bignum_digit_type digit_y = (*--scan_y);
 463       if (digit_x < digit_y)
 464         return (bignum_comparison_less);
 465       if (digit_x > digit_y)
 466         return (bignum_comparison_greater);
 467     }
 468   }
 469   return (bignum_comparison_equal);
 470 }
 471
 472 /* Addition */
 473
 474 /* Allocates memory */
 475 bignum* factor_vm::bignum_add_unsigned(bignum* x, bignum* y, int negative_p) {
 476   GC_BIGNUM(x);
 477   GC_BIGNUM(y);
 478
 479   if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
 480     std::swap(x, y);
 481   }
 482   {
 483     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 484
 485     bignum* r = (allot_bignum((x_length + 1), negative_p));
 486
 487     bignum_digit_type sum;
 488     bignum_digit_type carry = 0;
 489     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 490     bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
 491     {
 492       bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 493       bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
 494       while (scan_y < end_y) {
 495         sum = ((*scan_x++) + (*scan_y++) + carry);
 496         if (sum < BIGNUM_RADIX) {
 497           (*scan_r++) = sum;
 498           carry = 0;
 499         } else {
 500           (*scan_r++) = (sum - BIGNUM_RADIX);
 501           carry = 1;
 502         }
 503       }
 504     }
 505     {
 506       bignum_digit_type* end_x = ((BIGNUM_START_PTR(x)) + x_length);
 507       if (carry != 0)
 508         while (scan_x < end_x) {
 509           sum = ((*scan_x++) + 1);
 510           if (sum < BIGNUM_RADIX) {
 511             (*scan_r++) = sum;
 512             carry = 0;
 513             break;
 514           } else
 515             (*scan_r++) = (sum - BIGNUM_RADIX);
 516         }
 517       while (scan_x < end_x)
 518         (*scan_r++) = (*scan_x++);
 519     }
 520     if (carry != 0) {
 521       (*scan_r) = 1;
 522       return (r);
 523     }
 524     return (bignum_shorten_length(r, x_length));
 525   }
 526 }
 527
 528 /* Subtraction */
 529
 530 /* Allocates memory */
 531 bignum* factor_vm::bignum_subtract_unsigned(bignum* x, bignum* y) {
 532   GC_BIGNUM(x);
 533   GC_BIGNUM(y);
 534
 535   int negative_p = 0;
 536   switch (bignum_compare_unsigned(x, y)) {
 537     case bignum_comparison_equal:
 538       return (BIGNUM_ZERO());
 539     case bignum_comparison_less: {
 540       std::swap(x, y);
 541     }
 542       negative_p = 1;
 543       break;
 544     case bignum_comparison_greater:
 545       negative_p = 0;
 546       break;
 547   }
 548   {
 549     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 550
 551     bignum* r = (allot_bignum(x_length, negative_p));
 552
 553     bignum_digit_type difference;
 554     bignum_digit_type borrow = 0;
 555     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 556     bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
 557     {
 558       bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 559       bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
 560       while (scan_y < end_y) {
 561         difference = (((*scan_x++) - (*scan_y++)) - borrow);
 562         if (difference < 0) {
 563           (*scan_r++) = (difference + BIGNUM_RADIX);
 564           borrow = 1;
 565         } else {
 566           (*scan_r++) = difference;
 567           borrow = 0;
 568         }
 569       }
 570     }
 571     {
 572       bignum_digit_type* end_x = ((BIGNUM_START_PTR(x)) + x_length);
 573       if (borrow != 0)
 574         while (scan_x < end_x) {
 575           difference = ((*scan_x++) - borrow);
 576           if (difference < 0)
 577             (*scan_r++) = (difference + BIGNUM_RADIX);
 578           else {
 579             (*scan_r++) = difference;
 580             borrow = 0;
 581             break;
 582           }
 583         }
 584       BIGNUM_ASSERT(borrow == 0);
 585       while (scan_x < end_x)
 586         (*scan_r++) = (*scan_x++);
 587     }
 588     return (bignum_trim(r));
 589   }
 590 }
 591
 592 /* Multiplication
 593    Maximum value for product_low or product_high:
 594    ((R * R) + (R * (R - 2)) + (R - 1))
 595    Maximum value for carry: ((R * (R - 1)) + (R - 1))
 596    where R == BIGNUM_RADIX_ROOT */
 597
 598 /* Allocates memory */
 599 bignum* factor_vm::bignum_multiply_unsigned(bignum* x, bignum* y,
 600                                             int negative_p) {
 601   GC_BIGNUM(x);
 602   GC_BIGNUM(y);
 603
 604   if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
 605     std::swap(x, y);
 606   }
 607   {
 608     bignum_digit_type carry;
 609     bignum_digit_type y_digit_low;
 610     bignum_digit_type y_digit_high;
 611     bignum_digit_type x_digit_low;
 612     bignum_digit_type x_digit_high;
 613     bignum_digit_type product_low;
 614     bignum_digit_type* scan_r;
 615     bignum_digit_type* scan_y;
 616     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 617     bignum_length_type y_length = (BIGNUM_LENGTH(y));
 618
 619     bignum* r = (allot_bignum_zeroed((x_length + y_length), negative_p));
 620
 621     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 622     bignum_digit_type* end_x = (scan_x + x_length);
 623     bignum_digit_type* start_y = (BIGNUM_START_PTR(y));
 624     bignum_digit_type* end_y = (start_y + y_length);
 625     bignum_digit_type* start_r = (BIGNUM_START_PTR(r));
 626 #define x_digit x_digit_high
 627 #define y_digit y_digit_high
 628 #define product_high carry
 629     while (scan_x < end_x) {
 630       x_digit = (*scan_x++);
 631       x_digit_low = (HD_LOW(x_digit));
 632       x_digit_high = (HD_HIGH(x_digit));
 633       carry = 0;
 634       scan_y = start_y;
 635       scan_r = (start_r++);
 636       while (scan_y < end_y) {
 637         y_digit = (*scan_y++);
 638         y_digit_low = (HD_LOW(y_digit));
 639         y_digit_high = (HD_HIGH(y_digit));
 640         product_low =
 641             ((*scan_r) + (x_digit_low * y_digit_low) + (HD_LOW(carry)));
 642         product_high =
 643             ((x_digit_high * y_digit_low) + (x_digit_low * y_digit_high) +
 644              (HD_HIGH(product_low)) + (HD_HIGH(carry)));
 645         (*scan_r++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
 646         carry = ((x_digit_high * y_digit_high) + (HD_HIGH(product_high)));
 647       }
 648       (*scan_r) += carry;
 649     }
 650     return (bignum_trim(r));
 651 #undef x_digit
 652 #undef y_digit
 653 #undef product_high
 654   }
 655 }
 656
 657 /* Allocates memory */
 658 bignum* factor_vm::bignum_multiply_unsigned_small_factor(bignum* x,
 659                                                          bignum_digit_type y,
 660                                                          int negative_p) {
 661   GC_BIGNUM(x);
 662
 663   bignum_length_type length_x = (BIGNUM_LENGTH(x));
 664
 665   bignum* p = (allot_bignum((length_x + 1), negative_p));
 666
 667   bignum_destructive_copy(x, p);
 668   (BIGNUM_REF(p, length_x)) = 0;
 669   bignum_destructive_scale_up(p, y);
 670   return (bignum_trim(p));
 671 }
 672
 673 void factor_vm::bignum_destructive_add(bignum* bignum, bignum_digit_type n) {
 674   bignum_digit_type* scan = (BIGNUM_START_PTR(bignum));
 675   bignum_digit_type digit;
 676   digit = ((*scan) + n);
 677   if (digit < BIGNUM_RADIX) {
 678     (*scan) = digit;
 679     return;
 680   }
 681   (*scan++) = (digit - BIGNUM_RADIX);
 682   while (1) {
 683     digit = ((*scan) + 1);
 684     if (digit < BIGNUM_RADIX) {
 685       (*scan) = digit;
 686       return;
 687     }
 688     (*scan++) = (digit - BIGNUM_RADIX);
 689   }
 690 }
 691
 692 void factor_vm::bignum_destructive_scale_up(bignum* bignum,
 693                                             bignum_digit_type factor) {
 694   bignum_digit_type carry = 0;
 695   bignum_digit_type* scan = (BIGNUM_START_PTR(bignum));
 696   bignum_digit_type two_digits;
 697   bignum_digit_type product_low;
 698 #define product_high carry
 699   bignum_digit_type* end = (scan + (BIGNUM_LENGTH(bignum)));
 700   BIGNUM_ASSERT((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
 701   while (scan < end) {
 702     two_digits = (*scan);
 703     product_low = ((factor * (HD_LOW(two_digits))) + (HD_LOW(carry)));
 704     product_high = ((factor * (HD_HIGH(two_digits))) + (HD_HIGH(product_low)) +
 705                     (HD_HIGH(carry)));
 706     (*scan++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
 707     carry = (HD_HIGH(product_high));
 708   }
 709   /* A carry here would be an overflow, i.e. it would not fit.
 710            Hopefully the callers allocate enough space that this will
 711            never happen.
 712         */
 713   BIGNUM_ASSERT(carry == 0);
 714   return;
 715 #undef product_high
 716 }
 717
 718 /* Division */
 719
 720 /* For help understanding this algorithm, see:
 721    Knuth, Donald E., "The Art of Computer Programming",
 722    volume 2, "Seminumerical Algorithms"
 723    section 4.3.1, "Multiple-Precision Arithmetic". */
 724
 725 /* Allocates memory */
 726 void factor_vm::bignum_divide_unsigned_large_denominator(
 727     bignum* numerator, bignum* denominator, bignum** quotient,
 728     bignum** remainder, int q_negative_p, int r_negative_p) {
 729   GC_BIGNUM(numerator);
 730   GC_BIGNUM(denominator);
 731
 732   bignum_length_type length_n = ((BIGNUM_LENGTH(numerator)) + 1);
 733   bignum_length_type length_d = (BIGNUM_LENGTH(denominator));
 734
 735   bignum* q = ((quotient != ((bignum**)0))
 736                    ? (allot_bignum((length_n - length_d), q_negative_p))
 737                    : BIGNUM_OUT_OF_BAND);
 738   GC_BIGNUM(q);
 739
 740   bignum* u = (allot_bignum(length_n, r_negative_p));
 741   GC_BIGNUM(u);
 742
 743   int shift = 0;
 744   BIGNUM_ASSERT(length_d > 1);
 745   {
 746     bignum_digit_type v1 = (BIGNUM_REF((denominator), (length_d - 1)));
 747     while (v1 < (BIGNUM_RADIX / 2)) {
 748       v1 <<= 1;
 749       shift += 1;
 750     }
 751   }
 752   if (shift == 0) {
 753     bignum_destructive_copy(numerator, u);
 754     (BIGNUM_REF(u, (length_n - 1))) = 0;
 755     bignum_divide_unsigned_normalized(u, denominator, q);
 756   } else {
 757     bignum* v = (allot_bignum(length_d, 0));
 758
 759     bignum_destructive_normalization(numerator, u, shift);
 760     bignum_destructive_normalization(denominator, v, shift);
 761     bignum_divide_unsigned_normalized(u, v, q);
 762     if (remainder != ((bignum**)0))
 763       bignum_destructive_unnormalization(u, shift);
 764   }
 765
 766   if (q)
 767     q = bignum_trim(q);
 768
 769   u = bignum_trim(u);
 770
 771   if (quotient != ((bignum**)0))
 772     (*quotient) = q;
 773
 774   if (remainder != ((bignum**)0))
 775     (*remainder) = u;
 776
 777   return;
 778 }
 779
 780 void factor_vm::bignum_divide_unsigned_normalized(bignum* u, bignum* v,
 781                                                   bignum* q) {
 782   bignum_length_type u_length = (BIGNUM_LENGTH(u));
 783   bignum_length_type v_length = (BIGNUM_LENGTH(v));
 784   bignum_digit_type* u_start = (BIGNUM_START_PTR(u));
 785   bignum_digit_type* u_scan = (u_start + u_length);
 786   bignum_digit_type* u_scan_limit = (u_start + v_length);
 787   bignum_digit_type* u_scan_start = (u_scan - v_length);
 788   bignum_digit_type* v_start = (BIGNUM_START_PTR(v));
 789   bignum_digit_type* v_end = (v_start + v_length);
 790   bignum_digit_type* q_scan = NULL;
 791   bignum_digit_type v1 = (v_end[-1]);
 792   bignum_digit_type v2 = (v_end[-2]);
 793   bignum_digit_type ph; /* high half of double-digit product */
 794   bignum_digit_type pl; /* low half of double-digit product */
 795   bignum_digit_type guess;
 796   bignum_digit_type gh; /* high half-digit of guess */
 797   bignum_digit_type ch; /* high half of double-digit comparand */
 798   bignum_digit_type v2l = (HD_LOW(v2));
 799   bignum_digit_type v2h = (HD_HIGH(v2));
 800   bignum_digit_type cl; /* low half of double-digit comparand */
 801 #define gl ph           /* low half-digit of guess */
 802 #define uj pl
 803 #define qj ph
 804   bignum_digit_type gm; /* memory loc for reference parameter */
 805   if (q != BIGNUM_OUT_OF_BAND)
 806     q_scan = ((BIGNUM_START_PTR(q)) + (BIGNUM_LENGTH(q)));
 807   while (u_scan_limit < u_scan) {
 808     uj = (*--u_scan);
 809     if (uj != v1) {
 810       /* comparand =
 811                            (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
 812                            guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
 813       cl = (u_scan[-2]);
 814       ch = (bignum_digit_divide(uj, (u_scan[-1]), v1, (&gm)));
 815       guess = gm;
 816     } else {
 817       cl = (u_scan[-2]);
 818       ch = ((u_scan[-1]) + v1);
 819       guess = (BIGNUM_RADIX - 1);
 820     }
 821     while (1) {
 822       /* product = (guess * v2); */
 823       gl = (HD_LOW(guess));
 824       gh = (HD_HIGH(guess));
 825       pl = (v2l * gl);
 826       ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH(pl)));
 827       pl = (HD_CONS((HD_LOW(ph)), (HD_LOW(pl))));
 828       ph = ((v2h * gh) + (HD_HIGH(ph)));
 829       /* if (comparand >= product) */
 830       if ((ch > ph) || ((ch == ph) && (cl >= pl)))
 831         break;
 832       guess -= 1;
 833       /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
 834       ch += v1;
 835       /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
 836       if (ch >= BIGNUM_RADIX)
 837         break;
 838     }
 839     qj = (bignum_divide_subtract(v_start, v_end, guess, (--u_scan_start)));
 840     if (q != BIGNUM_OUT_OF_BAND)
 841       (*--q_scan) = qj;
 842   }
 843   return;
 844 #undef gl
 845 #undef uj
 846 #undef qj
 847 }
 848
 849 bignum_digit_type factor_vm::bignum_divide_subtract(
 850     bignum_digit_type* v_start, bignum_digit_type* v_end,
 851     bignum_digit_type guess, bignum_digit_type* u_start) {
 852   bignum_digit_type* v_scan = v_start;
 853   bignum_digit_type* u_scan = u_start;
 854   bignum_digit_type carry = 0;
 855   if (guess == 0)
 856     return (0);
 857   {
 858     bignum_digit_type gl = (HD_LOW(guess));
 859     bignum_digit_type gh = (HD_HIGH(guess));
 860     bignum_digit_type v;
 861     bignum_digit_type pl;
 862     bignum_digit_type vl;
 863 #define vh v
 864 #define ph carry
 865 #define diff pl
 866     while (v_scan < v_end) {
 867       v = (*v_scan++);
 868       vl = (HD_LOW(v));
 869       vh = (HD_HIGH(v));
 870       pl = ((vl * gl) + (HD_LOW(carry)));
 871       ph = ((vl * gh) + (vh * gl) + (HD_HIGH(pl)) + (HD_HIGH(carry)));
 872       diff = ((*u_scan) - (HD_CONS((HD_LOW(ph)), (HD_LOW(pl)))));
 873       if (diff < 0) {
 874         (*u_scan++) = (diff + BIGNUM_RADIX);
 875         carry = ((vh * gh) + (HD_HIGH(ph)) + 1);
 876       } else {
 877         (*u_scan++) = diff;
 878         carry = ((vh * gh) + (HD_HIGH(ph)));
 879       }
 880     }
 881     if (carry == 0)
 882       return (guess);
 883     diff = ((*u_scan) - carry);
 884     if (diff < 0)
 885       (*u_scan) = (diff + BIGNUM_RADIX);
 886     else {
 887       (*u_scan) = diff;
 888       return (guess);
 889     }
 890 #undef vh
 891 #undef ph
 892 #undef diff
 893   }
 894   /* Subtraction generated carry, implying guess is one too large.
 895            Add v back in to bring it back down. */
 896   v_scan = v_start;
 897   u_scan = u_start;
 898   carry = 0;
 899   while (v_scan < v_end) {
 900     bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
 901     if (sum < BIGNUM_RADIX) {
 902       (*u_scan++) = sum;
 903       carry = 0;
 904     } else {
 905       (*u_scan++) = (sum - BIGNUM_RADIX);
 906       carry = 1;
 907     }
 908   }
 909   if (carry == 1) {
 910     bignum_digit_type sum = ((*u_scan) + carry);
 911     (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
 912   }
 913   return (guess - 1);
 914 }
 915
 916 /* Allocates memory */
 917 void factor_vm::bignum_divide_unsigned_medium_denominator(
 918     bignum* numerator, bignum_digit_type denominator, bignum** quotient,
 919     bignum** remainder, int q_negative_p, int r_negative_p) {
 920   GC_BIGNUM(numerator);
 921
 922   bignum_length_type length_n = (BIGNUM_LENGTH(numerator));
 923   bignum_length_type length_q;
 924   bignum* q = NULL;
 925   GC_BIGNUM(q);
 926
 927   int shift = 0;
 928   /* Because `bignum_digit_divide' requires a normalized denominator. */
 929   while (denominator < (BIGNUM_RADIX / 2)) {
 930     denominator <<= 1;
 931     shift += 1;
 932   }
 933   if (shift == 0) {
 934     length_q = length_n;
 935
 936     q = (allot_bignum(length_q, q_negative_p));
 937     bignum_destructive_copy(numerator, q);
 938   } else {
 939     length_q = (length_n + 1);
 940
 941     q = (allot_bignum(length_q, q_negative_p));
 942     bignum_destructive_normalization(numerator, q, shift);
 943   }
 944   {
 945     bignum_digit_type r = 0;
 946     bignum_digit_type* start = (BIGNUM_START_PTR(q));
 947     bignum_digit_type* scan = (start + length_q);
 948     bignum_digit_type qj;
 949
 950     while (start < scan) {
 951       r = (bignum_digit_divide(r, (*--scan), denominator, (&qj)));
 952       (*scan) = qj;
 953     }
 954
 955     q = bignum_trim(q);
 956
 957     if (remainder != ((bignum**)0)) {
 958       if (shift != 0)
 959         r >>= shift;
 960
 961       (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
 962     }
 963
 964     if (quotient != ((bignum**)0))
 965       (*quotient) = q;
 966   }
 967   return;
 968 }
 969
 970 void factor_vm::bignum_destructive_normalization(bignum* source, bignum* target,
 971                                                  int shift_left) {
 972   bignum_digit_type digit;
 973   bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
 974   bignum_digit_type carry = 0;
 975   bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
 976   bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
 977   bignum_digit_type* end_target = (scan_target + (BIGNUM_LENGTH(target)));
 978   int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
 979   bignum_digit_type mask = (((cell)1 << shift_right) - 1);
 980   while (scan_source < end_source) {
 981     digit = (*scan_source++);
 982     (*scan_target++) = (((digit & mask) << shift_left) | carry);
 983     carry = (digit >> shift_right);
 984   }
 985   if (scan_target < end_target)
 986     (*scan_target) = carry;
 987   else
 988     BIGNUM_ASSERT(carry == 0);
 989   return;
 990 }
 991
 992 void factor_vm::bignum_destructive_unnormalization(bignum* bignum,
 993                                                    int shift_right) {
 994   bignum_digit_type* start = (BIGNUM_START_PTR(bignum));
 995   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bignum)));
 996   bignum_digit_type digit;
 997   bignum_digit_type carry = 0;
 998   int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
 999   bignum_digit_type mask = (((fixnum)1 << shift_right) - 1);
1000   while (start < scan) {
1001     digit = (*--scan);
1002     (*scan) = ((digit >> shift_right) | carry);
1003     carry = ((digit & mask) << shift_left);
1004   }
1005   BIGNUM_ASSERT(carry == 0);
1006   return;
1007 }
1008
1009 /* This is a reduced version of the division algorithm, applied to the
1010    case of dividing two bignum digits by one bignum digit.  It is
1011    assumed that the numerator, denominator are normalized. */
1012
1013 #define BDD_STEP(qn, j)                                          \
1014   {                                                              \
1015     uj = (u[j]);                                                 \
1016     if (uj != v1) {                                              \
1017       uj_uj1 = (HD_CONS(uj, (u[j + 1])));                        \
1018       guess = (uj_uj1 / v1);                                     \
1019       comparand = (HD_CONS((uj_uj1 % v1), (u[j + 2])));          \
1020     } else {                                                     \
1021       guess = (BIGNUM_RADIX_ROOT - 1);                           \
1022       comparand = (HD_CONS(((u[j + 1]) + v1), (u[j + 2])));      \
1023     }                                                            \
1024     while ((guess * v2) > comparand) {                           \
1025       guess -= 1;                                                \
1026       comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH);             \
1027       if (comparand >= BIGNUM_RADIX)                             \
1028         break;                                                   \
1029     }                                                            \
1030     qn = (bignum_digit_divide_subtract(v1, v2, guess, (&u[j]))); \
1031   }
1032
1033 bignum_digit_type factor_vm::bignum_digit_divide(
1034     bignum_digit_type uh, bignum_digit_type ul, bignum_digit_type v,
1035     bignum_digit_type* q) /* return value */
1036     {
1037   bignum_digit_type guess;
1038   bignum_digit_type comparand;
1039   bignum_digit_type v1 = (HD_HIGH(v));
1040   bignum_digit_type v2 = (HD_LOW(v));
1041   bignum_digit_type uj;
1042   bignum_digit_type uj_uj1;
1043   bignum_digit_type q1;
1044   bignum_digit_type q2;
1045   bignum_digit_type u[4];
1046   if (uh == 0) {
1047     if (ul < v) {
1048       (*q) = 0;
1049       return (ul);
1050     } else if (ul == v) {
1051       (*q) = 1;
1052       return (0);
1053     }
1054   }
1055   (u[0]) = (HD_HIGH(uh));
1056   (u[1]) = (HD_LOW(uh));
1057   (u[2]) = (HD_HIGH(ul));
1058   (u[3]) = (HD_LOW(ul));
1059   v1 = (HD_HIGH(v));
1060   v2 = (HD_LOW(v));
1061   BDD_STEP(q1, 0);
1062   BDD_STEP(q2, 1);
1063   (*q) = (HD_CONS(q1, q2));
1064   return (HD_CONS((u[2]), (u[3])));
1065 }
1066
1067 #undef BDD_STEP
1068
1069 #define BDDS_MULSUB(vn, un, carry_in)    \
1070   {                                      \
1071     product = ((vn * guess) + carry_in); \
1072     diff = (un - (HD_LOW(product)));     \
1073     if (diff < 0) {                      \
1074       un = (diff + BIGNUM_RADIX_ROOT);   \
1075       carry = ((HD_HIGH(product)) + 1);  \
1076     } else {                             \
1077       un = diff;                         \
1078       carry = (HD_HIGH(product));        \
1079     }                                    \
1080   }
1081
1082 #define BDDS_ADD(vn, un, carry_in)    \
1083   {                                   \
1084     sum = (vn + un + carry_in);       \
1085     if (sum < BIGNUM_RADIX_ROOT) {    \
1086       un = sum;                       \
1087       carry = 0;                      \
1088     } else {                          \
1089       un = (sum - BIGNUM_RADIX_ROOT); \
1090       carry = 1;                      \
1091     }                                 \
1092   }
1093
1094 bignum_digit_type factor_vm::bignum_digit_divide_subtract(
1095     bignum_digit_type v1, bignum_digit_type v2, bignum_digit_type guess,
1096     bignum_digit_type* u) {
1097   {
1098     bignum_digit_type product;
1099     bignum_digit_type diff;
1100     bignum_digit_type carry;
1101     BDDS_MULSUB(v2, (u[2]), 0);
1102     BDDS_MULSUB(v1, (u[1]), carry);
1103     if (carry == 0)
1104       return (guess);
1105     diff = ((u[0]) - carry);
1106     if (diff < 0)
1107       (u[0]) = (diff + BIGNUM_RADIX);
1108     else {
1109       (u[0]) = diff;
1110       return (guess);
1111     }
1112   }
1113   {
1114     bignum_digit_type sum;
1115     bignum_digit_type carry;
1116     BDDS_ADD(v2, (u[2]), 0);
1117     BDDS_ADD(v1, (u[1]), carry);
1118     if (carry == 1)
1119       (u[0]) += 1;
1120   }
1121   return (guess - 1);
1122 }
1123
1124 #undef BDDS_MULSUB
1125 #undef BDDS_ADD
1126
1127 /* Allocates memory */
1128 void factor_vm::bignum_divide_unsigned_small_denominator(
1129     bignum* numerator, bignum_digit_type denominator, bignum** quotient,
1130     bignum** remainder, int q_negative_p, int r_negative_p) {
1131   GC_BIGNUM(numerator);
1132
1133   bignum* q = (bignum_new_sign(numerator, q_negative_p));
1134   GC_BIGNUM(q);
1135
1136   bignum_digit_type r = (bignum_destructive_scale_down(q, denominator));
1137
1138   q = (bignum_trim(q));
1139
1140   if (remainder != ((bignum**)0))
1141     (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
1142
1143   (*quotient) = q;
1144
1145   return;
1146 }
1147
1148 /* Given (denominator > 1), it is fairly easy to show that
1149    (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
1150    that all digits are < BIGNUM_RADIX. */
1151
1152 bignum_digit_type factor_vm::bignum_destructive_scale_down(
1153     bignum* bignum, bignum_digit_type denominator) {
1154   bignum_digit_type numerator;
1155   bignum_digit_type remainder = 0;
1156   bignum_digit_type two_digits;
1157 #define quotient_high remainder
1158   bignum_digit_type* start = (BIGNUM_START_PTR(bignum));
1159   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bignum)));
1160   BIGNUM_ASSERT((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
1161   while (start < scan) {
1162     two_digits = (*--scan);
1163     numerator = (HD_CONS(remainder, (HD_HIGH(two_digits))));
1164     quotient_high = (numerator / denominator);
1165     numerator = (HD_CONS((numerator % denominator), (HD_LOW(two_digits))));
1166     (*scan) = (HD_CONS(quotient_high, (numerator / denominator)));
1167     remainder = (numerator % denominator);
1168   }
1169   return (remainder);
1170 #undef quotient_high
1171 }
1172
1173 /* Allocates memory */
1174 bignum* factor_vm::bignum_remainder_unsigned_small_denominator(
1175     bignum* n, bignum_digit_type d, int negative_p) {
1176   bignum_digit_type two_digits;
1177   bignum_digit_type* start = (BIGNUM_START_PTR(n));
1178   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(n)));
1179   bignum_digit_type r = 0;
1180   BIGNUM_ASSERT((d > 1) && (d < BIGNUM_RADIX_ROOT));
1181   while (start < scan) {
1182     two_digits = (*--scan);
1183     r = ((HD_CONS(((HD_CONS(r, (HD_HIGH(two_digits)))) % d),
1184                   (HD_LOW(two_digits)))) %
1185          d);
1186   }
1187   return (bignum_digit_to_bignum(r, negative_p));
1188 }
1189
1190 /* Allocates memory */
1191 bignum* factor_vm::bignum_digit_to_bignum(bignum_digit_type digit,
1192                                           int negative_p) {
1193   if (digit == 0)
1194     return (BIGNUM_ZERO());
1195   else {
1196     bignum* result = (allot_bignum(1, negative_p));
1197     (BIGNUM_REF(result, 0)) = digit;
1198     return (result);
1199   }
1200 }
1201
1202 /* Allocates memory */
1203 bignum* factor_vm::allot_bignum(bignum_length_type length, int negative_p) {
1204   BIGNUM_ASSERT((length >= 0) || (length < BIGNUM_RADIX));
1205   bignum* result = allot_uninitialized_array<bignum>(length + 1);
1206   BIGNUM_SET_NEGATIVE_P(result, negative_p);
1207   return (result);
1208 }
1209
1210 /* Allocates memory */
1211 bignum* factor_vm::allot_bignum_zeroed(bignum_length_type length,
1212                                        int negative_p) {
1213   bignum* result = allot_bignum(length, negative_p);
1214   bignum_digit_type* scan = (BIGNUM_START_PTR(result));
1215   bignum_digit_type* end = (scan + length);
1216   while (scan < end)
1217     (*scan++) = 0;
1218   return (result);
1219 }
1220
1221 /* can allocate if not in nursery or size is larger */
1222 /* Allocates memory conditionally */
1223 #define BIGNUM_REDUCE_LENGTH(source, length) \
1224   source = reallot_array(source, length + 1)
1225
1226 /* Allocates memory */
1227 bignum* factor_vm::bignum_shorten_length(bignum* bignum,
1228                                          bignum_length_type length) {
1229   bignum_length_type current_length = (BIGNUM_LENGTH(bignum));
1230   BIGNUM_ASSERT((length >= 0) || (length <= current_length));
1231   if (length < current_length) {
1232     GC_BIGNUM(bignum);
1233     BIGNUM_REDUCE_LENGTH(bignum, length);
1234     BIGNUM_SET_NEGATIVE_P(bignum, (length != 0) && (BIGNUM_NEGATIVE_P(bignum)));
1235   }
1236   return (bignum);
1237 }
1238
1239 /* Allocates memory */
1240 bignum* factor_vm::bignum_trim(bignum* bignum) {
1241   bignum_digit_type* start = (BIGNUM_START_PTR(bignum));
1242   bignum_digit_type* end = (start + (BIGNUM_LENGTH(bignum)));
1243   bignum_digit_type* scan = end;
1244   while ((start <= scan) && ((*--scan) == 0))
1245     ;
1246   scan += 1;
1247   if (scan < end) {
1248     GC_BIGNUM(bignum);
1249     bignum_length_type length = (scan - start);
1250     BIGNUM_REDUCE_LENGTH(bignum, length);
1251     BIGNUM_SET_NEGATIVE_P(bignum, (length != 0) && (BIGNUM_NEGATIVE_P(bignum)));
1252   }
1253   return (bignum);
1254 }
1255
1256 /* Copying */
1257
1258 /* Allocates memory */
1259 bignum* factor_vm::bignum_new_sign(bignum* x, int negative_p) {
1260   GC_BIGNUM(x);
1261   bignum* result = (allot_bignum((BIGNUM_LENGTH(x)), negative_p));
1262
1263   bignum_destructive_copy(x, result);
1264   return (result);
1265 }
1266
1267 /* Allocates memory */
1268 bignum* factor_vm::bignum_maybe_new_sign(bignum* x, int negative_p) {
1269   if ((BIGNUM_NEGATIVE_P(x)) ? negative_p : (!negative_p))
1270     return (x);
1271   else {
1272     GC_BIGNUM(x);
1273     bignum* result = (allot_bignum((BIGNUM_LENGTH(x)), negative_p));
1274     bignum_destructive_copy(x, result);
1275     return (result);
1276   }
1277 }
1278
1279 void factor_vm::bignum_destructive_copy(bignum* source, bignum* target) {
1280   bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
1281   bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1282   bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1283   while (scan_source < end_source)
1284     (*scan_target++) = (*scan_source++);
1285   return;
1286 }
1287
1288 /*
1289  * Added bitwise operations (and oddp).
1290  */
1291
1292 /* Allocates memory */
1293 bignum* factor_vm::bignum_bitwise_not(bignum* x) {
1294   GC_BIGNUM(x);
1295
1296   bignum_length_type size = BIGNUM_LENGTH(x);
1297   bignum_digit_type* scan_x, *end_x, *scan_y;
1298   bignum* y;
1299   int carry = 1;
1300
1301   if (BIGNUM_NEGATIVE_P(x)) {
1302     y = allot_bignum(size, 0);
1303     scan_x = BIGNUM_START_PTR(x);
1304     end_x = scan_x + size;
1305     scan_y = BIGNUM_START_PTR(y);
1306     while (scan_x < end_x) {
1307       if (*scan_x == 0) {
1308         *scan_y++ = BIGNUM_RADIX - 1;
1309         scan_x++;
1310       } else {
1311         *scan_y++ = *scan_x++ - 1;
1312         carry = 0;
1313         break;
1314       }
1315     }
1316   } else {
1317     y = allot_bignum(size, 1);
1318     scan_x = BIGNUM_START_PTR(x);
1319     end_x = scan_x + size;
1320     scan_y = BIGNUM_START_PTR(y);
1321     while (scan_x < end_x) {
1322       if (*scan_x == (BIGNUM_RADIX - 1)) {
1323         *scan_y++ = 0;
1324         scan_x++;
1325       } else {
1326         *scan_y++ = *scan_x++ + 1;
1327         carry = 0;
1328         break;
1329       }
1330     }
1331   }
1332
1333   while (scan_x < end_x) {
1334     *scan_y++ = *scan_x++;
1335   }
1336
1337   if (carry) {
1338     GC_BIGNUM(y);
1339     x = allot_bignum(size + 1, BIGNUM_NEGATIVE_P(y));
1340     bignum_destructive_copy(y, x);
1341     scan_x = BIGNUM_START_PTR(x);
1342     *(scan_x + size) = 1;
1343     return x;
1344   } else {
1345     return bignum_trim(y);
1346   }
1347 }
1348
1349 /* Allocates memory */
1350 bignum* factor_vm::bignum_arithmetic_shift(bignum* arg1, fixnum n) {
1351   if (BIGNUM_NEGATIVE_P(arg1) && n < 0)
1352     return bignum_bitwise_not(
1353         bignum_magnitude_ash(bignum_bitwise_not(arg1), n));
1354   else
1355     return bignum_magnitude_ash(arg1, n);
1356 }
1357
1358 #define AND_OP 0
1359 #define IOR_OP 1
1360 #define XOR_OP 2
1361
1362 /* Allocates memory */
1363 bignum* factor_vm::bignum_bitwise_and(bignum* arg1, bignum* arg2) {
1364   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1365               ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2)
1366               : bignum_posneg_bitwise_op(AND_OP, arg2, arg1)
1367               : (BIGNUM_NEGATIVE_P(arg2))
1368               ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2)
1369               : bignum_pospos_bitwise_op(AND_OP, arg1, arg2));
1370 }
1371
1372 /* Allocates memory */
1373 bignum* factor_vm::bignum_bitwise_ior(bignum* arg1, bignum* arg2) {
1374   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1375               ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2)
1376               : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1)
1377               : (BIGNUM_NEGATIVE_P(arg2))
1378               ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2)
1379               : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2));
1380 }
1381
1382 /* Allocates memory */
1383 bignum* factor_vm::bignum_bitwise_xor(bignum* arg1, bignum* arg2) {
1384   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1385               ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2)
1386               : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1)
1387               : (BIGNUM_NEGATIVE_P(arg2))
1388               ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2)
1389               : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2));
1390 }
1391
1392 /* Allocates memory */
1393 /* ash for the magnitude */
1394 /* assume arg1 is a big number, n is a long */
1395 bignum* factor_vm::bignum_magnitude_ash(bignum* arg1, fixnum n) {
1396   GC_BIGNUM(arg1);
1397
1398   bignum* result = NULL;
1399   bignum_digit_type* scan1;
1400   bignum_digit_type* scanr;
1401   bignum_digit_type* end;
1402
1403   fixnum digit_offset, bit_offset;
1404
1405   if (BIGNUM_ZERO_P(arg1))
1406     return (arg1);
1407
1408   if (n > 0) {
1409     digit_offset = n / BIGNUM_DIGIT_LENGTH;
1410     bit_offset = n % BIGNUM_DIGIT_LENGTH;
1411
1412     result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) + digit_offset + 1,
1413                                  BIGNUM_NEGATIVE_P(arg1));
1414
1415     scanr = BIGNUM_START_PTR(result) + digit_offset;
1416     scan1 = BIGNUM_START_PTR(arg1);
1417     end = scan1 + BIGNUM_LENGTH(arg1);
1418
1419     while (scan1 < end) {
1420       *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
1421       *scanr = *scanr & BIGNUM_DIGIT_MASK;
1422       scanr++;
1423       *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
1424       *scanr = *scanr & BIGNUM_DIGIT_MASK;
1425     }
1426   } else if (n < 0 && (-n >= (BIGNUM_LENGTH(arg1) * (bignum_length_type)
1427                               BIGNUM_DIGIT_LENGTH)))
1428     result = BIGNUM_ZERO();
1429
1430   else if (n < 0) {
1431     digit_offset = -n / BIGNUM_DIGIT_LENGTH;
1432     bit_offset = -n % BIGNUM_DIGIT_LENGTH;
1433
1434     result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) - digit_offset,
1435                                  BIGNUM_NEGATIVE_P(arg1));
1436
1437     scanr = BIGNUM_START_PTR(result);
1438     scan1 = BIGNUM_START_PTR(arg1) + digit_offset;
1439     end = scanr + BIGNUM_LENGTH(result) - 1;
1440
1441     while (scanr < end) {
1442       *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1443       *scanr = (*scanr | *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) &
1444                BIGNUM_DIGIT_MASK;
1445       scanr++;
1446     }
1447     *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1448   } else if (n == 0)
1449     result = arg1;
1450
1451   return (bignum_trim(result));
1452 }
1453
1454 /* Allocates memory */
1455 bignum* factor_vm::bignum_pospos_bitwise_op(int op, bignum* arg1,
1456                                             bignum* arg2) {
1457   GC_BIGNUM(arg1);
1458   GC_BIGNUM(arg2);
1459
1460   bignum* result;
1461   bignum_length_type max_length;
1462
1463   bignum_digit_type* scan1, *end1, digit1;
1464   bignum_digit_type* scan2, *end2, digit2;
1465   bignum_digit_type* scanr, *endr;
1466
1467   max_length =
1468       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1)
1469                                                   : BIGNUM_LENGTH(arg2);
1470
1471   result = allot_bignum(max_length, 0);
1472
1473   scanr = BIGNUM_START_PTR(result);
1474   scan1 = BIGNUM_START_PTR(arg1);
1475   scan2 = BIGNUM_START_PTR(arg2);
1476   endr = scanr + max_length;
1477   end1 = scan1 + BIGNUM_LENGTH(arg1);
1478   end2 = scan2 + BIGNUM_LENGTH(arg2);
1479
1480   while (scanr < endr) {
1481     digit1 = (scan1 < end1) ? *scan1++ : 0;
1482     digit2 = (scan2 < end2) ? *scan2++ : 0;
1483     *scanr++ =
1484         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1485                                                           : digit1 ^ digit2;
1486   }
1487   return bignum_trim(result);
1488 }
1489
1490 /* Allocates memory */
1491 bignum* factor_vm::bignum_posneg_bitwise_op(int op, bignum* arg1,
1492                                             bignum* arg2) {
1493   GC_BIGNUM(arg1);
1494   GC_BIGNUM(arg2);
1495
1496   bignum* result;
1497   bignum_length_type max_length;
1498
1499   bignum_digit_type* scan1, *end1, digit1;
1500   bignum_digit_type* scan2, *end2, digit2, carry2;
1501   bignum_digit_type* scanr, *endr;
1502
1503   char neg_p = op == IOR_OP || op == XOR_OP;
1504
1505   max_length =
1506       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1) ? BIGNUM_LENGTH(arg1)
1507                                                       : BIGNUM_LENGTH(arg2) + 1;
1508
1509   result = allot_bignum(max_length, neg_p);
1510
1511   scanr = BIGNUM_START_PTR(result);
1512   scan1 = BIGNUM_START_PTR(arg1);
1513   scan2 = BIGNUM_START_PTR(arg2);
1514   endr = scanr + max_length;
1515   end1 = scan1 + BIGNUM_LENGTH(arg1);
1516   end2 = scan2 + BIGNUM_LENGTH(arg2);
1517
1518   carry2 = 1;
1519
1520   while (scanr < endr) {
1521     digit1 = (scan1 < end1) ? *scan1++ : 0;
1522     digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1523
1524     if (digit2 < BIGNUM_RADIX)
1525       carry2 = 0;
1526     else {
1527       digit2 = (digit2 - BIGNUM_RADIX);
1528       carry2 = 1;
1529     }
1530
1531     *scanr++ =
1532         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1533                                                           : digit1 ^ digit2;
1534   }
1535
1536   if (neg_p)
1537     bignum_negate_magnitude(result);
1538
1539   return bignum_trim(result);
1540 }
1541
1542 /* Allocates memory */
1543 bignum* factor_vm::bignum_negneg_bitwise_op(int op, bignum* arg1,
1544                                             bignum* arg2) {
1545   GC_BIGNUM(arg1);
1546   GC_BIGNUM(arg2);
1547
1548   bignum* result;
1549   bignum_length_type max_length;
1550
1551   bignum_digit_type* scan1, *end1, digit1, carry1;
1552   bignum_digit_type* scan2, *end2, digit2, carry2;
1553   bignum_digit_type* scanr, *endr;
1554
1555   char neg_p = op == AND_OP || op == IOR_OP;
1556
1557   max_length =
1558       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) + 1
1559                                                   : BIGNUM_LENGTH(arg2) + 1;
1560
1561   result = allot_bignum(max_length, neg_p);
1562
1563   scanr = BIGNUM_START_PTR(result);
1564   scan1 = BIGNUM_START_PTR(arg1);
1565   scan2 = BIGNUM_START_PTR(arg2);
1566   endr = scanr + max_length;
1567   end1 = scan1 + BIGNUM_LENGTH(arg1);
1568   end2 = scan2 + BIGNUM_LENGTH(arg2);
1569
1570   carry1 = 1;
1571   carry2 = 1;
1572
1573   while (scanr < endr) {
1574     digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
1575     digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1576
1577     if (digit1 < BIGNUM_RADIX)
1578       carry1 = 0;
1579     else {
1580       digit1 = (digit1 - BIGNUM_RADIX);
1581       carry1 = 1;
1582     }
1583
1584     if (digit2 < BIGNUM_RADIX)
1585       carry2 = 0;
1586     else {
1587       digit2 = (digit2 - BIGNUM_RADIX);
1588       carry2 = 1;
1589     }
1590
1591     *scanr++ =
1592         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1593                                                           : digit1 ^ digit2;
1594   }
1595
1596   if (neg_p)
1597     bignum_negate_magnitude(result);
1598
1599   return bignum_trim(result);
1600 }
1601
1602 void factor_vm::bignum_negate_magnitude(bignum* arg) {
1603   bignum_digit_type* scan;
1604   bignum_digit_type* end;
1605   bignum_digit_type digit;
1606   bignum_digit_type carry;
1607
1608   scan = BIGNUM_START_PTR(arg);
1609   end = scan + BIGNUM_LENGTH(arg);
1610
1611   carry = 1;
1612
1613   while (scan < end) {
1614     digit = (~ * scan & BIGNUM_DIGIT_MASK) + carry;
1615
1616     if (digit < BIGNUM_RADIX)
1617       carry = 0;
1618     else {
1619       digit = (digit - BIGNUM_RADIX);
1620       carry = 1;
1621     }
1622
1623     *scan++ = digit;
1624   }
1625 }
1626
1627 /* Allocates memory */
1628 bignum* factor_vm::bignum_integer_length(bignum* x) {
1629   GC_BIGNUM(x);
1630
1631   bignum_length_type index = ((BIGNUM_LENGTH(x)) - 1);
1632   bignum_digit_type digit = (BIGNUM_REF(x, index));
1633   bignum_digit_type carry = 0;
1634   bignum* result;
1635
1636   while (digit > 1) {
1637     carry += 1;
1638     digit >>= 1;
1639   }
1640
1641   if (index < BIGNUM_RADIX_ROOT) {
1642      result = allot_bignum(1, 0);
1643      (BIGNUM_REF(result, 0)) = (index * BIGNUM_DIGIT_LENGTH) + carry;
1644   } else {
1645      result = allot_bignum(2, 0);
1646      (BIGNUM_REF(result, 0)) = index;
1647      (BIGNUM_REF(result, 1)) = 0;
1648      bignum_destructive_scale_up(result, BIGNUM_DIGIT_LENGTH);
1649      bignum_destructive_add(result, carry);
1650   }
1651   return (bignum_trim(result));
1652 }
1653
1654 /* Allocates memory */
1655 int factor_vm::bignum_logbitp(int shift, bignum* arg) {
1656   return ((BIGNUM_NEGATIVE_P(arg))
1657               ? !bignum_unsigned_logbitp(shift, bignum_bitwise_not(arg))
1658               : bignum_unsigned_logbitp(shift, arg));
1659 }
1660
1661 int factor_vm::bignum_unsigned_logbitp(int shift, bignum* bignum) {
1662   bignum_length_type len = (BIGNUM_LENGTH(bignum));
1663   int index = shift / BIGNUM_DIGIT_LENGTH;
1664   if (index >= len)
1665     return 0;
1666   bignum_digit_type digit = (BIGNUM_REF(bignum, index));
1667   int p = shift % BIGNUM_DIGIT_LENGTH;
1668   bignum_digit_type mask = ((fixnum)1) << p;
1669   return (digit & mask) ? 1 : 0;
1670 }
1671
1672 #ifdef _WIN64
1673 /* Allocates memory */
1674 bignum* factor_vm::bignum_gcd(bignum* a, bignum* b) {
1675   GC_BIGNUM(a);
1676   GC_BIGNUM(b);
1677   bignum* d;
1678   bignum_length_type size_a, size_b;
1679   bignum_digit_type* scan_a, *scan_b, *scan_d, *a_end, *b_end;
1680
1681   if (BIGNUM_NEGATIVE_P(a)) {
1682     size_a = BIGNUM_LENGTH(a);
1683     d = allot_bignum(size_a, 0);
1684     scan_d = BIGNUM_START_PTR(d);
1685     scan_a = BIGNUM_START_PTR(a);
1686     a_end = scan_a + size_a;
1687     while (scan_a < a_end)
1688       (*scan_d++) = (*scan_a++);
1689     a = d;
1690   }
1691
1692   if (BIGNUM_NEGATIVE_P(b)) {
1693     size_b = BIGNUM_LENGTH(b);
1694     d = allot_bignum(size_b, 0);
1695     scan_d = BIGNUM_START_PTR(d);
1696     scan_b = BIGNUM_START_PTR(b);
1697     b_end = scan_b + size_b;
1698     while (scan_b < b_end)
1699       (*scan_d++) = (*scan_b++);
1700     b = d;
1701   }
1702
1703   if (bignum_compare(a, b) == bignum_comparison_less) {
1704     std::swap(a, b);
1705   }
1706
1707   while (BIGNUM_LENGTH(b) != 0) {
1708     d = bignum_remainder(a, b);
1709     GC_BIGNUM(d);
1710     if (d == BIGNUM_OUT_OF_BAND) {
1711       return d;
1712     }
1713     a = b;
1714     b = d;
1715   }
1716
1717   return a;
1718 }
1719 #else
1720 /* Allocates memory */
1721 bignum* factor_vm::bignum_gcd(bignum* a, bignum* b) {
1722   GC_BIGNUM(a);
1723   GC_BIGNUM(b);
1724   bignum* c, *d, *e, *f;
1725   bignum_twodigit_type x, y, q, s, t, A, B, C, D;
1726   int nbits, k;
1727   bignum_length_type size_a, size_b, size_c;
1728   bignum_digit_type* scan_a, *scan_b, *scan_c, *scan_d;
1729   bignum_digit_type* a_end, *b_end, *c_end;
1730
1731   /* clone the bignums so we can modify them in-place */
1732   size_a = BIGNUM_LENGTH(a);
1733   c = allot_bignum(size_a, 0);
1734   scan_a = BIGNUM_START_PTR(a);
1735   a_end = scan_a + size_a;
1736   GC_BIGNUM(c);
1737   scan_c = BIGNUM_START_PTR(c);
1738   while (scan_a < a_end)
1739     (*scan_c++) = (*scan_a++);
1740   a = c;
1741   size_b = BIGNUM_LENGTH(b);
1742   d = allot_bignum(size_b, 0);
1743   scan_b = BIGNUM_START_PTR(b);
1744   b_end = scan_b + size_b;
1745   GC_BIGNUM(d);
1746   scan_d = BIGNUM_START_PTR(d);
1747   while (scan_b < b_end)
1748     (*scan_d++) = (*scan_b++);
1749   b = d;
1750
1751   /* Initial reduction: make sure that 0 <= b <= a. */
1752   if (bignum_compare(a, b) == bignum_comparison_less) {
1753     std::swap(a, b);
1754     std::swap(size_a, size_b);
1755   }
1756
1757   while (size_a > 1) {
1758     nbits = log2(BIGNUM_REF(a, size_a - 1));
1759     x = ((BIGNUM_REF(a, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)) |
1760          (BIGNUM_REF(a, size_a - 2) >> nbits));
1761     y = ((size_b >= size_a - 1 ? BIGNUM_REF(b, size_a - 2) >> nbits : 0) |
1762          (size_b >= size_a
1763               ? BIGNUM_REF(b, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)
1764               : 0));
1765
1766     /* inner loop of Lehmer's algorithm; */
1767     A = 1;
1768     B = 0;
1769     C = 0;
1770     D = 1;
1771     for (k = 0;; k++) {
1772       if (y - C == 0)
1773         break;
1774
1775       q = (x + (A - 1)) / (y - C);
1776
1777       s = B + (q * D);
1778       t = x - (q * y);
1779
1780       if (s > t)
1781         break;
1782
1783       x = y;
1784       y = t;
1785
1786       t = A + (q * C);
1787
1788       A = D;
1789       B = C;
1790       C = s;
1791       D = t;
1792     }
1793
1794     if (k == 0) {
1795       /* no progress; do a Euclidean step */
1796       if (size_b == 0) {
1797         return bignum_trim(a);
1798       }
1799       e = bignum_trim(a);
1800       GC_BIGNUM(e);
1801       f = bignum_trim(b);
1802       GC_BIGNUM(f);
1803       c = bignum_remainder(e, f);
1804       GC_BIGNUM(c);
1805       if (c == BIGNUM_OUT_OF_BAND) {
1806         return c;
1807       }
1808
1809       // copy 'b' to 'a'
1810       scan_a = BIGNUM_START_PTR(a);
1811       scan_b = BIGNUM_START_PTR(b);
1812       a_end = scan_a + size_a;
1813       b_end = scan_b + size_b;
1814       while (scan_b < b_end)
1815         *(scan_a++) = *(scan_b++);
1816       while (scan_a < a_end)
1817         *(scan_a++) = 0;
1818       size_a = size_b;
1819
1820       // copy 'c' to 'b'
1821       scan_b = BIGNUM_START_PTR(b);
1822       scan_c = BIGNUM_START_PTR(c);
1823       size_c = BIGNUM_LENGTH(c);
1824       c_end = scan_c + size_c;
1825       while (scan_c < c_end)
1826         *(scan_b++) = *(scan_c++);
1827       while (scan_b < b_end)
1828         *(scan_b++) = 0;
1829       size_b = size_c;
1830
1831       continue;
1832     }
1833
1834     /*
1835       a, b = A*b - B*a, D*a - C*b if k is odd
1836       a, b = A*a - B*b, D*b - C*a if k is even
1837     */
1838     scan_a = BIGNUM_START_PTR(a);
1839     scan_b = BIGNUM_START_PTR(b);
1840     scan_c = scan_a;
1841     scan_d = scan_b;
1842     a_end = scan_a + size_a;
1843     b_end = scan_b + size_b;
1844     s = 0;
1845     t = 0;
1846     if (k & 1) {
1847       while (scan_b < b_end) {
1848         s += (A * *scan_b) - (B * *scan_a);
1849         t += (D * *scan_a++) - (C * *scan_b++);
1850         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1851         *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1852         s >>= BIGNUM_DIGIT_LENGTH;
1853         t >>= BIGNUM_DIGIT_LENGTH;
1854       }
1855       while (scan_a < a_end) {
1856         s -= (B * *scan_a);
1857         t += (D * *scan_a++);
1858         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1859         //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1860         s >>= BIGNUM_DIGIT_LENGTH;
1861         t >>= BIGNUM_DIGIT_LENGTH;
1862       }
1863     } else {
1864       while (scan_b < b_end) {
1865         s += (A * *scan_a) - (B * *scan_b);
1866         t += (D * *scan_b++) - (C * *scan_a++);
1867         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1868         *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1869         s >>= BIGNUM_DIGIT_LENGTH;
1870         t >>= BIGNUM_DIGIT_LENGTH;
1871       }
1872       while (scan_a < a_end) {
1873         s += (A * *scan_a);
1874         t -= (C * *scan_a++);
1875         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1876         //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1877         s >>= BIGNUM_DIGIT_LENGTH;
1878         t >>= BIGNUM_DIGIT_LENGTH;
1879       }
1880     }
1881     BIGNUM_ASSERT(s == 0);
1882     BIGNUM_ASSERT(t == 0);
1883
1884     // update size_a and size_b to remove any zeros at end
1885     while (size_a > 0 && *(--scan_a) == 0)
1886       size_a--;
1887     while (size_b > 0 && *(--scan_b) == 0)
1888       size_b--;
1889
1890     BIGNUM_ASSERT(size_a >= size_b);
1891   }
1892
1893   /* a fits into a fixnum, so b must too */
1894   fixnum xx = bignum_to_fixnum(a);
1895   fixnum yy = bignum_to_fixnum(b);
1896   fixnum tt;
1897
1898   /* usual Euclidean algorithm for longs */
1899   while (yy != 0) {
1900     tt = yy;
1901     yy = xx % yy;
1902     xx = tt;
1903   }
1904
1905   return fixnum_to_bignum(xx);
1906 }
1907 #endif
1908
1909 }