vm/bignum.cpp

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