lib-src/numeric/bigint.cc

   1 //////////////////////////////////////////////////////////////////////////////
   2 // NOTICE:
   3 //
   4 // ADLib, Prop and their related set of tools and documentation are in the
   5 // public domain.   The author(s) of this software reserve no copyrights on
   6 // the source code and any code generated using the tools.  You are encouraged
   7 // to use ADLib and Prop to develop software, in both academic and commercial
   8 // settings, and are free to incorporate any part of ADLib and Prop into
   9 // your programs.
  10 //
  11 // Although you are under no obligation to do so, we strongly recommend that
  12 // you give away all software developed using our tools.
  13 //
  14 // We also ask that credit be given to us when ADLib and/or Prop are used in
  15 // your programs, and that this notice be preserved intact in all the source
  16 // code.
  17 //
  18 // This software is still under development and we welcome any suggestions
  19 // and help from the users.
  20 //
  21 // Allen Leung
  22 // 1994
  23 //////////////////////////////////////////////////////////////////////////////
  24
  25 //////////////////////////////////////////////////////////////////////////////
  26 // Variable length integers are represented using a sign/magnitude format.
  27 // Arithmetic is done using 16bit unsigned integers.  The machine is
  28 // assumed to implement unsigned 32 bit arithmetic as arithmetic modulo 2^32.
  29 //
  30 // This version is inspired by various ``bignum'' code, including G++'s
  31 // Integer class and SmallTalk's LargeInteger class.
  32 //////////////////////////////////////////////////////////////////////////////
  33
  34 #include <iostream>
  35 #include <stdlib.h>
  36 #include <string.h>
  37 #include <new>
  38 #include <ctype.h>
  39 #include <AD/generic/generic.h>
  40 #include <AD/numeric/bigint.h>
  41
  42 //////////////////////////////////////////////////////////////////////////////
  43 //  Make hidden types available
  44 //////////////////////////////////////////////////////////////////////////////
  45 typedef BigInt::a_BigInt a_BigInt;
  46 typedef a_BigInt::Digit  Digit;
  47 typedef a_BigInt::Unit   Unit;
  48
  49 //////////////////////////////////////////////////////////////////////////////
  50 //  Some manifest constants: assume 8 bits in a byte
  51 //////////////////////////////////////////////////////////////////////////////
  52 #define Bits_per_digit  (8 * sizeof(Digit))
  53 #define Max_digit       (1 << Bits_per_digit)
  54 #define Max_Unit        (~(Unit)0)
  55 #define LOW(n)          ((n) & (Max_digit - 1))
  56 #define HIGH(n)         ((n) >> Bits_per_digit)
  57
  58 //////////////////////////////////////////////////////////////////////////////
  59 //  Auxiliary routines
  60 //////////////////////////////////////////////////////////////////////////////
  61 inline int min(int a, int b)     { return a < b ? a : b; }
  62 inline int max(int a, int b)     { return a > b ? a : b; }
  63 inline unsigned long abs(long a) { return a > 0 ? a : -a; }
  64
  65 //////////////////////////////////////////////////////////////////////////////
  66 //  Default error handler
  67 //////////////////////////////////////////////////////////////////////////////
  68 static void bigint_error(const char message[])
  69 {  std::cerr << "BigInt error: " << message << '\n';
  70    exit(1);
  71 }
  72
  73 //////////////////////////////////////////////////////////////////////////////
  74 //  Error handler
  75 //////////////////////////////////////////////////////////////////////////////
  76 void (*BigInt::error_handler)(const char []) = bigint_error;
  77
  78 //////////////////////////////////////////////////////////////////////////////
  79 //  Representation for common values; these have capacity of 0 so that
  80 //  they won't be erroneously deallocated.
  81 //////////////////////////////////////////////////////////////////////////////
  82 static a_BigInt zero_rep = { 0,  0, 0, { 0 } };
  83 a_BigInt * BigInt::zero  = &zero_rep;
  84
  85 //////////////////////////////////////////////////////////////////////////////
  86 // Memory management
  87 //////////////////////////////////////////////////////////////////////////////
  88 inline void * a_BigInt::operator new(size_t size, int capacity)
  89 {  register int count = 2;
  90    while (count < capacity) count <<= 1;
  91    a_BigInt * x = (a_BigInt*)malloc(size + (count - 1) * sizeof(Digit));
  92    x->capacity = count;
  93    return x;
  94 }
  95
  96 void a_BigInt::operator delete(void * x)
  97 { if (((a_BigInt*)x)->capacity > 0) free(x); }
  98
  99 //////////////////////////////////////////////////////////////////////////////
 100 //  Normalize a number by eliminating leading zeros
 101 //////////////////////////////////////////////////////////////////////////////
 102 inline void normalize(a_BigInt * r)
 103 {  register int i;
 104    for (i = r->len-1; i >= 0; i--)
 105       if (r->digits[i] != 0) break;
 106    if ((r->len = i+1) == 0) r->sign = 0;
 107 }
 108
 109 //////////////////////////////////////////////////////////////////////////////
 110 // Delete a BigInt
 111 //////////////////////////////////////////////////////////////////////////////
 112 inline void kill(a_BigInt * r) { if (r && r->capacity > 0) free(r); }
 113
 114 //////////////////////////////////////////////////////////////////////////////
 115 // Make zero, one, or minus one
 116 //////////////////////////////////////////////////////////////////////////////
 117 inline a_BigInt * make_zero(a_BigInt * r)
 118 {  kill(r); return &zero_rep; }
 119
 120 //////////////////////////////////////////////////////////////////////////////
 121 //  Copy a bigint; reallocate if necessary
 122 //////////////////////////////////////////////////////////////////////////////
 123 a_BigInt * copy(a_BigInt * r, const a_BigInt * a)
 124 {  if (r == a) return r;
 125    if (r == 0) {
 126       r = new (a->len) BigInt::a_BigInt;
 127    } else {
 128       if ((r->capacity) < (a->len)) {
 129          kill(r); r = new (a->len) BigInt::a_BigInt;
 130       }
 131    }
 132    r->len = a->len; r->sign = a->sign;
 133    memcpy(r->digits,a->digits,a->len * sizeof(Digit));
 134    return r;
 135 }
 136
 137 //////////////////////////////////////////////////////////////////////////////
 138 //  Copy a bigint and negate; reallocate if necessary
 139 //////////////////////////////////////////////////////////////////////////////
 140 a_BigInt * neg(a_BigInt * r, const a_BigInt * a)
 141 {  if (r == 0 || (r->capacity) < (a->len)) {
 142       kill(r); r = new(a->len) a_BigInt;
 143    }
 144    r->len = a->len; r->sign = - a->sign;
 145    memcpy(r->digits,a->digits,a->len * sizeof(Digit));
 146    return r;
 147 }
 148
 149 //////////////////////////////////////////////////////////////////////////////
 150 //  Copy a long
 151 //////////////////////////////////////////////////////////////////////////////
 152 a_BigInt * copy(a_BigInt * r, long a, int sign)
 153 {  if (a == 0) return make_zero(r);
 154    if (r == 0 || r->capacity == 0) r = new (2) a_BigInt;
 155    Unit mag;
 156    if (a > 0) { r->sign = sign; mag = a; }
 157    else       { r->sign = -sign; mag = -a; }
 158    r->len       = mag >= Max_digit ? 2 : 1;
 159    r->digits[0] = LOW(mag);
 160    r->digits[1] = HIGH(mag);
 161    return r;
 162 }
 163
 164 //////////////////////////////////////////////////////////////////////////////
 165 //  Allocate a BigInt
 166 //////////////////////////////////////////////////////////////////////////////
 167 inline a_BigInt * alloc(a_BigInt * r, int len)
 168 {  if (r == 0 || (r->capacity) < len) return new(len) a_BigInt;
 169    else return r;
 170 }
 171
 172 //////////////////////////////////////////////////////////////////////////////
 173 //  Reallocate a BigInt
 174 //////////////////////////////////////////////////////////////////////////////
 175 inline a_BigInt * realloc(a_BigInt * r, int len)
 176 {  if (r == 0 || (r->capacity) < len) {
 177       a_BigInt * new_r = new(len) a_BigInt;
 178       register Digit * p, * q, * limit;
 179       for (p = r->digits, q = new_r->digits, limit = p + r->len; p < limit;)
 180          *q++ = *p++;
 181       return new_r;
 182    } else return r;
 183 }
 184
 185 //////////////////////////////////////////////////////////////////////////////
 186 //  Clear the upper digits
 187 //////////////////////////////////////////////////////////////////////////////
 188 inline void clear(register Digit * p, int start, int stop)
 189 {  register Digit * limit = p + stop;
 190    for (p += start; p < limit; ) *p++ = 0;
 191 }
 192
 193 //////////////////////////////////////////////////////////////////////////////
 194 //  Comparisions
 195 //////////////////////////////////////////////////////////////////////////////
 196 inline int ucmp(const a_BigInt * a, const a_BigInt * b)
 197 {  int r = a->len - b->len;
 198    if (r == 0) {
 199       register const Digit * x;
 200       register const Digit * y;
 201       register int len;
 202       for (x = a->digits + a->len - 1, y = b->digits + b->len - 1,
 203            len = a->len; len > 0; len--, x--, y--)
 204          if ((r = ((long)*x - (long)*y))) return r;
 205    }
 206    return r;
 207 }
 208
 209 //////////////////////////////////////////////////////////////////////////////
 210 //  Signed comparison between two BigInt's
 211 //////////////////////////////////////////////////////////////////////////////
 212 int cmp(const a_BigInt * a, const a_BigInt * b)
 213 {  if (a->sign > b->sign)            return 1;
 214    if (a->sign < b->sign)            return -1;
 215    if (a->sign == 0 && b->sign == 0) return 0;
 216    int r = a->len - b->len;
 217    if (r == 0) r = ucmp(a,b);
 218    return a->sign == 1 ? r : -r;
 219 }
 220
 221 //////////////////////////////////////////////////////////////////////////////
 222 //  Signed comparison between a BigInt and a long
 223 //////////////////////////////////////////////////////////////////////////////
 224 int cmp(const a_BigInt * a, long b)
 225 {  if (a->sign == 0  && b == 0) return 0;
 226    if (a->sign == 1  && b <= 0) return 1;
 227    if (a->sign == -1 && b >= 0) return -1;
 228    int b_sign;
 229    unsigned long b_mag;
 230    if (b > 0) { b_sign = 1;  b_mag = b; }
 231    else       { b_sign = -1; b_mag = -b; }
 232    int b_len = b_mag >= Max_digit ? 2 : 1;
 233    int r = a->len - b_len;
 234    if (r == 0) {
 235       if (b_len == 2)
 236          r = (long)(a->digits[1] * Max_digit + a->digits[0]) - (long)b_mag;
 237       else
 238          r = (long)a->digits[0] - (long)b_mag;
 239    }
 240    return a->sign ? r : -r;
 241 }
 242
 243 //////////////////////////////////////////////////////////////////////////////
 244 //  Compute the estimated length of an unsigned addition
 245 //  taking into account of possible carry.
 246 //////////////////////////////////////////////////////////////////////////////
 247 inline int uaddlen(const a_BigInt * a, const a_BigInt * b)
 248 {  if (a->len > b->len)
 249       if (a->digits[a->len-1] == Max_digit-1) return a->len + 1;
 250       else return a->len;
 251    if (a->len < b->len)
 252       if (b->digits[b->len-1] == Max_digit-1) return b->len + 1;
 253       else return b->len;
 254    if ((Unit)a->digits[a->len-1] + (Unit)b->digits[b->len-1] >= Max_digit-1)
 255       return a->len + 1;
 256    else
 257       return a->len;
 258 }
 259
 260 //////////////////////////////////////////////////////////////////////////////
 261 //  Unsigned addition
 262 //////////////////////////////////////////////////////////////////////////////
 263 inline void uadd(a_BigInt * r, const a_BigInt * a, const a_BigInt * b)
 264 {   register Unit sum;
 265     register const Digit * x, * y;
 266     register Digit * z;
 267     register int len = min(a->len,b->len);
 268     for (x=a->digits, y=b->digits, z=r->digits, sum = 0; len > 0; len--) {
 269        sum   +=  (Unit)*x++ + (Unit)*y++;
 270        *z++  =   LOW(sum);
 271        sum   =   HIGH(sum);
 272     }
 273     len = a->len - b->len;
 274     if (len < 0) { len = -len; x = y; }
 275     for ( ; len > 0; len--) {
 276        sum   +=  (Unit)*x++;
 277        *z++  =   LOW(sum);
 278        sum   =   HIGH(sum);
 279     }
 280     if (sum) *z = sum;
 281     r->len = z - r->digits;
 282 }
 283
 284 //////////////////////////////////////////////////////////////////////////////
 285 //  Unsigned substraction; assume abs(a) >= abs(b)
 286 //////////////////////////////////////////////////////////////////////////////
 287 inline void usub(a_BigInt * r, const a_BigInt * a, const a_BigInt * b)
 288 {   register Unit diff;
 289     register const Digit * x, * y;
 290     register Digit * z;
 291     register int len = b->len;
 292     for (x=a->digits, y=b->digits, z=r->digits, diff = 1; len > 0; len--) {
 293        diff += (Unit)*x++ - (Unit)*y++ + Max_digit - 1;
 294        *z++  = LOW(diff);
 295        diff  = HIGH(diff);
 296     }
 297     for (len = a->len - b->len; diff != 1 && len > 0; len--) {
 298        diff += (Unit)*x++ + Max_digit - 1;
 299        *z++  = LOW(diff);
 300        diff  = HIGH(diff);
 301     }
 302     for ( ;len > 0; len--) *z++ = *x++;
 303     r->len = a->len;
 304     normalize(r);
 305 }
 306
 307 ////////////////////////////////////////////////////////////////
 308 //  Performs signed addition/subtraction
 309 ////////////////////////////////////////////////////////////////
 310 a_BigInt * addsub
 311    (a_BigInt * r, const a_BigInt * a, const a_BigInt * b, Bool sub)
 312 {  a_BigInt * result;
 313
 314    ////////////////////////////////////////////////////////////////
 315    // Check for zero in one of its arguments
 316    ////////////////////////////////////////////////////////////////
 317    if (a->sign == 0 && b->sign == 0 || a == b && sub) return make_zero(r);
 318    if (a->sign == 0) return sub ? neg(r,b) : copy(r,b);
 319    if (b->sign == 0) return copy(r,a);
 320
 321    ////////////////////////////////////////////////////////////////
 322    // If the signs are different, we'll do a subtraction instead
 323    ////////////////////////////////////////////////////////////////
 324    if (sub ? (a->sign == b->sign) : (a->sign != b->sign)) {
 325       int relative = ucmp(a,b);
 326       if (relative == 0) return make_zero(r);
 327       result = alloc(r,max(a->len,b->len));
 328       if (relative < 0) { const a_BigInt * swap = a; a = b; b = swap; }
 329       usub(result,a,b);
 330       result->sign =
 331          sub ? (relative > 0 ? 1 : -1) :
 332                (relative > 0 ? a->sign : b->sign);
 333    } else {
 334       result = alloc(r,uaddlen(a,b));
 335       uadd(result,a,b);
 336       result->sign = a->sign;
 337    }
 338    if (result != r) kill(r);
 339    return result;
 340 }
 341
 342 //////////////////////////////////////////////////////////////////////////////
 343 //  Signed addition/subtraction with a long
 344 //////////////////////////////////////////////////////////////////////////////
 345 a_BigInt * addsub(a_BigInt * r, const a_BigInt * a, long b, Bool sub)
 346 {  struct {
 347       a_BigInt n;
 348       Digit    d;
 349    } c;
 350    if (b == 0) return copy(r,a);
 351    if (b > 0) { c.n.sign = 1; }
 352    else       { c.n.sign = -1; b = -b; }
 353    c.n.len       = 2;
 354    c.n.digits[0] = LOW(b);
 355    c.n.digits[1] = HIGH(b);
 356    return addsub(r,a,&c.n,sub);
 357 }
 358
 359 //////////////////////////////////////////////////////////////////////////////
 360 // Unsigned multiplication by a small constant
 361 //////////////////////////////////////////////////////////////////////////////
 362 static void umul(a_BigInt * r, const a_BigInt * a, Digit b)
 363 {
 364    if (r != a) {
 365       register const Digit * y, * y_limit;
 366       register Digit * x;
 367       x       = r->digits;
 368       y       = a->digits;
 369       y_limit = y + a->len;
 370       register Unit product;
 371       for ( product = 0; y < y_limit; ) {
 372          product += *y++ * b;
 373          *x++     = LOW(product);
 374          product  = HIGH(product);
 375       }
 376       if (product) *x++ = product;
 377       r->len = x - r->digits;
 378    } else {  // r == a!!!
 379       register Digit * x, * y_limit;
 380       register Digit * y;
 381       x = r->digits + r->len - 1;
 382       y_limit = r->digits + r->len;
 383       *y_limit = 0;
 384       register Unit product;
 385       for (; x >= a->digits; x--) {
 386          product  = *x * b;
 387          *x       = LOW(product);
 388          product  = HIGH(product);
 389          for (y = x+1; product && y <= y_limit; y++) {
 390             product += *y;
 391             *y       = LOW(product);
 392             product  = HIGH(product);
 393          }
 394       }
 395       r->len = a->len + 1;
 396       normalize(r);
 397    }
 398 }
 399
 400 //////////////////////////////////////////////////////////////////////////////
 401 // Unsigned division by a small constant
 402 //////////////////////////////////////////////////////////////////////////////
 403 inline Digit udiv(a_BigInt * r, const a_BigInt * a, Digit b)
 404 {   register Digit * x;
 405     register const Digit * y;
 406     x = r ? r->digits + a->len - 1 : 0;
 407     y = a->digits + a->len - 1;
 408     register Unit remainder;
 409     for (remainder = 0; y >= a->digits; ) {
 410        remainder += *y--;
 411        Digit quotient = remainder / b;
 412        if (x) *x-- = quotient;
 413        remainder = (remainder - b * quotient) << Bits_per_digit;
 414     }
 415     if (r) { r->len = a->len; normalize(r); }
 416     return HIGH(remainder);
 417 }
 418
 419 //////////////////////////////////////////////////////////////////////////////
 420 //  Signed multiplication
 421 //////////////////////////////////////////////////////////////////////////////
 422 a_BigInt * mul(a_BigInt * r, const a_BigInt * a, const a_BigInt * b)
 423 {
 424    //////////////////////////////////////////////////////////////
 425    // Check for zeros and ones first
 426    //////////////////////////////////////////////////////////////
 427    if (a->sign == 0 || b->sign == 0) return make_zero(r);
 428    int result_sign = a->sign == 1 ? b->sign : -b->sign;
 429    a_BigInt * result;
 430    if (a->len == 1 && a->digits[0] == 1)      result = copy(r,b);
 431    else if (b->len == 1 && b->digits[0] == 1) result = copy(r,a);
 432    else {
 433       int result_len = a->len + b->len;
 434
 435       //////////////////////////////////////////////////////////////
 436       //  Reorder the multiplicants so that abs(a) < abs(b)
 437       //  The larger multiplicant will be used as the inner
 438       //  loop argument if possible to minimize looping overhead.
 439       //  This is probably insignificant compared to the actual
 440       //  cost of doing the multiplications.
 441       //////////////////////////////////////////////////////////////
 442       if (a->len > b->len) { const a_BigInt * swap = a; a = b; b = swap; }
 443
 444       //////////////////////////////////////////////////////////////
 445       //  Check for duplicated arguments
 446       //////////////////////////////////////////////////////////////
 447       const Digit * x, * x_limit, * y, * y_start, * y_limit;
 448       Digit * z, * z_start;
 449
 450       if (a == r && b == r) {  // r = r * r
 451          //////////////////////////////////////////////////////////////
 452          // This is a special case where the arguments and
 453          // the result are all the same thing.
 454          // We compute the product of the most significant digit
 455          // of the product first and work our way down one by one.
 456          //////////////////////////////////////////////////////////////
 457          result  = realloc(r, result_len);
 458          x_limit = y_limit = a->digits + a->len - 1;
 459          z_start = result->digits + result_len - 2;
 460          clear(result->digits, a->len, result_len);
 461          for ( ; z_start >= result->digits; z_start--) {
 462             Unit sum = 0;
 463             for (x=x_limit, y=y_limit; x <= y_limit; x++, y--) {
 464                Unit product = sum;
 465                product  += (Unit)*x * (Unit)*y;
 466                sum       = LOW(product);
 467                product   = HIGH(product);
 468                for (z = z_start; product; ) {    // Propagate carry
 469                   product  += *++z;
 470                   *z        = LOW(product);
 471                   product   = HIGH(product);
 472                }
 473             }
 474             *z_start = sum;
 475             if (x_limit > a->digits) x_limit--; else y_limit--;
 476          }
 477       } else {
 478          //////////////////////////////////////////////////////////////
 479          // The argument in the inner loop must not be the result.
 480          // However, it is okay to have the outer argument to be the same as
 481          // the result.
 482          //////////////////////////////////////////////////////////////
 483          if (b == r) {     // r = a * r
 484             result  = realloc(r, result_len);
 485             x       = b->digits + b->len - 1;
 486             x_limit = b->digits;
 487             y_start = a->digits;
 488             y_limit = a->digits + a->len;
 489             z_start = result->digits + b->len;
 490             clear(result->digits, b->len, result_len);
 491          } else {                 // r = a * b   or  r = r * b
 492             if (a == r)
 493                result = realloc(r, result_len);
 494             else
 495                result = alloc(r, result_len);
 496             x       = a->digits + a->len - 1;
 497             x_limit = a->digits;
 498             y_start = b->digits;
 499             y_limit = b->digits + b->len;
 500             z_start = result->digits + a->len;
 501             clear(result->digits, a->len, result_len);
 502          }
 503
 504          //////////////////////////////////////////////////////////////
 505          // Now perform unsigned multiplication the traditional way
 506          //////////////////////////////////////////////////////////////
 507          for ( ;x >= x_limit; x--) {
 508             Unit  product;
 509             z = --z_start;
 510             Unit p = *x;
 511             *z = 0;
 512             if (p != 0) {  // Skip unnecessary mulitiplication if possible
 513                for (y = y_start, product = 0; y < y_limit; ) {
 514                   product  += (Unit)*z + (Unit)p * (Unit)*y++;
 515                   *z++      = LOW(product);
 516                   product   = HIGH(product);
 517                }
 518                while (product) { // propagate carry
 519                   product  += *z;
 520                   *z++      = LOW(product);
 521                   product   = HIGH(product);
 522                }
 523             }
 524          }
 525       }
 526       if (result != r) kill(r);
 527       result->len = result_len;
 528       normalize(result);
 529    }
 530    result->sign = result_sign;
 531    return result;
 532 }
 533
 534 //////////////////////////////////////////////////////////////////////////////
 535 //  Signed multiplication with a long
 536 //////////////////////////////////////////////////////////////////////////////
 537 a_BigInt * mul(a_BigInt * r, const a_BigInt * a, long b)
 538 {   if (b == 0)  return make_zero(r);
 539     if (b == 1)  return copy(r,a);
 540     if (b == -1) return neg(r,a);
 541     if (a->sign == 0) return make_zero(r);
 542     if (a->len == 1 && a->digits[0] == 1)  // a == 1 or a == -1
 543        return copy(r,b,a->sign);
 544
 545     struct {
 546        a_BigInt n;
 547        Digit    d;
 548     } c;
 549
 550     Unit mag;
 551     if (b > 0) { c.n.sign = 1; mag = b; }
 552     else       { c.n.sign = -1; mag = -b; }
 553     if ((c.n.digits[1] = HIGH(mag))) {
 554        c.n.len = 2;
 555        c.n.digits[0] = LOW(mag);
 556        return mul(r,a,&c.n);
 557     } else {
 558        a_BigInt * result = alloc(r,a->len+1);
 559        umul(result,a,mag);
 560        result->sign = a->sign == 1 ? c.n.sign : -c.n.sign;
 561        if (result != r) kill(r);
 562        return result;
 563     }
 564 }
 565
 566 //////////////////////////////////////////////////////////////////////////////
 567 //  Unsigned division and modulo arithmetic.
 568 //  Assume a > b.
 569 //////////////////////////////////////////////////////////////////////////////
 570 static a_BigInt * udiv_mod(a_BigInt * r, a_BigInt * a, const a_BigInt * b)
 571 {  Digit * x, * x_hi, * x_lo;
 572    const Digit * y, * y_hi, * y_lo;
 573    Digit * z;
 574
 575    Bool need_quotient = r != 0;
 576    /////////////////////////////////////////////////////////////
 577    //  Set up the indices
 578    /////////////////////////////////////////////////////////////
 579    z    = r->digits + a->len - b->len;
 580    x_hi = a->digits + a->len - 1;
 581    x_lo = a->digits + a->len - b->len;
 582    y_hi = b->digits + b->len - 1;
 583    y_lo = b->digits;
 584
 585    Digit m1, m2, m3;
 586    m1    = y_hi[0];   // The most significant digit of denominator
 587    m2    = b->len == 1 ? 0 : y_hi[-1];  // The next most significant digit
 588    m3    = m1 + 1;
 589
 590    Unit top = 0, quotient_min, quotient_max, quotient;
 591    for (  ;z >= r->digits; z--, x_hi--, x_lo--) {
 592       top = (top << Bits_per_digit) + *x_hi;
 593
 594       /////////////////////////////////////////////////////////////
 595       //  Compute the range of the quotient.  We'll have to very
 596       //  careful to avoid overflow here.
 597       /////////////////////////////////////////////////////////////
 598       quotient_max = (top + 1) / m1;
 599       quotient_min = top / m3;
 600
 601       /////////////////////////////////////////////////////////////
 602       //  Narrow down the quotient value using a binary search
 603       /////////////////////////////////////////////////////////////
 604       Digit x_low = x_hi == a->digits ? 0 : x_hi[-1];
 605       while ( quotient_min < quotient_max ) {
 606          quotient = (quotient_min + quotient_max) / 2;
 607          if (quotient == quotient_min) break;
 608          Digit low; Unit product;
 609          product = quotient * m2;
 610          low     = LOW(product);
 611          product = HIGH(product) + quotient * m1;
 612          if (top > product)       quotient_min = quotient;
 613          else if (top < product)  quotient_max = quotient;
 614          else if (x_low > low)    quotient_min = quotient;
 615          else if (x_low < low)    quotient_max = quotient;
 616          else break;
 617       }
 618
 619       /////////////////////////////////////////////////////////////
 620       //  Now we know the range of the quotient within 1
 621       /////////////////////////////////////////////////////////////
 622       quotient = (quotient_min + quotient_max + 1) / 2;
 623
 624       /////////////////////////////////////////////////////////////
 625       //  Compute a -= quotient * b
 626       /////////////////////////////////////////////////////////////
 627       if (quotient != 0) {
 628          Unit borrow;
 629          for (x = x_lo, y = y_lo, borrow = 0; x <= x_hi; x++, y++) {
 630             borrow  = *x - *y * quotient - borrow;
 631             *x      = LOW(borrow);
 632             borrow  = HIGH(borrow);
 633             if (borrow) borrow = Max_digit - borrow;
 634          }
 635          if (borrow) { // If we have missed by 1, add back one copy
 636             if (x < a->digits + a->len && *x >= borrow) *x -= borrow;
 637             else {
 638                quotient--;
 639                Unit carry;
 640                for (x = x_lo, y = y_lo, carry = 0; x <= x_hi; x++, y++) {
 641                   carry += *x + *y;
 642                   *x     = LOW(carry);
 643                   carry  = HIGH(carry);
 644                }
 645             }
 646          }
 647       }
 648
 649       if (need_quotient) *z = quotient;
 650       top = *x_hi;
 651    }
 652    if (need_quotient) { r->len = a->len - b->len + 1; normalize(r); }
 653    normalize(a);
 654    return r;
 655 }
 656
 657 //////////////////////////////////////////////////////////////////////////////
 658 //  Signed division
 659 //////////////////////////////////////////////////////////////////////////////
 660 a_BigInt * div(a_BigInt * r, const a_BigInt * a, const a_BigInt * b)
 661 {
 662    //////////////////////////////////////////////////////////////
 663    // Check for zeros and ones first
 664    //////////////////////////////////////////////////////////////
 665    if (a->sign == 0) return make_zero(r);
 666    if (b->sign == 0) { BigInt::error_handler("Division by zero");
 667                        return make_zero(r); }
 668    a_BigInt * result;
 669    int result_sign = a->sign ? b->sign : -b->sign;
 670    if (b->len == 1 && b->digits[0] == 1) result = copy(r,a);
 671    else if (b->len == 1) {
 672       //////////////////////////////////////////////////////////////
 673       // fast division for small values
 674       //////////////////////////////////////////////////////////////
 675       result = alloc(r,a->len);
 676       udiv(result,a,b->digits[0]);
 677       if (result != r) kill(r);
 678    } else {
 679       //////////////////////////////////////////////////////////////
 680       // If the magnitude of a is smaller than that of b, the result
 681       // is zero.  If the magnitude is the same, its either 1 or -1
 682       //////////////////////////////////////////////////////////////
 683       int mag_diff = ucmp(a,b);
 684       if (mag_diff == 0) return copy(r,1,result_sign);
 685       if (mag_diff < 0)  return make_zero(r);
 686
 687       //////////////////////////////////////////////////////////////
 688       //  Allocate space for the divisor and result
 689       //////////////////////////////////////////////////////////////
 690       int divisor_len = a->len - b->len + 1;
 691       if (r == a || r == b) {  // Allocate space for the result
 692          result = new (divisor_len) a_BigInt;
 693       } else {
 694          result = alloc(r,divisor_len);
 695       }
 696       a_BigInt * divisor = copy(0,a);
 697       //////////////////////////////////////////////////////////////
 698       // Now perform unsigned division.
 699       //////////////////////////////////////////////////////////////
 700       result = udiv_mod(result,divisor,b);
 701       delete divisor;
 702       if (result != r) kill(r);
 703    }
 704    result->sign = result_sign;
 705    return result;
 706 }
 707
 708 //////////////////////////////////////////////////////////////////////////////
 709 //  Signed division
 710 //////////////////////////////////////////////////////////////////////////////
 711 a_BigInt * div(a_BigInt * r, const a_BigInt * a, long b)
 712 {
 713    //////////////////////////////////////////////////////////////
 714    // Check for zeros and ones first
 715    //////////////////////////////////////////////////////////////
 716    if (a->sign == 0) return make_zero(r);
 717    if (b == 0) { BigInt::error_handler("Division by zero");
 718                  return make_zero(r); }
 719    a_BigInt * result;
 720    int result_sign = a->sign ? b > 0 : b < 0;
 721    if (b == 1 || b == -1) result = copy(r,a);
 722    else {
 723       result = alloc(r,a->len);
 724       Unit mag = abs(b);
 725       if (b >= Max_digit) {
 726          struct {
 727             a_BigInt n;
 728             Digit    d;
 729          } c;
 730          c.n.len = 2;
 731          c.n.digits[0] = LOW(mag);
 732          c.n.digits[1] = HIGH(mag);
 733          a_BigInt * divisor = copy(0,a);
 734          udiv_mod(result,divisor,&c.n);
 735          delete divisor;
 736       } else {
 737          udiv(result,a,mag);
 738       }
 739       if (result != r) kill(r);
 740    }
 741    result->sign = result_sign;
 742    return result;
 743 }
 744
 745 //////////////////////////////////////////////////////////////////////////////
 746 // Signed modulo
 747 //////////////////////////////////////////////////////////////////////////////
 748 a_BigInt * mod(a_BigInt * r, const a_BigInt * a, const a_BigInt * b)
 749 {
 750    //////////////////////////////////////////////////////////////
 751    // Check for zeros and ones first
 752    //////////////////////////////////////////////////////////////
 753    if (a->sign == 0) return make_zero(r);
 754    if (b->sign == 0) { BigInt::error_handler("Modulo by zero");
 755                        return make_zero(r); }
 756    if (b->len == 1 && b->digits[0] == 1) return make_zero(r);
 757    int result_sign = a->sign == 1 ? b->sign : -b->sign;
 758    if (b->len == 1) {
 759       Digit m = udiv(0,a,b->digits[0]);
 760       if (m == 0) return make_zero(r);
 761       else return copy(r,m,result_sign);
 762    }
 763    a_BigInt * result;
 764    int mag_diff = ucmp(a,b);
 765    if (mag_diff == 0) return make_zero(r);
 766    if (mag_diff < 0)  result = copy(r,a);
 767    else {
 768       //////////////////////////////////////////////////////////////
 769       //  Allocate space for the divisor and result
 770       //////////////////////////////////////////////////////////////
 771       if (r == a && r != b) result = copy(r,a);
 772       else                  result = copy(0,a);
 773       udiv_mod(0,result,b);
 774       if (result != r) kill(r);
 775    }
 776    result->sign = result_sign;
 777    return result;
 778 }
 779
 780 //////////////////////////////////////////////////////////////////////////////
 781 //  Signed modulo
 782 //////////////////////////////////////////////////////////////////////////////
 783 a_BigInt * mod(a_BigInt * r, const a_BigInt * a, long b)
 784 {
 785    //////////////////////////////////////////////////////////////
 786    // Check for zeros and ones first
 787    //////////////////////////////////////////////////////////////
 788    if (a->sign == 0) return make_zero(r);
 789    if (b == 0) { BigInt::error_handler("Modulo by zero");
 790                  return make_zero(r); }
 791    a_BigInt * result = 0;
 792    int result_sign = a->sign ? b > 0 : b < 0;
 793    if (b == 1 || b == -1) return make_zero(r);
 794    else {
 795       Unit mag = abs(b);
 796       if (b >= Max_digit) {
 797          struct {
 798             a_BigInt n;
 799             Digit    d;
 800          } c;
 801          if (r != a) result = copy(r,a);
 802          c.n.len = 2;
 803          c.n.digits[0] = LOW(mag);
 804          c.n.digits[1] = HIGH(mag);
 805          udiv_mod(0,result,&c.n);
 806          if (result != r) kill(r);
 807       } else {
 808          Digit mod = udiv(0,a,mag);
 809          if (mod == 0) return make_zero(r);
 810          else return copy(r,mod,result_sign);
 811       }
 812    }
 813    result->sign = result_sign;
 814    return result;
 815 }
 816
 817 //////////////////////////////////////////////////////////////////////////////
 818 //  Bit and/or/xor
 819 //////////////////////////////////////////////////////////////////////////////
 820 a_BigInt * andor(a_BigInt * r, const a_BigInt * a, const a_BigInt * b, char op)
 821 {  a_BigInt * result;
 822    if (a->len < b->len) { const a_BigInt * swap = a; a = b; b = swap; }
 823    int len = op == '&' ? b->len : a->len;
 824    result = alloc(r, len);
 825    result->sign = a->sign >= 0 ? b->sign : -b->sign;
 826    result->len = len;
 827    register const Digit * x, * y, * limit;
 828    register Digit * z;
 829    x = a->digits; y = b->digits; z = result->digits; limit = y + b->len;
 830    switch (op) {
 831       case '&':  while (y < limit) *z++ = *x++ & *y++; break;
 832       case '|':  while (y < limit) *z++ = *x++ | *y++;
 833                  for (limit = x + len; x < limit; ) *z++ = *x++;
 834                  break;
 835       case '^':  while (y < limit) *z++ = *x++ ^ *y++;
 836                  for (limit = x + len; x < limit; ) *z++ = *x++;
 837                  break;
 838    }
 839    normalize(result);
 840    if (result != r) kill(r);
 841    return result;
 842 }
 843
 844 //////////////////////////////////////////////////////////////////////////////
 845 //  Bit and/or/xor with a long
 846 //////////////////////////////////////////////////////////////////////////////
 847 a_BigInt * andor(a_BigInt * r, const a_BigInt * a, long b, char op)
 848 {  struct {
 849       a_BigInt n;
 850       Digit    d;
 851    } c;
 852    c.n.sign      = 0;
 853    c.n.digits[0] = b % Max_digit;
 854    c.n.digits[1] = b / Max_digit;
 855    c.n.len       = 2;
 856    return andor(r,a,&c.n,op);
 857 }
 858
 859 //////////////////////////////////////////////////////////////////////////////
 860 //  Bit one's complement
 861 //////////////////////////////////////////////////////////////////////////////
 862 a_BigInt * not_(a_BigInt * r, const a_BigInt * a)
 863 {  if (a->sign == 0) return make_zero(r);
 864    a_BigInt * result = alloc(r,a->len);
 865    register const Digit * x, * limit;
 866    register Digit * y;
 867    for (x = a->digits, y = result->digits, limit = x + a->len; x < limit; )
 868       *y++ = ~ *x++;
 869    result->len  = a->len;
 870    result->sign = a->sign;
 871    normalize(result);
 872    if (result != r) kill(r);
 873    return result;
 874 }
 875
 876 //////////////////////////////////////////////////////////////////////////////
 877 //  Left/right shift
 878 //////////////////////////////////////////////////////////////////////////////
 879 a_BigInt * shift(a_BigInt * r, const a_BigInt * a, const a_BigInt * b, int sign)
 880 {  if (b->sign == 0) return copy(r,a);   // a shift 0 == a
 881    int len = b->len == 1 ? b->digits[0] : b->digits[0] + b->digits[1] * Max_digit;
 882    if (b->sign == -1) len = -len;
 883    return shift(r,a,len,sign);
 884 }
 885
 886 //////////////////////////////////////////////////////////////////////////////
 887 //  Left/right shift with a long
 888 //////////////////////////////////////////////////////////////////////////////
 889 a_BigInt * shift(a_BigInt * r, const a_BigInt * a, long len, int sign)
 890 {  if (a->sign == 0) return make_zero(r);
 891    if (sign < 0) len = -len;
 892    int result_len = a->len + (len + Bits_per_digit + 1) / Bits_per_digit;
 893    a_BigInt * result = alloc(r,result_len);
 894    if (result != r) kill(r);
 895
 896
 897    return result;
 898 }
 899
 900 //////////////////////////////////////////////////////////////////////////////
 901 //  Type conversion to a double
 902 //////////////////////////////////////////////////////////////////////////////
 903 BigInt::operator double () const
 904 {  register double r = 0.0;
 905    register Digit * p, * limit;
 906    for (p = D->digits + D->len - 1, limit = D->digits; p >= limit; )
 907       r = r * Max_digit + *p--;
 908    if (D->sign == -1) r = -r;
 909    return r;
 910 }
 911
 912 //////////////////////////////////////////////////////////////////////////////
 913 //  Return the length (bit number of highest order) in the number
 914 //////////////////////////////////////////////////////////////////////////////
 915 int BigInt::length() const
 916 {  int l;
 917    for (l = D->len * sizeof(Digit) * 8; l > 0; l--)
 918       if ((*this)[l] != 0) break;
 919    return l;
 920 }
 921
 922 //////////////////////////////////////////////////////////////////////////////
 923 //  Miscellaneous
 924 //////////////////////////////////////////////////////////////////////////////
 925 void div_mod(const BigInt&, const BigInt&, BigInt&, BigInt&)
 926 {
 927   /* unimplemented */
 928 }
 929
 930 BigInt pow (const BigInt& a, const BigInt& b)
 931 {  if (a.sign() == 0) return a;
 932    if (b.sign() == 0) return 1;
 933    if (b.sign() < 0) { BigInt::error_handler("negative pow() unsupported");
 934                        return 0; }
 935    BigInt n;
 936    int i;
 937    for (n = 1, i = b.length(); i > 0; i--) {
 938       n *= n; if (b[i]) n *= a;
 939    }
 940    return n;
 941 }
 942
 943
 944 BigInt pow_mod (const BigInt& a, const BigInt& b, const BigInt& c)
 945 {  if (a.sign() == 0) return a;
 946    if (b.sign() == 0) return 1;
 947    if (b.sign() < 0) { BigInt::error_handler("negative pow_mod() unsupported");
 948                        return 0; }
 949    BigInt n;
 950    int i;
 951    for (n = 1, i = b.length(); i > 0; i--) {
 952       n *= n; if (b[i]) n *= a; n %= c;
 953    }
 954    return n;
 955 }
 956
 957 BigInt gcd (BigInt a, BigInt b)
 958 {  for (;;) {
 959       if (a == (int)1) return b;
 960       if (b == (int)1) return a;
 961       if (a > b) a %= b; else b %= a;
 962    }
 963 }
 964
 965 //////////////////////////////////////////////////////////////////////////////
 966 //  Hashing
 967 //////////////////////////////////////////////////////////////////////////////
 968 unsigned int hash(const BigInt& a)
 969 {  unsigned int h = 0;
 970    for (int i = a.D->len - 1; i >= 0; i--)
 971       h += h + a.D->digits[i];
 972    return h;
 973 }
 974
 975 //////////////////////////////////////////////////////////////////////////////
 976 //  String conversion
 977 //////////////////////////////////////////////////////////////////////////////
 978 char * BigInt::makeString(char buf[], size_t length, unsigned int base) const
 979 {  char * cursor = buf;
 980    char * limit  = buf + length;
 981
 982    if (D->sign == 0) {       // Check for zero
 983       if (length >= 2) {
 984          buf[0] = '0'; buf[1] = '\0'; return buf+2;
 985       } else return buf;
 986    }
 987
 988    a_BigInt * a = copy(0,D);
 989    while (a->len > 0) {
 990       int digit = udiv(a,a,base);
 991       digit += digit >= 10 ? 'a' - 10 : '0';
 992       if (cursor == limit) break;
 993       *cursor++ = digit;
 994    }
 995    kill(a);
 996
 997    if (D->sign == -1 && cursor < limit) {    // Insert sign
 998       *cursor++ = '-';
 999    }
1000
1001    //////////////////////////////////////////////////////
1002    // Now swap characters within the buffer
1003    //////////////////////////////////////////////////////
1004    register char * p, * q;
1005    for (p = buf, q = cursor - 1; p < q; p++, q--) {
1006       char c = *p; *p = *q; *q = c;
1007    }
1008    if (cursor < limit) *cursor = '\0';
1009    return cursor;
1010 }
1011
1012 //////////////////////////////////////////////////////////////////////////////
1013 // Print a BigInt
1014 //////////////////////////////////////////////////////////////////////////////
1015 std::ostream& operator << (std::ostream& out, const BigInt& n)
1016 {  char buffer[256];
1017    char * buf;
1018    size_t len = n.D->len * 5 + 3;
1019    if (len > sizeof(buffer)) buf = new char [len];
1020    else buf = buffer;
1021    n.makeString(buf,len,10);
1022    out << buf;
1023    if (buf != buffer) delete [] buf;
1024    return out;
1025 }
1026
1027 //////////////////////////////////////////////////////////////////////////////
1028 // Construct a BigInt
1029 //////////////////////////////////////////////////////////////////////////////
1030 BigInt::BigInt(const char * number) : D(0)
1031 {
1032    operator = (number);
1033 }
1034
1035 //////////////////////////////////////////////////////////////////////////////
1036 //  Make a BigInt from a string
1037 //////////////////////////////////////////////////////////////////////////////
1038 BigInt& BigInt::operator = (const char * number)
1039 {  parseString(number);
1040    return *this;
1041 }
1042
1043 //////////////////////////////////////////////////////////////////////////////
1044 //  Make a BigInt from a string
1045 //////////////////////////////////////////////////////////////////////////////
1046 int BigInt::parseString(const char * number, unsigned int base)
1047 {  BigInt n = 0;
1048    Bool negative = false;
1049    const char * p;
1050    char c;
1051
1052    for (p = number; (c = *p); p++)
1053    {  if (c == '-') negative = ! negative;
1054       else break;
1055    }
1056
1057    for (  ; (c = *p); p++)
1058    {  unsigned int digit;
1059       if (isdigit(c)) digit = c - '0';
1060       else if (islower(c)) digit = (c - 'a') + 10;
1061       else if (isupper(c)) digit = (c - 'A') + 10;
1062       else if (isspace(c)) continue;
1063       else break;
1064       if (digit >= base) break;
1065       n = n * 10 + (long)digit;
1066    }
1067    if (negative) n = -n;
1068    *this = n;
1069    return number - p;
1070 }