vm/bignum.cpp

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