vm/bignum.cpp

   1 // Copyright (C) 1989-94 Massachusetts Institute of Technology
   2 // Portions copyright (C) 2004-2008 Slava Pestov
   3
   4 // This material was developed by the Scheme project at the Massachusetts
   5 // Institute of Technology, Department of Electrical Engineering and
   6 // Computer Science.  Permission to copy and modify this software, to
   7 // redistribute either the original software or a modified version, and
   8 // to use this software for any purpose is granted, subject to the
   9 // following restrictions and understandings.
  10
  11 // 1. Any copy made of this software must include this copyright notice
  12 // in full.
  13
  14 // 2. Users of this software agree to make their best efforts (a) to
  15 // return to the MIT Scheme project any improvements or extensions that
  16 // they make, so that these may be included in future releases; and (b)
  17 // to inform MIT of noteworthy uses of this software.
  18
  19 // 3. All materials developed as a consequence of the use of this
  20 // software shall duly acknowledge such use, in accordance with the usual
  21 // standards of acknowledging credit in academic research.
  22
  23 // 4. MIT has made no warrantee or representation that the operation of
  24 // this software will be error-free, and MIT is under no obligation to
  25 // provide any services, by way of maintenance, update, or otherwise.
  26
  27 // 5. In conjunction with products arising from the use of this material,
  28 // there shall be no use of the name of the Massachusetts Institute of
  29 // Technology nor of any adaptation thereof in any advertising,
  30 // promotional, or sales literature without prior written consent from
  31 // MIT in each case.
  32
  33 // Changes for Scheme 48:
  34 // *  - Converted to ANSI.
  35 // *  - Added bitwise operations.
  36 // *  - Added s48 to the beginning of all externally visible names.
  37 // *  - Cached the bignum representations of -1, 0, and 1.
  38
  39 // Changes for Factor:
  40 // *  - Adapt bignumint.h for Factor memory manager
  41 // *  - Add more bignum <-> C type conversions
  42 // *  - Remove unused functions
  43 // *  - Add local variable GC root recording
  44 // *  - Remove s48 prefix from function names
  45 // *  - Various fixes for Win64
  46 // *  - Port to C++
  47 // *  - Added bignum_gcd implementation
  48
  49 #include "master.hpp"
  50
  51 namespace factor {
  52
  53 // Exports
  54
  55 int factor_vm::bignum_equal_p(bignum* x, bignum* y) {
  56   return ((BIGNUM_ZERO_P(x))
  57               ? (BIGNUM_ZERO_P(y))
  58               : ((!(BIGNUM_ZERO_P(y))) &&
  59                  ((BIGNUM_NEGATIVE_P(x)) ? (BIGNUM_NEGATIVE_P(y))
  60                                          : (!(BIGNUM_NEGATIVE_P(y)))) &&
  61                  (bignum_equal_p_unsigned(x, y))));
  62 }
  63
  64 enum bignum_comparison factor_vm::bignum_compare(bignum* x, bignum* y) {
  65   return ((BIGNUM_ZERO_P(x)) ? ((BIGNUM_ZERO_P(y)) ? BIGNUM_COMPARISON_EQUAL
  66                                                    : (BIGNUM_NEGATIVE_P(y))
  67                                     ? BIGNUM_COMPARISON_GREATER
  68                                     : BIGNUM_COMPARISON_LESS)
  69                              : (BIGNUM_ZERO_P(y))
  70               ? ((BIGNUM_NEGATIVE_P(x)) ? BIGNUM_COMPARISON_LESS
  71                                         : BIGNUM_COMPARISON_GREATER)
  72               : (BIGNUM_NEGATIVE_P(x))
  73               ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_compare_unsigned(y, x))
  74                                         : (BIGNUM_COMPARISON_LESS))
  75               : ((BIGNUM_NEGATIVE_P(y)) ? (BIGNUM_COMPARISON_GREATER)
  76                                         : (bignum_compare_unsigned(x, y))));
  77 }
  78
  79 // Allocates memory
  80 bignum* factor_vm::bignum_add(bignum* x, bignum* y) {
  81   return (
  82       (BIGNUM_ZERO_P(x)) ? (y) : (BIGNUM_ZERO_P(y))
  83           ? (x)
  84           : ((BIGNUM_NEGATIVE_P(x))
  85                  ? ((BIGNUM_NEGATIVE_P(y)) ? (bignum_add_unsigned(x, y, 1))
  86                                            : (bignum_subtract_unsigned(y, x)))
  87                  : ((BIGNUM_NEGATIVE_P(y)) ? (bignum_subtract_unsigned(x, y))
  88                                            : (bignum_add_unsigned(x, y, 0)))));
  89 }
  90
  91 // Allocates memory
  92 bignum* factor_vm::bignum_subtract(bignum* x, bignum* y) {
  93   return ((BIGNUM_ZERO_P(x))
  94               ? ((BIGNUM_ZERO_P(y)) ? (y) : (bignum_new_sign(
  95                                                 y, (!(BIGNUM_NEGATIVE_P(y))))))
  96               : ((BIGNUM_ZERO_P(y))
  97                      ? (x)
  98                      : ((BIGNUM_NEGATIVE_P(x))
  99                             ? ((BIGNUM_NEGATIVE_P(y))
 100                                    ? (bignum_subtract_unsigned(y, x))
 101                                    : (bignum_add_unsigned(x, y, 1)))
 102                             : ((BIGNUM_NEGATIVE_P(y))
 103                                    ? (bignum_add_unsigned(x, y, 0))
 104                                    : (bignum_subtract_unsigned(x, y))))));
 105 }
 106
 107 #ifdef _WIN64
 108 bignum *factor_vm::bignum_square(bignum* x_)
 109 {
 110     return bignum_multiply(x_, x_);
 111 }
 112 #else
 113 // Allocates memory
 114 bignum *factor_vm::bignum_square(bignum* x_)
 115 {
 116     data_root<bignum> x(x_, this);
 117
 118     bignum_length_type length = (BIGNUM_LENGTH (x));
 119     bignum * z = (allot_bignum_zeroed ((length + length), 0));
 120
 121     bignum_digit_type * scan_z = BIGNUM_START_PTR (z);
 122     bignum_digit_type * scan_x = BIGNUM_START_PTR (x);
 123     bignum_digit_type * end_x = scan_x + length;
 124
 125     for (int i = 0; i < length; ++i) {
 126         bignum_twodigit_type carry;
 127         bignum_twodigit_type f = BIGNUM_REF (x, i);
 128         bignum_digit_type *pz = scan_z + (i << 1);
 129         bignum_digit_type *px = scan_x + i + 1;
 130
 131         carry = *pz + f * f;
 132         *pz++ = carry & BIGNUM_DIGIT_MASK;
 133         carry >>= BIGNUM_DIGIT_LENGTH;
 134         BIGNUM_ASSERT (carry <= BIGNUM_DIGIT_MASK);
 135
 136         f <<= 1;
 137         while (px < end_x)
 138         {
 139             carry += *pz + *px++ * f;
 140             *pz++ = carry & BIGNUM_DIGIT_MASK;
 141             carry >>= BIGNUM_DIGIT_LENGTH;
 142             BIGNUM_ASSERT (carry <= (BIGNUM_DIGIT_MASK << 1));
 143         }
 144         if (carry) {
 145             carry += *pz;
 146             *pz++ = carry & BIGNUM_DIGIT_MASK;
 147             carry >>= BIGNUM_DIGIT_LENGTH;
 148         }
 149         if (carry)
 150             *pz += carry & BIGNUM_DIGIT_MASK;
 151         BIGNUM_ASSERT ((carry >> BIGNUM_DIGIT_LENGTH) == 0);
 152     }
 153     return (bignum_trim (z));
 154 }
 155 #endif
 156
 157 // Allocates memory
 158 bignum* factor_vm::bignum_multiply(bignum* x, bignum* y) {
 159
 160 #ifndef _WIN64
 161   if (x == y) {
 162     return bignum_square(x);
 163   }
 164 #endif
 165
 166   bignum_length_type x_length = (BIGNUM_LENGTH(x));
 167   bignum_length_type y_length = (BIGNUM_LENGTH(y));
 168   int negative_p = ((BIGNUM_NEGATIVE_P(x)) ? (!(BIGNUM_NEGATIVE_P(y)))
 169                                            : (BIGNUM_NEGATIVE_P(y)));
 170   if (BIGNUM_ZERO_P(x))
 171     return (x);
 172   if (BIGNUM_ZERO_P(y))
 173     return (y);
 174   if (x_length == 1) {
 175     bignum_digit_type digit = (BIGNUM_REF(x, 0));
 176     if (digit == 1)
 177       return (bignum_maybe_new_sign(y, negative_p));
 178     if (digit < BIGNUM_RADIX_ROOT)
 179       return (bignum_multiply_unsigned_small_factor(y, digit, negative_p));
 180   }
 181   if (y_length == 1) {
 182     bignum_digit_type digit = (BIGNUM_REF(y, 0));
 183     if (digit == 1)
 184       return (bignum_maybe_new_sign(x, negative_p));
 185     if (digit < BIGNUM_RADIX_ROOT)
 186       return (bignum_multiply_unsigned_small_factor(x, digit, negative_p));
 187   }
 188   return (bignum_multiply_unsigned(x, y, negative_p));
 189 }
 190
 191 // Allocates memory
 192 void factor_vm::bignum_divide(bignum* numerator, bignum* denominator,
 193                               bignum** quotient, bignum** remainder) {
 194   if (BIGNUM_ZERO_P(denominator)) {
 195     divide_by_zero_error();
 196     (*quotient) = denominator;
 197     (*remainder) = denominator;
 198     return;
 199   }
 200   if (BIGNUM_ZERO_P(numerator)) {
 201     (*quotient) = numerator;
 202     (*remainder) = numerator;
 203   } else {
 204     int r_negative_p = (BIGNUM_NEGATIVE_P(numerator));
 205     int q_negative_p =
 206         ((BIGNUM_NEGATIVE_P(denominator)) ? (!r_negative_p) : r_negative_p);
 207     switch (bignum_compare_unsigned(numerator, denominator)) {
 208       case BIGNUM_COMPARISON_EQUAL: {
 209         (*quotient) = (BIGNUM_ONE(q_negative_p));
 210         (*remainder) = (BIGNUM_ZERO());
 211         break;
 212       }
 213       case BIGNUM_COMPARISON_LESS: {
 214         (*quotient) = (BIGNUM_ZERO());
 215         (*remainder) = numerator;
 216         break;
 217       }
 218       case BIGNUM_COMPARISON_GREATER: {
 219         if ((BIGNUM_LENGTH(denominator)) == 1) {
 220           bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
 221           if (digit == 1) {
 222             (*quotient) = (bignum_maybe_new_sign(numerator, q_negative_p));
 223             (*remainder) = (BIGNUM_ZERO());
 224             break;
 225           } else if (digit < BIGNUM_RADIX_ROOT) {
 226             bignum_divide_unsigned_small_denominator(numerator, digit, quotient,
 227                                                      remainder, q_negative_p,
 228                                                      r_negative_p);
 229             break;
 230           } else {
 231             bignum_divide_unsigned_medium_denominator(
 232                 numerator, digit, quotient, remainder, q_negative_p,
 233                 r_negative_p);
 234             break;
 235           }
 236         }
 237         bignum_divide_unsigned_large_denominator(
 238             numerator, denominator, quotient, remainder, q_negative_p,
 239             r_negative_p);
 240         break;
 241       }
 242     }
 243   }
 244 }
 245
 246 // Allocates memory
 247 bignum* factor_vm::bignum_quotient(bignum* numerator, bignum* denominator) {
 248   if (BIGNUM_ZERO_P(denominator)) {
 249     divide_by_zero_error();
 250     return (BIGNUM_OUT_OF_BAND);
 251   }
 252   if (BIGNUM_ZERO_P(numerator))
 253     return numerator;
 254   {
 255     int q_negative_p =
 256         ((BIGNUM_NEGATIVE_P(denominator)) ? (!(BIGNUM_NEGATIVE_P(numerator)))
 257                                           : (BIGNUM_NEGATIVE_P(numerator)));
 258     switch (bignum_compare_unsigned(numerator, denominator)) {
 259       case BIGNUM_COMPARISON_EQUAL:
 260         return (BIGNUM_ONE(q_negative_p));
 261       case BIGNUM_COMPARISON_LESS:
 262         return (BIGNUM_ZERO());
 263       case BIGNUM_COMPARISON_GREATER:
 264       default: // to appease gcc -Wall
 265                {
 266         bignum* quotient;
 267         if ((BIGNUM_LENGTH(denominator)) == 1) {
 268           bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
 269           if (digit == 1)
 270             return (bignum_maybe_new_sign(numerator, q_negative_p));
 271           if (digit < BIGNUM_RADIX_ROOT)
 272             bignum_divide_unsigned_small_denominator(
 273                 numerator, digit, (&quotient), ((bignum**)0), q_negative_p, 0);
 274           else
 275             bignum_divide_unsigned_medium_denominator(
 276                 numerator, digit, (&quotient), ((bignum**)0), q_negative_p, 0);
 277         } else
 278           bignum_divide_unsigned_large_denominator(
 279               numerator, denominator, (&quotient), ((bignum**)0), q_negative_p,
 280               0);
 281         return (quotient);
 282       }
 283     }
 284   }
 285 }
 286
 287 // Allocates memory
 288 bignum* factor_vm::bignum_remainder(bignum* numerator, bignum* denominator) {
 289   if (BIGNUM_ZERO_P(denominator)) {
 290     divide_by_zero_error();
 291     return (BIGNUM_OUT_OF_BAND);
 292   }
 293   if (BIGNUM_ZERO_P(numerator))
 294     return numerator;
 295   switch (bignum_compare_unsigned(numerator, denominator)) {
 296     case BIGNUM_COMPARISON_EQUAL:
 297       return (BIGNUM_ZERO());
 298     case BIGNUM_COMPARISON_LESS:
 299       return numerator;
 300     case BIGNUM_COMPARISON_GREATER:
 301     default: // to appease gcc -Wall
 302              {
 303       bignum* remainder;
 304       if ((BIGNUM_LENGTH(denominator)) == 1) {
 305         bignum_digit_type digit = (BIGNUM_REF(denominator, 0));
 306         if (digit == 1)
 307           return (BIGNUM_ZERO());
 308         if (digit < BIGNUM_RADIX_ROOT)
 309           return (bignum_remainder_unsigned_small_denominator(
 310               numerator, digit, (BIGNUM_NEGATIVE_P(numerator))));
 311         bignum_divide_unsigned_medium_denominator(
 312             numerator, digit, ((bignum**)0), (&remainder), 0,
 313             (BIGNUM_NEGATIVE_P(numerator)));
 314       } else
 315         bignum_divide_unsigned_large_denominator(
 316             numerator, denominator, ((bignum**)0), (&remainder), 0,
 317             (BIGNUM_NEGATIVE_P(numerator)));
 318       return (remainder);
 319     }
 320   }
 321 }
 322
 323 // cell_to_bignum, fixnum_to_bignum, long_long_to_bignum, ulong_long_to_bignum
 324
 325 // Allocates memory
 326 #define FOO_TO_BIGNUM_SIGNED(name, type, utype)                       \
 327   bignum* factor_vm::name##_to_bignum(type n) {                       \
 328     int negative_p;                                                   \
 329     /* Special cases win when these small constants are cached. */    \
 330     if (n == 0)                                                       \
 331       return (BIGNUM_ZERO());                                         \
 332     if (n == 1)                                                       \
 333       return (BIGNUM_ONE(0));                                         \
 334     if (n < (type) 0 && n == (type) -1)                               \
 335       return (BIGNUM_ONE(1));                                         \
 336     {                                                                 \
 337       utype accumulator =                                             \
 338           ((negative_p = n < (type) 0) ? -n : n);                     \
 339       if (accumulator < BIGNUM_RADIX)                                 \
 340       {                                                               \
 341         bignum* result = allot_bignum(1, negative_p);                 \
 342         bignum_digit_type* scan = (BIGNUM_START_PTR(result));         \
 343         *scan = (accumulator & BIGNUM_DIGIT_MASK);                    \
 344         return (result);                                              \
 345       } else {                                                        \
 346         bignum_digit_type result_digits[BIGNUM_DIGITS_FOR(type)];     \
 347         bignum_digit_type* end_digits = result_digits;                \
 348         do {                                                          \
 349           (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);        \
 350           accumulator >>= BIGNUM_DIGIT_LENGTH;                        \
 351         } while (accumulator != 0);                                   \
 352         bignum* result =                                              \
 353            (allot_bignum((end_digits - result_digits), negative_p));  \
 354         bignum_digit_type* scan_digits = result_digits;               \
 355         bignum_digit_type* scan_result = (BIGNUM_START_PTR(result));  \
 356         while (scan_digits < end_digits)                              \
 357           (*scan_result++) = (*scan_digits++);                        \
 358         return (result);                                              \
 359       }                                                               \
 360     }                                                                 \
 361   }
 362
 363 // Allocates memory
 364 #define FOO_TO_BIGNUM_UNSIGNED(name, type, utype)                     \
 365   bignum* factor_vm::name##_to_bignum(type n) {                       \
 366     /* Special cases win when these small constants are cached. */    \
 367     if (n == 0)                                                       \
 368       return (BIGNUM_ZERO());                                         \
 369     if (n == 1)                                                       \
 370       return (BIGNUM_ONE(0));                                         \
 371     {                                                                 \
 372       utype accumulator = n;                                          \
 373       if (accumulator < BIGNUM_RADIX)                                 \
 374       {                                                               \
 375         bignum* result = allot_bignum(1, false);                      \
 376         bignum_digit_type* scan = (BIGNUM_START_PTR(result));         \
 377         *scan = (accumulator & BIGNUM_DIGIT_MASK);                    \
 378         return (result);                                              \
 379       } else {                                                        \
 380         bignum_digit_type result_digits[BIGNUM_DIGITS_FOR(type)];     \
 381         bignum_digit_type* end_digits = result_digits;                \
 382         do {                                                          \
 383           (*end_digits++) = (accumulator & BIGNUM_DIGIT_MASK);        \
 384           accumulator >>= BIGNUM_DIGIT_LENGTH;                        \
 385         } while (accumulator != 0);                                   \
 386         bignum* result =                                              \
 387            (allot_bignum((end_digits - result_digits), false));       \
 388         bignum_digit_type* scan_digits = result_digits;               \
 389         bignum_digit_type* scan_result = (BIGNUM_START_PTR(result));  \
 390         while (scan_digits < end_digits)                              \
 391           (*scan_result++) = (*scan_digits++);                        \
 392         return (result);                                              \
 393       }                                                               \
 394     }                                                                 \
 395   }
 396
 397 FOO_TO_BIGNUM_SIGNED(fixnum, fixnum, cell)
 398 FOO_TO_BIGNUM_UNSIGNED(cell, cell, cell)
 399 FOO_TO_BIGNUM_SIGNED(long_long, int64_t, uint64_t)
 400 FOO_TO_BIGNUM_UNSIGNED(ulong_long, uint64_t, uint64_t)
 401
 402 // cannot allocate memory
 403 // bignum_to_cell, fixnum_to_cell, long_long_to_cell, ulong_long_to_cell
 404 #define BIGNUM_TO_FOO(name, type, stype, utype)                            \
 405   type bignum_to_##name(bignum* bn) {                                      \
 406     if (BIGNUM_ZERO_P(bn))                                                 \
 407       return (0);                                                          \
 408     {                                                                      \
 409       utype accumulator = 0;                                               \
 410       bignum_digit_type* start = (BIGNUM_START_PTR(bn));                   \
 411       bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));             \
 412       while (start < scan)                                                 \
 413         accumulator = ((accumulator << BIGNUM_DIGIT_LENGTH) + (*--scan));  \
 414       return ((BIGNUM_NEGATIVE_P(bn)) ? ((type)(-(stype) accumulator))     \
 415                                       : accumulator);                      \
 416     }                                                                      \
 417   }
 418
 419 BIGNUM_TO_FOO(cell, cell, fixnum, cell)
 420 BIGNUM_TO_FOO(fixnum, fixnum, fixnum, cell)
 421 BIGNUM_TO_FOO(long_long, int64_t, int64_t, uint64_t)
 422 BIGNUM_TO_FOO(ulong_long, uint64_t, int64_t, uint64_t)
 423
 424 bool bignum_fits_fixnum_p(bignum* bn) {
 425   fixnum len = BIGNUM_LENGTH(bn);
 426   if (len == 0)
 427     return true;
 428   if (len > 1)
 429     return false;
 430   bignum_digit_type dig = BIGNUM_START_PTR(bn)[0];
 431   return (BIGNUM_NEGATIVE_P(bn) && dig <= -fixnum_min) ||
 432       (!BIGNUM_NEGATIVE_P(bn) && dig <= fixnum_max);
 433 }
 434
 435 cell bignum_maybe_to_fixnum(bignum* bn) {
 436   if (bignum_fits_fixnum_p(bn))
 437     return tag_fixnum(bignum_to_fixnum(bn));
 438   return tag<bignum>(bn);
 439 }
 440
 441 // cannot allocate memory
 442 fixnum factor_vm::bignum_to_fixnum_strict(bignum* bn) {
 443
 444   if (!bignum_fits_fixnum_p(bn)) {
 445      general_error(ERROR_OUT_OF_FIXNUM_RANGE, tag<bignum>(bn), false_object);
 446   }
 447   fixnum fix = bignum_to_fixnum(bn);
 448   FACTOR_ASSERT(fix <= fixnum_max && fix >= fixnum_min);
 449   return fix;
 450 }
 451
 452 #define DTB_WRITE_DIGIT(factor)                \
 453   {                                            \
 454     significand *= (factor);                   \
 455     digit = ((bignum_digit_type) significand); \
 456     (*--scan) = digit;                         \
 457     significand -= ((double)digit);            \
 458   }
 459
 460 #define inf std::numeric_limits<double>::infinity()
 461
 462 // Allocates memory
 463 bignum* factor_vm::double_to_bignum(double x) {
 464   if (x == inf || x == -inf || x != x)
 465     return (BIGNUM_ZERO());
 466   int exponent;
 467   double significand = (frexp(x, (&exponent)));
 468   if (exponent <= 0)
 469     return (BIGNUM_ZERO());
 470   if (exponent == 1)
 471     return (BIGNUM_ONE(x < 0));
 472   if (significand < 0)
 473     significand = (-significand);
 474   {
 475     bignum_length_type length = (BIGNUM_BITS_TO_DIGITS(exponent));
 476     bignum* result = (allot_bignum(length, (x < 0)));
 477     bignum_digit_type* start = (BIGNUM_START_PTR(result));
 478     bignum_digit_type* scan = (start + length);
 479     bignum_digit_type digit;
 480     int odd_bits = (exponent % BIGNUM_DIGIT_LENGTH);
 481     if (odd_bits > 0)
 482       DTB_WRITE_DIGIT((fixnum)1 << odd_bits);
 483     while (start < scan) {
 484       if (significand == 0) {
 485         while (start < scan)
 486           (*--scan) = 0;
 487         break;
 488       }
 489       DTB_WRITE_DIGIT(BIGNUM_RADIX);
 490     }
 491     return (result);
 492   }
 493 }
 494
 495 #undef DTB_WRITE_DIGIT
 496
 497 // Comparisons
 498
 499 int factor_vm::bignum_equal_p_unsigned(bignum* x, bignum* y) {
 500   bignum_length_type length = (BIGNUM_LENGTH(x));
 501   if (length != (BIGNUM_LENGTH(y)))
 502     return (0);
 503   else {
 504     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 505     bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 506     bignum_digit_type* end_x = (scan_x + length);
 507     while (scan_x < end_x)
 508       if ((*scan_x++) != (*scan_y++))
 509         return (0);
 510     return (1);
 511   }
 512 }
 513
 514 enum bignum_comparison factor_vm::bignum_compare_unsigned(bignum* x,
 515                                                           bignum* y) {
 516   bignum_length_type x_length = (BIGNUM_LENGTH(x));
 517   bignum_length_type y_length = (BIGNUM_LENGTH(y));
 518   if (x_length < y_length)
 519     return BIGNUM_COMPARISON_LESS;
 520   if (x_length > y_length)
 521     return BIGNUM_COMPARISON_GREATER;
 522   {
 523     bignum_digit_type* start_x = (BIGNUM_START_PTR(x));
 524     bignum_digit_type* scan_x = (start_x + x_length);
 525     bignum_digit_type* scan_y = ((BIGNUM_START_PTR(y)) + y_length);
 526     while (start_x < scan_x) {
 527       bignum_digit_type digit_x = (*--scan_x);
 528       bignum_digit_type digit_y = (*--scan_y);
 529       if (digit_x < digit_y)
 530         return BIGNUM_COMPARISON_LESS;
 531       if (digit_x > digit_y)
 532         return BIGNUM_COMPARISON_GREATER;
 533     }
 534   }
 535   return BIGNUM_COMPARISON_EQUAL;
 536 }
 537
 538 // Addition
 539
 540 // Allocates memory
 541 bignum* factor_vm::bignum_add_unsigned(bignum* x_, bignum* y_, int negative_p) {
 542
 543   data_root<bignum> x(x_, this);
 544   data_root<bignum> y(y_, this);
 545   if ((BIGNUM_LENGTH(y)) > (BIGNUM_LENGTH(x))) {
 546     swap(x, y);
 547   }
 548   {
 549     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 550
 551     bignum* r = (allot_bignum((x_length + 1), negative_p));
 552
 553     bignum_digit_type sum;
 554     bignum_digit_type carry = 0;
 555     bignum_digit_type* scan_x = (BIGNUM_START_PTR(x));
 556     bignum_digit_type* scan_r = (BIGNUM_START_PTR(r));
 557     {
 558       bignum_digit_type* scan_y = (BIGNUM_START_PTR(y));
 559       bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
 560       while (scan_y < end_y) {
 561         sum = ((*scan_x++) + (*scan_y++) + carry);
 562         if (sum < BIGNUM_RADIX) {
 563           (*scan_r++) = sum;
 564           carry = 0;
 565         } else {
 566           (*scan_r++) = (sum - BIGNUM_RADIX);
 567           carry = 1;
 568         }
 569       }
 570     }
 571     {
 572       bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length;
 573       if (carry != 0)
 574         while (scan_x < end_x) {
 575           sum = ((*scan_x++) + 1);
 576           if (sum < BIGNUM_RADIX) {
 577             (*scan_r++) = sum;
 578             carry = 0;
 579             break;
 580           } else
 581             (*scan_r++) = (sum - BIGNUM_RADIX);
 582         }
 583       while (scan_x < end_x)
 584         (*scan_r++) = (*scan_x++);
 585     }
 586     if (carry != 0) {
 587       (*scan_r) = 1;
 588       return (r);
 589     }
 590     return (bignum_shorten_length(r, x_length));
 591   }
 592 }
 593
 594 // Subtraction
 595
 596 // Allocates memory
 597 bignum* factor_vm::bignum_subtract_unsigned(bignum* x_, bignum* y_) {
 598
 599   data_root<bignum> x(x_, this);
 600   data_root<bignum> y(y_, this);
 601
 602   int negative_p = 0;
 603   switch (bignum_compare_unsigned(x.untagged(), y.untagged())) {
 604     case BIGNUM_COMPARISON_EQUAL:
 605       return (BIGNUM_ZERO());
 606     case BIGNUM_COMPARISON_LESS:
 607       swap(x, y);
 608       negative_p = 1;
 609       break;
 610     case BIGNUM_COMPARISON_GREATER:
 611       negative_p = 0;
 612       break;
 613   }
 614   {
 615     bignum_length_type x_length = (BIGNUM_LENGTH(x));
 616
 617     bignum* r = (allot_bignum(x_length, negative_p));
 618
 619     bignum_digit_type difference;
 620     bignum_digit_type borrow = 0;
 621     bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
 622     bignum_digit_type* scan_r = BIGNUM_START_PTR(r);
 623     {
 624       bignum_digit_type* scan_y = BIGNUM_START_PTR(y);
 625       bignum_digit_type* end_y = (scan_y + (BIGNUM_LENGTH(y)));
 626       while (scan_y < end_y) {
 627         difference = (((*scan_x++) - (*scan_y++)) - borrow);
 628         if (difference < 0) {
 629           (*scan_r++) = (difference + BIGNUM_RADIX);
 630           borrow = 1;
 631         } else {
 632           (*scan_r++) = difference;
 633           borrow = 0;
 634         }
 635       }
 636     }
 637     {
 638       bignum_digit_type* end_x = BIGNUM_START_PTR(x) + x_length;
 639       if (borrow != 0)
 640         while (scan_x < end_x) {
 641           difference = ((*scan_x++) - borrow);
 642           if (difference < 0)
 643             (*scan_r++) = (difference + BIGNUM_RADIX);
 644           else {
 645             (*scan_r++) = difference;
 646             borrow = 0;
 647             break;
 648           }
 649         }
 650       BIGNUM_ASSERT(borrow == 0);
 651       while (scan_x < end_x)
 652         (*scan_r++) = (*scan_x++);
 653     }
 654     return (bignum_trim(r));
 655   }
 656 }
 657
 658 // Multiplication
 659 // Maximum value for product_low or product_high:
 660 // ((R * R) + (R * (R - 2)) + (R - 1))
 661 // Maximum value for carry: ((R * (R - 1)) + (R - 1))
 662 // where R == BIGNUM_RADIX_ROOT
 663
 664 // Allocates memory
 665 bignum* factor_vm::bignum_multiply_unsigned(bignum* x_, bignum* y_,
 666                                             int negative_p) {
 667
 668   data_root<bignum> x(x_, this);
 669   data_root<bignum> y(y_, this);
 670
 671   if (BIGNUM_LENGTH(y) > BIGNUM_LENGTH(x)) {
 672     swap(x, y);
 673   }
 674   {
 675     bignum_digit_type carry;
 676     bignum_digit_type y_digit_low;
 677     bignum_digit_type y_digit_high;
 678     bignum_digit_type x_digit_low;
 679     bignum_digit_type x_digit_high;
 680     bignum_digit_type product_low;
 681     bignum_digit_type* scan_r;
 682     bignum_digit_type* scan_y;
 683     bignum_length_type x_length = BIGNUM_LENGTH(x);
 684     bignum_length_type y_length = BIGNUM_LENGTH(y);
 685
 686     bignum* r = (allot_bignum_zeroed((x_length + y_length), negative_p));
 687
 688     bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
 689     bignum_digit_type* end_x = (scan_x + x_length);
 690     bignum_digit_type* start_y = BIGNUM_START_PTR(y);
 691     bignum_digit_type* end_y = (start_y + y_length);
 692     bignum_digit_type* start_r = (BIGNUM_START_PTR(r));
 693 #define x_digit x_digit_high
 694 #define y_digit y_digit_high
 695 #define product_high carry
 696     while (scan_x < end_x) {
 697       x_digit = (*scan_x++);
 698       x_digit_low = (HD_LOW(x_digit));
 699       x_digit_high = (HD_HIGH(x_digit));
 700       carry = 0;
 701       scan_y = start_y;
 702       scan_r = (start_r++);
 703       while (scan_y < end_y) {
 704         y_digit = (*scan_y++);
 705         y_digit_low = (HD_LOW(y_digit));
 706         y_digit_high = (HD_HIGH(y_digit));
 707         product_low =
 708             ((*scan_r) + (x_digit_low * y_digit_low) + (HD_LOW(carry)));
 709         product_high =
 710             ((x_digit_high * y_digit_low) + (x_digit_low * y_digit_high) +
 711              (HD_HIGH(product_low)) + (HD_HIGH(carry)));
 712         (*scan_r++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
 713         carry = ((x_digit_high * y_digit_high) + (HD_HIGH(product_high)));
 714       }
 715       (*scan_r) += carry;
 716     }
 717     return (bignum_trim(r));
 718 #undef x_digit
 719 #undef y_digit
 720 #undef product_high
 721   }
 722 }
 723
 724 // Allocates memory
 725 bignum* factor_vm::bignum_multiply_unsigned_small_factor(bignum* x_,
 726                                                          bignum_digit_type y,
 727                                                          int negative_p) {
 728   data_root<bignum> x(x_, this);
 729
 730   bignum_length_type length_x = (BIGNUM_LENGTH(x));
 731
 732   bignum* p = (allot_bignum((length_x + 1), negative_p));
 733
 734   bignum_destructive_copy(x.untagged(), p);
 735   (BIGNUM_REF(p, length_x)) = 0;
 736   bignum_destructive_scale_up(p, y);
 737   return (bignum_trim(p));
 738 }
 739
 740 void factor_vm::bignum_destructive_add(bignum* bn, bignum_digit_type n) {
 741   bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
 742   bignum_digit_type digit;
 743   digit = ((*scan) + n);
 744   if (digit < BIGNUM_RADIX) {
 745     (*scan) = digit;
 746     return;
 747   }
 748   (*scan++) = (digit - BIGNUM_RADIX);
 749   while (1) {
 750     digit = ((*scan) + 1);
 751     if (digit < BIGNUM_RADIX) {
 752       (*scan) = digit;
 753       return;
 754     }
 755     (*scan++) = (digit - BIGNUM_RADIX);
 756   }
 757 }
 758
 759 void factor_vm::bignum_destructive_scale_up(bignum* bn,
 760                                             bignum_digit_type factor) {
 761   bignum_digit_type carry = 0;
 762   bignum_digit_type* scan = (BIGNUM_START_PTR(bn));
 763   bignum_digit_type two_digits;
 764   bignum_digit_type product_low;
 765 #define product_high carry
 766   bignum_digit_type* end = (scan + (BIGNUM_LENGTH(bn)));
 767   BIGNUM_ASSERT((factor > 1) && (factor < BIGNUM_RADIX_ROOT));
 768   while (scan < end) {
 769     two_digits = (*scan);
 770     product_low = ((factor * (HD_LOW(two_digits))) + (HD_LOW(carry)));
 771     product_high = ((factor * (HD_HIGH(two_digits))) + (HD_HIGH(product_low)) +
 772                     (HD_HIGH(carry)));
 773     (*scan++) = (HD_CONS((HD_LOW(product_high)), (HD_LOW(product_low))));
 774     carry = (HD_HIGH(product_high));
 775   }
 776   // A carry here would be an overflow, i.e. it would not fit.
 777   // Hopefully the callers allocate enough space that this will
 778   // never happen.
 779   BIGNUM_ASSERT(carry == 0);
 780   return;
 781 #undef product_high
 782 }
 783
 784 // Division
 785
 786 // For help understanding this algorithm, see:
 787 // Knuth, Donald E., "The Art of Computer Programming",
 788 // volume 2, "Seminumerical Algorithms"
 789 // section 4.3.1, "Multiple-Precision Arithmetic".
 790
 791 // Allocates memory
 792 void factor_vm::bignum_divide_unsigned_large_denominator(
 793     bignum* numerator_, bignum* denominator_,
 794     bignum** quotient, bignum** remainder,
 795     int q_negative_p, int r_negative_p) {
 796
 797   data_root<bignum> numerator(numerator_, this);
 798   data_root<bignum> denominator(denominator_, this);
 799
 800   bignum_length_type length_n = BIGNUM_LENGTH(numerator) + 1;
 801   bignum_length_type length_d = BIGNUM_LENGTH(denominator);
 802
 803   data_root<bignum> u(allot_bignum(length_n, r_negative_p), this);
 804
 805   int shift = 0;
 806   BIGNUM_ASSERT(length_d > 1);
 807   {
 808     bignum_digit_type v1 = BIGNUM_REF(denominator.untagged(), length_d - 1);
 809     while (v1 < (BIGNUM_RADIX / 2)) {
 810       v1 <<= 1;
 811       shift += 1;
 812     }
 813   }
 814
 815   if (quotient != NULL) {
 816     bignum *q_ = allot_bignum(length_n - length_d, q_negative_p);
 817     data_root<bignum> q(q_, this);
 818
 819     if (shift == 0) {
 820       bignum_destructive_copy(numerator.untagged(), u.untagged());
 821       (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0;
 822       bignum_divide_unsigned_normalized(u.untagged(),
 823                                         denominator.untagged(),
 824                                         q.untagged());
 825     } else {
 826       bignum* v = allot_bignum(length_d, 0);
 827       bignum_destructive_normalization(numerator.untagged(),
 828                                        u.untagged(),
 829                                        shift);
 830       bignum_destructive_normalization(denominator.untagged(), v, shift);
 831       bignum_divide_unsigned_normalized(u.untagged(), v, q.untagged());
 832       if (remainder != NULL)
 833         bignum_destructive_unnormalization(u.untagged(), shift);
 834     }
 835
 836     q.set_untagged(bignum_trim(q.untagged()));
 837     *quotient = q.untagged();
 838   } else {
 839
 840     if (shift == 0) {
 841         bignum_destructive_copy(numerator.untagged(), u.untagged());
 842         (BIGNUM_REF(u.untagged(), (length_n - 1))) = 0;
 843         bignum_divide_unsigned_normalized(u.untagged(),
 844                                           denominator.untagged(),
 845                                           NULL);
 846       } else {
 847         bignum* v = allot_bignum(length_d, 0);
 848         bignum_destructive_normalization(numerator.untagged(),
 849                                          u.untagged(),
 850                                          shift);
 851         bignum_destructive_normalization(denominator.untagged(),
 852                                          v,
 853                                          shift);
 854         bignum_divide_unsigned_normalized(u.untagged(), v, NULL);
 855         if (remainder != NULL)
 856           bignum_destructive_unnormalization(u.untagged(), shift);
 857       }
 858   }
 859
 860   u.set_untagged(bignum_trim(u.untagged()));
 861   if (remainder != NULL)
 862     *remainder = u.untagged();
 863 }
 864
 865 void factor_vm::bignum_divide_unsigned_normalized(bignum* u, bignum* v,
 866                                                   bignum* q) {
 867   bignum_length_type u_length = (BIGNUM_LENGTH(u));
 868   bignum_length_type v_length = (BIGNUM_LENGTH(v));
 869   bignum_digit_type* u_start = (BIGNUM_START_PTR(u));
 870   bignum_digit_type* u_scan = (u_start + u_length);
 871   bignum_digit_type* u_scan_limit = (u_start + v_length);
 872   bignum_digit_type* u_scan_start = (u_scan - v_length);
 873   bignum_digit_type* v_start = (BIGNUM_START_PTR(v));
 874   bignum_digit_type* v_end = (v_start + v_length);
 875   bignum_digit_type* q_scan = NULL;
 876   bignum_digit_type v1 = (v_end[-1]);
 877   bignum_digit_type v2 = (v_end[-2]);
 878   bignum_digit_type ph; // high half of double-digit product
 879   bignum_digit_type pl; // low half of double-digit product
 880   bignum_digit_type guess;
 881   bignum_digit_type gh; // high half-digit of guess
 882   bignum_digit_type ch; // high half of double-digit comparand
 883   bignum_digit_type v2l = (HD_LOW(v2));
 884   bignum_digit_type v2h = (HD_HIGH(v2));
 885   bignum_digit_type cl; // low half of double-digit comparand
 886 #define gl ph           // low half-digit of guess
 887 #define uj pl
 888 #define qj ph
 889   bignum_digit_type gm; // memory loc for reference parameter
 890   if (q != BIGNUM_OUT_OF_BAND)
 891     q_scan = ((BIGNUM_START_PTR(q)) + (BIGNUM_LENGTH(q)));
 892   while (u_scan_limit < u_scan) {
 893     uj = (*--u_scan);
 894     if (uj != v1) {
 895       // comparand =
 896       // (((((uj * BIGNUM_RADIX) + uj1) % v1) * BIGNUM_RADIX) + uj2);
 897       // guess = (((uj * BIGNUM_RADIX) + uj1) / v1);
 898       cl = (u_scan[-2]);
 899       ch = (bignum_digit_divide(uj, (u_scan[-1]), v1, (&gm)));
 900       guess = gm;
 901     } else {
 902       cl = (u_scan[-2]);
 903       ch = ((u_scan[-1]) + v1);
 904       guess = (BIGNUM_RADIX - 1);
 905     }
 906     while (1) {
 907       // product = (guess * v2);
 908       gl = (HD_LOW(guess));
 909       gh = (HD_HIGH(guess));
 910       pl = (v2l * gl);
 911       ph = ((v2l * gh) + (v2h * gl) + (HD_HIGH(pl)));
 912       pl = (HD_CONS((HD_LOW(ph)), (HD_LOW(pl))));
 913       ph = ((v2h * gh) + (HD_HIGH(ph)));
 914       // if (comparand >= product)
 915       if ((ch > ph) || ((ch == ph) && (cl >= pl)))
 916         break;
 917       guess -= 1;
 918       // comparand += (v1 << BIGNUM_DIGIT_LENGTH)
 919       ch += v1;
 920       // if (comparand >= (BIGNUM_RADIX * BIGNUM_RADIX))
 921       if (ch >= BIGNUM_RADIX)
 922         break;
 923     }
 924     qj = (bignum_divide_subtract(v_start, v_end, guess, (--u_scan_start)));
 925     if (q != BIGNUM_OUT_OF_BAND)
 926       (*--q_scan) = qj;
 927   }
 928   return;
 929 #undef gl
 930 #undef uj
 931 #undef qj
 932 }
 933
 934 bignum_digit_type factor_vm::bignum_divide_subtract(
 935     bignum_digit_type* v_start, bignum_digit_type* v_end,
 936     bignum_digit_type guess, bignum_digit_type* u_start) {
 937   bignum_digit_type* v_scan = v_start;
 938   bignum_digit_type* u_scan = u_start;
 939   bignum_digit_type carry = 0;
 940   if (guess == 0)
 941     return (0);
 942   {
 943     bignum_digit_type gl = (HD_LOW(guess));
 944     bignum_digit_type gh = (HD_HIGH(guess));
 945     bignum_digit_type v;
 946     bignum_digit_type pl;
 947     bignum_digit_type vl;
 948 #define vh v
 949 #define ph carry
 950 #define diff pl
 951     while (v_scan < v_end) {
 952       v = (*v_scan++);
 953       vl = (HD_LOW(v));
 954       vh = (HD_HIGH(v));
 955       pl = ((vl * gl) + (HD_LOW(carry)));
 956       ph = ((vl * gh) + (vh * gl) + (HD_HIGH(pl)) + (HD_HIGH(carry)));
 957       diff = ((*u_scan) - (HD_CONS((HD_LOW(ph)), (HD_LOW(pl)))));
 958       if (diff < 0) {
 959         (*u_scan++) = (diff + BIGNUM_RADIX);
 960         carry = ((vh * gh) + (HD_HIGH(ph)) + 1);
 961       } else {
 962         (*u_scan++) = diff;
 963         carry = ((vh * gh) + (HD_HIGH(ph)));
 964       }
 965     }
 966     if (carry == 0)
 967       return (guess);
 968     diff = ((*u_scan) - carry);
 969     if (diff < 0)
 970       (*u_scan) = (diff + BIGNUM_RADIX);
 971     else {
 972       (*u_scan) = diff;
 973       return (guess);
 974     }
 975 #undef vh
 976 #undef ph
 977 #undef diff
 978   }
 979   // Subtraction generated carry, implying guess is one too large.
 980   // Add v back in to bring it back down.
 981   v_scan = v_start;
 982   u_scan = u_start;
 983   carry = 0;
 984   while (v_scan < v_end) {
 985     bignum_digit_type sum = ((*v_scan++) + (*u_scan) + carry);
 986     if (sum < BIGNUM_RADIX) {
 987       (*u_scan++) = sum;
 988       carry = 0;
 989     } else {
 990       (*u_scan++) = (sum - BIGNUM_RADIX);
 991       carry = 1;
 992     }
 993   }
 994   if (carry == 1) {
 995     bignum_digit_type sum = ((*u_scan) + carry);
 996     (*u_scan) = ((sum < BIGNUM_RADIX) ? sum : (sum - BIGNUM_RADIX));
 997   }
 998   return (guess - 1);
 999 }
1000
1001 // Allocates memory
1002 void factor_vm::bignum_divide_unsigned_medium_denominator(
1003     bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
1004     bignum** remainder, int q_negative_p, int r_negative_p) {
1005
1006   data_root<bignum> numerator(numerator_, this);
1007
1008   bignum_length_type length_n = (BIGNUM_LENGTH(numerator));
1009
1010   int shift = 0;
1011   // Because `bignum_digit_divide' requires a normalized denominator.
1012   while (denominator < (BIGNUM_RADIX / 2)) {
1013     denominator <<= 1;
1014     shift += 1;
1015   }
1016
1017   bignum_length_type length_q = (shift == 0) ? length_n : length_n + 1;
1018   data_root<bignum> q(allot_bignum(length_q, q_negative_p), this);
1019   if (shift == 0) {
1020     bignum_destructive_copy(numerator.untagged(), q.untagged());
1021   } else {
1022     bignum_destructive_normalization(numerator.untagged(), q.untagged(), shift);
1023   }
1024   {
1025     bignum_digit_type r = 0;
1026     bignum_digit_type* start = (BIGNUM_START_PTR(q));
1027     bignum_digit_type* scan = (start + length_q);
1028     bignum_digit_type qj;
1029
1030     while (start < scan) {
1031       r = (bignum_digit_divide(r, (*--scan), denominator, (&qj)));
1032       (*scan) = qj;
1033     }
1034
1035     q.set_untagged(bignum_trim(q.untagged()));
1036
1037     if (remainder != ((bignum**)0)) {
1038       if (shift != 0)
1039         r >>= shift;
1040
1041       (*remainder) = (bignum_digit_to_bignum(r, r_negative_p));
1042     }
1043
1044     if (quotient != ((bignum**)0))
1045       (*quotient) = q.untagged();
1046   }
1047   return;
1048 }
1049
1050 void factor_vm::bignum_destructive_normalization(bignum* source, bignum* target,
1051                                                  int shift_left) {
1052   bignum_digit_type digit;
1053   bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
1054   bignum_digit_type carry = 0;
1055   bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1056   bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1057   bignum_digit_type* end_target = (scan_target + (BIGNUM_LENGTH(target)));
1058   int shift_right = (BIGNUM_DIGIT_LENGTH - shift_left);
1059   bignum_digit_type mask = (((cell)1 << shift_right) - 1);
1060   while (scan_source < end_source) {
1061     digit = (*scan_source++);
1062     (*scan_target++) = (((digit & mask) << shift_left) | carry);
1063     carry = (digit >> shift_right);
1064   }
1065   if (scan_target < end_target)
1066     (*scan_target) = carry;
1067   else
1068     BIGNUM_ASSERT(carry == 0);
1069   return;
1070 }
1071
1072 void factor_vm::bignum_destructive_unnormalization(bignum* bn,
1073                                                    int shift_right) {
1074   bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1075   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
1076   bignum_digit_type digit;
1077   bignum_digit_type carry = 0;
1078   int shift_left = (BIGNUM_DIGIT_LENGTH - shift_right);
1079   bignum_digit_type mask = (((fixnum)1 << shift_right) - 1);
1080   while (start < scan) {
1081     digit = (*--scan);
1082     (*scan) = ((digit >> shift_right) | carry);
1083     carry = ((digit & mask) << shift_left);
1084   }
1085   BIGNUM_ASSERT(carry == 0);
1086   return;
1087 }
1088
1089 // This is a reduced version of the division algorithm, applied to the
1090 // case of dividing two bignum digits by one bignum digit. It is
1091 // assumed that the numerator, denominator are normalized.
1092
1093 #define BDD_STEP(qn, j)                                          \
1094   {                                                              \
1095     uj = (u[j]);                                                 \
1096     if (uj != v1) {                                              \
1097       uj_uj1 = (HD_CONS(uj, (u[j + 1])));                        \
1098       guess = (uj_uj1 / v1);                                     \
1099       comparand = (HD_CONS((uj_uj1 % v1), (u[j + 2])));          \
1100     } else {                                                     \
1101       guess = (BIGNUM_RADIX_ROOT - 1);                           \
1102       comparand = (HD_CONS(((u[j + 1]) + v1), (u[j + 2])));      \
1103     }                                                            \
1104     while ((guess * v2) > comparand) {                           \
1105       guess -= 1;                                                \
1106       comparand += (v1 << BIGNUM_HALF_DIGIT_LENGTH);             \
1107       if (comparand >= BIGNUM_RADIX)                             \
1108         break;                                                   \
1109     }                                                            \
1110     qn = (bignum_digit_divide_subtract(v1, v2, guess, (&u[j]))); \
1111   }
1112
1113 bignum_digit_type factor_vm::bignum_digit_divide(
1114     bignum_digit_type uh, bignum_digit_type ul, bignum_digit_type v,
1115     bignum_digit_type* q) // return value
1116     {
1117   bignum_digit_type guess;
1118   bignum_digit_type comparand;
1119   bignum_digit_type v1 = (HD_HIGH(v));
1120   bignum_digit_type v2 = (HD_LOW(v));
1121   bignum_digit_type uj;
1122   bignum_digit_type uj_uj1;
1123   bignum_digit_type q1;
1124   bignum_digit_type q2;
1125   bignum_digit_type u[4];
1126   if (uh == 0) {
1127     if (ul < v) {
1128       (*q) = 0;
1129       return (ul);
1130     } else if (ul == v) {
1131       (*q) = 1;
1132       return (0);
1133     }
1134   }
1135   (u[0]) = (HD_HIGH(uh));
1136   (u[1]) = (HD_LOW(uh));
1137   (u[2]) = (HD_HIGH(ul));
1138   (u[3]) = (HD_LOW(ul));
1139   v1 = (HD_HIGH(v));
1140   v2 = (HD_LOW(v));
1141   BDD_STEP(q1, 0);
1142   BDD_STEP(q2, 1);
1143   (*q) = (HD_CONS(q1, q2));
1144   return (HD_CONS((u[2]), (u[3])));
1145 }
1146
1147 #undef BDD_STEP
1148
1149 #define BDDS_MULSUB(vn, un, carry_in)    \
1150   {                                      \
1151     product = ((vn * guess) + carry_in); \
1152     diff = (un - (HD_LOW(product)));     \
1153     if (diff < 0) {                      \
1154       un = (diff + BIGNUM_RADIX_ROOT);   \
1155       carry = ((HD_HIGH(product)) + 1);  \
1156     } else {                             \
1157       un = diff;                         \
1158       carry = (HD_HIGH(product));        \
1159     }                                    \
1160   }
1161
1162 #define BDDS_ADD(vn, un, carry_in)    \
1163   {                                   \
1164     sum = (vn + un + carry_in);       \
1165     if (sum < BIGNUM_RADIX_ROOT) {    \
1166       un = sum;                       \
1167       carry = 0;                      \
1168     } else {                          \
1169       un = (sum - BIGNUM_RADIX_ROOT); \
1170       carry = 1;                      \
1171     }                                 \
1172   }
1173
1174 bignum_digit_type factor_vm::bignum_digit_divide_subtract(
1175     bignum_digit_type v1, bignum_digit_type v2, bignum_digit_type guess,
1176     bignum_digit_type* u) {
1177   {
1178     bignum_digit_type product;
1179     bignum_digit_type diff;
1180     bignum_digit_type carry;
1181     BDDS_MULSUB(v2, (u[2]), 0);
1182     BDDS_MULSUB(v1, (u[1]), carry);
1183     if (carry == 0)
1184       return (guess);
1185     diff = ((u[0]) - carry);
1186     if (diff < 0)
1187       (u[0]) = (diff + BIGNUM_RADIX);
1188     else {
1189       (u[0]) = diff;
1190       return (guess);
1191     }
1192   }
1193   {
1194     bignum_digit_type sum;
1195     bignum_digit_type carry;
1196     BDDS_ADD(v2, (u[2]), 0);
1197     BDDS_ADD(v1, (u[1]), carry);
1198     if (carry == 1)
1199       (u[0]) += 1;
1200   }
1201   return (guess - 1);
1202 }
1203
1204 #undef BDDS_MULSUB
1205 #undef BDDS_ADD
1206
1207 // Allocates memory
1208 void factor_vm::bignum_divide_unsigned_small_denominator(
1209     bignum* numerator_, bignum_digit_type denominator, bignum** quotient,
1210     bignum** remainder, int q_negative_p, int r_negative_p) {
1211   data_root<bignum> numerator(numerator_, this);
1212
1213   bignum* q_ = bignum_new_sign(numerator.untagged(), q_negative_p);
1214   data_root<bignum> q(q_, this);
1215
1216   bignum_digit_type r = bignum_destructive_scale_down(q.untagged(), denominator);
1217
1218   q.set_untagged(bignum_trim(q.untagged()));
1219
1220   if (remainder != ((bignum**)0))
1221     (*remainder) = bignum_digit_to_bignum(r, r_negative_p);
1222
1223   (*quotient) = q.untagged();
1224
1225   return;
1226 }
1227
1228 // Given (denominator > 1), it is fairly easy to show that
1229 // (quotient_high < BIGNUM_RADIX_ROOT), after which it is easy to see
1230 // that all digits are < BIGNUM_RADIX.
1231
1232 bignum_digit_type factor_vm::bignum_destructive_scale_down(
1233     bignum* bn, bignum_digit_type denominator) {
1234   bignum_digit_type numerator;
1235   bignum_digit_type remainder = 0;
1236   bignum_digit_type two_digits;
1237 #define quotient_high remainder
1238   bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1239   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(bn)));
1240   BIGNUM_ASSERT((denominator > 1) && (denominator < BIGNUM_RADIX_ROOT));
1241   while (start < scan) {
1242     two_digits = (*--scan);
1243     numerator = (HD_CONS(remainder, (HD_HIGH(two_digits))));
1244     quotient_high = (numerator / denominator);
1245     numerator = (HD_CONS((numerator % denominator), (HD_LOW(two_digits))));
1246     (*scan) = (HD_CONS(quotient_high, (numerator / denominator)));
1247     remainder = (numerator % denominator);
1248   }
1249   return (remainder);
1250 #undef quotient_high
1251 }
1252
1253 // Allocates memory
1254 bignum* factor_vm::bignum_remainder_unsigned_small_denominator(
1255     bignum* n, bignum_digit_type d, int negative_p) {
1256   bignum_digit_type two_digits;
1257   bignum_digit_type* start = (BIGNUM_START_PTR(n));
1258   bignum_digit_type* scan = (start + (BIGNUM_LENGTH(n)));
1259   bignum_digit_type r = 0;
1260   BIGNUM_ASSERT((d > 1) && (d < BIGNUM_RADIX_ROOT));
1261   while (start < scan) {
1262     two_digits = (*--scan);
1263     r = ((HD_CONS(((HD_CONS(r, (HD_HIGH(two_digits)))) % d),
1264                   (HD_LOW(two_digits)))) %
1265          d);
1266   }
1267   return (bignum_digit_to_bignum(r, negative_p));
1268 }
1269
1270 // Allocates memory
1271 bignum* factor_vm::bignum_digit_to_bignum(bignum_digit_type digit,
1272                                           int negative_p) {
1273   if (digit == 0)
1274     return (BIGNUM_ZERO());
1275   else {
1276     bignum* result = (allot_bignum(1, negative_p));
1277     (BIGNUM_REF(result, 0)) = digit;
1278     return (result);
1279   }
1280 }
1281
1282 // Allocates memory
1283 bignum* factor_vm::allot_bignum(bignum_length_type length, int negative_p) {
1284   BIGNUM_ASSERT((length >= 0) || (length < BIGNUM_RADIX));
1285   bignum* result = allot_uninitialized_array<bignum>(length + 1);
1286   BIGNUM_SET_NEGATIVE_P(result, negative_p);
1287   return (result);
1288 }
1289
1290 // Allocates memory
1291 bignum* factor_vm::allot_bignum_zeroed(bignum_length_type length,
1292                                        int negative_p) {
1293   bignum* result = allot_bignum(length, negative_p);
1294   bignum_digit_type* scan = (BIGNUM_START_PTR(result));
1295   bignum_digit_type* end = (scan + length);
1296   while (scan < end)
1297     (*scan++) = 0;
1298   return (result);
1299 }
1300
1301 // Allocates memory
1302 bignum* factor_vm::bignum_shorten_length(bignum* bn,
1303                                          bignum_length_type length) {
1304   bignum_length_type current_length = (BIGNUM_LENGTH(bn));
1305   BIGNUM_ASSERT((length >= 0) || (length <= current_length));
1306   if (length < current_length) {
1307     bn = reallot_array(bn, length + 1);
1308     BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
1309   }
1310   return (bn);
1311 }
1312
1313 // Allocates memory
1314 bignum* factor_vm::bignum_trim(bignum* bn) {
1315   bignum_digit_type* start = (BIGNUM_START_PTR(bn));
1316   bignum_digit_type* end = (start + (BIGNUM_LENGTH(bn)));
1317   bignum_digit_type* scan = end;
1318   while ((start <= scan) && ((*--scan) == 0))
1319     ;
1320   scan += 1;
1321   if (scan < end) {
1322     bignum_length_type length = (scan - start);
1323     bn = reallot_array(bn, length + 1);
1324     BIGNUM_SET_NEGATIVE_P(bn, (length != 0) && (BIGNUM_NEGATIVE_P(bn)));
1325   }
1326   return (bn);
1327 }
1328
1329 // Copying
1330
1331 // Allocates memory
1332 bignum* factor_vm::bignum_new_sign(bignum* x_, int negative_p) {
1333   data_root<bignum> x(x_, this);
1334   bignum* result = allot_bignum(BIGNUM_LENGTH(x), negative_p);
1335   bignum_destructive_copy(x.untagged(), result);
1336   return result;
1337 }
1338
1339 // Allocates memory
1340 bignum* factor_vm::bignum_maybe_new_sign(bignum* x_, int negative_p) {
1341   if ((BIGNUM_NEGATIVE_P(x_)) ? negative_p : (!negative_p))
1342     return x_;
1343   else {
1344     return bignum_new_sign(x_, negative_p);
1345   }
1346 }
1347
1348 void factor_vm::bignum_destructive_copy(bignum* source, bignum* target) {
1349   bignum_digit_type* scan_source = (BIGNUM_START_PTR(source));
1350   bignum_digit_type* end_source = (scan_source + (BIGNUM_LENGTH(source)));
1351   bignum_digit_type* scan_target = (BIGNUM_START_PTR(target));
1352   while (scan_source < end_source)
1353     (*scan_target++) = (*scan_source++);
1354   return;
1355 }
1356
1357 // * Added bitwise operations (and oddp).
1358
1359 // Allocates memory
1360 bignum* factor_vm::bignum_bitwise_not(bignum* x_) {
1361
1362   int carry = 1;
1363   bignum_length_type size = BIGNUM_LENGTH(x_);
1364   int is_negative = BIGNUM_NEGATIVE_P(x_);
1365   data_root<bignum> x(x_, this);
1366   data_root<bignum> y(allot_bignum(size, is_negative ? 0 : 1), this);
1367
1368   bignum_digit_type* scan_x = BIGNUM_START_PTR(x);
1369   bignum_digit_type* end_x = scan_x + size;
1370   bignum_digit_type* scan_y = BIGNUM_START_PTR(y);
1371
1372   if (is_negative) {
1373     while (scan_x < end_x) {
1374       if (*scan_x == 0) {
1375         *scan_y++ = BIGNUM_RADIX - 1;
1376         scan_x++;
1377       } else {
1378         *scan_y++ = *scan_x++ - 1;
1379         carry = 0;
1380         break;
1381       }
1382     }
1383   } else {
1384     while (scan_x < end_x) {
1385       if (*scan_x == (BIGNUM_RADIX - 1)) {
1386         *scan_y++ = 0;
1387         scan_x++;
1388       } else {
1389         *scan_y++ = *scan_x++ + 1;
1390         carry = 0;
1391         break;
1392       }
1393     }
1394   }
1395
1396   while (scan_x < end_x) {
1397     *scan_y++ = *scan_x++;
1398   }
1399
1400   if (carry) {
1401     bignum* ret = allot_bignum(size + 1, BIGNUM_NEGATIVE_P(y));
1402     bignum_destructive_copy(y.untagged(), ret);
1403     bignum_digit_type* ret_start = BIGNUM_START_PTR(ret);
1404     *(ret_start + size) = 1;
1405     return ret;
1406   } else {
1407     return bignum_trim(y.untagged());
1408   }
1409 }
1410
1411 // Allocates memory
1412 bignum* factor_vm::bignum_arithmetic_shift(bignum* arg1, fixnum n) {
1413   if (BIGNUM_NEGATIVE_P(arg1) && n < 0)
1414     return bignum_bitwise_not(
1415         bignum_magnitude_ash(bignum_bitwise_not(arg1), n));
1416   else
1417     return bignum_magnitude_ash(arg1, n);
1418 }
1419
1420 #define AND_OP 0
1421 #define IOR_OP 1
1422 #define XOR_OP 2
1423
1424 // Allocates memory
1425 bignum* factor_vm::bignum_bitwise_and(bignum* arg1, bignum* arg2) {
1426   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1427               ? bignum_negneg_bitwise_op(AND_OP, arg1, arg2)
1428               : bignum_posneg_bitwise_op(AND_OP, arg2, arg1)
1429               : (BIGNUM_NEGATIVE_P(arg2))
1430               ? bignum_posneg_bitwise_op(AND_OP, arg1, arg2)
1431               : bignum_pospos_bitwise_op(AND_OP, arg1, arg2));
1432 }
1433
1434 // Allocates memory
1435 bignum* factor_vm::bignum_bitwise_ior(bignum* arg1, bignum* arg2) {
1436   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1437               ? bignum_negneg_bitwise_op(IOR_OP, arg1, arg2)
1438               : bignum_posneg_bitwise_op(IOR_OP, arg2, arg1)
1439               : (BIGNUM_NEGATIVE_P(arg2))
1440               ? bignum_posneg_bitwise_op(IOR_OP, arg1, arg2)
1441               : bignum_pospos_bitwise_op(IOR_OP, arg1, arg2));
1442 }
1443
1444 // Allocates memory
1445 bignum* factor_vm::bignum_bitwise_xor(bignum* arg1, bignum* arg2) {
1446   return ((BIGNUM_NEGATIVE_P(arg1)) ? (BIGNUM_NEGATIVE_P(arg2))
1447               ? bignum_negneg_bitwise_op(XOR_OP, arg1, arg2)
1448               : bignum_posneg_bitwise_op(XOR_OP, arg2, arg1)
1449               : (BIGNUM_NEGATIVE_P(arg2))
1450               ? bignum_posneg_bitwise_op(XOR_OP, arg1, arg2)
1451               : bignum_pospos_bitwise_op(XOR_OP, arg1, arg2));
1452 }
1453
1454 // Allocates memory
1455 // ash for the magnitude
1456 // assume arg1 is a big number, n is a long
1457 bignum* factor_vm::bignum_magnitude_ash(bignum* arg1_, fixnum n) {
1458
1459   data_root<bignum> arg1(arg1_, this);
1460
1461   bignum* result = NULL;
1462   bignum_digit_type* scan1;
1463   bignum_digit_type* scanr;
1464   bignum_digit_type* end;
1465
1466   fixnum digit_offset, bit_offset;
1467
1468   if (BIGNUM_ZERO_P(arg1))
1469     return arg1.untagged();
1470
1471   if (n > 0) {
1472     digit_offset = n / BIGNUM_DIGIT_LENGTH;
1473     bit_offset = n % BIGNUM_DIGIT_LENGTH;
1474
1475     result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) + digit_offset + 1,
1476                                  BIGNUM_NEGATIVE_P(arg1));
1477
1478     scanr = BIGNUM_START_PTR(result) + digit_offset;
1479     scan1 = BIGNUM_START_PTR(arg1);
1480     end = scan1 + BIGNUM_LENGTH(arg1);
1481
1482     while (scan1 < end) {
1483       *scanr = *scanr | (*scan1 & BIGNUM_DIGIT_MASK) << bit_offset;
1484       *scanr = *scanr & BIGNUM_DIGIT_MASK;
1485       scanr++;
1486       *scanr = *scan1++ >> (BIGNUM_DIGIT_LENGTH - bit_offset);
1487       *scanr = *scanr & BIGNUM_DIGIT_MASK;
1488     }
1489   } else if (n < 0 && (-n >= (BIGNUM_LENGTH(arg1) * (bignum_length_type)
1490                               BIGNUM_DIGIT_LENGTH))) {
1491     result = BIGNUM_ZERO();
1492   } else if (n < 0) {
1493     digit_offset = -n / BIGNUM_DIGIT_LENGTH;
1494     bit_offset = -n % BIGNUM_DIGIT_LENGTH;
1495
1496     result = allot_bignum_zeroed(BIGNUM_LENGTH(arg1) - digit_offset,
1497                                  BIGNUM_NEGATIVE_P(arg1));
1498
1499     scanr = BIGNUM_START_PTR(result);
1500     scan1 = BIGNUM_START_PTR(arg1) + digit_offset;
1501     end = scanr + BIGNUM_LENGTH(result) - 1;
1502
1503     while (scanr < end) {
1504       *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1505       *scanr = (*scanr | *scan1 << (BIGNUM_DIGIT_LENGTH - bit_offset)) &
1506                BIGNUM_DIGIT_MASK;
1507       scanr++;
1508     }
1509     *scanr = (*scan1++ & BIGNUM_DIGIT_MASK) >> bit_offset;
1510   } else if (n == 0) {
1511     result = arg1.untagged();
1512   }
1513
1514   return bignum_trim(result);
1515 }
1516
1517 // Allocates memory
1518 bignum* factor_vm::bignum_pospos_bitwise_op(int op, bignum* arg1_,
1519                                             bignum* arg2_) {
1520   data_root<bignum> arg1(arg1_, this);
1521   data_root<bignum> arg2(arg2_, this);
1522
1523   bignum_length_type max_length;
1524
1525   bignum_digit_type* scan1, *end1, digit1;
1526   bignum_digit_type* scan2, *end2, digit2;
1527   bignum_digit_type* scanr, *endr;
1528
1529   max_length =
1530       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1)
1531                                                   : BIGNUM_LENGTH(arg2);
1532
1533   bignum* result = allot_bignum(max_length, 0);
1534
1535   scanr = BIGNUM_START_PTR(result);
1536   scan1 = BIGNUM_START_PTR(arg1);
1537   scan2 = BIGNUM_START_PTR(arg2);
1538   endr = scanr + max_length;
1539   end1 = scan1 + BIGNUM_LENGTH(arg1);
1540   end2 = scan2 + BIGNUM_LENGTH(arg2);
1541
1542   while (scanr < endr) {
1543     digit1 = (scan1 < end1) ? *scan1++ : 0;
1544     digit2 = (scan2 < end2) ? *scan2++ : 0;
1545     *scanr++ =
1546         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1547                                                           : digit1 ^ digit2;
1548   }
1549   return bignum_trim(result);
1550 }
1551
1552 // Allocates memory
1553 bignum* factor_vm::bignum_posneg_bitwise_op(int op, bignum* arg1_,
1554                                             bignum* arg2_) {
1555   data_root<bignum> arg1(arg1_, this);
1556   data_root<bignum> arg2(arg2_, this);
1557
1558   bignum_length_type max_length;
1559
1560   bignum_digit_type* scan1, *end1, digit1;
1561   bignum_digit_type* scan2, *end2, digit2, carry2;
1562   bignum_digit_type* scanr, *endr;
1563
1564   char neg_p = op == IOR_OP || op == XOR_OP;
1565
1566   max_length =
1567       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2) + 1) ? BIGNUM_LENGTH(arg1)
1568                                                       : BIGNUM_LENGTH(arg2) + 1;
1569
1570   bignum* result = allot_bignum(max_length, neg_p);
1571
1572   scanr = BIGNUM_START_PTR(result);
1573   scan1 = BIGNUM_START_PTR(arg1);
1574   scan2 = BIGNUM_START_PTR(arg2);
1575   endr = scanr + max_length;
1576   end1 = scan1 + BIGNUM_LENGTH(arg1);
1577   end2 = scan2 + BIGNUM_LENGTH(arg2);
1578
1579   carry2 = 1;
1580
1581   while (scanr < endr) {
1582     digit1 = (scan1 < end1) ? *scan1++ : 0;
1583     digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1584
1585     if (digit2 < BIGNUM_RADIX)
1586       carry2 = 0;
1587     else {
1588       digit2 = (digit2 - BIGNUM_RADIX);
1589       carry2 = 1;
1590     }
1591
1592     *scanr++ =
1593         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1594                                                           : digit1 ^ digit2;
1595   }
1596
1597   if (neg_p)
1598     bignum_negate_magnitude(result);
1599
1600   return bignum_trim(result);
1601 }
1602
1603 // Allocates memory
1604 bignum* factor_vm::bignum_negneg_bitwise_op(int op, bignum* arg1_,
1605                                             bignum* arg2_) {
1606   data_root<bignum> arg1(arg1_, this);
1607   data_root<bignum> arg2(arg2_, this);
1608
1609   bignum_length_type max_length;
1610
1611   bignum_digit_type* scan1, *end1, digit1, carry1;
1612   bignum_digit_type* scan2, *end2, digit2, carry2;
1613   bignum_digit_type* scanr, *endr;
1614
1615   char neg_p = op == AND_OP || op == IOR_OP;
1616
1617   max_length =
1618       (BIGNUM_LENGTH(arg1) > BIGNUM_LENGTH(arg2)) ? BIGNUM_LENGTH(arg1) + 1
1619                                                   : BIGNUM_LENGTH(arg2) + 1;
1620
1621   bignum* result = allot_bignum(max_length, neg_p);
1622
1623   scanr = BIGNUM_START_PTR(result);
1624   scan1 = BIGNUM_START_PTR(arg1);
1625   scan2 = BIGNUM_START_PTR(arg2);
1626   endr = scanr + max_length;
1627   end1 = scan1 + BIGNUM_LENGTH(arg1);
1628   end2 = scan2 + BIGNUM_LENGTH(arg2);
1629
1630   carry1 = 1;
1631   carry2 = 1;
1632
1633   while (scanr < endr) {
1634     digit1 = (~((scan1 < end1) ? *scan1++ : 0) & BIGNUM_DIGIT_MASK) + carry1;
1635     digit2 = (~((scan2 < end2) ? *scan2++ : 0) & BIGNUM_DIGIT_MASK) + carry2;
1636
1637     if (digit1 < BIGNUM_RADIX)
1638       carry1 = 0;
1639     else {
1640       digit1 = (digit1 - BIGNUM_RADIX);
1641       carry1 = 1;
1642     }
1643
1644     if (digit2 < BIGNUM_RADIX)
1645       carry2 = 0;
1646     else {
1647       digit2 = (digit2 - BIGNUM_RADIX);
1648       carry2 = 1;
1649     }
1650
1651     *scanr++ =
1652         (op == AND_OP) ? digit1 & digit2 : (op == IOR_OP) ? digit1 | digit2
1653                                                           : digit1 ^ digit2;
1654   }
1655
1656   if (neg_p)
1657     bignum_negate_magnitude(result);
1658
1659   return bignum_trim(result);
1660 }
1661
1662 void factor_vm::bignum_negate_magnitude(bignum* arg) {
1663   bignum_digit_type* scan;
1664   bignum_digit_type* end;
1665   bignum_digit_type digit;
1666   bignum_digit_type carry;
1667
1668   scan = BIGNUM_START_PTR(arg);
1669   end = scan + BIGNUM_LENGTH(arg);
1670
1671   carry = 1;
1672
1673   while (scan < end) {
1674     digit = (~ * scan & BIGNUM_DIGIT_MASK) + carry;
1675
1676     if (digit < BIGNUM_RADIX)
1677       carry = 0;
1678     else {
1679       digit = (digit - BIGNUM_RADIX);
1680       carry = 1;
1681     }
1682
1683     *scan++ = digit;
1684   }
1685 }
1686
1687 // Allocates memory
1688 bignum* factor_vm::bignum_integer_length(bignum* x_) {
1689   data_root<bignum> x(x_, this);
1690   bignum_length_type index = ((BIGNUM_LENGTH(x)) - 1);
1691   bignum_digit_type digit = (BIGNUM_REF(x, index));
1692   bignum_digit_type carry = 0;
1693   bignum* result;
1694
1695   while (digit > 1) {
1696     carry += 1;
1697     digit >>= 1;
1698   }
1699
1700   if (index < BIGNUM_RADIX_ROOT) {
1701      result = allot_bignum(1, 0);
1702      (BIGNUM_REF(result, 0)) = (index * BIGNUM_DIGIT_LENGTH) + carry;
1703   } else {
1704      result = allot_bignum(2, 0);
1705      (BIGNUM_REF(result, 0)) = index;
1706      (BIGNUM_REF(result, 1)) = 0;
1707      bignum_destructive_scale_up(result, BIGNUM_DIGIT_LENGTH);
1708      bignum_destructive_add(result, carry);
1709   }
1710   return (bignum_trim(result));
1711 }
1712
1713 // Allocates memory
1714 int factor_vm::bignum_logbitp(int shift, bignum* arg) {
1715   return ((BIGNUM_NEGATIVE_P(arg))
1716               ? !bignum_unsigned_logbitp(shift, bignum_bitwise_not(arg))
1717               : bignum_unsigned_logbitp(shift, arg));
1718 }
1719
1720 int factor_vm::bignum_unsigned_logbitp(int shift, bignum* bn) {
1721   bignum_length_type len = (BIGNUM_LENGTH(bn));
1722   int index = shift / BIGNUM_DIGIT_LENGTH;
1723   if (index >= len)
1724     return 0;
1725   bignum_digit_type digit = (BIGNUM_REF(bn, index));
1726   int p = shift % BIGNUM_DIGIT_LENGTH;
1727   bignum_digit_type mask = ((fixnum)1) << p;
1728   return (digit & mask) ? 1 : 0;
1729 }
1730
1731 #ifdef _WIN64
1732 // Allocates memory.
1733 bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {
1734
1735   data_root<bignum> a(a_, this);
1736   data_root<bignum> b(b_, this);
1737
1738   // Copies of a and b with that are both positive.
1739   data_root<bignum> ac(bignum_maybe_new_sign(a.untagged(), 0), this);
1740   data_root<bignum> bc(bignum_maybe_new_sign(b.untagged(), 0), this);
1741
1742   if (bignum_compare(ac.untagged(), bc.untagged()) == BIGNUM_COMPARISON_LESS) {
1743     swap(ac, bc);
1744   }
1745
1746   while (BIGNUM_LENGTH(bc) != 0) {
1747     data_root<bignum> d(bignum_remainder(ac.untagged(), bc.untagged()), this);
1748     if (d.untagged() == BIGNUM_OUT_OF_BAND) {
1749       return d.untagged();
1750     }
1751     ac = bc;
1752     bc = d;
1753   }
1754   return ac.untagged();
1755 }
1756 #else
1757 // Allocates memory
1758 bignum* factor_vm::bignum_gcd(bignum* a_, bignum* b_) {
1759   data_root<bignum> a(a_, this);
1760   data_root<bignum> b(b_, this);
1761   bignum_twodigit_type x, y, q, s, t, A, B, C, D;
1762   unsigned long nbits;
1763   int k;
1764   bignum_length_type size_a, size_b, size_c;
1765   bignum_digit_type* scan_a, *scan_b, *scan_c, *scan_d;
1766   bignum_digit_type* a_end, *b_end, *c_end;
1767
1768   // clone the bignums so we can modify them in-place
1769   size_a = BIGNUM_LENGTH(a);
1770   data_root<bignum> c(allot_bignum(size_a, 0), this);
1771   // c = allot_bignum(size_a, 0);
1772   scan_a = BIGNUM_START_PTR(a);
1773   a_end = scan_a + size_a;
1774   scan_c = BIGNUM_START_PTR(c);
1775   while (scan_a < a_end)
1776     (*scan_c++) = (*scan_a++);
1777   a = c;
1778   size_b = BIGNUM_LENGTH(b);
1779   data_root<bignum> d(allot_bignum(size_b, 0), this);
1780   scan_b = BIGNUM_START_PTR(b);
1781   b_end = scan_b + size_b;
1782   scan_d = BIGNUM_START_PTR(d);
1783   while (scan_b < b_end)
1784     (*scan_d++) = (*scan_b++);
1785   b = d;
1786
1787   // Initial reduction: make sure that 0 <= b <= a.
1788   if (bignum_compare(a.untagged(), b.untagged()) == BIGNUM_COMPARISON_LESS) {
1789     swap(a, b);
1790     std::swap(size_a, size_b);
1791   }
1792
1793   while (size_a > 1) {
1794     nbits = log2(BIGNUM_REF(a, size_a - 1));
1795     x = ((BIGNUM_REF(a, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)) |
1796          (BIGNUM_REF(a, size_a - 2) >> nbits));
1797     y = ((size_b >= size_a - 1 ? BIGNUM_REF(b, size_a - 2) >> nbits : 0) |
1798          (size_b >= size_a
1799               ? BIGNUM_REF(b, size_a - 1) << (BIGNUM_DIGIT_LENGTH - nbits)
1800               : 0));
1801
1802     // inner loop of Lehmer's algorithm;
1803     A = 1;
1804     B = 0;
1805     C = 0;
1806     D = 1;
1807     for (k = 0;; k++) {
1808       if (y - C == 0)
1809         break;
1810
1811       q = (x + (A - 1)) / (y - C);
1812
1813       s = B + (q * D);
1814       t = x - (q * y);
1815
1816       if (s > t)
1817         break;
1818
1819       x = y;
1820       y = t;
1821
1822       t = A + (q * C);
1823
1824       A = D;
1825       B = C;
1826       C = s;
1827       D = t;
1828     }
1829
1830     if (k == 0) {
1831       // no progress; do a Euclidean step
1832       if (size_b == 0) {
1833         return bignum_trim(a.untagged());
1834       }
1835       data_root<bignum> e(bignum_trim(a.untagged()), this);
1836       data_root<bignum> f(bignum_trim(b.untagged()), this);
1837       c.set_untagged(bignum_remainder(e.untagged(), f.untagged()));
1838       if (c.untagged() == BIGNUM_OUT_OF_BAND) {
1839         return c.untagged();
1840       }
1841
1842       // copy 'b' to 'a'
1843       scan_a = BIGNUM_START_PTR(a);
1844       scan_b = BIGNUM_START_PTR(b);
1845       a_end = scan_a + size_a;
1846       b_end = scan_b + size_b;
1847       while (scan_b < b_end)
1848         *(scan_a++) = *(scan_b++);
1849       while (scan_a < a_end)
1850         *(scan_a++) = 0;
1851       size_a = size_b;
1852
1853       // copy 'c' to 'b'
1854       scan_b = BIGNUM_START_PTR(b);
1855       scan_c = BIGNUM_START_PTR(c);
1856       size_c = BIGNUM_LENGTH(c);
1857       c_end = scan_c + size_c;
1858       while (scan_c < c_end)
1859         *(scan_b++) = *(scan_c++);
1860       while (scan_b < b_end)
1861         *(scan_b++) = 0;
1862       size_b = size_c;
1863
1864       continue;
1865     }
1866
1867     // a, b = A*b - B*a, D*a - C*b if k is odd
1868     // a, b = A*a - B*b, D*b - C*a if k is even
1869
1870     scan_a = BIGNUM_START_PTR(a);
1871     scan_b = BIGNUM_START_PTR(b);
1872     scan_c = scan_a;
1873     scan_d = scan_b;
1874     a_end = scan_a + size_a;
1875     b_end = scan_b + size_b;
1876     s = 0;
1877     t = 0;
1878     if (k & 1) {
1879       while (scan_b < b_end) {
1880         s += (A * *scan_b) - (B * *scan_a);
1881         t += (D * *scan_a++) - (C * *scan_b++);
1882         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1883         *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1884         s >>= BIGNUM_DIGIT_LENGTH;
1885         t >>= BIGNUM_DIGIT_LENGTH;
1886       }
1887       while (scan_a < a_end) {
1888         s -= (B * *scan_a);
1889         t += (D * *scan_a++);
1890         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1891         //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1892         s >>= BIGNUM_DIGIT_LENGTH;
1893         t >>= BIGNUM_DIGIT_LENGTH;
1894       }
1895     } else {
1896       while (scan_b < b_end) {
1897         s += (A * *scan_a) - (B * *scan_b);
1898         t += (D * *scan_b++) - (C * *scan_a++);
1899         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1900         *scan_d++ = (bignum_digit_type)(t & BIGNUM_DIGIT_MASK);
1901         s >>= BIGNUM_DIGIT_LENGTH;
1902         t >>= BIGNUM_DIGIT_LENGTH;
1903       }
1904       while (scan_a < a_end) {
1905         s += (A * *scan_a);
1906         t -= (C * *scan_a++);
1907         *scan_c++ = (bignum_digit_type)(s & BIGNUM_DIGIT_MASK);
1908         //*scan_d++ = (bignum_digit_type) (t & BIGNUM_DIGIT_MASK);
1909         s >>= BIGNUM_DIGIT_LENGTH;
1910         t >>= BIGNUM_DIGIT_LENGTH;
1911       }
1912     }
1913     BIGNUM_ASSERT(s == 0);
1914     BIGNUM_ASSERT(t == 0);
1915
1916     // update size_a and size_b to remove any zeros at end
1917     while (size_a > 0 && *(--scan_a) == 0)
1918       size_a--;
1919     while (size_b > 0 && *(--scan_b) == 0)
1920       size_b--;
1921
1922     BIGNUM_ASSERT(size_a >= size_b);
1923   }
1924
1925   // a fits into a fixnum, so b must too
1926   fixnum xx = bignum_to_fixnum(a.untagged());
1927   fixnum yy = bignum_to_fixnum(b.untagged());
1928   fixnum tt;
1929
1930   // usual Euclidean algorithm for longs
1931   while (yy != 0) {
1932     tt = yy;
1933     yy = xx % yy;
1934     xx = tt;
1935   }
1936
1937   return fixnum_to_bignum(xx);
1938 }
1939 #endif
1940
1941 }