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