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