vm/bignum.cpp

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