| 1 | /* |
| 2 | * Matt McCutchen's Big Integer Library |
| 3 | */ |
| 4 | |
| 5 | #include "BigUnsigned.hh" |
| 6 | |
| 7 | // The "management" routines that used to be here are now in NumberlikeArray.hh. |
| 8 | |
| 9 | /* |
| 10 | * The steps for construction of a BigUnsigned |
| 11 | * from an integral value x are as follows: |
| 12 | * 1. If x is zero, create an empty BigUnsigned and stop. |
| 13 | * 2. If x is negative, throw an exception. |
| 14 | * 3. Allocate a one-block number array. |
| 15 | * 4. If x is of a signed type, convert x to the unsigned |
| 16 | * type of the same length. |
| 17 | * 5. Expand x to a Blk, and store it in the number array. |
| 18 | * |
| 19 | * Since 2005.01.06, NumberlikeArray uses `NULL' rather |
| 20 | * than a real array if one of zero length is needed. |
| 21 | * These constructors implicitly call NumberlikeArray's |
| 22 | * default constructor, which sets `blk = NULL, cap = len = 0'. |
| 23 | * So if the input number is zero, they can just return. |
| 24 | * See remarks in `NumberlikeArray.hh'. |
| 25 | */ |
| 26 | |
| 27 | BigUnsigned::BigUnsigned(unsigned long x) { |
| 28 | if (x == 0) |
| 29 | ; // NumberlikeArray already did all the work |
| 30 | else { |
| 31 | cap = 1; |
| 32 | blk = new Blk[1]; |
| 33 | len = 1; |
| 34 | blk[0] = Blk(x); |
| 35 | } |
| 36 | } |
| 37 | |
| 38 | BigUnsigned::BigUnsigned(long x) { |
| 39 | if (x == 0) |
| 40 | ; |
| 41 | else if (x > 0) { |
| 42 | cap = 1; |
| 43 | blk = new Blk[1]; |
| 44 | len = 1; |
| 45 | blk[0] = Blk(x); |
| 46 | } else |
| 47 | throw "BigUnsigned::BigUnsigned(long): Cannot construct a BigUnsigned from a negative number"; |
| 48 | } |
| 49 | |
| 50 | BigUnsigned::BigUnsigned(unsigned int x) { |
| 51 | if (x == 0) |
| 52 | ; |
| 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 | ; |
| 64 | else if (x > 0) { |
| 65 | cap = 1; |
| 66 | blk = new Blk[1]; |
| 67 | len = 1; |
| 68 | blk[0] = Blk(x); |
| 69 | } else |
| 70 | throw "BigUnsigned::BigUnsigned(int): Cannot construct a BigUnsigned from a negative number"; |
| 71 | } |
| 72 | |
| 73 | BigUnsigned::BigUnsigned(unsigned short x) { |
| 74 | if (x == 0) |
| 75 | ; |
| 76 | else { |
| 77 | cap = 1; |
| 78 | blk = new Blk[1]; |
| 79 | len = 1; |
| 80 | blk[0] = Blk(x); |
| 81 | } |
| 82 | } |
| 83 | |
| 84 | BigUnsigned::BigUnsigned(short x) { |
| 85 | if (x == 0) |
| 86 | ; |
| 87 | else if (x > 0) { |
| 88 | cap = 1; |
| 89 | blk = new Blk[1]; |
| 90 | len = 1; |
| 91 | blk[0] = Blk(x); |
| 92 | } else |
| 93 | throw "BigUnsigned::BigUnsigned(short): Cannot construct a BigUnsigned from a negative number"; |
| 94 | } |
| 95 | |
| 96 | // CONVERTERS |
| 97 | /* |
| 98 | * The steps for conversion of a BigUnsigned to an |
| 99 | * integral type are as follows: |
| 100 | * 1. If the BigUnsigned is zero, return zero. |
| 101 | * 2. If it is more than one block long or its lowest |
| 102 | * block has bits set out of the range of the target |
| 103 | * type, throw an exception. |
| 104 | * 3. Otherwise, convert the lowest block to the |
| 105 | * target type and return it. |
| 106 | */ |
| 107 | |
| 108 | namespace { |
| 109 | // These masks are used to test whether a Blk has bits |
| 110 | // set out of the range of a smaller integral type. Note |
| 111 | // that this range is not considered to include the sign bit. |
| 112 | const BigUnsigned::Blk lMask = ~0 >> 1; |
| 113 | const BigUnsigned::Blk uiMask = (unsigned int)(~0); |
| 114 | const BigUnsigned::Blk iMask = uiMask >> 1; |
| 115 | const BigUnsigned::Blk usMask = (unsigned short)(~0); |
| 116 | const BigUnsigned::Blk sMask = usMask >> 1; |
| 117 | } |
| 118 | |
| 119 | BigUnsigned::operator unsigned long() const { |
| 120 | if (len == 0) |
| 121 | return 0; |
| 122 | else if (len == 1) |
| 123 | return (unsigned long) blk[0]; |
| 124 | else |
| 125 | throw "BigUnsigned::operator unsigned long: Value is too big for an unsigned long"; |
| 126 | } |
| 127 | |
| 128 | BigUnsigned::operator long() const { |
| 129 | if (len == 0) |
| 130 | return 0; |
| 131 | else if (len == 1 && (blk[0] & lMask) == blk[0]) |
| 132 | return (long) blk[0]; |
| 133 | else |
| 134 | throw "BigUnsigned::operator long: Value is too big for a long"; |
| 135 | } |
| 136 | |
| 137 | BigUnsigned::operator unsigned int() const { |
| 138 | if (len == 0) |
| 139 | return 0; |
| 140 | else if (len == 1 && (blk[0] & uiMask) == blk[0]) |
| 141 | return (unsigned int) blk[0]; |
| 142 | else |
| 143 | throw "BigUnsigned::operator unsigned int: Value is too big for an unsigned int"; |
| 144 | } |
| 145 | |
| 146 | BigUnsigned::operator int() const { |
| 147 | if (len == 0) |
| 148 | return 0; |
| 149 | else if (len == 1 && (blk[0] & iMask) == blk[0]) |
| 150 | return (int) blk[0]; |
| 151 | else |
| 152 | throw "BigUnsigned::operator int: Value is too big for an int"; |
| 153 | } |
| 154 | |
| 155 | BigUnsigned::operator unsigned short() const { |
| 156 | if (len == 0) |
| 157 | return 0; |
| 158 | else if (len == 1 && (blk[0] & usMask) == blk[0]) |
| 159 | return (unsigned short) blk[0]; |
| 160 | else |
| 161 | throw "BigUnsigned::operator unsigned short: Value is too big for an unsigned short"; |
| 162 | } |
| 163 | |
| 164 | BigUnsigned::operator short() const { |
| 165 | if (len == 0) |
| 166 | return 0; |
| 167 | else if (len == 1 && (blk[0] & sMask) == blk[0]) |
| 168 | return (short) blk[0]; |
| 169 | else |
| 170 | throw "BigUnsigned::operator short: Value is too big for a short"; |
| 171 | } |
| 172 | |
| 173 | // COMPARISON |
| 174 | BigUnsigned::CmpRes BigUnsigned::compareTo(const BigUnsigned &x) const { |
| 175 | // A bigger length implies a bigger number. |
| 176 | if (len < x.len) |
| 177 | return less; |
| 178 | else if (len > x.len) |
| 179 | return greater; |
| 180 | else { |
| 181 | // Compare blocks one by one from left to right. |
| 182 | Index i = len; |
| 183 | while (i > 0) { |
| 184 | i--; |
| 185 | if (blk[i] == x.blk[i]) |
| 186 | continue; |
| 187 | else if (blk[i] > x.blk[i]) |
| 188 | return greater; |
| 189 | else |
| 190 | return less; |
| 191 | } |
| 192 | // If no blocks differed, the numbers are equal. |
| 193 | return equal; |
| 194 | } |
| 195 | } |
| 196 | |
| 197 | // PUT-HERE OPERATIONS |
| 198 | |
| 199 | /* |
| 200 | * Below are implementations of the four basic arithmetic operations |
| 201 | * for `BigUnsigned's. Their purpose is to use a mechanism that can |
| 202 | * calculate the sum, difference, product, and quotient/remainder of |
| 203 | * two individual blocks in order to calculate the sum, difference, |
| 204 | * product, and quotient/remainder of two multi-block BigUnsigned |
| 205 | * numbers. |
| 206 | * |
| 207 | * As alluded to in the comment before class `BigUnsigned', |
| 208 | * these algorithms bear a remarkable similarity (in purpose, if |
| 209 | * not in implementation) to the way humans operate on big numbers. |
| 210 | * The built-in `+', `-', `*', `/' and `%' operators are analogous |
| 211 | * to elementary-school ``math facts'' and ``times tables''; the |
| 212 | * four routines below are analogous to ``long division'' and its |
| 213 | * relatives. (Only a computer can ``memorize'' a times table with |
| 214 | * 18446744073709551616 entries! (For 32-bit blocks.)) |
| 215 | * |
| 216 | * The discovery of these four algorithms, called the ``classical |
| 217 | * algorithms'', marked the beginning of the study of computer science. |
| 218 | * See Section 4.3.1 of Knuth's ``The Art of Computer Programming''. |
| 219 | */ |
| 220 | |
| 221 | /* |
| 222 | * On most calls to put-here operations, it's safe to read the inputs little by |
| 223 | * little and write the outputs little by little. However, if one of the |
| 224 | * inputs is coming from the same variable into which the output is to be |
| 225 | * stored (an "aliased" call), we risk overwriting the input before we read it. |
| 226 | * In this case, we first compute the result into a temporary BigUnsigned |
| 227 | * variable and then copy it into the requested output variable *this. |
| 228 | * Each put-here operation uses the DOTR_ALIASED macro (Do The Right Thing on |
| 229 | * aliased calls) to generate code for this check. |
| 230 | * |
| 231 | * I adopted this approach on 2007.02.13 (see Assignment Operators in |
| 232 | * BigUnsigned.hh). Before then, put-here operations rejected aliased calls |
| 233 | * with an exception. I think doing the right thing is better. |
| 234 | * |
| 235 | * Some of the put-here operations can probably handle aliased calls safely |
| 236 | * without the extra copy because (for example) they process blocks strictly |
| 237 | * right-to-left. At some point I might determine which ones don't need the |
| 238 | * copy, but my reasoning would need to be verified very carefully. For now |
| 239 | * I'll leave in the copy. |
| 240 | */ |
| 241 | #define DOTR_ALIASED(cond, op) \ |
| 242 | if (cond) { \ |
| 243 | BigUnsigned tmpThis; \ |
| 244 | tmpThis.op; \ |
| 245 | *this = tmpThis; \ |
| 246 | return; \ |
| 247 | } |
| 248 | |
| 249 | // Addition |
| 250 | void BigUnsigned::add(const BigUnsigned &a, const BigUnsigned &b) { |
| 251 | DOTR_ALIASED(this == &a || this == &b, add(a, b)); |
| 252 | // If one argument is zero, copy the other. |
| 253 | if (a.len == 0) { |
| 254 | operator =(b); |
| 255 | return; |
| 256 | } else if (b.len == 0) { |
| 257 | operator =(a); |
| 258 | return; |
| 259 | } |
| 260 | // Some variables... |
| 261 | // Carries in and out of an addition stage |
| 262 | bool carryIn, carryOut; |
| 263 | Blk temp; |
| 264 | Index i; |
| 265 | // a2 points to the longer input, b2 points to the shorter |
| 266 | const BigUnsigned *a2, *b2; |
| 267 | if (a.len >= b.len) { |
| 268 | a2 = &a; |
| 269 | b2 = &b; |
| 270 | } else { |
| 271 | a2 = &b; |
| 272 | b2 = &a; |
| 273 | } |
| 274 | // Set prelimiary length and make room in this BigUnsigned |
| 275 | len = a2->len + 1; |
| 276 | allocate(len); |
| 277 | // For each block index that is present in both inputs... |
| 278 | for (i = 0, carryIn = false; i < b2->len; i++) { |
| 279 | // Add input blocks |
| 280 | temp = a2->blk[i] + b2->blk[i]; |
| 281 | // If a rollover occurred, the result is less than either input. |
| 282 | // This test is used many times in the BigUnsigned code. |
| 283 | carryOut = (temp < a2->blk[i]); |
| 284 | // If a carry was input, handle it |
| 285 | if (carryIn) { |
| 286 | temp++; |
| 287 | carryOut |= (temp == 0); |
| 288 | } |
| 289 | blk[i] = temp; // Save the addition result |
| 290 | carryIn = carryOut; // Pass the carry along |
| 291 | } |
| 292 | // If there is a carry left over, increase blocks until |
| 293 | // one does not roll over. |
| 294 | for (; i < a2->len && carryIn; i++) { |
| 295 | temp = a2->blk[i] + 1; |
| 296 | carryIn = (temp == 0); |
| 297 | blk[i] = temp; |
| 298 | } |
| 299 | // If the carry was resolved but the larger number |
| 300 | // still has blocks, copy them over. |
| 301 | for (; i < a2->len; i++) |
| 302 | blk[i] = a2->blk[i]; |
| 303 | // Set the extra block if there's still a carry, decrease length otherwise |
| 304 | if (carryIn) |
| 305 | blk[i] = 1; |
| 306 | else |
| 307 | len--; |
| 308 | } |
| 309 | |
| 310 | // Subtraction |
| 311 | void BigUnsigned::subtract(const BigUnsigned &a, const BigUnsigned &b) { |
| 312 | DOTR_ALIASED(this == &a || this == &b, subtract(a, b)); |
| 313 | // If b is zero, copy a. If a is shorter than b, the result is negative. |
| 314 | if (b.len == 0) { |
| 315 | operator =(a); |
| 316 | return; |
| 317 | } else if (a.len < b.len) |
| 318 | throw "BigUnsigned::subtract: Negative result in unsigned calculation"; |
| 319 | // Some variables... |
| 320 | bool borrowIn, borrowOut; |
| 321 | Blk temp; |
| 322 | Index i; |
| 323 | // Set preliminary length and make room |
| 324 | len = a.len; |
| 325 | allocate(len); |
| 326 | // For each block index that is present in both inputs... |
| 327 | for (i = 0, borrowIn = false; i < b.len; i++) { |
| 328 | temp = a.blk[i] - b.blk[i]; |
| 329 | // If a reverse rollover occurred, the result is greater than the block from a. |
| 330 | borrowOut = (temp > a.blk[i]); |
| 331 | // Handle an incoming borrow |
| 332 | if (borrowIn) { |
| 333 | borrowOut |= (temp == 0); |
| 334 | temp--; |
| 335 | } |
| 336 | blk[i] = temp; // Save the subtraction result |
| 337 | borrowIn = borrowOut; // Pass the borrow along |
| 338 | } |
| 339 | // If there is a borrow left over, decrease blocks until |
| 340 | // one does not reverse rollover. |
| 341 | for (; i < a.len && borrowIn; i++) { |
| 342 | borrowIn = (a.blk[i] == 0); |
| 343 | blk[i] = a.blk[i] - 1; |
| 344 | } |
| 345 | // If there's still a borrow, the result is negative. |
| 346 | // Throw an exception, but zero out this object first just in case. |
| 347 | if (borrowIn) { |
| 348 | len = 0; |
| 349 | throw "BigUnsigned::subtract: Negative result in unsigned calculation"; |
| 350 | } else // Copy over the rest of the blocks |
| 351 | for (; i < a.len; i++) |
| 352 | blk[i] = a.blk[i]; |
| 353 | // Zap leading zeros |
| 354 | zapLeadingZeros(); |
| 355 | } |
| 356 | |
| 357 | /* |
| 358 | * About the multiplication and division algorithms: |
| 359 | * |
| 360 | * I searched unsucessfully for fast built-in operations like the `b_0' |
| 361 | * and `c_0' Knuth describes in Section 4.3.1 of ``The Art of Computer |
| 362 | * Programming'' (replace `place' by `Blk'): |
| 363 | * |
| 364 | * ``b_0[:] multiplication of a one-place integer by another one-place |
| 365 | * integer, giving a two-place answer; |
| 366 | * |
| 367 | * ``c_0[:] division of a two-place integer by a one-place integer, |
| 368 | * provided that the quotient is a one-place integer, and yielding |
| 369 | * also a one-place remainder.'' |
| 370 | * |
| 371 | * I also missed his note that ``[b]y adjusting the word size, if |
| 372 | * necessary, nearly all computers will have these three operations |
| 373 | * available'', so I gave up on trying to use algorithms similar to his. |
| 374 | * A future version of the library might include such algorithms; I |
| 375 | * would welcome contributions from others for this. |
| 376 | * |
| 377 | * I eventually decided to use bit-shifting algorithms. To multiply `a' |
| 378 | * and `b', we zero out the result. Then, for each `1' bit in `a', we |
| 379 | * shift `b' left the appropriate amount and add it to the result. |
| 380 | * Similarly, to divide `a' by `b', we shift `b' left varying amounts, |
| 381 | * repeatedly trying to subtract it from `a'. When we succeed, we note |
| 382 | * the fact by setting a bit in the quotient. While these algorithms |
| 383 | * have the same O(n^2) time complexity as Knuth's, the ``constant factor'' |
| 384 | * is likely to be larger. |
| 385 | * |
| 386 | * Because I used these algorithms, which require single-block addition |
| 387 | * and subtraction rather than single-block multiplication and division, |
| 388 | * the innermost loops of all four routines are very similar. Study one |
| 389 | * of them and all will become clear. |
| 390 | */ |
| 391 | |
| 392 | /* |
| 393 | * This is a little inline function used by both the multiplication |
| 394 | * routine and the division routine. |
| 395 | * |
| 396 | * `getShiftedBlock' returns the `x'th block of `num << y'. |
| 397 | * `y' may be anything from 0 to N - 1, and `x' may be anything from |
| 398 | * 0 to `num.len'. |
| 399 | * |
| 400 | * Two things contribute to this block: |
| 401 | * |
| 402 | * (1) The `N - y' low bits of `num.blk[x]', shifted `y' bits left. |
| 403 | * |
| 404 | * (2) The `y' high bits of `num.blk[x-1]', shifted `N - y' bits right. |
| 405 | * |
| 406 | * But we must be careful if `x == 0' or `x == num.len', in |
| 407 | * which case we should use 0 instead of (2) or (1), respectively. |
| 408 | * |
| 409 | * If `y == 0', then (2) contributes 0, as it should. However, |
| 410 | * in some computer environments, for a reason I cannot understand, |
| 411 | * `a >> b' means `a >> (b % N)'. This means `num.blk[x-1] >> (N - y)' |
| 412 | * will return `num.blk[x-1]' instead of the desired 0 when `y == 0'; |
| 413 | * the test `y == 0' handles this case specially. |
| 414 | */ |
| 415 | inline BigUnsigned::Blk getShiftedBlock(const BigUnsigned &num, |
| 416 | BigUnsigned::Index x, unsigned int y) { |
| 417 | BigUnsigned::Blk part1 = (x == 0 || y == 0) ? 0 : (num.blk[x - 1] >> (BigUnsigned::N - y)); |
| 418 | BigUnsigned::Blk part2 = (x == num.len) ? 0 : (num.blk[x] << y); |
| 419 | return part1 | part2; |
| 420 | } |
| 421 | |
| 422 | // Multiplication |
| 423 | void BigUnsigned::multiply(const BigUnsigned &a, const BigUnsigned &b) { |
| 424 | DOTR_ALIASED(this == &a || this == &b, multiply(a, b)); |
| 425 | // If either a or b is zero, set to zero. |
| 426 | if (a.len == 0 || b.len == 0) { |
| 427 | len = 0; |
| 428 | return; |
| 429 | } |
| 430 | /* |
| 431 | * Overall method: |
| 432 | * |
| 433 | * Set this = 0. |
| 434 | * For each 1-bit of `a' (say the `i2'th bit of block `i'): |
| 435 | * Add `b << (i blocks and i2 bits)' to *this. |
| 436 | */ |
| 437 | // Variables for the calculation |
| 438 | Index i, j, k; |
| 439 | unsigned int i2; |
| 440 | Blk temp; |
| 441 | bool carryIn, carryOut; |
| 442 | // Set preliminary length and make room |
| 443 | len = a.len + b.len; |
| 444 | allocate(len); |
| 445 | // Zero out this object |
| 446 | for (i = 0; i < len; i++) |
| 447 | blk[i] = 0; |
| 448 | // For each block of the first number... |
| 449 | for (i = 0; i < a.len; i++) { |
| 450 | // For each 1-bit of that block... |
| 451 | for (i2 = 0; i2 < N; i2++) { |
| 452 | if ((a.blk[i] & (Blk(1) << i2)) == 0) |
| 453 | continue; |
| 454 | /* |
| 455 | * Add b to this, shifted left i blocks and i2 bits. |
| 456 | * j is the index in b, and k = i + j is the index in this. |
| 457 | * |
| 458 | * `getShiftedBlock', a short inline function defined above, |
| 459 | * is now used for the bit handling. It replaces the more |
| 460 | * complex `bHigh' code, in which each run of the loop dealt |
| 461 | * immediately with the low bits and saved the high bits to |
| 462 | * be picked up next time. The last run of the loop used to |
| 463 | * leave leftover high bits, which were handled separately. |
| 464 | * Instead, this loop runs an additional time with j == b.len. |
| 465 | * These changes were made on 2005.01.11. |
| 466 | */ |
| 467 | for (j = 0, k = i, carryIn = false; j <= b.len; j++, k++) { |
| 468 | /* |
| 469 | * The body of this loop is very similar to the body of the first loop |
| 470 | * in `add', except that this loop does a `+=' instead of a `+'. |
| 471 | */ |
| 472 | temp = blk[k] + getShiftedBlock(b, j, i2); |
| 473 | carryOut = (temp < blk[k]); |
| 474 | if (carryIn) { |
| 475 | temp++; |
| 476 | carryOut |= (temp == 0); |
| 477 | } |
| 478 | blk[k] = temp; |
| 479 | carryIn = carryOut; |
| 480 | } |
| 481 | // No more extra iteration to deal with `bHigh'. |
| 482 | // Roll-over a carry as necessary. |
| 483 | for (; carryIn; k++) { |
| 484 | blk[k]++; |
| 485 | carryIn = (blk[k] == 0); |
| 486 | } |
| 487 | } |
| 488 | } |
| 489 | // Zap possible leading zero |
| 490 | if (blk[len - 1] == 0) |
| 491 | len--; |
| 492 | } |
| 493 | |
| 494 | /* |
| 495 | * DIVISION WITH REMAINDER |
| 496 | * The functionality of divide, modulo, and %= is included in this one monstrous call, |
| 497 | * which deserves some explanation. |
| 498 | * |
| 499 | * The division *this / b is performed. |
| 500 | * Afterwards, q has the quotient, and *this has the remainder. |
| 501 | * Thus, a call is like q = *this / b, *this %= b. |
| 502 | * |
| 503 | * This seemingly bizarre pattern of inputs and outputs has a justification. The |
| 504 | * ``put-here operations'' are supposed to be fast. Therefore, they accept inputs |
| 505 | * and provide outputs in the most convenient places so that no value ever needs |
| 506 | * to be copied in its entirety. That way, the client can perform exactly the |
| 507 | * copying it needs depending on where the inputs are and where it wants the output. |
| 508 | * A better name for this function might be "modWithQuotient", but I would rather |
| 509 | * not change the name now. |
| 510 | */ |
| 511 | void BigUnsigned::divideWithRemainder(const BigUnsigned &b, BigUnsigned &q) { |
| 512 | /* |
| 513 | * Defending against aliased calls is a bit tricky because we are |
| 514 | * writing to both *this and q. |
| 515 | * |
| 516 | * It would be silly to try to write quotient and remainder to the |
| 517 | * same variable. Rule that out right away. |
| 518 | */ |
| 519 | if (this == &q) |
| 520 | throw "BigUnsigned::divideWithRemainder: Cannot write quotient and remainder into the same variable"; |
| 521 | /* |
| 522 | * Now *this and q are separate, so the only concern is that b might be |
| 523 | * aliased to one of them. If so, use a temporary copy of b. |
| 524 | */ |
| 525 | if (this == &b || &q == &b) { |
| 526 | BigUnsigned tmpB(b); |
| 527 | divideWithRemainder(tmpB, q); |
| 528 | return; |
| 529 | } |
| 530 | |
| 531 | /* |
| 532 | * Note that the mathematical definition of mod (I'm trusting Knuth) is somewhat |
| 533 | * different from the way the normal C++ % operator behaves in the case of division by 0. |
| 534 | * This function does it Knuth's way. |
| 535 | * |
| 536 | * We let a / 0 == 0 (it doesn't matter) and a % 0 == a, no exceptions thrown. |
| 537 | * This allows us to preserve both Knuth's demand that a mod 0 == a |
| 538 | * and the useful property that (a / b) * b + (a % b) == a. |
| 539 | */ |
| 540 | if (b.len == 0) { |
| 541 | q.len = 0; |
| 542 | return; |
| 543 | } |
| 544 | |
| 545 | /* |
| 546 | * If *this.len < b.len, then *this < b, and we can be sure that b doesn't go into |
| 547 | * *this at all. The quotient is 0 and *this is already the remainder (so leave it alone). |
| 548 | */ |
| 549 | if (len < b.len) { |
| 550 | q.len = 0; |
| 551 | return; |
| 552 | } |
| 553 | |
| 554 | /* |
| 555 | * At this point we know *this > b > 0. (Whew!) |
| 556 | */ |
| 557 | |
| 558 | /* |
| 559 | * Overall method: |
| 560 | * |
| 561 | * For each appropriate i and i2, decreasing: |
| 562 | * Try to subtract (b << (i blocks and i2 bits)) from *this. |
| 563 | * (`work2' holds the result of this subtraction.) |
| 564 | * If the result is nonnegative: |
| 565 | * Turn on bit i2 of block i of the quotient q. |
| 566 | * Save the result of the subtraction back into *this. |
| 567 | * Otherwise: |
| 568 | * Bit i2 of block i remains off, and *this is unchanged. |
| 569 | * |
| 570 | * Eventually q will contain the entire quotient, and *this will |
| 571 | * be left with the remainder. |
| 572 | * |
| 573 | * We use work2 to temporarily store the result of a subtraction. |
| 574 | * work2[x] corresponds to blk[x], not blk[x+i], since 2005.01.11. |
| 575 | * If the subtraction is successful, we copy work2 back to blk. |
| 576 | * (There's no `work1'. In a previous version, when division was |
| 577 | * coded for a read-only dividend, `work1' played the role of |
| 578 | * the here-modifiable `*this' and got the remainder.) |
| 579 | * |
| 580 | * We never touch the i lowest blocks of either blk or work2 because |
| 581 | * they are unaffected by the subtraction: we are subtracting |
| 582 | * (b << (i blocks and i2 bits)), which ends in at least `i' zero blocks. |
| 583 | */ |
| 584 | // Variables for the calculation |
| 585 | Index i, j, k; |
| 586 | unsigned int i2; |
| 587 | Blk temp; |
| 588 | bool borrowIn, borrowOut; |
| 589 | |
| 590 | /* |
| 591 | * Make sure we have an extra zero block just past the value. |
| 592 | * |
| 593 | * When we attempt a subtraction, we might shift `b' so |
| 594 | * its first block begins a few bits left of the dividend, |
| 595 | * and then we'll try to compare these extra bits with |
| 596 | * a nonexistent block to the left of the dividend. The |
| 597 | * extra zero block ensures sensible behavior; we need |
| 598 | * an extra block in `work2' for exactly the same reason. |
| 599 | * |
| 600 | * See below `divideWithRemainder' for the interesting and |
| 601 | * amusing story of this section of code. |
| 602 | */ |
| 603 | Index origLen = len; // Save real length. |
| 604 | // 2006.05.03: Copy the number and then change the length! |
| 605 | allocateAndCopy(len + 1); // Get the space. |
| 606 | len++; // Increase the length. |
| 607 | blk[origLen] = 0; // Zero the extra block. |
| 608 | |
| 609 | // work2 holds part of the result of a subtraction; see above. |
| 610 | Blk *work2 = new Blk[len]; |
| 611 | |
| 612 | // Set preliminary length for quotient and make room |
| 613 | q.len = origLen - b.len + 1; |
| 614 | q.allocate(q.len); |
| 615 | // Zero out the quotient |
| 616 | for (i = 0; i < q.len; i++) |
| 617 | q.blk[i] = 0; |
| 618 | |
| 619 | // For each possible left-shift of b in blocks... |
| 620 | i = q.len; |
| 621 | while (i > 0) { |
| 622 | i--; |
| 623 | // For each possible left-shift of b in bits... |
| 624 | // (Remember, N is the number of bits in a Blk.) |
| 625 | q.blk[i] = 0; |
| 626 | i2 = N; |
| 627 | while (i2 > 0) { |
| 628 | i2--; |
| 629 | /* |
| 630 | * Subtract b, shifted left i blocks and i2 bits, from *this, |
| 631 | * and store the answer in work2. In the for loop, `k == i + j'. |
| 632 | * |
| 633 | * Compare this to the middle section of `multiply'. They |
| 634 | * are in many ways analogous. See especially the discussion |
| 635 | * of `getShiftedBlock'. |
| 636 | */ |
| 637 | for (j = 0, k = i, borrowIn = false; j <= b.len; j++, k++) { |
| 638 | temp = blk[k] - getShiftedBlock(b, j, i2); |
| 639 | borrowOut = (temp > blk[k]); |
| 640 | if (borrowIn) { |
| 641 | borrowOut |= (temp == 0); |
| 642 | temp--; |
| 643 | } |
| 644 | // Since 2005.01.11, indices of `work2' directly match those of `blk', so use `k'. |
| 645 | work2[k] = temp; |
| 646 | borrowIn = borrowOut; |
| 647 | } |
| 648 | // No more extra iteration to deal with `bHigh'. |
| 649 | // Roll-over a borrow as necessary. |
| 650 | for (; k < origLen && borrowIn; k++) { |
| 651 | borrowIn = (blk[k] == 0); |
| 652 | work2[k] = blk[k] - 1; |
| 653 | } |
| 654 | /* |
| 655 | * If the subtraction was performed successfully (!borrowIn), |
| 656 | * set bit i2 in block i of the quotient. |
| 657 | * |
| 658 | * Then, copy the portion of work2 filled by the subtraction |
| 659 | * back to *this. This portion starts with block i and ends-- |
| 660 | * where? Not necessarily at block `i + b.len'! Well, we |
| 661 | * increased k every time we saved a block into work2, so |
| 662 | * the region of work2 we copy is just [i, k). |
| 663 | */ |
| 664 | if (!borrowIn) { |
| 665 | q.blk[i] |= (Blk(1) << i2); |
| 666 | while (k > i) { |
| 667 | k--; |
| 668 | blk[k] = work2[k]; |
| 669 | } |
| 670 | } |
| 671 | } |
| 672 | } |
| 673 | // Zap possible leading zero in quotient |
| 674 | if (q.blk[q.len - 1] == 0) |
| 675 | q.len--; |
| 676 | // Zap any/all leading zeros in remainder |
| 677 | zapLeadingZeros(); |
| 678 | // Deallocate temporary array. |
| 679 | // (Thanks to Brad Spencer for noticing my accidental omission of this!) |
| 680 | delete [] work2; |
| 681 | |
| 682 | } |
| 683 | /* |
| 684 | * The out-of-bounds accesses story: |
| 685 | * |
| 686 | * On 2005.01.06 or 2005.01.07 (depending on your time zone), |
| 687 | * Milan Tomic reported out-of-bounds memory accesses in |
| 688 | * the Big Integer Library. To investigate the problem, I |
| 689 | * added code to bounds-check every access to the `blk' array |
| 690 | * of a `NumberlikeArray'. |
| 691 | * |
| 692 | * This gave me warnings that fell into two categories of false |
| 693 | * positives. The bounds checker was based on length, not |
| 694 | * capacity, and in two places I had accessed memory that I knew |
| 695 | * was inside the capacity but that wasn't inside the length: |
| 696 | * |
| 697 | * (1) The extra zero block at the left of `*this'. Earlier |
| 698 | * versions said `allocateAndCopy(len + 1); blk[len] = 0;' |
| 699 | * but did not increment `len'. |
| 700 | * |
| 701 | * (2) The entire digit array in the conversion constructor |
| 702 | * ``BigUnsignedInABase(BigUnsigned)''. It was allocated with |
| 703 | * a conservatively high capacity, but the length wasn't set |
| 704 | * until the end of the constructor. |
| 705 | * |
| 706 | * To simplify matters, I changed both sections of code so that |
| 707 | * all accesses occurred within the length. The messages went |
| 708 | * away, and I told Milan that I couldn't reproduce the problem, |
| 709 | * sending a development snapshot of the bounds-checked code. |
| 710 | * |
| 711 | * Then, on 2005.01.09-10, he told me his debugger still found |
| 712 | * problems, specifically at the line `delete [] work2'. |
| 713 | * It was `work2', not `blk', that was causing the problems; |
| 714 | * this possibility had not occurred to me at all. In fact, |
| 715 | * the problem was that `work2' needed an extra block just |
| 716 | * like `*this'. Go ahead and laugh at me for finding (1) |
| 717 | * without seeing what was actually causing the trouble. :-) |
| 718 | * |
| 719 | * The 2005.01.11 version fixes this problem. I hope this is |
| 720 | * the last of my memory-related bloopers. So this is what |
| 721 | * starts happening to your C++ code if you use Java too much! |
| 722 | */ |
| 723 | |
| 724 | // Bitwise and |
| 725 | void BigUnsigned::bitAnd(const BigUnsigned &a, const BigUnsigned &b) { |
| 726 | DOTR_ALIASED(this == &a || this == &b, bitAnd(a, b)); |
| 727 | len = (a.len >= b.len) ? b.len : a.len; |
| 728 | allocate(len); |
| 729 | Index i; |
| 730 | for (i = 0; i < len; i++) |
| 731 | blk[i] = a.blk[i] & b.blk[i]; |
| 732 | zapLeadingZeros(); |
| 733 | } |
| 734 | |
| 735 | // Bitwise or |
| 736 | void BigUnsigned::bitOr(const BigUnsigned &a, const BigUnsigned &b) { |
| 737 | DOTR_ALIASED(this == &a || this == &b, bitOr(a, b)); |
| 738 | Index i; |
| 739 | const BigUnsigned *a2, *b2; |
| 740 | if (a.len >= b.len) { |
| 741 | a2 = &a; |
| 742 | b2 = &b; |
| 743 | } else { |
| 744 | a2 = &b; |
| 745 | b2 = &a; |
| 746 | } |
| 747 | allocate(a2->len); |
| 748 | for (i = 0; i < b2->len; i++) |
| 749 | blk[i] = a2->blk[i] | b2->blk[i]; |
| 750 | for (; i < a2->len; i++) |
| 751 | blk[i] = a2->blk[i]; |
| 752 | len = a2->len; |
| 753 | } |
| 754 | |
| 755 | // Bitwise xor |
| 756 | void BigUnsigned::bitXor(const BigUnsigned &a, const BigUnsigned &b) { |
| 757 | DOTR_ALIASED(this == &a || this == &b, bitXor(a, b)); |
| 758 | Index i; |
| 759 | const BigUnsigned *a2, *b2; |
| 760 | if (a.len >= b.len) { |
| 761 | a2 = &a; |
| 762 | b2 = &b; |
| 763 | } else { |
| 764 | a2 = &b; |
| 765 | b2 = &a; |
| 766 | } |
| 767 | allocate(a2->len); |
| 768 | for (i = 0; i < b2->len; i++) |
| 769 | blk[i] = a2->blk[i] ^ b2->blk[i]; |
| 770 | for (; i < a2->len; i++) |
| 771 | blk[i] = a2->blk[i]; |
| 772 | len = a2->len; |
| 773 | zapLeadingZeros(); |
| 774 | } |
| 775 | |
| 776 | // INCREMENT/DECREMENT OPERATORS |
| 777 | |
| 778 | // Prefix increment |
| 779 | void BigUnsigned::operator ++() { |
| 780 | Index i; |
| 781 | bool carry = true; |
| 782 | for (i = 0; i < len && carry; i++) { |
| 783 | blk[i]++; |
| 784 | carry = (blk[i] == 0); |
| 785 | } |
| 786 | if (carry) { |
| 787 | // Matt fixed a bug 2004.12.24: next 2 lines used to say allocateAndCopy(len + 1) |
| 788 | // Matt fixed another bug 2006.04.24: |
| 789 | // old number only has len blocks, so copy before increasing length |
| 790 | allocateAndCopy(len + 1); |
| 791 | len++; |
| 792 | blk[i] = 1; |
| 793 | } |
| 794 | } |
| 795 | |
| 796 | // Postfix increment: same as prefix |
| 797 | void BigUnsigned::operator ++(int) { |
| 798 | operator ++(); |
| 799 | } |
| 800 | |
| 801 | // Prefix decrement |
| 802 | void BigUnsigned::operator --() { |
| 803 | if (len == 0) |
| 804 | throw "BigUnsigned::operator --(): Cannot decrement an unsigned zero"; |
| 805 | Index i; |
| 806 | bool borrow = true; |
| 807 | for (i = 0; borrow; i++) { |
| 808 | borrow = (blk[i] == 0); |
| 809 | blk[i]--; |
| 810 | } |
| 811 | // Zap possible leading zero (there can only be one) |
| 812 | if (blk[len - 1] == 0) |
| 813 | len--; |
| 814 | } |
| 815 | |
| 816 | // Postfix decrement: same as prefix |
| 817 | void BigUnsigned::operator --(int) { |
| 818 | operator --(); |
| 819 | } |