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