vm/bignum.cpp

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