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