vm/bignum.cpp

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