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