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