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 /* does allocate memory */
 387 fixnum factor_vm::bignum_to_fixnum_strict(bignum* bignum_in) {
 388   fixnum fix = bignum_to_fixnum(bignum_in);
 389   bignum* bignum_out = fixnum_to_bignum(fix);
 390   GC_BIGNUM(bignum_out);
 391   if (bignum_compare(bignum_in, bignum_out) != bignum_comparison_equal) {
 392     general_error(ERROR_OUT_OF_FIXNUM_RANGE, tag<bignum>(bignum_in), false_object);
 393   }
 394   return fix;
 395 }
 396
 397 #define DTB_WRITE_DIGIT(factor)                \
 398   {                                            \
 399     significand *= (factor);                   \
 400     digit = ((bignum_digit_type) significand); \
 401     (*--scan) = digit;                         \
 402     significand -= ((double)digit);            \
 403   }
 404
 405 #define inf std::numeric_limits<double>::infinity()
 406
 407 /* Allocates memory */
 408 bignum* factor_vm::double_to_bignum(double x) {
 409   if (x == inf || x == -inf || x != x)
 410     return (BIGNUM_ZERO());
 411   int exponent;
 412   double significand = (frexp(x, (&exponent)));
 413   if (exponent <= 0)
 414     return (BIGNUM_ZERO());
 415   if (exponent == 1)
 416     return (BIGNUM_ONE(x < 0));
 417   if (significand < 0)
 418     significand = (-significand);
 419   {
 420     bignum_length_type length = (BIGNUM_BITS_TO_DIGITS(exponent));
 421     bignum* result = (allot_bignum(length, (x < 0)));
 422     bignum_digit_type* start = (BIGNUM_START_PTR(result));
 423     bignum_digit_type* scan = (start + length);
 424     bignum_digit_type digit;
 425     int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
 426     if (odd_bits > 0)
 427       DTB_WRITE_DIGIT((fixnum)1 << odd_bits);
 428     while (start < scan) {
 429       if (significand == 0) {
 430         while (start < scan)
 431           (*--scan) = 0;
 432         break;
 433       }
 434       DTB_WRITE_DIGIT(BIGNUM_RADIX);
 435     }
 436     return (result);
 437   }
 438 }
 439
 440 #undef DTB_WRITE_DIGIT
 441
 442 /* Comparisons */
 443
 444 int factor_vm::bignum_equal_p_unsigned(bignum* x, bignum* y) {
 445   bignum_length_type length = (BIGNUM_LENGTH(x));
 446   if (length != (BIGNUM_LENGTH(y)))
 447     return (0);
 448   else {
 449     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 450     bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 451     bignum_digit_type* end_x = (scan_x + length);
 452     while (scan_x < end_x)
 453       if ((*scan_x++) != (*scan_y++))
 454         return (0);
 455     return (1);
 456   }
 457 }
 458
 459 enum bignum_comparison factor_vm::bignum_compare_unsigned(bignum* x,
 460                                                           bignum* y) {
 461   bignum_length_type x_length = (BIGNUM_LENGTH(x));
 462   bignum_length_type y_length = (BIGNUM_LENGTH(y));
 463   if (x_length < y_length)
 464     return (bignum_comparison_less);
 465   if (x_length > y_length)
 466     return (bignum_comparison_greater);
 467   {
 468     bignum_digit_type* start_x = (BIGNUM_START_PTR(x));
 469     bignum_digit_type* scan_x = (start_x + x_length);
 470     bignum_digit_type* scan_y = ((BIGNUM_START_PTR(y)) + y_length);
 471     while (start_x < scan_x) {
 472       bignum_digit_type digit_x = (*--scan_x);
 473       bignum_digit_type digit_y = (*--scan_y);
 474       if (digit_x < digit_y)
 475         return (bignum_comparison_less);
 476       if (digit_x > digit_y)
 477         return (bignum_comparison_greater);
 478     }
 479   }
 480   return (bignum_comparison_equal);
 481 }
 482
 483 /* Addition */
 484
 485 /* Allocates memory */
 486 bignum* factor_vm::bignum_add_unsigned(bignum* x, bignum* y, int negative_p) {
 487   GC_BIGNUM(x);
 488   GC_BIGNUM(y);
 489
 490   if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
 491     std::swap(x, y);
 492   }
 493   {
 494     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 495
 496     bignum* r = (allot_bignum((x_length + 1), negative_p));
 497
 498     bignum_digit_type sum;
 499     bignum_digit_type carry = 0;
 500     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 501     bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
 502     {
 503       bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 504       bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
 505       while (scan_y < end_y) {
 506         sum = ((*scan_x++) + (*scan_y++) + carry);
 507         if (sum < BIGNUM_RADIX) {
 508           (*scan_r++) = sum;
 509           carry = 0;
 510         } else {
 511           (*scan_r++) = (sum - BIGNUM_RADIX);
 512           carry = 1;
 513         }
 514       }
 515     }
 516     {
 517       bignum_digit_type* end_x = ((BIGNUM_START_PTR(x)) + x_length);
 518       if (carry != 0)
 519         while (scan_x < end_x) {
 520           sum = ((*scan_x++) + 1);
 521           if (sum < BIGNUM_RADIX) {
 522             (*scan_r++) = sum;
 523             carry = 0;
 524             break;
 525           } else
 526             (*scan_r++) = (sum - BIGNUM_RADIX);
 527         }
 528       while (scan_x < end_x)
 529         (*scan_r++) = (*scan_x++);
 530     }
 531     if (carry != 0) {
 532       (*scan_r) = 1;
 533       return (r);
 534     }
 535     return (bignum_shorten_length(r, x_length));
 536   }
 537 }
 538
 539 /* Subtraction */
 540
 541 /* Allocates memory */
 542 bignum* factor_vm::bignum_subtract_unsigned(bignum* x, bignum* y) {
 543   GC_BIGNUM(x);
 544   GC_BIGNUM(y);
 545
 546   int negative_p = 0;
 547   switch (bignum_compare_unsigned(x, y)) {
 548     case bignum_comparison_equal:
 549       return (BIGNUM_ZERO());
 550     case bignum_comparison_less: {
 551       std::swap(x, y);
 552     }
 553       negative_p = 1;
 554       break;
 555     case bignum_comparison_greater:
 556       negative_p = 0;
 557       break;
 558   }
 559   {
 560     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 561
 562     bignum* r = (allot_bignum(x_length, negative_p));
 563
 564     bignum_digit_type difference;
 565     bignum_digit_type borrow = 0;
 566     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 567     bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
 568     {
 569       bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 570       bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
 571       while (scan_y < end_y) {
 572         difference = (((*scan_x++) - (*scan_y++)) - borrow);
 573         if (difference < 0) {
 574           (*scan_r++) = (difference + BIGNUM_RADIX);
 575           borrow = 1;
 576         } else {
 577           (*scan_r++) = difference;
 578           borrow = 0;
 579         }
 580       }
 581     }
 582     {
 583       bignum_digit_type* end_x = ((BIGNUM_START_PTR(x)) + x_length);
 584       if (borrow != 0)
 585         while (scan_x < end_x) {
 586           difference = ((*scan_x++) - borrow);
 587           if (difference < 0)
 588             (*scan_r++) = (difference + BIGNUM_RADIX);
 589           else {
 590             (*scan_r++) = difference;
 591             borrow = 0;
 592             break;
 593           }
 594         }
 595       BIGNUM_ASSERT(borrow == 0);
 596       while (scan_x < end_x)
 597         (*scan_r++) = (*scan_x++);
 598     }
 599     return (bignum_trim(r));
 600   }
 601 }
 602
 603 /* Multiplication
 604    Maximum value for product_low or product_high:
 605    ((R * R) + (R * (R - 2)) + (R - 1))
 606    Maximum value for carry: ((R * (R - 1)) + (R - 1))
 607    where R == BIGNUM_RADIX_ROOT */
 608
 609 /* Allocates memory */
 610 bignum* factor_vm::bignum_multiply_unsigned(bignum* x, bignum* y,
 611                                             int negative_p) {
 612   GC_BIGNUM(x);
 613   GC_BIGNUM(y);
 614
 615   if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
 616     std::swap(x, y);
 617   }
 618   {
 619     bignum_digit_type carry;
 620     bignum_digit_type y_digit_low;
 621     bignum_digit_type y_digit_high;
 622     bignum_digit_type x_digit_low;
 623     bignum_digit_type x_digit_high;
 624     bignum_digit_type product_low;
 625     bignum_digit_type* scan_r;
 626     bignum_digit_type* scan_y;
 627     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 628     bignum_length_type y_length = (BIGNUM_LENGTH(y));
 629
 630     bignum* r = (allot_bignum_zeroed((x_length + y_length), negative_p));
 631
 632     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 633     bignum_digit_type* end_x = (scan_x + x_length);
 634     bignum_digit_type* start_y = (BIGNUM_START_PTR(y));
 635     bignum_digit_type* end_y = (start_y + y_length);
 636     bignum_digit_type* start_r = (BIGNUM_START_PTR(r));
 637 #define x_digit x_digit_high
 638 #define y_digit y_digit_high
 639 #define product_high carry
 640     while (scan_x < end_x) {
 641       x_digit = (*scan_x++);
 642       x_digit_low = (HD_LOW(x_digit));
 643       x_digit_high = (HD_HIGH(x_digit));
 644       carry = 0;
 645       scan_y = start_y;
 646       scan_r = (start_r++);
 647       while (scan_y < end_y) {
 648         y_digit = (*scan_y++);
 649         y_digit_low = (HD_LOW(y_digit));
 650         y_digit_high = (HD_HIGH(y_digit));
 651         product_low =
 652             ((*scan_r) + (x_digit_low * y_digit_low) + (HD_LOW(carry)));
 653         product_high =
 654             ((x_digit_high * y_digit_low) + (x_digit_low * y_digit_high) +
 655              (HD_HIGH(product_low)) + (HD_HIGH(carry)));
 656         (*scan_r++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
 657         carry = ((x_digit_high * y_digit_high) + (HD_HIGH(product_high)));
 658       }
 659       (*scan_r) += carry;
 660     }
 661     return (bignum_trim(r));
 662 #undef x_digit
 663 #undef y_digit
 664 #undef product_high
 665   }
 666 }
 667
 668 /* Allocates memory */
 669 bignum* factor_vm::bignum_multiply_unsigned_small_factor(bignum* x,
 670                                                          bignum_digit_type y,
 671                                                          int negative_p) {
 672   GC_BIGNUM(x);
 673
 674   bignum_length_type length_x = (BIGNUM_LENGTH(x));
 675
 676   bignum* p = (allot_bignum((length_x + 1), negative_p));
 677
 678   bignum_destructive_copy(x, p);
 679   (BIGNUM_REF(p, length_x)) = 0;
 680   bignum_destructive_scale_up(p, y);
 681   return (bignum_trim(p));
 682 }
 683
 684 void factor_vm::bignum_destructive_add(bignum* bignum, bignum_digit_type n) {
 685   bignum_digit_type* scan = (BIGNUM_START_PTR(bignum));
 686   bignum_digit_type digit;
 687   digit = ((*scan) + n);
 688   if (digit < BIGNUM_RADIX) {
 689     (*scan) = digit;
 690     return;
 691   }
 692   (*scan++) = (digit - BIGNUM_RADIX);
 693   while (1) {
 694     digit = ((*scan) + 1);
 695     if (digit < BIGNUM_RADIX) {
 696       (*scan) = digit;
 697       return;
 698     }
 699     (*scan++) = (digit - BIGNUM_RADIX);
 700   }
 701 }
 702
 703 void factor_vm::bignum_destructive_scale_up(bignum* bignum,
 704                                             bignum_digit_type factor) {
 705   bignum_digit_type carry = 0;
 706   bignum_digit_type* scan = (BIGNUM_START_PTR(bignum));
 707   bignum_digit_type two_digits;
 708   bignum_digit_type product_low;
 709 #define product_high carry
 710   bignum_digit_type* end = (scan + (BIGNUM_LENGTH(bignum)));
 711   BIGNUM_ASSERT((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
 712   while (scan < end) {
 713     two_digits = (*scan);
 714     product_low = ((factor * (HD_LOW(two_digits))) + (HD_LOW(carry)));
 715     product_high = ((factor * (HD_HIGH(two_digits))) + (HD_HIGH(product_low)) +
 716                     (HD_HIGH(carry)));
 717     (*scan++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
 718     carry = (HD_HIGH(product_high));
 719   }
 720   /* A carry here would be an overflow, i.e. it would not fit.
 721            Hopefully the callers allocate enough space that this will
 722            never happen.
 723         */
 724   BIGNUM_ASSERT(carry == 0);
 725   return;
 726 #undef product_high
 727 }
 728
 729 /* Division */
 730
 731 /* For help understanding this algorithm, see:
 732    Knuth, Donald E., "The Art of Computer Programming",
 733    volume 2, "Seminumerical Algorithms"
 734    section 4.3.1, "Multiple-Precision Arithmetic". */
 735
 736 /* Allocates memory */
 737 void factor_vm::bignum_divide_unsigned_large_denominator(
 738     bignum* numerator, bignum* denominator, bignum** quotient,
 739     bignum** remainder, int q_negative_p, int r_negative_p) {
 740   GC_BIGNUM(numerator);
 741   GC_BIGNUM(denominator);
 742
 743   bignum_length_type length_n = ((BIGNUM_LENGTH(numerator)) + 1);
 744   bignum_length_type length_d = (BIGNUM_LENGTH(denominator));
 745
 746   bignum* q = ((quotient != ((bignum**)0))
 747                    ? (allot_bignum((length_n - length_d), q_negative_p))
 748                    : BIGNUM_OUT_OF_BAND);
 749   GC_BIGNUM(q);
 750
 751   bignum* u = (allot_bignum(length_n, r_negative_p));
 752   GC_BIGNUM(u);
 753
 754   int shift = 0;
 755   BIGNUM_ASSERT(length_d > 1);
 756   {
 757     bignum_digit_type v1 = (BIGNUM_REF((denominator), (length_d - 1)));
 758     while (v1 < (BIGNUM_RADIX / 2)) {
 759       v1 <<= 1;
 760       shift += 1;
 761     }
 762   }
 763   if (shift == 0) {
 764     bignum_destructive_copy(numerator, u);
 765     (BIGNUM_REF(u, (length_n - 1))) = 0;
 766     bignum_divide_unsigned_normalized(u, denominator, q);
 767   } else {
 768     bignum* v = (allot_bignum(length_d, 0));
 769
 770     bignum_destructive_normalization(numerator, u, shift);
 771     bignum_destructive_normalization(denominator, v, shift);
 772     bignum_divide_unsigned_normalized(u, v, q);
 773     if (remainder != ((bignum**)0))
 774       bignum_destructive_unnormalization(u, shift);
 775   }
 776
 777   if (q)
 778     q = bignum_trim(q);
 779
 780   u = bignum_trim(u);
 781
 782   if (quotient != ((bignum**)0))
 783     (*quotient) = q;
 784
 785   if (remainder != ((bignum**)0))
 786     (*remainder) = u;
 787
 788   return;
 789 }
 790
 791 void factor_vm::bignum_divide_unsigned_normalized(bignum* u, bignum* v,
 792                                                   bignum* q) {
 793   bignum_length_type u_length = (BIGNUM_LENGTH(u));
 794   bignum_length_type v_length = (BIGNUM_LENGTH(v));
 795   bignum_digit_type* u_start = (BIGNUM_START_PTR(u));
 796   bignum_digit_type* u_scan = (u_start + u_length);
 797   bignum_digit_type* u_scan_limit = (u_start + v_length);
 798   bignum_digit_type* u_scan_start = (u_scan - v_length);
 799   bignum_digit_type* v_start = (BIGNUM_START_PTR(v));
 800   bignum_digit_type* v_end = (v_start + v_length);
 801   bignum_digit_type* q_scan = NULL;
 802   bignum_digit_type v1 = (v_end[-1]);
 803   bignum_digit_type v2 = (v_end[-2]);
 804   bignum_digit_type ph; /* high half of double-digit product */
 805   bignum_digit_type pl; /* low half of double-digit product */
 806   bignum_digit_type guess;
 807   bignum_digit_type gh; /* high half-digit of guess */
 808   bignum_digit_type ch; /* high half of double-digit comparand */
 809   bignum_digit_type v2l = (HD_LOW(v2));
 810   bignum_digit_type v2h = (HD_HIGH(v2));
 811   bignum_digit_type cl; /* low half of double-digit comparand */
 812 #define gl ph           /* low half-digit of guess */
 813 #define uj pl
 814 #define qj ph
 815   bignum_digit_type gm; /* memory loc for reference parameter */
 816   if (q != BIGNUM_OUT_OF_BAND)
 817     q_scan = ((BIGNUM_START_PTR(q)) + (BIGNUM_LENGTH(q)));
 818   while (u_scan_limit < u_scan) {
 819     uj = (*--u_scan);
 820     if (uj != v1) {
 821       /* comparand =
 822                            (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
 823                            guess = (((uj * BIGNUM_RADIX) + uj1) / v1); */
 824       cl = (u_scan[-2]);
 825       ch = (bignum_digit_divide(uj, (u_scan[-1]), v1, (&gm)));
 826       guess = gm;
 827     } else {
 828       cl = (u_scan[-2]);
 829       ch = ((u_scan[-1]) + v1);
 830       guess = (BIGNUM_RADIX - 1);
 831     }
 832     while (1) {
 833       /* product = (guess * v2); */
 834       gl = (HD_LOW(guess));
 835       gh = (HD_HIGH(guess));
 836       pl = (v2l * gl);
 837       ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH(pl)));
 838       pl = (HD_CONS((HD_LOW(ph)), (HD_LOW(pl))));
 839       ph = ((v2h * gh) + (HD_HIGH(ph)));
 840       /* if (comparand >= product) */
 841       if ((ch > ph) || ((ch == ph) && (cl >= pl)))
 842         break;
 843       guess -= 1;
 844       /* comparand += (v1 << BIGNUM_DIGIT_LENGTH) */
 845       ch += v1;
 846       /* if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX)) */
 847       if (ch >= BIGNUM_RADIX)
 848         break;
 849     }
 850     qj = (bignum_divide_subtract(v_start, v_end, guess, (--u_scan_start)));
 851     if (q != BIGNUM_OUT_OF_BAND)
 852       (*--q_scan) = qj;
 853   }
 854   return;
 855 #undef gl
 856 #undef uj
 857 #undef qj
 858 }
 859
 860 bignum_digit_type factor_vm::bignum_divide_subtract(
 861     bignum_digit_type* v_start, bignum_digit_type* v_end,
 862     bignum_digit_type guess, bignum_digit_type* u_start) {
 863   bignum_digit_type* v_scan = v_start;
 864   bignum_digit_type* u_scan = u_start;
 865   bignum_digit_type carry = 0;
 866   if (guess == 0)
 867     return (0);
 868   {
 869     bignum_digit_type gl = (HD_LOW(guess));
 870     bignum_digit_type gh = (HD_HIGH(guess));
 871     bignum_digit_type v;
 872     bignum_digit_type pl;
 873     bignum_digit_type vl;
 874 #define vh v
 875 #define ph carry
 876 #define diff pl
 877     while (v_scan < v_end) {
 878       v = (*v_scan++);
 879       vl = (HD_LOW(v));
 880       vh = (HD_HIGH(v));
 881       pl = ((vl * gl) + (HD_LOW(carry)));
 882       ph = ((vl * gh) + (vh * gl) + (HD_HIGH(pl)) + (HD_HIGH(carry)));
 883       diff = ((*u_scan) - (HD_CONS((HD_LOW(ph)), (HD_LOW(pl)))));
 884       if (diff < 0) {
 885         (*u_scan++) = (diff + BIGNUM_RADIX);
 886         carry = ((vh * gh) + (HD_HIGH(ph)) + 1);
 887       } else {
 888         (*u_scan++) = diff;
 889         carry = ((vh * gh) + (HD_HIGH(ph)));
 890       }
 891     }
 892     if (carry == 0)
 893       return (guess);
 894     diff = ((*u_scan) - carry);
 895     if (diff < 0)
 896       (*u_scan) = (diff + BIGNUM_RADIX);
 897     else {
 898       (*u_scan) = diff;
 899       return (guess);
 900     }
 901 #undef vh
 902 #undef ph
 903 #undef diff
 904   }
 905   /* Subtraction generated carry, implying guess is one too large.
 906            Add v back in to bring it back down. */
 907   v_scan = v_start;
 908   u_scan = u_start;
 909   carry = 0;
 910   while (v_scan < v_end) {
 911     bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
 912     if (sum < BIGNUM_RADIX) {
 913       (*u_scan++) = sum;
 914       carry = 0;
 915     } else {
 916       (*u_scan++) = (sum - BIGNUM_RADIX);
 917       carry = 1;
 918     }
 919   }
 920   if (carry == 1) {
 921     bignum_digit_type sum = ((*u_scan) + carry);
 922     (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
 923   }
 924   return (guess - 1);
 925 }
 926
 927 /* Allocates memory */
 928 void factor_vm::bignum_divide_unsigned_medium_denominator(
 929     bignum* numerator, bignum_digit_type denominator, bignum** quotient,
 930     bignum** remainder, int q_negative_p, int r_negative_p) {
 931   GC_BIGNUM(numerator);
 932
 933   bignum_length_type length_n = (BIGNUM_LENGTH(numerator));
 934   bignum_length_type length_q;
 935   bignum* q = NULL;
 936   GC_BIGNUM(q);
 937
 938   int shift = 0;
 939   /* Because `bignum_digit_divide' requires a normalized denominator. */
 940   while (denominator < (BIGNUM_RADIX / 2)) {
 941     denominator <<= 1;
 942     shift += 1;
 943   }
 944   if (shift == 0) {
 945     length_q = length_n;
 946
 947     q = (allot_bignum(length_q, q_negative_p));
 948     bignum_destructive_copy(numerator, q);
 949   } else {
 950     length_q = (length_n + 1);
 951
 952     q = (allot_bignum(length_q, q_negative_p));
 953     bignum_destructive_normalization(numerator, q, shift);
 954   }
 955   {
 956     bignum_digit_type r = 0;
 957     bignum_digit_type* start = (BIGNUM_START_PTR(q));
 958     bignum_digit_type* scan = (start + length_q);
 959     bignum_digit_type qj;
 960
 961     while (start < scan) {
 962       r = (bignum_digit_divide(r, (*--scan), denominator, (&qj)));
 963       (*scan) = qj;
 964     }
 965
 966     q = bignum_trim(q);
 967
 968     if (remainder != ((bignum**)0)) {
 969       if (shift != 0)
 970         r >>= shift;
 971
 972       (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
 973     }
 974
 975     if (quotient != ((bignum**)0))
 976       (*quotient) = q;
 977   }
 978   return;
 979 }
 980
 981 void factor_vm::bignum_destructive_normalization(bignum* source, bignum* target,
 982                                                  int shift_left) {
 983   bignum_digit_type digit;
 984   bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
 985   bignum_digit_type carry = 0;
 986   bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
 987   bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
 988   bignum_digit_type* end_target = (scan_target + (BIGNUM_LENGTH(target)));
 989   int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
 990   bignum_digit_type mask = (((cell)1 << shift_right) - 1);
 991   while (scan_source < end_source) {
 992     digit = (*scan_source++);
 993     (*scan_target++) = (((digit & mask) << shift_left) | carry);
 994     carry = (digit >> shift_right);
 995   }
 996   if (scan_target < end_target)
 997     (*scan_target) = carry;
 998   else
 999     BIGNUM_ASSERT(carry == 0);
1000   return;
1001 }
1002
1003 void factor_vm::bignum_destructive_unnormalization(bignum* bignum,
1004                                                    int shift_right) {
1005   bignum_digit_type* start = (BIGNUM_START_PTR(bignum));
1006   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bignum)));
1007   bignum_digit_type digit;
1008   bignum_digit_type carry = 0;
1009   int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
1010   bignum_digit_type mask = (((fixnum)1 << shift_right) - 1);
1011   while (start < scan) {
1012     digit = (*--scan);
1013     (*scan) = ((digit >> shift_right) | carry);
1014     carry = ((digit & mask) << shift_left);
1015   }
1016   BIGNUM_ASSERT(carry == 0);
1017   return;
1018 }
1019
1020 /* This is a reduced version of the division algorithm, applied to the
1021    case of dividing two bignum digits by one bignum digit.  It is
1022    assumed that the numerator, denominator are normalized. */
1023
1024 #define BDD_STEP(qn, j)                                          \
1025   {                                                              \
1026     uj = (u[j]);                                                 \
1027     if (uj != v1) {                                              \
1028       uj_uj1 = (HD_CONS(uj, (u[j + 1])));                        \
1029       guess = (uj_uj1 / v1);                                     \
1030       comparand = (HD_CONS((uj_uj1 % v1), (u[j + 2])));          \
1031     } else {                                                     \
1032       guess = (BIGNUM_RADIX_ROOT - 1);                           \
1033       comparand = (HD_CONS(((u[j + 1]) + v1), (u[j + 2])));      \
1034     }                                                            \
1035     while ((guess * v2) > comparand) {                           \
1036       guess -= 1;                                                \
1037       comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH);             \
1038       if (comparand >= BIGNUM_RADIX)                             \
1039         break;                                                   \
1040     }                                                            \
1041     qn = (bignum_digit_divide_subtract(v1, v2, guess, (&u[j]))); \
1042   }
1043
1044 bignum_digit_type factor_vm::bignum_digit_divide(
1045     bignum_digit_type uh, bignum_digit_type ul, bignum_digit_type v,
1046     bignum_digit_type* q) /* return value */
1047     {
1048   bignum_digit_type guess;
1049   bignum_digit_type comparand;
1050   bignum_digit_type v1 = (HD_HIGH(v));
1051   bignum_digit_type v2 = (HD_LOW(v));
1052   bignum_digit_type uj;
1053   bignum_digit_type uj_uj1;
1054   bignum_digit_type q1;
1055   bignum_digit_type q2;
1056   bignum_digit_type u[4];
1057   if (uh == 0) {
1058     if (ul < v) {
1059       (*q) = 0;
1060       return (ul);
1061     } else if (ul == v) {
1062       (*q) = 1;
1063       return (0);
1064     }
1065   }
1066   (u[0]) = (HD_HIGH(uh));
1067   (u[1]) = (HD_LOW(uh));
1068   (u[2]) = (HD_HIGH(ul));
1069   (u[3]) = (HD_LOW(ul));
1070   v1 = (HD_HIGH(v));
1071   v2 = (HD_LOW(v));
1072   BDD_STEP(q1, 0);
1073   BDD_STEP(q2, 1);
1074   (*q) = (HD_CONS(q1, q2));
1075   return (HD_CONS((u[2]), (u[3])));
1076 }
1077
1078 #undef BDD_STEP
1079
1080 #define BDDS_MULSUB(vn, un, carry_in)    \
1081   {                                      \
1082     product = ((vn * guess) + carry_in); \
1083     diff = (un - (HD_LOW(product)));     \
1084     if (diff < 0) {                      \
1085       un = (diff + BIGNUM_RADIX_ROOT);   \
1086       carry = ((HD_HIGH(product)) + 1);  \
1087     } else {                             \
1088       un = diff;                         \
1089       carry = (HD_HIGH(product));        \
1090     }                                    \
1091   }
1092
1093 #define BDDS_ADD(vn, un, carry_in)    \
1094   {                                   \
1095     sum = (vn + un + carry_in);       \
1096     if (sum < BIGNUM_RADIX_ROOT) {    \
1097       un = sum;                       \
1098       carry = 0;                      \
1099     } else {                          \
1100       un = (sum - BIGNUM_RADIX_ROOT); \
1101       carry = 1;                      \
1102     }                                 \
1103   }
1104
1105 bignum_digit_type factor_vm::bignum_digit_divide_subtract(
1106     bignum_digit_type v1, bignum_digit_type v2, bignum_digit_type guess,
1107     bignum_digit_type* u) {
1108   {
1109     bignum_digit_type product;
1110     bignum_digit_type diff;
1111     bignum_digit_type carry;
1112     BDDS_MULSUB(v2, (u[2]), 0);
1113     BDDS_MULSUB(v1, (u[1]), carry);
1114     if (carry == 0)
1115       return (guess);
1116     diff = ((u[0]) - carry);
1117     if (diff < 0)
1118       (u[0]) = (diff + BIGNUM_RADIX);
1119     else {
1120       (u[0]) = diff;
1121       return (guess);
1122     }
1123   }
1124   {
1125     bignum_digit_type sum;
1126     bignum_digit_type carry;
1127     BDDS_ADD(v2, (u[2]), 0);
1128     BDDS_ADD(v1, (u[1]), carry);
1129     if (carry == 1)
1130       (u[0]) += 1;
1131   }
1132   return (guess - 1);
1133 }
1134
1135 #undef BDDS_MULSUB
1136 #undef BDDS_ADD
1137
1138 /* Allocates memory */
1139 void factor_vm::bignum_divide_unsigned_small_denominator(
1140     bignum* numerator, bignum_digit_type denominator, bignum** quotient,
1141     bignum** remainder, int q_negative_p, int r_negative_p) {
1142   GC_BIGNUM(numerator);
1143
1144   bignum* q = (bignum_new_sign(numerator, q_negative_p));
1145   GC_BIGNUM(q);
1146
1147   bignum_digit_type r = (bignum_destructive_scale_down(q, denominator));
1148
1149   q = (bignum_trim(q));
1150
1151   if (remainder != ((bignum**)0))
1152     (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
1153
1154   (*quotient) = q;
1155
1156   return;
1157 }
1158
1159 /* Given (denominator > 1), it is fairly easy to show that
1160    (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
1161    that all digits are < BIGNUM_RADIX. */
1162
1163 bignum_digit_type factor_vm::bignum_destructive_scale_down(
1164     bignum* bignum, bignum_digit_type denominator) {
1165   bignum_digit_type numerator;
1166   bignum_digit_type remainder = 0;
1167   bignum_digit_type two_digits;
1168 #define quotient_high remainder
1169   bignum_digit_type* start = (BIGNUM_START_PTR(bignum));
1170   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bignum)));
1171   BIGNUM_ASSERT((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
1172   while (start < scan) {
1173     two_digits = (*--scan);
1174     numerator = (HD_CONS(remainder, (HD_HIGH(two_digits))));
1175     quotient_high = (numerator / denominator);
1176     numerator = (HD_CONS((numerator % denominator), (HD_LOW(two_digits))));
1177     (*scan) = (HD_CONS(quotient_high, (numerator / denominator)));
1178     remainder = (numerator % denominator);
1179   }
1180   return (remainder);
1181 #undef quotient_high
1182 }
1183
1184 /* Allocates memory */
1185 bignum* factor_vm::bignum_remainder_unsigned_small_denominator(
1186     bignum* n, bignum_digit_type d, int negative_p) {
1187   bignum_digit_type two_digits;
1188   bignum_digit_type* start = (BIGNUM_START_PTR(n));
1189   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(n)));
1190   bignum_digit_type r = 0;
1191   BIGNUM_ASSERT((d > 1) && (d < BIGNUM_RADIX_ROOT));
1192   while (start < scan) {
1193     two_digits = (*--scan);
1194     r = ((HD_CONS(((HD_CONS(r, (HD_HIGH(two_digits)))) % d),
1195                   (HD_LOW(two_digits)))) %
1196          d);
1197   }
1198   return (bignum_digit_to_bignum(r, negative_p));
1199 }
1200
1201 /* Allocates memory */
1202 bignum* factor_vm::bignum_digit_to_bignum(bignum_digit_type digit,
1203                                           int negative_p) {
1204   if (digit == 0)
1205     return (BIGNUM_ZERO());
1206   else {
1207     bignum* result = (allot_bignum(1, negative_p));
1208     (BIGNUM_REF(result, 0)) = digit;
1209     return (result);
1210   }
1211 }
1212
1213 /* Allocates memory */
1214 bignum* factor_vm::allot_bignum(bignum_length_type length, int negative_p) {
1215   BIGNUM_ASSERT((length >= 0) || (length < BIGNUM_RADIX));
1216   bignum* result = allot_uninitialized_array<bignum>(length + 1);
1217   BIGNUM_SET_NEGATIVE_P(result, negative_p);
1218   return (result);
1219 }
1220
1221 /* Allocates memory */
1222 bignum* factor_vm::allot_bignum_zeroed(bignum_length_type length,
1223                                        int negative_p) {
1224   bignum* result = allot_bignum(length, negative_p);
1225   bignum_digit_type* scan = (BIGNUM_START_PTR(result));
1226   bignum_digit_type* end = (scan + length);
1227   while (scan < end)
1228     (*scan++) = 0;
1229   return (result);
1230 }
1231
1232 /* can allocate if not in nursery or size is larger */
1233 /* Allocates memory conditionally */
1234 #define BIGNUM_REDUCE_LENGTH(source, length) \
1235   source = reallot_array(source, length + 1)
1236
1237 /* Allocates memory */
1238 bignum* factor_vm::bignum_shorten_length(bignum* bignum,
1239                                          bignum_length_type length) {
1240   bignum_length_type current_length = (BIGNUM_LENGTH(bignum));
1241   BIGNUM_ASSERT((length >= 0) || (length <= current_length));
1242   if (length < current_length) {
1243     GC_BIGNUM(bignum);
1244     BIGNUM_REDUCE_LENGTH(bignum, length);
1245     BIGNUM_SET_NEGATIVE_P(bignum, (length != 0) && (BIGNUM_NEGATIVE_P(bignum)));
1246   }
1247   return (bignum);
1248 }
1249
1250 /* Allocates memory */
1251 bignum* factor_vm::bignum_trim(bignum* bignum) {
1252   bignum_digit_type* start = (BIGNUM_START_PTR(bignum));
1253   bignum_digit_type* end = (start + (BIGNUM_LENGTH(bignum)));
1254   bignum_digit_type* scan = end;
1255   while ((start <= scan) && ((*--scan) == 0))
1256     ;
1257   scan += 1;
1258   if (scan < end) {
1259     GC_BIGNUM(bignum);
1260     bignum_length_type length = (scan - start);
1261     BIGNUM_REDUCE_LENGTH(bignum, length);
1262     BIGNUM_SET_NEGATIVE_P(bignum, (length != 0) && (BIGNUM_NEGATIVE_P(bignum)));
1263   }
1264   return (bignum);
1265 }
1266
1267 /* Copying */
1268
1269 /* Allocates memory */
1270 bignum* factor_vm::bignum_new_sign(bignum* x, int negative_p) {
1271   GC_BIGNUM(x);
1272   bignum* result = (allot_bignum((BIGNUM_LENGTH(x)), negative_p));
1273
1274   bignum_destructive_copy(x, result);
1275   return (result);
1276 }
1277
1278 /* Allocates memory */
1279 bignum* factor_vm::bignum_maybe_new_sign(bignum* x, int negative_p) {
1280   if ((BIGNUM_NEGATIVE_P(x)) ? negative_p : (!negative_p))
1281     return (x);
1282   else {
1283     GC_BIGNUM(x);
1284     bignum* result = (allot_bignum((BIGNUM_LENGTH(x)), negative_p));
1285     bignum_destructive_copy(x, result);
1286     return (result);
1287   }
1288 }
1289
1290 void factor_vm::bignum_destructive_copy(bignum* source, bignum* target) {
1291   bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
1292   bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1293   bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1294   while (scan_source < end_source)
1295     (*scan_target++) = (*scan_source++);
1296   return;
1297 }
1298
1299 /*
1300  * Added bitwise operations (and oddp).
1301  */
1302
1303 /* Allocates memory */
1304 bignum* factor_vm::bignum_bitwise_not(bignum* x) {
1305   GC_BIGNUM(x);
1306
1307   bignum_length_type size = BIGNUM_LENGTH(x);
1308   bignum_digit_type* scan_x, *end_x, *scan_y;
1309   bignum* y;
1310   int carry = 1;
1311
1312   if (BIGNUM_NEGATIVE_P(x)) {
1313     y = allot_bignum(size, 0);
1314     scan_x = BIGNUM_START_PTR(x);
1315     end_x = scan_x + size;
1316     scan_y = BIGNUM_START_PTR(y);
1317     while (scan_x < end_x) {
1318       if (*scan_x == 0) {
1319         *scan_y++ = BIGNUM_RADIX - 1;
1320         scan_x++;
1321       } else {
1322         *scan_y++ = *scan_x++ - 1;
1323         carry = 0;
1324         break;
1325       }
1326     }
1327   } else {
1328     y = allot_bignum(size, 1);
1329     scan_x = BIGNUM_START_PTR(x);
1330     end_x = scan_x + size;
1331     scan_y = BIGNUM_START_PTR(y);
1332     while (scan_x < end_x) {
1333       if (*scan_x == (BIGNUM_RADIX - 1)) {
1334         *scan_y++ = 0;
1335         scan_x++;
1336       } else {
1337         *scan_y++ = *scan_x++ + 1;
1338         carry = 0;
1339         break;
1340       }
1341     }
1342   }
1343
1344   while (scan_x < end_x) {
1345     *scan_y++ = *scan_x++;
1346   }
1347
1348   if (carry) {
1349     GC_BIGNUM(y);
1350     x = allot_bignum(size + 1, BIGNUM_NEGATIVE_P(y));
1351     bignum_destructive_copy(y, x);
1352     scan_x = BIGNUM_START_PTR(x);
1353     *(scan_x + size) = 1;
1354     return x;
1355   } else {
1356     return bignum_trim(y);
1357   }
1358 }
1359
1360 /* Allocates memory */
1361 bignum* factor_vm::bignum_arithmetic_shift(bignum* arg1, fixnum n) {
1362   if (BIGNUM_NEGATIVE_P(arg1) && n < 0)
1363     return bignum_bitwise_not(
1364         bignum_magnitude_ash(bignum_bitwise_not(arg1), n));
1365   else
1366     return bignum_magnitude_ash(arg1, n);
1367 }
1368
1369 #define AND_OP 0
1370 #define IOR_OP 1
1371 #define XOR_OP 2
1372
1373 /* Allocates memory */
1374 bignum* factor_vm::bignum_bitwise_and(bignum* arg1, bignum* arg2) {
1375   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1376               ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2)
1377               : bignum_posneg_bitwise_op(AND_OP, arg2, arg1)
1378               : (BIGNUM_NEGATIVE_P(arg2))
1379               ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2)
1380               : bignum_pospos_bitwise_op(AND_OP, arg1, arg2));
1381 }
1382
1383 /* Allocates memory */
1384 bignum* factor_vm::bignum_bitwise_ior(bignum* arg1, bignum* arg2) {
1385   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1386               ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2)
1387               : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1)
1388               : (BIGNUM_NEGATIVE_P(arg2))
1389               ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2)
1390               : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2));
1391 }
1392
1393 /* Allocates memory */
1394 bignum* factor_vm::bignum_bitwise_xor(bignum* arg1, bignum* arg2) {
1395   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1396               ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2)
1397               : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1)
1398               : (BIGNUM_NEGATIVE_P(arg2))
1399               ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2)
1400               : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2));
1401 }
1402
1403 /* Allocates memory */
1404 /* ash for the magnitude */
1405 /* assume arg1 is a big number, n is a long */
1406 bignum* factor_vm::bignum_magnitude_ash(bignum* arg1, fixnum n) {
1407   GC_BIGNUM(arg1);
1408
1409   bignum* result = NULL;
1410   bignum_digit_type* scan1;
1411   bignum_digit_type* scanr;
1412   bignum_digit_type* end;
1413
1414   fixnum digit_offset, bit_offset;
1415
1416   if (BIGNUM_ZERO_P(arg1))
1417     return (arg1);
1418
1419   if (n > 0) {
1420     digit_offset = n / BIGNUM_DIGIT_LENGTH;
1421     bit_offset = n % BIGNUM_DIGIT_LENGTH;
1422
1423     result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) + digit_offset + 1,
1424                                  BIGNUM_NEGATIVE_P(arg1));
1425
1426     scanr = BIGNUM_START_PTR(result) + digit_offset;
1427     scan1 = BIGNUM_START_PTR(arg1);
1428     end = scan1 + BIGNUM_LENGTH(arg1);
1429
1430     while (scan1 < end) {
1431       *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
1432       *scanr = *scanr & BIGNUM_DIGIT_MASK;
1433       scanr++;
1434       *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
1435       *scanr = *scanr & BIGNUM_DIGIT_MASK;
1436     }
1437   } else if (n < 0 && (-n >= (BIGNUM_LENGTH(arg1) * (bignum_length_type)
1438                               BIGNUM_DIGIT_LENGTH)))
1439     result = BIGNUM_ZERO();
1440
1441   else if (n < 0) {
1442     digit_offset = -n / BIGNUM_DIGIT_LENGTH;
1443     bit_offset = -n % BIGNUM_DIGIT_LENGTH;
1444
1445     result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) - digit_offset,
1446                                  BIGNUM_NEGATIVE_P(arg1));
1447
1448     scanr = BIGNUM_START_PTR(result);
1449     scan1 = BIGNUM_START_PTR(arg1) + digit_offset;
1450     end = scanr + BIGNUM_LENGTH(result) - 1;
1451
1452     while (scanr < end) {
1453       *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1454       *scanr = (*scanr | *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) &
1455                BIGNUM_DIGIT_MASK;
1456       scanr++;
1457     }
1458     *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1459   } else if (n == 0)
1460     result = arg1;
1461
1462   return (bignum_trim(result));
1463 }
1464
1465 /* Allocates memory */
1466 bignum* factor_vm::bignum_pospos_bitwise_op(int op, bignum* arg1,
1467                                             bignum* arg2) {
1468   GC_BIGNUM(arg1);
1469   GC_BIGNUM(arg2);
1470
1471   bignum* result;
1472   bignum_length_type max_length;
1473
1474   bignum_digit_type* scan1, *end1, digit1;
1475   bignum_digit_type* scan2, *end2, digit2;
1476   bignum_digit_type* scanr, *endr;
1477
1478   max_length =
1479       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1)
1480                                                   : BIGNUM_LENGTH(arg2);
1481
1482   result = allot_bignum(max_length, 0);
1483
1484   scanr = BIGNUM_START_PTR(result);
1485   scan1 = BIGNUM_START_PTR(arg1);
1486   scan2 = BIGNUM_START_PTR(arg2);
1487   endr = scanr + max_length;
1488   end1 = scan1 + BIGNUM_LENGTH(arg1);
1489   end2 = scan2 + BIGNUM_LENGTH(arg2);
1490
1491   while (scanr < endr) {
1492     digit1 = (scan1 < end1) ? *scan1++ : 0;
1493     digit2 = (scan2 < end2) ? *scan2++ : 0;
1494     *scanr++ =
1495         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1496                                                           : digit1 ^ digit2;
1497   }
1498   return bignum_trim(result);
1499 }
1500
1501 /* Allocates memory */
1502 bignum* factor_vm::bignum_posneg_bitwise_op(int op, bignum* arg1,
1503                                             bignum* arg2) {
1504   GC_BIGNUM(arg1);
1505   GC_BIGNUM(arg2);
1506
1507   bignum* result;
1508   bignum_length_type max_length;
1509
1510   bignum_digit_type* scan1, *end1, digit1;
1511   bignum_digit_type* scan2, *end2, digit2, carry2;
1512   bignum_digit_type* scanr, *endr;
1513
1514   char neg_p = op == IOR_OP || op == XOR_OP;
1515
1516   max_length =
1517       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1) ? BIGNUM_LENGTH(arg1)
1518                                                       : BIGNUM_LENGTH(arg2) + 1;
1519
1520   result = allot_bignum(max_length, neg_p);
1521
1522   scanr = BIGNUM_START_PTR(result);
1523   scan1 = BIGNUM_START_PTR(arg1);
1524   scan2 = BIGNUM_START_PTR(arg2);
1525   endr = scanr + max_length;
1526   end1 = scan1 + BIGNUM_LENGTH(arg1);
1527   end2 = scan2 + BIGNUM_LENGTH(arg2);
1528
1529   carry2 = 1;
1530
1531   while (scanr < endr) {
1532     digit1 = (scan1 < end1) ? *scan1++ : 0;
1533     digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1534
1535     if (digit2 < BIGNUM_RADIX)
1536       carry2 = 0;
1537     else {
1538       digit2 = (digit2 - BIGNUM_RADIX);
1539       carry2 = 1;
1540     }
1541
1542     *scanr++ =
1543         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1544                                                           : digit1 ^ digit2;
1545   }
1546
1547   if (neg_p)
1548     bignum_negate_magnitude(result);
1549
1550   return bignum_trim(result);
1551 }
1552
1553 /* Allocates memory */
1554 bignum* factor_vm::bignum_negneg_bitwise_op(int op, bignum* arg1,
1555                                             bignum* arg2) {
1556   GC_BIGNUM(arg1);
1557   GC_BIGNUM(arg2);
1558
1559   bignum* result;
1560   bignum_length_type max_length;
1561
1562   bignum_digit_type* scan1, *end1, digit1, carry1;
1563   bignum_digit_type* scan2, *end2, digit2, carry2;
1564   bignum_digit_type* scanr, *endr;
1565
1566   char neg_p = op == AND_OP || op == IOR_OP;
1567
1568   max_length =
1569       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) + 1
1570                                                   : BIGNUM_LENGTH(arg2) + 1;
1571
1572   result = allot_bignum(max_length, neg_p);
1573
1574   scanr = BIGNUM_START_PTR(result);
1575   scan1 = BIGNUM_START_PTR(arg1);
1576   scan2 = BIGNUM_START_PTR(arg2);
1577   endr = scanr + max_length;
1578   end1 = scan1 + BIGNUM_LENGTH(arg1);
1579   end2 = scan2 + BIGNUM_LENGTH(arg2);
1580
1581   carry1 = 1;
1582   carry2 = 1;
1583
1584   while (scanr < endr) {
1585     digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
1586     digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1587
1588     if (digit1 < BIGNUM_RADIX)
1589       carry1 = 0;
1590     else {
1591       digit1 = (digit1 - BIGNUM_RADIX);
1592       carry1 = 1;
1593     }
1594
1595     if (digit2 < BIGNUM_RADIX)
1596       carry2 = 0;
1597     else {
1598       digit2 = (digit2 - BIGNUM_RADIX);
1599       carry2 = 1;
1600     }
1601
1602     *scanr++ =
1603         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1604                                                           : digit1 ^ digit2;
1605   }
1606
1607   if (neg_p)
1608     bignum_negate_magnitude(result);
1609
1610   return bignum_trim(result);
1611 }
1612
1613 void factor_vm::bignum_negate_magnitude(bignum* arg) {
1614   bignum_digit_type* scan;
1615   bignum_digit_type* end;
1616   bignum_digit_type digit;
1617   bignum_digit_type carry;
1618
1619   scan = BIGNUM_START_PTR(arg);
1620   end = scan + BIGNUM_LENGTH(arg);
1621
1622   carry = 1;
1623
1624   while (scan < end) {
1625     digit = (~ * scan & BIGNUM_DIGIT_MASK) + carry;
1626
1627     if (digit < BIGNUM_RADIX)
1628       carry = 0;
1629     else {
1630       digit = (digit - BIGNUM_RADIX);
1631       carry = 1;
1632     }
1633
1634     *scan++ = digit;
1635   }
1636 }
1637
1638 /* Allocates memory */
1639 bignum* factor_vm::bignum_integer_length(bignum* x) {
1640   GC_BIGNUM(x);
1641
1642   bignum_length_type index = ((BIGNUM_LENGTH(x)) - 1);
1643   bignum_digit_type digit = (BIGNUM_REF(x, index));
1644   bignum_digit_type carry = 0;
1645   bignum* result;
1646
1647   while (digit > 1) {
1648     carry += 1;
1649     digit >>= 1;
1650   }
1651
1652   if (index < BIGNUM_RADIX_ROOT) {
1653      result = allot_bignum(1, 0);
1654      (BIGNUM_REF(result, 0)) = (index * BIGNUM_DIGIT_LENGTH) + carry;
1655   } else {
1656      result = allot_bignum(2, 0);
1657      (BIGNUM_REF(result, 0)) = index;
1658      (BIGNUM_REF(result, 1)) = 0;
1659      bignum_destructive_scale_up(result, BIGNUM_DIGIT_LENGTH);
1660      bignum_destructive_add(result, carry);
1661   }
1662   return (bignum_trim(result));
1663 }
1664
1665 /* Allocates memory */
1666 int factor_vm::bignum_logbitp(int shift, bignum* arg) {
1667   return ((BIGNUM_NEGATIVE_P(arg))
1668               ? !bignum_unsigned_logbitp(shift, bignum_bitwise_not(arg))
1669               : bignum_unsigned_logbitp(shift, arg));
1670 }
1671
1672 int factor_vm::bignum_unsigned_logbitp(int shift, bignum* bignum) {
1673   bignum_length_type len = (BIGNUM_LENGTH(bignum));
1674   int index = shift / BIGNUM_DIGIT_LENGTH;
1675   if (index >= len)
1676     return 0;
1677   bignum_digit_type digit = (BIGNUM_REF(bignum, index));
1678   int p = shift % BIGNUM_DIGIT_LENGTH;
1679   bignum_digit_type mask = ((fixnum)1) << p;
1680   return (digit & mask) ? 1 : 0;
1681 }
1682
1683 #ifdef _WIN64
1684 /* Allocates memory */
1685 bignum* factor_vm::bignum_gcd(bignum* a, bignum* b) {
1686   GC_BIGNUM(a);
1687   GC_BIGNUM(b);
1688   bignum* d;
1689   bignum_length_type size_a, size_b;
1690   bignum_digit_type* scan_a, *scan_b, *scan_d, *a_end, *b_end;
1691
1692   if (BIGNUM_NEGATIVE_P(a)) {
1693     size_a = BIGNUM_LENGTH(a);
1694     d = allot_bignum(size_a, 0);
1695     scan_d = BIGNUM_START_PTR(d);
1696     scan_a = BIGNUM_START_PTR(a);
1697     a_end = scan_a + size_a;
1698     while (scan_a < a_end)
1699       (*scan_d++) = (*scan_a++);
1700     a = d;
1701   }
1702
1703   if (BIGNUM_NEGATIVE_P(b)) {
1704     size_b = BIGNUM_LENGTH(b);
1705     d = allot_bignum(size_b, 0);
1706     scan_d = BIGNUM_START_PTR(d);
1707     scan_b = BIGNUM_START_PTR(b);
1708     b_end = scan_b + size_b;
1709     while (scan_b < b_end)
1710       (*scan_d++) = (*scan_b++);
1711     b = d;
1712   }
1713
1714   if (bignum_compare(a, b) == bignum_comparison_less) {
1715     std::swap(a, b);
1716   }
1717
1718   while (BIGNUM_LENGTH(b) != 0) {
1719     d = bignum_remainder(a, b);
1720     GC_BIGNUM(d);
1721     if (d == BIGNUM_OUT_OF_BAND) {
1722       return d;
1723     }
1724     a = b;
1725     b = d;
1726   }
1727
1728   return a;
1729 }
1730 #else
1731 /* Allocates memory */
1732 bignum* factor_vm::bignum_gcd(bignum* a, bignum* b) {
1733   GC_BIGNUM(a);
1734   GC_BIGNUM(b);
1735   bignum* c, *d, *e, *f;
1736   bignum_twodigit_type x, y, q, s, t, A, B, C, D;
1737   int nbits, k;
1738   bignum_length_type size_a, size_b, size_c;
1739   bignum_digit_type* scan_a, *scan_b, *scan_c, *scan_d;
1740   bignum_digit_type* a_end, *b_end, *c_end;
1741
1742   /* clone the bignums so we can modify them in-place */
1743   size_a = BIGNUM_LENGTH(a);
1744   c = allot_bignum(size_a, 0);
1745   scan_a = BIGNUM_START_PTR(a);
1746   a_end = scan_a + size_a;
1747   GC_BIGNUM(c);
1748   scan_c = BIGNUM_START_PTR(c);
1749   while (scan_a < a_end)
1750     (*scan_c++) = (*scan_a++);
1751   a = c;
1752   size_b = BIGNUM_LENGTH(b);
1753   d = allot_bignum(size_b, 0);
1754   scan_b = BIGNUM_START_PTR(b);
1755   b_end = scan_b + size_b;
1756   GC_BIGNUM(d);
1757   scan_d = BIGNUM_START_PTR(d);
1758   while (scan_b < b_end)
1759     (*scan_d++) = (*scan_b++);
1760   b = d;
1761
1762   /* Initial reduction: make sure that 0 <= b <= a. */
1763   if (bignum_compare(a, b) == bignum_comparison_less) {
1764     std::swap(a, b);
1765     std::swap(size_a, size_b);
1766   }
1767
1768   while (size_a > 1) {
1769     nbits = log2(BIGNUM_REF(a, size_a - 1));
1770     x = ((BIGNUM_REF(a, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)) |
1771          (BIGNUM_REF(a, size_a - 2) >> nbits));
1772     y = ((size_b >= size_a - 1 ? BIGNUM_REF(b, size_a - 2) >> nbits : 0) |
1773          (size_b >= size_a
1774               ? BIGNUM_REF(b, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)
1775               : 0));
1776
1777     /* inner loop of Lehmer's algorithm; */
1778     A = 1;
1779     B = 0;
1780     C = 0;
1781     D = 1;
1782     for (k = 0;; k++) {
1783       if (y - C == 0)
1784         break;
1785
1786       q = (x + (A - 1)) / (y - C);
1787
1788       s = B + (q * D);
1789       t = x - (q * y);
1790
1791       if (s > t)
1792         break;
1793
1794       x = y;
1795       y = t;
1796
1797       t = A + (q * C);
1798
1799       A = D;
1800       B = C;
1801       C = s;
1802       D = t;
1803     }
1804
1805     if (k == 0) {
1806       /* no progress; do a Euclidean step */
1807       if (size_b == 0) {
1808         return bignum_trim(a);
1809       }
1810       e = bignum_trim(a);
1811       GC_BIGNUM(e);
1812       f = bignum_trim(b);
1813       GC_BIGNUM(f);
1814       c = bignum_remainder(e, f);
1815       GC_BIGNUM(c);
1816       if (c == BIGNUM_OUT_OF_BAND) {
1817         return c;
1818       }
1819
1820       // copy 'b' to 'a'
1821       scan_a = BIGNUM_START_PTR(a);
1822       scan_b = BIGNUM_START_PTR(b);
1823       a_end = scan_a + size_a;
1824       b_end = scan_b + size_b;
1825       while (scan_b < b_end)
1826         *(scan_a++) = *(scan_b++);
1827       while (scan_a < a_end)
1828         *(scan_a++) = 0;
1829       size_a = size_b;
1830
1831       // copy 'c' to 'b'
1832       scan_b = BIGNUM_START_PTR(b);
1833       scan_c = BIGNUM_START_PTR(c);
1834       size_c = BIGNUM_LENGTH(c);
1835       c_end = scan_c + size_c;
1836       while (scan_c < c_end)
1837         *(scan_b++) = *(scan_c++);
1838       while (scan_b < b_end)
1839         *(scan_b++) = 0;
1840       size_b = size_c;
1841
1842       continue;
1843     }
1844
1845     /*
1846       a, b = A*b - B*a, D*a - C*b if k is odd
1847       a, b = A*a - B*b, D*b - C*a if k is even
1848     */
1849     scan_a = BIGNUM_START_PTR(a);
1850     scan_b = BIGNUM_START_PTR(b);
1851     scan_c = scan_a;
1852     scan_d = scan_b;
1853     a_end = scan_a + size_a;
1854     b_end = scan_b + size_b;
1855     s = 0;
1856     t = 0;
1857     if (k & 1) {
1858       while (scan_b < b_end) {
1859         s += (A * *scan_b) - (B * *scan_a);
1860         t += (D * *scan_a++) - (C * *scan_b++);
1861         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1862         *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1863         s >>= BIGNUM_DIGIT_LENGTH;
1864         t >>= BIGNUM_DIGIT_LENGTH;
1865       }
1866       while (scan_a < a_end) {
1867         s -= (B * *scan_a);
1868         t += (D * *scan_a++);
1869         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1870         //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1871         s >>= BIGNUM_DIGIT_LENGTH;
1872         t >>= BIGNUM_DIGIT_LENGTH;
1873       }
1874     } else {
1875       while (scan_b < b_end) {
1876         s += (A * *scan_a) - (B * *scan_b);
1877         t += (D * *scan_b++) - (C * *scan_a++);
1878         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1879         *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1880         s >>= BIGNUM_DIGIT_LENGTH;
1881         t >>= BIGNUM_DIGIT_LENGTH;
1882       }
1883       while (scan_a < a_end) {
1884         s += (A * *scan_a);
1885         t -= (C * *scan_a++);
1886         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1887         //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1888         s >>= BIGNUM_DIGIT_LENGTH;
1889         t >>= BIGNUM_DIGIT_LENGTH;
1890       }
1891     }
1892     BIGNUM_ASSERT(s == 0);
1893     BIGNUM_ASSERT(t == 0);
1894
1895     // update size_a and size_b to remove any zeros at end
1896     while (size_a > 0 && *(--scan_a) == 0)
1897       size_a--;
1898     while (size_b > 0 && *(--scan_b) == 0)
1899       size_b--;
1900
1901     BIGNUM_ASSERT(size_a >= size_b);
1902   }
1903
1904   /* a fits into a fixnum, so b must too */
1905   fixnum xx = bignum_to_fixnum(a);
1906   fixnum yy = bignum_to_fixnum(b);
1907   fixnum tt;
1908
1909   /* usual Euclidean algorithm for longs */
1910   while (yy != 0) {
1911     tt = yy;
1912     yy = xx % yy;
1913     xx = tt;
1914   }
1915
1916   return fixnum_to_bignum(xx);
1917 }
1918 #endif
1919
1920 }