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