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