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