vm/bignum.cpp

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