vm/bignum.c

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