BigUnsigned.cc

   1 /*
   2 * Matt McCutchen's Big Integer Library
   3 * http://mysite.verizon.net/mccutchen/bigint/
   4 */
   5
   6 #include "BigUnsigned.hh"
   7
   8 // The "management" routines that used to be here are now in NumberlikeArray.cpp.
   9
  10 /*
  11 * The steps for construction of a BigUnsigned
  12 * from an integral value x are as follows:
  13 * 1. If x is zero, create an empty BigUnsigned and stop.
  14 * 2. If x is negative, throw an exception.
  15 * 3. Allocate a one-block number array.
  16 * 4. If x is of a signed type, convert x to the unsigned
  17 *    type of the same length.
  18 * 5. Expand x to a Blk, and store it in the number array.
  19 */
  20
  21 BigUnsigned::BigUnsigned(unsigned long x) {
  22         if (x == 0) {
  23                 cap = 0;
  24                 blk = new Blk[0];
  25                 len = 0;
  26         } else {
  27                 cap = 1;
  28                 blk = new Blk[1];
  29                 len = 1;
  30                 blk[0] = Blk(x);
  31         }
  32 }
  33
  34 BigUnsigned::BigUnsigned(long x) {
  35         if (x == 0) {
  36                 cap = 0;
  37                 blk = new Blk[0];
  38                 len = 0;
  39         } else if (x > 0) {
  40                 cap = 1;
  41                 blk = new Blk[1];
  42                 len = 1;
  43                 blk[0] = Blk(x);
  44         } else
  45     throw "BigUnsigned::BigUnsigned(long): Cannot construct a BigUnsigned from a negative number";
  46 }
  47
  48 BigUnsigned::BigUnsigned(unsigned int x) {
  49         if (x == 0) {
  50                 cap = 0;
  51                 blk = new Blk[0];
  52                 len = 0;
  53         } else {
  54                 cap = 1;
  55                 blk = new Blk[1];
  56                 len = 1;
  57                 blk[0] = Blk(x);
  58         }
  59 }
  60
  61 BigUnsigned::BigUnsigned(int x) {
  62         if (x == 0) {
  63                 cap = 0;
  64                 blk = new Blk[0];
  65                 len = 0;
  66         } else if (x > 0) {
  67                 cap = 1;
  68                 blk = new Blk[1];
  69                 len = 1;
  70                 blk[0] = Blk(x);
  71         } else
  72     throw "BigUnsigned::BigUnsigned(int): Cannot construct a BigUnsigned from a negative number";
  73 }
  74
  75 BigUnsigned::BigUnsigned(unsigned short x) {
  76         if (x == 0) {
  77                 cap = 0;
  78                 blk = new Blk[0];
  79                 len = 0;
  80         } else {
  81                 cap = 1;
  82                 blk = new Blk[1];
  83                 len = 1;
  84                 blk[0] = Blk(x);
  85         }
  86 }
  87
  88 BigUnsigned::BigUnsigned(short x) {
  89         if (x == 0) {
  90                 cap = 0;
  91                 blk = new Blk[0];
  92                 len = 0;
  93         } else if (x > 0) {
  94                 cap = 1;
  95                 blk = new Blk[1];
  96                 len = 1;
  97                 blk[0] = Blk(x);
  98         } else
  99     throw "BigUnsigned::BigUnsigned(short): Cannot construct a BigUnsigned from a negative number";
 100 }
 101
 102 // CONVERTERS
 103 /*
 104 * The steps for conversion of a BigUnsigned to an
 105 * integral type are as follows:
 106 * 1. If the BigUnsigned is zero, return zero.
 107 * 2. If it is more than one block long or its lowest
 108 *    block has bits set out of the range of the target
 109 *    type, throw an exception.
 110 * 3. Otherwise, convert the lowest block to the
 111 *    target type and return it.
 112 */
 113
 114 namespace {
 115         // These masks are used to test whether a Blk has bits
 116         // set out of the range of a smaller integral type.  Note
 117         // that this range is not considered to include the sign bit.
 118         const BigUnsigned::Blk  lMask = ~0 >> 1;
 119         const BigUnsigned::Blk uiMask = (unsigned int)(~0);
 120         const BigUnsigned::Blk  iMask = uiMask >> 1;
 121         const BigUnsigned::Blk usMask = (unsigned short)(~0);
 122         const BigUnsigned::Blk  sMask = usMask >> 1;
 123 }
 124
 125 BigUnsigned::operator unsigned long() const {
 126         if (len == 0)
 127                 return 0;
 128         else if (len == 1)
 129                 return (unsigned long) blk[0];
 130         else
 131                 throw "BigUnsigned::operator unsigned long: Value is too big for an unsigned long";
 132 }
 133
 134 BigUnsigned::operator long() const {
 135         if (len == 0)
 136                 return 0;
 137         else if (len == 1 && (blk[0] & lMask) == blk[0])
 138                 return (long) blk[0];
 139         else
 140                 throw "BigUnsigned::operator long: Value is too big for a long";
 141 }
 142
 143 BigUnsigned::operator unsigned int() const {
 144         if (len == 0)
 145                 return 0;
 146         else if (len == 1 && (blk[0] & uiMask) == blk[0])
 147                 return (unsigned int) blk[0];
 148         else
 149                 throw "BigUnsigned::operator unsigned int: Value is too big for an unsigned int";
 150 }
 151
 152 BigUnsigned::operator int() const {
 153         if (len == 0)
 154                 return 0;
 155         else if (len == 1 && (blk[0] & iMask) == blk[0])
 156                 return (int) blk[0];
 157         else
 158                 throw "BigUnsigned::operator int: Value is too big for an int";
 159 }
 160
 161 BigUnsigned::operator unsigned short() const {
 162         if (len == 0)
 163                 return 0;
 164         else if (len == 1 && (blk[0] & usMask) == blk[0])
 165                 return (unsigned short) blk[0];
 166         else
 167                 throw "BigUnsigned::operator unsigned short: Value is too big for an unsigned short";
 168 }
 169
 170 BigUnsigned::operator short() const {
 171         if (len == 0)
 172                 return 0;
 173         else if (len == 1 && (blk[0] & sMask) == blk[0])
 174                 return (short) blk[0];
 175         else
 176                 throw "BigUnsigned::operator short: Value is too big for a short";
 177 }
 178
 179 // COMPARISON
 180 BigUnsigned::CmpRes BigUnsigned::compareTo(const BigUnsigned &x) const {
 181         // A bigger length implies a bigger number.
 182         if (len < x.len)
 183                 return less;
 184         else if (len > x.len)
 185                 return greater;
 186         else {
 187                 // Compare blocks one by one from left to right.
 188                 Index i = len;
 189                 while (i > 0) {
 190                         i--;
 191                         if (blk[i] == x.blk[i])
 192                                 continue;
 193                         else if (blk[i] > x.blk[i])
 194                                 return greater;
 195                         else
 196                                 return less;
 197                 }
 198                 // If no blocks differed, the numbers are equal.
 199                 return equal;
 200         }
 201 }
 202
 203 // PUT-HERE OPERATIONS
 204
 205 // Addition
 206 void BigUnsigned::add(const BigUnsigned &a, const BigUnsigned &b) {
 207         // Block unsafe calls
 208         if (this == &a || this == &b)
 209                 throw "BigUnsigned::add: One of the arguments is the invoked object";
 210         // If one argument is zero, copy the other.
 211         if (a.len == 0) {
 212                 operator =(b);
 213                 return;
 214         } else if (b.len == 0) {
 215                 operator =(a);
 216                 return;
 217         }
 218         // Carries in and out of an addition stage
 219         bool carryIn, carryOut;
 220         Blk temp;
 221         Index i;
 222         // a2 points to the longer input, b2 points to the shorter
 223         const BigUnsigned *a2, *b2;
 224         if (a.len >= b.len) {
 225                 a2 = &a;
 226                 b2 = &b;
 227         } else {
 228                 a2 = &b;
 229                 b2 = &a;
 230         }
 231         // Set prelimiary length and make room in this BigUnsigned
 232         len = a2->len + 1;
 233         allocate(len);
 234         // For each block index that is present in both inputs...
 235         for (i = 0, carryIn = false; i < b2->len; i++) {
 236                 // Add input blocks
 237                 temp = a2->blk[i] + b2->blk[i];
 238                 // If a rollover occurred, the result is less than either input.
 239                 // This test is used many times in the BigUnsigned code.
 240                 carryOut = (temp < a2->blk[i]);
 241                 // If a carry was input, handle it
 242                 if (carryIn) {
 243                         temp++;
 244                         carryOut |= (temp == 0);
 245                 }
 246                 blk[i] = temp; // Save the addition result
 247                 carryIn = carryOut; // Pass the carry along
 248         }
 249         // If there is a carry left over, increase blocks until
 250         // one does not roll over.
 251         for (; i < a2->len && carryIn; i++) {
 252                 temp = a2->blk[i] + 1;
 253                 carryIn = (temp == 0);
 254                 blk[i] = temp;
 255         }
 256         // If the carry was resolved but the larger number
 257         // still has blocks, copy them over.
 258         for (; i < a2->len; i++)
 259                 blk[i] = a2->blk[i];
 260         // Set the extra block if there's still a carry, decrease length otherwise
 261         if (carryIn)
 262                 blk[i] = 1;
 263         else
 264                 len--;
 265 }
 266
 267 // Subtraction
 268 void BigUnsigned::subtract(const BigUnsigned &a, const BigUnsigned &b) {
 269         // Block unsafe calls
 270         if (this == &a || this == &b)
 271                 throw "BigUnsigned::subtract: One of the arguments is the invoked object";
 272         // If b is zero, copy a.  If a is shorter than b, the result is negative.
 273         if (b.len == 0) {
 274                 operator =(a);
 275                 return;
 276         } else if (a.len < b.len)
 277         throw "BigUnsigned::subtract: Negative result in unsigned calculation";
 278         bool borrowIn, borrowOut;
 279         Blk temp;
 280         Index i;
 281         // Set preliminary length and make room
 282         len = a.len;
 283         allocate(len);
 284         // For each block index that is present in both inputs...
 285         for (i = 0, borrowIn = false; i < b.len; i++) {
 286                 temp = a.blk[i] - b.blk[i];
 287                 // If a reverse rollover occurred, the result is greater than the block from a.
 288                 borrowOut = (temp > a.blk[i]);
 289                 // Handle an incoming borrow
 290                 if (borrowIn) {
 291                         borrowOut |= (temp == 0);
 292                         temp--;
 293                 }
 294                 blk[i] = temp; // Save the subtraction result
 295                 borrowIn = borrowOut; // Pass the borrow along
 296         }
 297         // If there is a borrow left over, decrease blocks until
 298         // one does not reverse rollover.
 299         for (; i < a.len && borrowIn; i++) {
 300                 borrowIn = (a.blk[i] == 0);
 301                 blk[i] = a.blk[i] - 1;
 302         }
 303         // If there's still a borrow, the result is negative.
 304         // Throw an exception, but zero out this object first just in case.
 305         if (borrowIn) {
 306                 len = 0;
 307                 throw "BigUnsigned::subtract: Negative result in unsigned calculation";
 308         } else // Copy over the rest of the blocks
 309         for (; i < a.len; i++)
 310                 blk[i] = a.blk[i];
 311         // Zap leading zeros
 312         zapLeadingZeros();
 313 }
 314
 315 // Multiplication
 316 void BigUnsigned::multiply(const BigUnsigned &a, const BigUnsigned &b) {
 317         // Block unsafe calls
 318         if (this == &a || this == &b)
 319                 throw "BigUnsigned::multiply: One of the arguments is the invoked object";
 320         // If either a or b is zero, set to zero.
 321         if (a.len == 0 || b.len == 0) {
 322                 len = 0;
 323                 return;
 324         }
 325         // Overall method: this = 0, then for each 1-bit of a, add b
 326         // to this shifted the appropriate amount.
 327         // Variables for the calculation
 328         Index i, j, k;
 329         unsigned int i2;
 330         Blk aBlk, bHigh, temp;
 331         bool carryIn, carryOut;
 332         // Set preliminary length and make room
 333         len = a.len + b.len;
 334         allocate(len);
 335         // Zero out this object
 336         for (i = 0; i < len; i++)
 337                 blk[i] = 0;
 338         // For each block of the first number...
 339         for (i = 0; i < a.len; i++) {
 340                 // For each 1-bit of that block...
 341                 for (i2 = 0, aBlk = a.blk[i]; aBlk != 0; i2++, aBlk >>= 1) {
 342                         if ((aBlk & 1) == 0)
 343                                 continue;
 344                         /* Add b to this, shifted left i blocks and i2 bits.
 345                         * j is the index in b, and k = i + j is the index in this.
 346                         * The low bits of b.blk[j] are shifted and added to blk[k].
 347                         * bHigh is used to carry the high bits to the next addition. */
 348                         bHigh = 0;
 349                         for (j = 0, k = i, carryIn = false; j < b.len; j++, k++) {
 350                                 temp = blk[k] + ((b.blk[j] << i2) | bHigh);
 351                                 carryOut = (temp < blk[k]);
 352                                 if (carryIn) {
 353                                         temp++;
 354                                         carryOut |= (temp == 0);
 355                                 }
 356                                 blk[k] = temp;
 357                                 carryIn = carryOut;
 358                                 bHigh = (i2 == 0) ? 0 : b.blk[j] >> (8 * sizeof(Blk) - i2);
 359                         }
 360                         temp = blk[k] + bHigh;
 361                         carryOut = (temp < blk[k]);
 362                         if (carryIn) {
 363                                 temp++;
 364                                 carryOut |= (temp == 0);
 365                         }
 366                         blk[k] = temp;
 367                         carryIn = carryOut;
 368                         k++; // Added by Matt 2004.12.23: Move to the next block.  It belongs here (and there was a corresponding line in the division routine), but I'm not certain whether it ever matters.
 369                         for (; carryIn; k++) {
 370                                 blk[k]++;
 371                                 carryIn = (blk[k] == 0);
 372                         }
 373                 }
 374         }
 375         // Zap possible leading zero
 376         if (blk[len - 1] == 0)
 377                 len--;
 378 }
 379
 380 /*
 381 * DIVISION WITH REMAINDER
 382 * The functionality of divide, modulo, and %= is included in this one monstrous call,
 383 * which deserves some explanation.
 384 *
 385 * The division *this / b is performed.
 386 * Afterwards, q has the quotient, and *this has the remainder.
 387 * Thus, a call is like q = *this / b, *this %= b.
 388 *
 389 * This seemingly bizarre pattern of inputs and outputs has a justification.  The
 390 * ``put-here operations'' are supposed to be fast.  Therefore, they accept inputs
 391 * and provide outputs in the most convenient places so that no value ever needs
 392 * to be copied in its entirety.  That way, the client can perform exactly the
 393 * copying it needs depending on where the inputs are and where it wants the output.
 394 */
 395 void BigUnsigned::divideWithRemainder(const BigUnsigned &b, BigUnsigned &q) {
 396         // Block unsafe calls
 397         if (this == &b || &q == &b || this == &q)
 398                 throw "BigUnsigned::divideWithRemainder: Some two objects involved are the same";
 399
 400         /*
 401         * Note that the mathematical definition of mod (I'm trusting Knuth) is somewhat
 402         * different from the way the normal C++ % operator behaves in the case of division by 0.
 403         * This function does it Knuth's way.
 404         *
 405         * We let a / 0 == 0 (it doesn't matter) and a % 0 == a, no exceptions thrown.
 406         * This allows us to preserve both Knuth's demand that a mod 0 == a
 407         * and the useful property that (a / b) * b + (a % b) == a.
 408         */
 409         if (b.len == 0) {
 410                 q.len = 0;
 411                 return;
 412         }
 413
 414         /*
 415         * If *this.len < b.len, then *this < b, and we can be sure that b doesn't go into
 416         * *this at all.  The quotient is 0 and *this is already the remainder (so leave it alone).
 417         */
 418         if (len < b.len) {
 419                 q.len = 0;
 420                 return;
 421         }
 422
 423         /*
 424         * At this point we know *this > b > 0.  (Whew!)
 425         */
 426
 427         /* DEBUG *
 428         std::cout << "divideWithRemainder starting\n"
 429         << "length of dividend: " << len
 430         << "\nlast block of dividend: " << getBlock(0)
 431         << "\nlength of divisor: " << b.len
 432         << "\nlast block of divisor: " << b.getBlock(0)
 433         << std::endl; */
 434
 435         /*
 436         * Overall method: Subtract b, shifted varying amounts to
 437         * the left, from this, setting the bit in the quotient q
 438         * whenever the subtraction succeeds.  Eventually q will contain the entire
 439         * quotient, and this will be left with the remainder.
 440         *
 441         * We use work2 to temporarily store the result of a subtraction.
 442         * But we don't even compute the i lowest blocks of the result,
 443         * because they are unaffected (we shift left i places).
 444         * */
 445         // Variables for the calculation
 446         Index i, j, k;
 447         unsigned int i2;
 448         Blk bHigh, temp;
 449         bool borrowIn, borrowOut;
 450
 451         // Make sure we have an extra zero block just past the value,
 452         // but don't increase the logical length.  A shifted subtraction
 453         // (for example, subtracting 1 << 2 from 4) might stick into
 454         // this block.
 455         allocateAndCopy(len + 1);
 456         blk[len] = 0;
 457
 458         // work2 holds part of the result of a subtraction.
 459         // (There's no work1.  The name work2 is from a previous version.)
 460         Blk *work2 = new Blk[len];
 461
 462         // Set preliminary length for quotient and make room
 463         q.len = len - b.len + 1;
 464         q.allocate(q.len);
 465         // Zero out the quotient
 466         for (i = 0; i < q.len; i++)
 467                 q.blk[i] = 0;
 468
 469         // For each possible left-shift of b in blocks...
 470         i = q.len;
 471         while (i > 0) {
 472                 i--;
 473                 // For each possible left-shift of b in bits...
 474                 q.blk[i] = 0;
 475                 i2 = 8 * sizeof(Blk);
 476                 while (i2 > 0) {
 477                         i2--;
 478                         /*
 479                         * Subtract b, shifted left i blocks and i2 bits, from this.
 480                         * and store the answer in work2.
 481                         *
 482                         * Compare this to the middle section of `multiply'.  They
 483                         * are in many ways analogous.
 484                         */
 485                         bHigh = 0;
 486                         for (j = 0, k = i, borrowIn = false; j < b.len; j++, k++) {
 487                                 temp = blk[k] - ((b.blk[j] << i2) | bHigh);
 488                                 borrowOut = (temp > blk[k]);
 489                                 if (borrowIn) {
 490                                         borrowOut |= (temp == 0);
 491                                         temp--;
 492                                 }
 493                                 work2[j] = temp;
 494                                 borrowIn = borrowOut;
 495                                 bHigh = (i2 == 0) ? 0 : b.blk[j] >> (8 * sizeof(Blk) - i2);
 496                         }
 497                         temp = blk[k] - bHigh;
 498                         borrowOut = (temp > blk[k]);
 499                         if (borrowIn) {
 500                                 borrowOut |= (temp == 0);
 501                                 temp--;
 502                         }
 503                         work2[j] = temp;
 504                         borrowIn = borrowOut;
 505                         j++;
 506                         k++;
 507                         for (; k < len && borrowIn; j++, k++) {
 508                                 borrowIn = (blk[k] == 0);
 509                                 work2[j] = blk[k] - 1;
 510                         }
 511                         /* If the subtraction was performed successfully (!borrowIn), set bit i2
 512                         * in block i of the quotient, and copy the changed portion of
 513                         * work2 back to this. Otherwise, reset that bit and move on. */
 514                         if (!borrowIn) {
 515                                 q.blk[i] |= (1 << i2);
 516                                 while (j > 0) {
 517                                         j--;
 518                                         k--;
 519                                         blk[k] = work2[j];
 520                                 }
 521                         }
 522                 }
 523         }
 524         // Zap possible leading zero in quotient
 525         if (q.blk[q.len - 1] == 0)
 526                 q.len--;
 527         // Zap any/all leading zeros in remainder
 528         zapLeadingZeros();
 529         // Deallocate temporary array.
 530         // (Thanks to Brad Spencer for noticing my accidental omission of this!)
 531         delete [] work2;
 532
 533         /* DEBUG *
 534         std::cout << "divideWithRemainder complete\n"
 535         << "length of quotient: " << q.len
 536         << "\nlast block of quotient: " << q.getBlock(0)
 537         << "\nlength of remainder: " << len
 538         << "\nlast block of remainder: " << getBlock(0)
 539         << std::endl; */
 540 }
 541
 542 // Bitwise and
 543 void BigUnsigned::bitAnd(const BigUnsigned &a, const BigUnsigned &b) {
 544         // Block unsafe calls
 545         if (this == &a || this == &b)
 546                 throw "BigUnsigned::bitAnd: One of the arguments is the invoked object";
 547         len = (a.len >= b.len) ? b.len : a.len;
 548         allocate(len);
 549         Index i;
 550         for (i = 0; i < len; i++)
 551                 blk[i] = a.blk[i] & b.blk[i];
 552         zapLeadingZeros();
 553 }
 554
 555 // Bitwise or
 556 void BigUnsigned::bitOr(const BigUnsigned &a, const BigUnsigned &b) {
 557         // Block unsafe calls
 558         if (this == &a || this == &b)
 559                 throw "BigUnsigned::bitOr: One of the arguments is the invoked object";
 560         Index i;
 561         const BigUnsigned *a2, *b2;
 562         if (a.len >= b.len) {
 563                 a2 = &a;
 564                 b2 = &b;
 565         } else {
 566                 a2 = &b;
 567                 b2 = &a;
 568         }
 569         allocate(a2->len);
 570         for (i = 0; i < b2->len; i++)
 571                 blk[i] = a2->blk[i] | b2->blk[i];
 572         for (; i < a2->len; i++)
 573                 blk[i] = a2->blk[i];
 574         len = a2->len;
 575 }
 576
 577 // Bitwise xor
 578 void BigUnsigned::bitXor(const BigUnsigned &a, const BigUnsigned &b) {
 579         // Block unsafe calls
 580         if (this == &a || this == &b)
 581                 throw "BigUnsigned::bitXor: One of the arguments is the invoked object";
 582         Index i;
 583         const BigUnsigned *a2, *b2;
 584         if (a.len >= b.len) {
 585                 a2 = &a;
 586                 b2 = &b;
 587         } else {
 588                 a2 = &b;
 589                 b2 = &a;
 590         }
 591         allocate(b2->len);
 592         for (i = 0; i < b2->len; i++)
 593                 blk[i] = a2->blk[i] ^ b2->blk[i];
 594         for (; i < a2->len; i++)
 595                 blk[i] = a2->blk[i];
 596         len = a2->len;
 597         zapLeadingZeros();
 598 }
 599
 600 // INCREMENT/DECREMENT OPERATORS
 601
 602 // Prefix increment
 603 void BigUnsigned::operator ++() {
 604         Index i;
 605         bool carry = true;
 606         for (i = 0; i < len && carry; i++) {
 607                 blk[i]++;
 608                 carry = (blk[i] == 0);
 609         }
 610         if (carry) {
 611                 // Matt fixed a bug 2004.12.24: next 2 lines used to say allocateAndCopy(len + 1)
 612                 len++;
 613                 allocateAndCopy(len);
 614                 blk[i] = 1;
 615         }
 616 }
 617
 618 // Postfix increment: same as prefix
 619 void BigUnsigned::operator ++(int) {
 620         operator ++();
 621 }
 622
 623 // Prefix decrement
 624 void BigUnsigned::operator --() {
 625         if (len == 0)
 626                 throw "BigUnsigned::operator --(): Cannot decrement an unsigned zero";
 627         Index i;
 628         bool borrow = true;
 629         for (i = 0; borrow; i++) {
 630                 borrow = (blk[i] == 0);
 631                 blk[i]--;
 632         }
 633         // Zap possible leading zero (there can only be one)
 634         if (blk[len - 1] == 0)
 635                 len--;
 636 }
 637
 638 // Postfix decrement: same as prefix
 639 void BigUnsigned::operator --(int) {
 640         operator --();
 641 }