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