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