1 #include "BigUnsigned.hh"
3 // Memory management definitions have moved to the bottom of NumberlikeArray.hh.
5 // The templates used by these constructors and converters are at the bottom of
8 BigUnsigned::BigUnsigned(unsigned long x) { initFromPrimitive (x); }
9 BigUnsigned::BigUnsigned(unsigned int x) { initFromPrimitive (x); }
10 BigUnsigned::BigUnsigned(unsigned short x) { initFromPrimitive (x); }
11 BigUnsigned::BigUnsigned( long x) { initFromSignedPrimitive(x); }
12 BigUnsigned::BigUnsigned( int x) { initFromSignedPrimitive(x); }
13 BigUnsigned::BigUnsigned( short x) { initFromSignedPrimitive(x); }
15 unsigned long BigUnsigned::toUnsignedLong () const { return convertToPrimitive <unsigned long >(); }
16 unsigned int BigUnsigned::toUnsignedInt () const { return convertToPrimitive <unsigned int >(); }
17 unsigned short BigUnsigned::toUnsignedShort() const { return convertToPrimitive <unsigned short>(); }
18 long BigUnsigned::toLong () const { return convertToSignedPrimitive< long >(); }
19 int BigUnsigned::toInt () const { return convertToSignedPrimitive< int >(); }
20 short BigUnsigned::toShort () const { return convertToSignedPrimitive< short>(); }
22 // BIT/BLOCK ACCESSORS
24 void BigUnsigned::setBlock(Index i, Blk newBlock) {
30 // If i >= len, no effect.
33 // The nonzero block extends the number.
35 // Zero any added blocks that we aren't setting.
36 for (Index j = len; j < i; j++)
44 /* Evidently the compiler wants BigUnsigned:: on the return type because, at
45 * that point, it hasn't yet parsed the BigUnsigned:: on the name to get the
47 BigUnsigned::Index BigUnsigned::bitLength() const {
51 Blk leftmostBlock = getBlock(len - 1);
52 Index leftmostBlockLen = 0;
53 while (leftmostBlock != 0) {
57 return leftmostBlockLen + (len - 1) * N;
61 void BigUnsigned::setBit(Index bi, bool newBit) {
62 Index blockI = bi / N;
63 Blk block = getBlock(blockI), mask = 1 << (bi % N);
64 block = newBit ? (block | mask) : (block & ~mask);
65 setBlock(blockI, block);
69 BigUnsigned::CmpRes BigUnsigned::compareTo(const BigUnsigned &x) const {
70 // A bigger length implies a bigger number.
76 // Compare blocks one by one from left to right.
80 if (blk[i] == x.blk[i])
82 else if (blk[i] > x.blk[i])
87 // If no blocks differed, the numbers are equal.
92 // COPY-LESS OPERATIONS
95 * On most calls to copy-less operations, it's safe to read the inputs little by
96 * little and write the outputs little by little. However, if one of the
97 * inputs is coming from the same variable into which the output is to be
98 * stored (an "aliased" call), we risk overwriting the input before we read it.
99 * In this case, we first compute the result into a temporary BigUnsigned
100 * variable and then copy it into the requested output variable *this.
101 * Each put-here operation uses the DTRT_ALIASED macro (Do The Right Thing on
102 * aliased calls) to generate code for this check.
104 * I adopted this approach on 2007.02.13 (see Assignment Operators in
105 * BigUnsigned.hh). Before then, put-here operations rejected aliased calls
106 * with an exception. I think doing the right thing is better.
108 * Some of the put-here operations can probably handle aliased calls safely
109 * without the extra copy because (for example) they process blocks strictly
110 * right-to-left. At some point I might determine which ones don't need the
111 * copy, but my reasoning would need to be verified very carefully. For now
112 * I'll leave in the copy.
114 #define DTRT_ALIASED(cond, op) \
116 BigUnsigned tmpThis; \
124 void BigUnsigned::add(const BigUnsigned &a, const BigUnsigned &b) {
125 DTRT_ALIASED(this == &a || this == &b, add(a, b));
126 // If one argument is zero, copy the other.
130 } else if (b.len == 0) {
135 // Carries in and out of an addition stage
136 bool carryIn, carryOut;
139 // a2 points to the longer input, b2 points to the shorter
140 const BigUnsigned *a2, *b2;
141 if (a.len >= b.len) {
148 // Set prelimiary length and make room in this BigUnsigned
151 // For each block index that is present in both inputs...
152 for (i = 0, carryIn = false; i < b2->len; i++) {
154 temp = a2->blk[i] + b2->blk[i];
155 // If a rollover occurred, the result is less than either input.
156 // This test is used many times in the BigUnsigned code.
157 carryOut = (temp < a2->blk[i]);
158 // If a carry was input, handle it
161 carryOut |= (temp == 0);
163 blk[i] = temp; // Save the addition result
164 carryIn = carryOut; // Pass the carry along
166 // If there is a carry left over, increase blocks until
167 // one does not roll over.
168 for (; i < a2->len && carryIn; i++) {
169 temp = a2->blk[i] + 1;
170 carryIn = (temp == 0);
173 // If the carry was resolved but the larger number
174 // still has blocks, copy them over.
175 for (; i < a2->len; i++)
177 // Set the extra block if there's still a carry, decrease length otherwise
184 void BigUnsigned::subtract(const BigUnsigned &a, const BigUnsigned &b) {
185 DTRT_ALIASED(this == &a || this == &b, subtract(a, b));
187 // If b is zero, copy a.
190 } else if (a.len < b.len)
191 // If a is shorter than b, the result is negative.
192 throw "BigUnsigned::subtract: "
193 "Negative result in unsigned calculation";
195 bool borrowIn, borrowOut;
198 // Set preliminary length and make room
201 // For each block index that is present in both inputs...
202 for (i = 0, borrowIn = false; i < b.len; i++) {
203 temp = a.blk[i] - b.blk[i];
204 // If a reverse rollover occurred,
205 // the result is greater than the block from a.
206 borrowOut = (temp > a.blk[i]);
207 // Handle an incoming borrow
209 borrowOut |= (temp == 0);
212 blk[i] = temp; // Save the subtraction result
213 borrowIn = borrowOut; // Pass the borrow along
215 // If there is a borrow left over, decrease blocks until
216 // one does not reverse rollover.
217 for (; i < a.len && borrowIn; i++) {
218 borrowIn = (a.blk[i] == 0);
219 blk[i] = a.blk[i] - 1;
221 /* If there's still a borrow, the result is negative.
222 * Throw an exception, but zero out this object so as to leave it in a
223 * predictable state. */
226 throw "BigUnsigned::subtract: Negative result in unsigned calculation";
228 // Copy over the rest of the blocks
229 for (; i < a.len; i++)
236 * About the multiplication and division algorithms:
238 * I searched unsucessfully for fast C++ built-in operations like the `b_0'
239 * and `c_0' Knuth describes in Section 4.3.1 of ``The Art of Computer
240 * Programming'' (replace `place' by `Blk'):
242 * ``b_0[:] multiplication of a one-place integer by another one-place
243 * integer, giving a two-place answer;
245 * ``c_0[:] division of a two-place integer by a one-place integer,
246 * provided that the quotient is a one-place integer, and yielding
247 * also a one-place remainder.''
249 * I also missed his note that ``[b]y adjusting the word size, if
250 * necessary, nearly all computers will have these three operations
251 * available'', so I gave up on trying to use algorithms similar to his.
252 * A future version of the library might include such algorithms; I
253 * would welcome contributions from others for this.
255 * I eventually decided to use bit-shifting algorithms. To multiply `a'
256 * and `b', we zero out the result. Then, for each `1' bit in `a', we
257 * shift `b' left the appropriate amount and add it to the result.
258 * Similarly, to divide `a' by `b', we shift `b' left varying amounts,
259 * repeatedly trying to subtract it from `a'. When we succeed, we note
260 * the fact by setting a bit in the quotient. While these algorithms
261 * have the same O(n^2) time complexity as Knuth's, the ``constant factor''
262 * is likely to be larger.
264 * Because I used these algorithms, which require single-block addition
265 * and subtraction rather than single-block multiplication and division,
266 * the innermost loops of all four routines are very similar. Study one
267 * of them and all will become clear.
271 * This is a little inline function used by both the multiplication
272 * routine and the division routine.
274 * `getShiftedBlock' returns the `x'th block of `num << y'.
275 * `y' may be anything from 0 to N - 1, and `x' may be anything from
278 * Two things contribute to this block:
280 * (1) The `N - y' low bits of `num.blk[x]', shifted `y' bits left.
282 * (2) The `y' high bits of `num.blk[x-1]', shifted `N - y' bits right.
284 * But we must be careful if `x == 0' or `x == num.len', in
285 * which case we should use 0 instead of (2) or (1), respectively.
287 * If `y == 0', then (2) contributes 0, as it should. However,
288 * in some computer environments, for a reason I cannot understand,
289 * `a >> b' means `a >> (b % N)'. This means `num.blk[x-1] >> (N - y)'
290 * will return `num.blk[x-1]' instead of the desired 0 when `y == 0';
291 * the test `y == 0' handles this case specially.
293 inline BigUnsigned::Blk getShiftedBlock(const BigUnsigned &num,
294 BigUnsigned::Index x, unsigned int y) {
295 BigUnsigned::Blk part1 = (x == 0 || y == 0) ? 0 : (num.blk[x - 1] >> (BigUnsigned::N - y));
296 BigUnsigned::Blk part2 = (x == num.len) ? 0 : (num.blk[x] << y);
297 return part1 | part2;
300 void BigUnsigned::multiply(const BigUnsigned &a, const BigUnsigned &b) {
301 DTRT_ALIASED(this == &a || this == &b, multiply(a, b));
302 // If either a or b is zero, set to zero.
303 if (a.len == 0 || b.len == 0) {
311 * For each 1-bit of `a' (say the `i2'th bit of block `i'):
312 * Add `b << (i blocks and i2 bits)' to *this.
314 // Variables for the calculation
318 bool carryIn, carryOut;
319 // Set preliminary length and make room
322 // Zero out this object
323 for (i = 0; i < len; i++)
325 // For each block of the first number...
326 for (i = 0; i < a.len; i++) {
327 // For each 1-bit of that block...
328 for (i2 = 0; i2 < N; i2++) {
329 if ((a.blk[i] & (Blk(1) << i2)) == 0)
332 * Add b to this, shifted left i blocks and i2 bits.
333 * j is the index in b, and k = i + j is the index in this.
335 * `getShiftedBlock', a short inline function defined above,
336 * is now used for the bit handling. It replaces the more
337 * complex `bHigh' code, in which each run of the loop dealt
338 * immediately with the low bits and saved the high bits to
339 * be picked up next time. The last run of the loop used to
340 * leave leftover high bits, which were handled separately.
341 * Instead, this loop runs an additional time with j == b.len.
342 * These changes were made on 2005.01.11.
344 for (j = 0, k = i, carryIn = false; j <= b.len; j++, k++) {
346 * The body of this loop is very similar to the body of the first loop
347 * in `add', except that this loop does a `+=' instead of a `+'.
349 temp = blk[k] + getShiftedBlock(b, j, i2);
350 carryOut = (temp < blk[k]);
353 carryOut |= (temp == 0);
358 // No more extra iteration to deal with `bHigh'.
359 // Roll-over a carry as necessary.
360 for (; carryIn; k++) {
362 carryIn = (blk[k] == 0);
366 // Zap possible leading zero
367 if (blk[len - 1] == 0)
372 * DIVISION WITH REMAINDER
373 * This monstrous function mods *this by the given divisor b while storing the
374 * quotient in the given object q; at the end, *this contains the remainder.
375 * The seemingly bizarre pattern of inputs and outputs was chosen so that the
376 * function copies as little as possible (since it is implemented by repeated
377 * subtraction of multiples of b from *this).
379 * "modWithQuotient" might be a better name for this function, but I would
380 * rather not change the name now.
382 void BigUnsigned::divideWithRemainder(const BigUnsigned &b, BigUnsigned &q) {
383 /* Defending against aliased calls is more complex than usual because we
384 * are writing to both *this and q.
386 * It would be silly to try to write quotient and remainder to the
387 * same variable. Rule that out right away. */
389 throw "BigUnsigned::divideWithRemainder: Cannot write quotient and remainder into the same variable";
390 /* Now *this and q are separate, so the only concern is that b might be
391 * aliased to one of them. If so, use a temporary copy of b. */
392 if (this == &b || &q == &b) {
394 divideWithRemainder(tmpB, q);
399 * Knuth's definition of mod (which this function uses) is somewhat
400 * different from the C++ definition of % in case of division by 0.
402 * We let a / 0 == 0 (it doesn't matter much) and a % 0 == a, no
403 * exceptions thrown. This allows us to preserve both Knuth's demand
404 * that a mod 0 == a and the useful property that
405 * (a / b) * b + (a % b) == a.
413 * If *this.len < b.len, then *this < b, and we can be sure that b doesn't go into
414 * *this at all. The quotient is 0 and *this is already the remainder (so leave it alone).
421 // At this point we know (*this).len >= b.len > 0. (Whew!)
426 * For each appropriate i and i2, decreasing:
427 * Subtract (b << (i blocks and i2 bits)) from *this, storing the
428 * result in subtractBuf.
429 * If the subtraction succeeds with a nonnegative result:
430 * Turn on bit i2 of block i of the quotient q.
431 * Copy subtractBuf back into *this.
432 * Otherwise bit i2 of block i remains off, and *this is unchanged.
434 * Eventually q will contain the entire quotient, and *this will
435 * be left with the remainder.
437 * subtractBuf[x] corresponds to blk[x], not blk[x+i], since 2005.01.11.
438 * But on a single iteration, we don't touch the i lowest blocks of blk
439 * (and don't use those of subtractBuf) because these blocks are
440 * unaffected by the subtraction: we are subtracting
441 * (b << (i blocks and i2 bits)), which ends in at least `i' zero
443 // Variables for the calculation
447 bool borrowIn, borrowOut;
450 * Make sure we have an extra zero block just past the value.
452 * When we attempt a subtraction, we might shift `b' so
453 * its first block begins a few bits left of the dividend,
454 * and then we'll try to compare these extra bits with
455 * a nonexistent block to the left of the dividend. The
456 * extra zero block ensures sensible behavior; we need
457 * an extra block in `subtractBuf' for exactly the same reason.
459 Index origLen = len; // Save real length.
460 /* To avoid an out-of-bounds access in case of reallocation, allocate
461 * first and then increment the logical length. */
462 allocateAndCopy(len + 1);
464 blk[origLen] = 0; // Zero the added block.
466 // subtractBuf holds part of the result of a subtraction; see above.
467 Blk *subtractBuf = new Blk[len];
469 // Set preliminary length for quotient and make room
470 q.len = origLen - b.len + 1;
472 // Zero out the quotient
473 for (i = 0; i < q.len; i++)
476 // For each possible left-shift of b in blocks...
480 // For each possible left-shift of b in bits...
481 // (Remember, N is the number of bits in a Blk.)
487 * Subtract b, shifted left i blocks and i2 bits, from *this,
488 * and store the answer in subtractBuf. In the for loop, `k == i + j'.
490 * Compare this to the middle section of `multiply'. They
491 * are in many ways analogous. See especially the discussion
492 * of `getShiftedBlock'.
494 for (j = 0, k = i, borrowIn = false; j <= b.len; j++, k++) {
495 temp = blk[k] - getShiftedBlock(b, j, i2);
496 borrowOut = (temp > blk[k]);
498 borrowOut |= (temp == 0);
501 // Since 2005.01.11, indices of `subtractBuf' directly match those of `blk', so use `k'.
502 subtractBuf[k] = temp;
503 borrowIn = borrowOut;
505 // No more extra iteration to deal with `bHigh'.
506 // Roll-over a borrow as necessary.
507 for (; k < origLen && borrowIn; k++) {
508 borrowIn = (blk[k] == 0);
509 subtractBuf[k] = blk[k] - 1;
512 * If the subtraction was performed successfully (!borrowIn),
513 * set bit i2 in block i of the quotient.
515 * Then, copy the portion of subtractBuf filled by the subtraction
516 * back to *this. This portion starts with block i and ends--
517 * where? Not necessarily at block `i + b.len'! Well, we
518 * increased k every time we saved a block into subtractBuf, so
519 * the region of subtractBuf we copy is just [i, k).
522 q.blk[i] |= (Blk(1) << i2);
525 blk[k] = subtractBuf[k];
530 // Zap possible leading zero in quotient
531 if (q.blk[q.len - 1] == 0)
533 // Zap any/all leading zeros in remainder
535 // Deallocate subtractBuf.
536 // (Thanks to Brad Spencer for noticing my accidental omission of this!)
537 delete [] subtractBuf;
541 * These are straightforward blockwise operations except that they differ in
542 * the output length and the necessity of zapLeadingZeros. */
544 void BigUnsigned::bitAnd(const BigUnsigned &a, const BigUnsigned &b) {
545 DTRT_ALIASED(this == &a || this == &b, bitAnd(a, b));
546 // The bitwise & can't be longer than either operand.
547 len = (a.len >= b.len) ? b.len : a.len;
550 for (i = 0; i < len; i++)
551 blk[i] = a.blk[i] & b.blk[i];
555 void BigUnsigned::bitOr(const BigUnsigned &a, const BigUnsigned &b) {
556 DTRT_ALIASED(this == &a || this == &b, bitOr(a, b));
558 const BigUnsigned *a2, *b2;
559 if (a.len >= b.len) {
567 for (i = 0; i < b2->len; i++)
568 blk[i] = a2->blk[i] | b2->blk[i];
569 for (; i < a2->len; i++)
572 // Doesn't need zapLeadingZeros.
575 void BigUnsigned::bitXor(const BigUnsigned &a, const BigUnsigned &b) {
576 DTRT_ALIASED(this == &a || this == &b, bitXor(a, b));
578 const BigUnsigned *a2, *b2;
579 if (a.len >= b.len) {
587 for (i = 0; i < b2->len; i++)
588 blk[i] = a2->blk[i] ^ b2->blk[i];
589 for (; i < a2->len; i++)
595 void BigUnsigned::bitShiftLeft(const BigUnsigned &a, int b) {
596 DTRT_ALIASED(this == &a, bitShiftLeft(a, b));
599 throw "BigUnsigned::bitShiftLeft: "
600 "Pathological shift amount not implemented";
602 bitShiftRight(a, -b);
606 Index shiftBlocks = b / N;
607 unsigned int shiftBits = b % N;
608 // + 1: room for high bits nudged left into another block
609 len = a.len + shiftBlocks + 1;
612 for (i = 0; i < shiftBlocks; i++)
614 for (j = 0, i = shiftBlocks; j <= a.len; j++, i++)
615 blk[i] = getShiftedBlock(a, j, shiftBits);
616 // Zap possible leading zero
617 if (blk[len - 1] == 0)
621 void BigUnsigned::bitShiftRight(const BigUnsigned &a, int b) {
622 DTRT_ALIASED(this == &a, bitShiftRight(a, b));
625 throw "BigUnsigned::bitShiftRight: "
626 "Pathological shift amount not implemented";
632 // This calculation is wacky, but expressing the shift as a left bit shift
633 // within each block lets us use getShiftedBlock.
634 Index rightShiftBlocks = (b + N - 1) / N;
635 unsigned int leftShiftBits = N * rightShiftBlocks - b;
636 // Now (N * rightShiftBlocks - leftShiftBits) == b
637 // and 0 <= leftShiftBits < N.
638 if (rightShiftBlocks >= a.len + 1) {
639 // All of a is guaranteed to be shifted off, even considering the left
644 // Now we're allocating a positive amount.
645 // + 1: room for high bits nudged left into another block
646 len = a.len + 1 - rightShiftBlocks;
649 for (j = rightShiftBlocks, i = 0; j <= a.len; j++, i++)
650 blk[i] = getShiftedBlock(a, j, leftShiftBits);
651 // Zap possible leading zero
652 if (blk[len - 1] == 0)
656 // INCREMENT/DECREMENT OPERATORS
659 void BigUnsigned::operator ++() {
662 for (i = 0; i < len && carry; i++) {
664 carry = (blk[i] == 0);
667 // Allocate and then increase length, as in divideWithRemainder
668 allocateAndCopy(len + 1);
674 // Postfix increment: same as prefix
675 void BigUnsigned::operator ++(int) {
680 void BigUnsigned::operator --() {
682 throw "BigUnsigned::operator --(): Cannot decrement an unsigned zero";
685 for (i = 0; borrow; i++) {
686 borrow = (blk[i] == 0);
689 // Zap possible leading zero (there can only be one)
690 if (blk[len - 1] == 0)
694 // Postfix decrement: same as prefix
695 void BigUnsigned::operator --(int) {