1 /*
2 * Copyright (c) 1999, 2013, Oracle and/or its affiliates. All rights reserved.
3 * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4 *
5 * This code is free software; you can redistribute it and/or modify it
6 * under the terms of the GNU General Public License version 2 only, as
7 * published by the Free Software Foundation. Oracle designates this
8 * particular file as subject to the "Classpath" exception as provided
9 * by Oracle in the LICENSE file that accompanied this code.
10 *
11 * This code is distributed in the hope that it will be useful, but WITHOUT
12 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
14 * version 2 for more details (a copy is included in the LICENSE file that
15 * accompanied this code).
16 *
17 * You should have received a copy of the GNU General Public License version
18 * 2 along with this work; if not, write to the Free Software Foundation,
19 * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20 *
21 * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22 * or visit www.oracle.com if you need additional information or have any
23 * questions.
24 */
25
26 package java.math;
27
28 /**
29 * A class used to represent multiprecision integers that makes efficient
30 * use of allocated space by allowing a number to occupy only part of
31 * an array so that the arrays do not have to be reallocated as often.
32 * When performing an operation with many iterations the array used to
33 * hold a number is only reallocated when necessary and does not have to
34 * be the same size as the number it represents. A mutable number allows
35 * calculations to occur on the same number without having to create
36 * a new number for every step of the calculation as occurs with
37 * BigIntegers.
38 *
39 * @see BigInteger
40 * @author Michael McCloskey
41 * @author Timothy Buktu
42 * @since 1.3
43 */
44
45 import static java.math.BigDecimal.INFLATED;
46 import static java.math.BigInteger.LONG_MASK;
47 import java.util.Arrays;
48
49 class MutableBigInteger {
50 /**
51 * Holds the magnitude of this MutableBigInteger in big endian order.
52 * The magnitude may start at an offset into the value array, and it may
53 * end before the length of the value array.
54 */
55 int[] value;
56
57 /**
58 * The number of ints of the value array that are currently used
59 * to hold the magnitude of this MutableBigInteger. The magnitude starts
60 * at an offset and offset + intLen may be less than value.length.
61 */
62 int intLen;
63
64 /**
65 * The offset into the value array where the magnitude of this
66 * MutableBigInteger begins.
67 */
68 int offset = 0;
69
70 // Constants
71 /**
72 * MutableBigInteger with one element value array with the value 1. Used by
73 * BigDecimal divideAndRound to increment the quotient. Use this constant
74 * only when the method is not going to modify this object.
75 */
76 static final MutableBigInteger ONE = new MutableBigInteger(1);
77
78 /**
79 * The minimum {@code intLen} for cancelling powers of two before
80 * dividing.
81 * If the number of ints is less than this threshold,
82 * {@code divideKnuth} does not eliminate common powers of two from
83 * the dividend and divisor.
84 */
85 static final int KNUTH_POW2_THRESH_LEN = 6;
86
87 /**
88 * The minimum number of trailing zero ints for cancelling powers of two
89 * before dividing.
90 * If the dividend and divisor don't share at least this many zero ints
91 * at the end, {@code divideKnuth} does not eliminate common powers
92 * of two from the dividend and divisor.
93 */
94 static final int KNUTH_POW2_THRESH_ZEROS = 3;
95
96 // Constructors
97
98 /**
99 * The default constructor. An empty MutableBigInteger is created with
100 * a one word capacity.
101 */
102 MutableBigInteger() {
103 value = new int[1];
104 intLen = 0;
105 }
106
107 /**
108 * Construct a new MutableBigInteger with a magnitude specified by
109 * the int val.
110 */
111 MutableBigInteger(int val) {
112 value = new int[1];
113 intLen = 1;
114 value[0] = val;
115 }
116
117 /**
118 * Construct a new MutableBigInteger with the specified value array
119 * up to the length of the array supplied.
120 */
121 MutableBigInteger(int[] val) {
122 value = val;
123 intLen = val.length;
124 }
125
126 /**
127 * Construct a new MutableBigInteger with a magnitude equal to the
128 * specified BigInteger.
129 */
130 MutableBigInteger(BigInteger b) {
131 intLen = b.mag.length;
132 value = Arrays.copyOf(b.mag, intLen);
133 }
134
135 /**
136 * Construct a new MutableBigInteger with a magnitude equal to the
137 * specified MutableBigInteger.
138 */
139 MutableBigInteger(MutableBigInteger val) {
140 intLen = val.intLen;
141 value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
142 }
143
144 /**
145 * Makes this number an {@code n}-int number all of whose bits are ones.
146 * Used by Burnikel-Ziegler division.
147 * @param n number of ints in the {@code value} array
148 * @return a number equal to {@code ((1<<(32*n)))-1}
149 */
150 private void ones(int n) {
151 if (n > value.length)
152 value = new int[n];
153 Arrays.fill(value, -1);
154 offset = 0;
155 intLen = n;
156 }
157
158 /**
159 * Internal helper method to return the magnitude array. The caller is not
160 * supposed to modify the returned array.
161 */
162 private int[] getMagnitudeArray() {
163 if (offset > 0 || value.length != intLen)
164 return Arrays.copyOfRange(value, offset, offset + intLen);
165 return value;
166 }
167
168 /**
169 * Convert this MutableBigInteger to a long value. The caller has to make
170 * sure this MutableBigInteger can be fit into long.
171 */
172 private long toLong() {
173 assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
174 if (intLen == 0)
175 return 0;
176 long d = value[offset] & LONG_MASK;
177 return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
178 }
179
180 /**
181 * Convert this MutableBigInteger to a BigInteger object.
182 */
183 BigInteger toBigInteger(int sign) {
184 if (intLen == 0 || sign == 0)
185 return BigInteger.ZERO;
186 return new BigInteger(getMagnitudeArray(), sign);
187 }
188
189 /**
190 * Converts this number to a nonnegative {@code BigInteger}.
191 */
192 BigInteger toBigInteger() {
193 normalize();
194 return toBigInteger(isZero() ? 0 : 1);
195 }
196
197 /**
198 * Convert this MutableBigInteger to BigDecimal object with the specified sign
199 * and scale.
200 */
201 BigDecimal toBigDecimal(int sign, int scale) {
202 if (intLen == 0 || sign == 0)
203 return BigDecimal.zeroValueOf(scale);
204 int[] mag = getMagnitudeArray();
205 int len = mag.length;
206 int d = mag[0];
207 // If this MutableBigInteger can't be fit into long, we need to
208 // make a BigInteger object for the resultant BigDecimal object.
209 if (len > 2 || (d < 0 && len == 2))
210 return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
211 long v = (len == 2) ?
212 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
213 d & LONG_MASK;
214 return BigDecimal.valueOf(sign == -1 ? -v : v, scale);
215 }
216
217 /**
218 * This is for internal use in converting from a MutableBigInteger
219 * object into a long value given a specified sign.
220 * returns INFLATED if value is not fit into long
221 */
222 long toCompactValue(int sign) {
223 if (intLen == 0 || sign == 0)
224 return 0L;
225 int[] mag = getMagnitudeArray();
226 int len = mag.length;
227 int d = mag[0];
228 // If this MutableBigInteger can not be fitted into long, we need to
229 // make a BigInteger object for the resultant BigDecimal object.
230 if (len > 2 || (d < 0 && len == 2))
231 return INFLATED;
232 long v = (len == 2) ?
233 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
234 d & LONG_MASK;
235 return sign == -1 ? -v : v;
236 }
237
238 /**
239 * Clear out a MutableBigInteger for reuse.
240 */
241 void clear() {
242 offset = intLen = 0;
243 for (int index=0, n=value.length; index < n; index++)
244 value[index] = 0;
245 }
246
247 /**
248 * Set a MutableBigInteger to zero, removing its offset.
249 */
250 void reset() {
251 offset = intLen = 0;
252 }
253
254 /**
255 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
256 * as this MutableBigInteger is numerically less than, equal to, or
257 * greater than <tt>b</tt>.
258 */
259 final int compare(MutableBigInteger b) {
260 int blen = b.intLen;
261 if (intLen < blen)
262 return -1;
263 if (intLen > blen)
264 return 1;
265
266 // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
267 // comparison.
268 int[] bval = b.value;
269 for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
270 int b1 = value[i] + 0x80000000;
271 int b2 = bval[j] + 0x80000000;
272 if (b1 < b2)
273 return -1;
274 if (b1 > b2)
275 return 1;
276 }
277 return 0;
278 }
279
280 /**
281 * Returns a value equal to what {@code b.leftShift(32*ints); return compare(b);}
282 * would return, but doesn't change the value of {@code b}.
283 */
284 private int compareShifted(MutableBigInteger b, int ints) {
285 int blen = b.intLen;
286 int alen = intLen - ints;
287 if (alen < blen)
288 return -1;
289 if (alen > blen)
290 return 1;
291
292 // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
293 // comparison.
294 int[] bval = b.value;
295 for (int i = offset, j = b.offset; i < alen + offset; i++, j++) {
296 int b1 = value[i] + 0x80000000;
297 int b2 = bval[j] + 0x80000000;
298 if (b1 < b2)
299 return -1;
300 if (b1 > b2)
301 return 1;
302 }
303 return 0;
304 }
305
306 /**
307 * Compare this against half of a MutableBigInteger object (Needed for
308 * remainder tests).
309 * Assumes no leading unnecessary zeros, which holds for results
310 * from divide().
311 */
312 final int compareHalf(MutableBigInteger b) {
313 int blen = b.intLen;
314 int len = intLen;
315 if (len <= 0)
316 return blen <= 0 ? 0 : -1;
317 if (len > blen)
318 return 1;
319 if (len < blen - 1)
320 return -1;
321 int[] bval = b.value;
322 int bstart = 0;
323 int carry = 0;
324 // Only 2 cases left:len == blen or len == blen - 1
325 if (len != blen) { // len == blen - 1
326 if (bval[bstart] == 1) {
327 ++bstart;
328 carry = 0x80000000;
329 } else
330 return -1;
331 }
332 // compare values with right-shifted values of b,
333 // carrying shifted-out bits across words
334 int[] val = value;
335 for (int i = offset, j = bstart; i < len + offset;) {
336 int bv = bval[j++];
337 long hb = ((bv >>> 1) + carry) & LONG_MASK;
338 long v = val[i++] & LONG_MASK;
339 if (v != hb)
340 return v < hb ? -1 : 1;
341 carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
342 }
343 return carry == 0 ? 0 : -1;
344 }
345
346 /**
347 * Return the index of the lowest set bit in this MutableBigInteger. If the
348 * magnitude of this MutableBigInteger is zero, -1 is returned.
349 */
350 private final int getLowestSetBit() {
351 if (intLen == 0)
352 return -1;
353 int j, b;
354 for (j=intLen-1; (j > 0) && (value[j+offset] == 0); j--)
355 ;
356 b = value[j+offset];
357 if (b == 0)
358 return -1;
359 return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
360 }
361
362 /**
363 * Return the int in use in this MutableBigInteger at the specified
364 * index. This method is not used because it is not inlined on all
365 * platforms.
366 */
367 private final int getInt(int index) {
368 return value[offset+index];
369 }
370
371 /**
372 * Return a long which is equal to the unsigned value of the int in
373 * use in this MutableBigInteger at the specified index. This method is
374 * not used because it is not inlined on all platforms.
375 */
376 private final long getLong(int index) {
377 return value[offset+index] & LONG_MASK;
378 }
379
380 /**
381 * Ensure that the MutableBigInteger is in normal form, specifically
382 * making sure that there are no leading zeros, and that if the
383 * magnitude is zero, then intLen is zero.
384 */
385 final void normalize() {
386 if (intLen == 0) {
387 offset = 0;
388 return;
389 }
390
391 int index = offset;
392 if (value[index] != 0)
393 return;
394
395 int indexBound = index+intLen;
396 do {
397 index++;
398 } while(index < indexBound && value[index] == 0);
399
400 int numZeros = index - offset;
401 intLen -= numZeros;
402 offset = (intLen == 0 ? 0 : offset+numZeros);
403 }
404
405 /**
406 * If this MutableBigInteger cannot hold len words, increase the size
407 * of the value array to len words.
408 */
409 private final void ensureCapacity(int len) {
410 if (value.length < len) {
411 value = new int[len];
412 offset = 0;
413 intLen = len;
414 }
415 }
416
417 /**
418 * Convert this MutableBigInteger into an int array with no leading
419 * zeros, of a length that is equal to this MutableBigInteger's intLen.
420 */
421 int[] toIntArray() {
422 int[] result = new int[intLen];
423 for(int i=0; i < intLen; i++)
424 result[i] = value[offset+i];
425 return result;
426 }
427
428 /**
429 * Sets the int at index+offset in this MutableBigInteger to val.
430 * This does not get inlined on all platforms so it is not used
431 * as often as originally intended.
432 */
433 void setInt(int index, int val) {
434 value[offset + index] = val;
435 }
436
437 /**
438 * Sets this MutableBigInteger's value array to the specified array.
439 * The intLen is set to the specified length.
440 */
441 void setValue(int[] val, int length) {
442 value = val;
443 intLen = length;
444 offset = 0;
445 }
446
447 /**
448 * Sets this MutableBigInteger's value array to a copy of the specified
449 * array. The intLen is set to the length of the new array.
450 */
451 void copyValue(MutableBigInteger src) {
452 int len = src.intLen;
453 if (value.length < len)
454 value = new int[len];
455 System.arraycopy(src.value, src.offset, value, 0, len);
456 intLen = len;
457 offset = 0;
458 }
459
460 /**
461 * Sets this MutableBigInteger's value array to a copy of the specified
462 * array. The intLen is set to the length of the specified array.
463 */
464 void copyValue(int[] val) {
465 int len = val.length;
466 if (value.length < len)
467 value = new int[len];
468 System.arraycopy(val, 0, value, 0, len);
469 intLen = len;
470 offset = 0;
471 }
472
473 /**
474 * Returns true iff this MutableBigInteger has a value of one.
475 */
476 boolean isOne() {
477 return (intLen == 1) && (value[offset] == 1);
478 }
479
480 /**
481 * Returns true iff this MutableBigInteger has a value of zero.
482 */
483 boolean isZero() {
484 return (intLen == 0);
485 }
486
487 /**
488 * Returns true iff this MutableBigInteger is even.
489 */
490 boolean isEven() {
491 return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
492 }
493
494 /**
495 * Returns true iff this MutableBigInteger is odd.
496 */
497 boolean isOdd() {
498 return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
499 }
500
501 /**
502 * Returns true iff this MutableBigInteger is in normal form. A
503 * MutableBigInteger is in normal form if it has no leading zeros
504 * after the offset, and intLen + offset <= value.length.
505 */
506 boolean isNormal() {
507 if (intLen + offset > value.length)
508 return false;
509 if (intLen == 0)
510 return true;
511 return (value[offset] != 0);
512 }
513
514 /**
515 * Returns a String representation of this MutableBigInteger in radix 10.
516 */
517 public String toString() {
518 BigInteger b = toBigInteger(1);
519 return b.toString();
520 }
521
522 /**
523 * Like {@link #rightShift(int)} but {@code n} can be greater than the length of the number.
524 */
525 void safeRightShift(int n) {
526 if (n/32 >= intLen) {
527 reset();
528 } else {
529 rightShift(n);
530 }
531 }
532
533 /**
534 * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
535 * in normal form.
536 */
537 void rightShift(int n) {
538 if (intLen == 0)
539 return;
540 int nInts = n >>> 5;
541 int nBits = n & 0x1F;
542 this.intLen -= nInts;
543 if (nBits == 0)
544 return;
545 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
546 if (nBits >= bitsInHighWord) {
547 this.primitiveLeftShift(32 - nBits);
548 this.intLen--;
549 } else {
550 primitiveRightShift(nBits);
551 }
552 }
553
554 /**
555 * Like {@link #leftShift(int)} but {@code n} can be zero.
556 */
557 void safeLeftShift(int n) {
558 if (n > 0) {
559 leftShift(n);
560 }
561 }
562
563 /**
564 * Left shift this MutableBigInteger n bits.
565 */
566 void leftShift(int n) {
567 /*
568 * If there is enough storage space in this MutableBigInteger already
569 * the available space will be used. Space to the right of the used
570 * ints in the value array is faster to utilize, so the extra space
571 * will be taken from the right if possible.
572 */
573 if (intLen == 0)
574 return;
575 int nInts = n >>> 5;
576 int nBits = n&0x1F;
577 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
578
579 // If shift can be done without moving words, do so
580 if (n <= (32-bitsInHighWord)) {
581 primitiveLeftShift(nBits);
582 return;
583 }
584
585 int newLen = intLen + nInts +1;
586 if (nBits <= (32-bitsInHighWord))
587 newLen--;
588 if (value.length < newLen) {
589 // The array must grow
590 int[] result = new int[newLen];
591 for (int i=0; i < intLen; i++)
592 result[i] = value[offset+i];
593 setValue(result, newLen);
594 } else if (value.length - offset >= newLen) {
595 // Use space on right
596 for(int i=0; i < newLen - intLen; i++)
597 value[offset+intLen+i] = 0;
598 } else {
599 // Must use space on left
600 for (int i=0; i < intLen; i++)
601 value[i] = value[offset+i];
602 for (int i=intLen; i < newLen; i++)
603 value[i] = 0;
604 offset = 0;
605 }
606 intLen = newLen;
607 if (nBits == 0)
608 return;
609 if (nBits <= (32-bitsInHighWord))
610 primitiveLeftShift(nBits);
611 else
612 primitiveRightShift(32 -nBits);
613 }
614
615 /**
616 * A primitive used for division. This method adds in one multiple of the
617 * divisor a back to the dividend result at a specified offset. It is used
618 * when qhat was estimated too large, and must be adjusted.
619 */
620 private int divadd(int[] a, int[] result, int offset) {
621 long carry = 0;
622
623 for (int j=a.length-1; j >= 0; j--) {
624 long sum = (a[j] & LONG_MASK) +
625 (result[j+offset] & LONG_MASK) + carry;
626 result[j+offset] = (int)sum;
627 carry = sum >>> 32;
628 }
629 return (int)carry;
630 }
631
632 /**
633 * This method is used for division. It multiplies an n word input a by one
634 * word input x, and subtracts the n word product from q. This is needed
635 * when subtracting qhat*divisor from dividend.
636 */
637 private int mulsub(int[] q, int[] a, int x, int len, int offset) {
638 long xLong = x & LONG_MASK;
639 long carry = 0;
640 offset += len;
641
642 for (int j=len-1; j >= 0; j--) {
643 long product = (a[j] & LONG_MASK) * xLong + carry;
644 long difference = q[offset] - product;
645 q[offset--] = (int)difference;
646 carry = (product >>> 32)
647 + (((difference & LONG_MASK) >
648 (((~(int)product) & LONG_MASK))) ? 1:0);
649 }
650 return (int)carry;
651 }
652
653 /**
654 * The method is the same as mulsun, except the fact that q array is not
655 * updated, the only result of the method is borrow flag.
656 */
657 private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
658 long xLong = x & LONG_MASK;
659 long carry = 0;
660 offset += len;
661 for (int j=len-1; j >= 0; j--) {
662 long product = (a[j] & LONG_MASK) * xLong + carry;
663 long difference = q[offset--] - product;
664 carry = (product >>> 32)
665 + (((difference & LONG_MASK) >
666 (((~(int)product) & LONG_MASK))) ? 1:0);
667 }
668 return (int)carry;
669 }
670
671 /**
672 * Right shift this MutableBigInteger n bits, where n is
673 * less than 32.
674 * Assumes that intLen > 0, n > 0 for speed
675 */
676 private final void primitiveRightShift(int n) {
677 int[] val = value;
678 int n2 = 32 - n;
679 for (int i=offset+intLen-1, c=val[i]; i > offset; i--) {
680 int b = c;
681 c = val[i-1];
682 val[i] = (c << n2) | (b >>> n);
683 }
684 val[offset] >>>= n;
685 }
686
687 /**
688 * Left shift this MutableBigInteger n bits, where n is
689 * less than 32.
690 * Assumes that intLen > 0, n > 0 for speed
691 */
692 private final void primitiveLeftShift(int n) {
693 int[] val = value;
694 int n2 = 32 - n;
695 for (int i=offset, c=val[i], m=i+intLen-1; i < m; i++) {
696 int b = c;
697 c = val[i+1];
698 val[i] = (b << n) | (c >>> n2);
699 }
700 val[offset+intLen-1] <<= n;
701 }
702
703 /**
704 * Returns a {@code BigInteger} equal to the {@code n}
705 * low ints of this number.
706 */
707 private BigInteger getLower(int n) {
708 if (isZero()) {
709 return BigInteger.ZERO;
710 } else if (intLen < n) {
711 return toBigInteger(1);
712 } else {
713 // strip zeros
714 int len = n;
715 while (len > 0 && value[offset+intLen-len] == 0)
716 len--;
717 int sign = len > 0 ? 1 : 0;
718 return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
719 }
720 }
721
722 /**
723 * Discards all ints whose index is greater than {@code n}.
724 */
725 private void keepLower(int n) {
726 if (intLen >= n) {
727 offset += intLen - n;
728 intLen = n;
729 }
730 }
731
732 /**
733 * Adds the contents of two MutableBigInteger objects.The result
734 * is placed within this MutableBigInteger.
735 * The contents of the addend are not changed.
736 */
737 void add(MutableBigInteger addend) {
738 int x = intLen;
739 int y = addend.intLen;
740 int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
741 int[] result = (value.length < resultLen ? new int[resultLen] : value);
742
743 int rstart = result.length-1;
744 long sum;
745 long carry = 0;
746
747 // Add common parts of both numbers
748 while(x > 0 && y > 0) {
749 x--; y--;
750 sum = (value[x+offset] & LONG_MASK) +
751 (addend.value[y+addend.offset] & LONG_MASK) + carry;
752 result[rstart--] = (int)sum;
753 carry = sum >>> 32;
754 }
755
756 // Add remainder of the longer number
757 while(x > 0) {
758 x--;
759 if (carry == 0 && result == value && rstart == (x + offset))
760 return;
761 sum = (value[x+offset] & LONG_MASK) + carry;
762 result[rstart--] = (int)sum;
763 carry = sum >>> 32;
764 }
765 while(y > 0) {
766 y--;
767 sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
768 result[rstart--] = (int)sum;
769 carry = sum >>> 32;
770 }
771
772 if (carry > 0) { // Result must grow in length
773 resultLen++;
774 if (result.length < resultLen) {
775 int temp[] = new int[resultLen];
776 // Result one word longer from carry-out; copy low-order
777 // bits into new result.
778 System.arraycopy(result, 0, temp, 1, result.length);
779 temp[0] = 1;
780 result = temp;
781 } else {
782 result[rstart--] = 1;
783 }
784 }
785
786 value = result;
787 intLen = resultLen;
788 offset = result.length - resultLen;
789 }
790
791 /**
792 * Adds the value of {@code addend} shifted {@code n} ints to the left.
793 * Has the same effect as {@code addend.leftShift(32*ints); add(addend);}
794 * but doesn't change the value of {@code addend}.
795 */
796 void addShifted(MutableBigInteger addend, int n) {
797 if (addend.isZero()) {
798 return;
799 }
800
801 int x = intLen;
802 int y = addend.intLen + n;
803 int resultLen = (intLen > y ? intLen : y);
804 int[] result = (value.length < resultLen ? new int[resultLen] : value);
805
806 int rstart = result.length-1;
807 long sum;
808 long carry = 0;
809
810 // Add common parts of both numbers
811 while (x > 0 && y > 0) {
812 x--; y--;
813 int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
814 sum = (value[x+offset] & LONG_MASK) +
815 (bval & LONG_MASK) + carry;
816 result[rstart--] = (int)sum;
817 carry = sum >>> 32;
818 }
819
820 // Add remainder of the longer number
821 while (x > 0) {
822 x--;
823 if (carry == 0 && result == value && rstart == (x + offset)) {
824 return;
825 }
826 sum = (value[x+offset] & LONG_MASK) + carry;
827 result[rstart--] = (int)sum;
828 carry = sum >>> 32;
829 }
830 while (y > 0) {
831 y--;
832 int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
833 sum = (bval & LONG_MASK) + carry;
834 result[rstart--] = (int)sum;
835 carry = sum >>> 32;
836 }
837
838 if (carry > 0) { // Result must grow in length
839 resultLen++;
840 if (result.length < resultLen) {
841 int temp[] = new int[resultLen];
842 // Result one word longer from carry-out; copy low-order
843 // bits into new result.
844 System.arraycopy(result, 0, temp, 1, result.length);
845 temp[0] = 1;
846 result = temp;
847 } else {
848 result[rstart--] = 1;
849 }
850 }
851
852 value = result;
853 intLen = resultLen;
854 offset = result.length - resultLen;
855 }
856
857 /**
858 * Like {@link #addShifted(MutableBigInteger, int)} but {@code this.intLen} must
859 * not be greater than {@code n}. In other words, concatenates {@code this}
860 * and {@code addend}.
861 */
862 void addDisjoint(MutableBigInteger addend, int n) {
863 if (addend.isZero())
864 return;
865
866 int x = intLen;
867 int y = addend.intLen + n;
868 int resultLen = (intLen > y ? intLen : y);
869 int[] result;
870 if (value.length < resultLen)
871 result = new int[resultLen];
872 else {
873 result = value;
874 Arrays.fill(value, offset+intLen, value.length, 0);
875 }
876
877 int rstart = result.length-1;
878
879 // copy from this if needed
880 System.arraycopy(value, offset, result, rstart+1-x, x);
881 y -= x;
882 rstart -= x;
883
884 int len = Math.min(y, addend.value.length-addend.offset);
885 System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len);
886
887 // zero the gap
888 for (int i=rstart+1-y+len; i < rstart+1; i++)
889 result[i] = 0;
890
891 value = result;
892 intLen = resultLen;
893 offset = result.length - resultLen;
894 }
895
896 /**
897 * Adds the low {@code n} ints of {@code addend}.
898 */
899 void addLower(MutableBigInteger addend, int n) {
900 MutableBigInteger a = new MutableBigInteger(addend);
901 if (a.offset + a.intLen >= n) {
902 a.offset = a.offset + a.intLen - n;
903 a.intLen = n;
904 }
905 a.normalize();
906 add(a);
907 }
908
909 /**
910 * Subtracts the smaller of this and b from the larger and places the
911 * result into this MutableBigInteger.
912 */
913 int subtract(MutableBigInteger b) {
914 MutableBigInteger a = this;
915
916 int[] result = value;
917 int sign = a.compare(b);
918
919 if (sign == 0) {
920 reset();
921 return 0;
922 }
923 if (sign < 0) {
924 MutableBigInteger tmp = a;
925 a = b;
926 b = tmp;
927 }
928
929 int resultLen = a.intLen;
930 if (result.length < resultLen)
931 result = new int[resultLen];
932
933 long diff = 0;
934 int x = a.intLen;
935 int y = b.intLen;
936 int rstart = result.length - 1;
937
938 // Subtract common parts of both numbers
939 while (y > 0) {
940 x--; y--;
941
942 diff = (a.value[x+a.offset] & LONG_MASK) -
943 (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
944 result[rstart--] = (int)diff;
945 }
946 // Subtract remainder of longer number
947 while (x > 0) {
948 x--;
949 diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
950 result[rstart--] = (int)diff;
951 }
952
953 value = result;
954 intLen = resultLen;
955 offset = value.length - resultLen;
956 normalize();
957 return sign;
958 }
959
960 /**
961 * Subtracts the smaller of a and b from the larger and places the result
962 * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
963 * operation was performed.
964 */
965 private int difference(MutableBigInteger b) {
966 MutableBigInteger a = this;
967 int sign = a.compare(b);
968 if (sign == 0)
969 return 0;
970 if (sign < 0) {
971 MutableBigInteger tmp = a;
972 a = b;
973 b = tmp;
974 }
975
976 long diff = 0;
977 int x = a.intLen;
978 int y = b.intLen;
979
980 // Subtract common parts of both numbers
981 while (y > 0) {
982 x--; y--;
983 diff = (a.value[a.offset+ x] & LONG_MASK) -
984 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
985 a.value[a.offset+x] = (int)diff;
986 }
987 // Subtract remainder of longer number
988 while (x > 0) {
989 x--;
990 diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
991 a.value[a.offset+x] = (int)diff;
992 }
993
994 a.normalize();
995 return sign;
996 }
997
998 /**
999 * Multiply the contents of two MutableBigInteger objects. The result is
1000 * placed into MutableBigInteger z. The contents of y are not changed.
1001 */
1002 void multiply(MutableBigInteger y, MutableBigInteger z) {
1003 int xLen = intLen;
1004 int yLen = y.intLen;
1005 int newLen = xLen + yLen;
1006
1007 // Put z into an appropriate state to receive product
1008 if (z.value.length < newLen)
1009 z.value = new int[newLen];
1010 z.offset = 0;
1011 z.intLen = newLen;
1012
1013 // The first iteration is hoisted out of the loop to avoid extra add
1014 long carry = 0;
1015 for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
1016 long product = (y.value[j+y.offset] & LONG_MASK) *
1017 (value[xLen-1+offset] & LONG_MASK) + carry;
1018 z.value[k] = (int)product;
1019 carry = product >>> 32;
1020 }
1021 z.value[xLen-1] = (int)carry;
1022
1023 // Perform the multiplication word by word
1024 for (int i = xLen-2; i >= 0; i--) {
1025 carry = 0;
1026 for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
1027 long product = (y.value[j+y.offset] & LONG_MASK) *
1028 (value[i+offset] & LONG_MASK) +
1029 (z.value[k] & LONG_MASK) + carry;
1030 z.value[k] = (int)product;
1031 carry = product >>> 32;
1032 }
1033 z.value[i] = (int)carry;
1034 }
1035
1036 // Remove leading zeros from product
1037 z.normalize();
1038 }
1039
1040 /**
1041 * Multiply the contents of this MutableBigInteger by the word y. The
1042 * result is placed into z.
1043 */
1044 void mul(int y, MutableBigInteger z) {
1045 if (y == 1) {
1046 z.copyValue(this);
1047 return;
1048 }
1049
1050 if (y == 0) {
1051 z.clear();
1052 return;
1053 }
1054
1055 // Perform the multiplication word by word
1056 long ylong = y & LONG_MASK;
1057 int[] zval = (z.value.length < intLen+1 ? new int[intLen + 1]
1058 : z.value);
1059 long carry = 0;
1060 for (int i = intLen-1; i >= 0; i--) {
1061 long product = ylong * (value[i+offset] & LONG_MASK) + carry;
1062 zval[i+1] = (int)product;
1063 carry = product >>> 32;
1064 }
1065
1066 if (carry == 0) {
1067 z.offset = 1;
1068 z.intLen = intLen;
1069 } else {
1070 z.offset = 0;
1071 z.intLen = intLen + 1;
1072 zval[0] = (int)carry;
1073 }
1074 z.value = zval;
1075 }
1076
1077 /**
1078 * This method is used for division of an n word dividend by a one word
1079 * divisor. The quotient is placed into quotient. The one word divisor is
1080 * specified by divisor.
1081 *
1082 * @return the remainder of the division is returned.
1083 *
1084 */
1085 int divideOneWord(int divisor, MutableBigInteger quotient) {
1086 long divisorLong = divisor & LONG_MASK;
1087
1088 // Special case of one word dividend
1089 if (intLen == 1) {
1090 long dividendValue = value[offset] & LONG_MASK;
1091 int q = (int) (dividendValue / divisorLong);
1092 int r = (int) (dividendValue - q * divisorLong);
1093 quotient.value[0] = q;
1094 quotient.intLen = (q == 0) ? 0 : 1;
1095 quotient.offset = 0;
1096 return r;
1097 }
1098
1099 if (quotient.value.length < intLen)
1100 quotient.value = new int[intLen];
1101 quotient.offset = 0;
1102 quotient.intLen = intLen;
1103
1104 // Normalize the divisor
1105 int shift = Integer.numberOfLeadingZeros(divisor);
1106
1107 int rem = value[offset];
1108 long remLong = rem & LONG_MASK;
1109 if (remLong < divisorLong) {
1110 quotient.value[0] = 0;
1111 } else {
1112 quotient.value[0] = (int)(remLong / divisorLong);
1113 rem = (int) (remLong - (quotient.value[0] * divisorLong));
1114 remLong = rem & LONG_MASK;
1115 }
1116 int xlen = intLen;
1117 while (--xlen > 0) {
1118 long dividendEstimate = (remLong << 32) |
1119 (value[offset + intLen - xlen] & LONG_MASK);
1120 int q;
1121 if (dividendEstimate >= 0) {
1122 q = (int) (dividendEstimate / divisorLong);
1123 rem = (int) (dividendEstimate - q * divisorLong);
1124 } else {
1125 long tmp = divWord(dividendEstimate, divisor);
1126 q = (int) (tmp & LONG_MASK);
1127 rem = (int) (tmp >>> 32);
1128 }
1129 quotient.value[intLen - xlen] = q;
1130 remLong = rem & LONG_MASK;
1131 }
1132
1133 quotient.normalize();
1134 // Unnormalize
1135 if (shift > 0)
1136 return rem % divisor;
1137 else
1138 return rem;
1139 }
1140
1141 /**
1142 * Calculates the quotient of this div b and places the quotient in the
1143 * provided MutableBigInteger objects and the remainder object is returned.
1144 *
1145 */
1146 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
1147 return divide(b,quotient,true);
1148 }
1149
1150 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1151 if (intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD ||
1152 b.intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
1153 return divideKnuth(b, quotient, needRemainder);
1154 } else {
1155 return divideAndRemainderBurnikelZiegler(b, quotient);
1156 }
1157 }
1158
1159 /**
1160 * @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
1161 */
1162 MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) {
1163 return divideKnuth(b,quotient,true);
1164 }
1165
1166 /**
1167 * Calculates the quotient of this div b and places the quotient in the
1168 * provided MutableBigInteger objects and the remainder object is returned.
1169 *
1170 * Uses Algorithm D in Knuth section 4.3.1.
1171 * Many optimizations to that algorithm have been adapted from the Colin
1172 * Plumb C library.
1173 * It special cases one word divisors for speed. The content of b is not
1174 * changed.
1175 *
1176 */
1177 MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1178 if (b.intLen == 0)
1179 throw new ArithmeticException("BigInteger divide by zero");
1180
1181 // Dividend is zero
1182 if (intLen == 0) {
1183 quotient.intLen = quotient.offset = 0;
1184 return needRemainder ? new MutableBigInteger() : null;
1185 }
1186
1187 int cmp = compare(b);
1188 // Dividend less than divisor
1189 if (cmp < 0) {
1190 quotient.intLen = quotient.offset = 0;
1191 return needRemainder ? new MutableBigInteger(this) : null;
1192 }
1193 // Dividend equal to divisor
1194 if (cmp == 0) {
1195 quotient.value[0] = quotient.intLen = 1;
1196 quotient.offset = 0;
1197 return needRemainder ? new MutableBigInteger() : null;
1198 }
1199
1200 quotient.clear();
1201 // Special case one word divisor
1202 if (b.intLen == 1) {
1203 int r = divideOneWord(b.value[b.offset], quotient);
1204 if(needRemainder) {
1205 if (r == 0)
1206 return new MutableBigInteger();
1207 return new MutableBigInteger(r);
1208 } else {
1209 return null;
1210 }
1211 }
1212
1213 // Cancel common powers of two if we're above the KNUTH_POW2_* thresholds
1214 if (intLen >= KNUTH_POW2_THRESH_LEN) {
1215 int trailingZeroBits = Math.min(getLowestSetBit(), b.getLowestSetBit());
1216 if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS*32) {
1217 MutableBigInteger a = new MutableBigInteger(this);
1218 b = new MutableBigInteger(b);
1219 a.rightShift(trailingZeroBits);
1220 b.rightShift(trailingZeroBits);
1221 MutableBigInteger r = a.divideKnuth(b, quotient);
1222 r.leftShift(trailingZeroBits);
1223 return r;
1224 }
1225 }
1226
1227 return divideMagnitude(b, quotient, needRemainder);
1228 }
1229
1230 /**
1231 * Computes {@code this/b} and {@code this%b} using the
1232 * <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
1233 * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
1234 * The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are
1235 * multiples of 32 bits.<br/>
1236 * {@code this} and {@code b} must be nonnegative.
1237 * @param b the divisor
1238 * @param quotient output parameter for {@code this/b}
1239 * @return the remainder
1240 */
1241 MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) {
1242 int r = intLen;
1243 int s = b.intLen;
1244
1245 if (r < s) {
1246 return this;
1247 } else {
1248 // Unlike Knuth division, we don't check for common powers of two here because
1249 // BZ already runs faster if both numbers contain powers of two and cancelling them has no
1250 // additional benefit.
1251
1252 // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
1253 int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
1254
1255 int j = (s+m-1) / m; // step 2a: j = ceil(s/m)
1256 int n = j * m; // step 2b: block length in 32-bit units
1257 int n32 = 32 * n; // block length in bits
1258 int sigma = Math.max(0, n32 - b.bitLength()); // step 3: sigma = max{T | (2^T)*B < beta^n}
1259 MutableBigInteger bShifted = new MutableBigInteger(b);
1260 bShifted.safeLeftShift(sigma); // step 4a: shift b so its length is a multiple of n
1261 safeLeftShift(sigma); // step 4b: shift this by the same amount
1262
1263 // step 5: t is the number of blocks needed to accommodate this plus one additional bit
1264 int t = (bitLength()+n32) / n32;
1265 if (t < 2) {
1266 t = 2;
1267 }
1268
1269 // step 6: conceptually split this into blocks a[t-1], ..., a[0]
1270 MutableBigInteger a1 = getBlock(t-1, t, n); // the most significant block of this
1271
1272 // step 7: z[t-2] = [a[t-1], a[t-2]]
1273 MutableBigInteger z = getBlock(t-2, t, n); // the second to most significant block
1274 z.addDisjoint(a1, n); // z[t-2]
1275
1276 // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
1277 MutableBigInteger qi = new MutableBigInteger();
1278 MutableBigInteger ri;
1279 quotient.offset = quotient.intLen = 0;
1280 for (int i=t-2; i > 0; i--) {
1281 // step 8a: compute (qi,ri) such that z=b*qi+ri
1282 ri = z.divide2n1n(bShifted, qi);
1283
1284 // step 8b: z = [ri, a[i-1]]
1285 z = getBlock(i-1, t, n); // a[i-1]
1286 z.addDisjoint(ri, n);
1287 quotient.addShifted(qi, i*n); // update q (part of step 9)
1288 }
1289 // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
1290 ri = z.divide2n1n(bShifted, qi);
1291 quotient.add(qi);
1292
1293 ri.rightShift(sigma); // step 9: this and b were shifted, so shift back
1294 return ri;
1295 }
1296 }
1297
1298 /**
1299 * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
1300 * It divides a 2n-digit number by a n-digit number.<br/>
1301 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
1302 * <br/>
1303 * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
1304 * @param b a positive number such that {@code b.bitLength()} is even
1305 * @param quotient output parameter for {@code this/b}
1306 * @return {@code this%b}
1307 */
1308 private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
1309 int n = b.intLen;
1310
1311 // step 1: base case
1312 if (n%2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
1313 return divideKnuth(b, quotient);
1314 }
1315
1316 // step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less
1317 MutableBigInteger aUpper = new MutableBigInteger(this);
1318 aUpper.safeRightShift(32*(n/2)); // aUpper = [a1,a2,a3]
1319 keepLower(n/2); // this = a4
1320
1321 // step 3: q1=aUpper/b, r1=aUpper%b
1322 MutableBigInteger q1 = new MutableBigInteger();
1323 MutableBigInteger r1 = aUpper.divide3n2n(b, q1);
1324
1325 // step 4: quotient=[r1,this]/b, r2=[r1,this]%b
1326 addDisjoint(r1, n/2); // this = [r1,this]
1327 MutableBigInteger r2 = divide3n2n(b, quotient);
1328
1329 // step 5: let quotient=[q1,quotient] and return r2
1330 quotient.addDisjoint(q1, n/2);
1331 return r2;
1332 }
1333
1334 /**
1335 * This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.
1336 * It divides a 3n-digit number by a 2n-digit number.<br/>
1337 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.<br/>
1338 * <br/>
1339 * {@code this} must be a nonnegative number such that {@code 2*this.bitLength() <= 3*b.bitLength()}
1340 * @param quotient output parameter for {@code this/b}
1341 * @return {@code this%b}
1342 */
1343 private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) {
1344 int n = b.intLen / 2; // half the length of b in ints
1345
1346 // step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2]
1347 MutableBigInteger a12 = new MutableBigInteger(this);
1348 a12.safeRightShift(32*n);
1349
1350 // step 2: view b as [b1,b2] where each bi is n ints or less
1351 MutableBigInteger b1 = new MutableBigInteger(b);
1352 b1.safeRightShift(n * 32);
1353 BigInteger b2 = b.getLower(n);
1354
1355 MutableBigInteger r;
1356 MutableBigInteger d;
1357 if (compareShifted(b, n) < 0) {
1358 // step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1
1359 r = a12.divide2n1n(b1, quotient);
1360
1361 // step 4: d=quotient*b2
1362 d = new MutableBigInteger(quotient.toBigInteger().multiply(b2));
1363 } else {
1364 // step 3b: if a1>=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1
1365 quotient.ones(n);
1366 a12.add(b1);
1367 b1.leftShift(32*n);
1368 a12.subtract(b1);
1369 r = a12;
1370
1371 // step 4: d=quotient*b2=(b2 << 32*n) - b2
1372 d = new MutableBigInteger(b2);
1373 d.leftShift(32 * n);
1374 d.subtract(new MutableBigInteger(b2));
1375 }
1376
1377 // step 5: r = r*beta^n + a3 - d (paper says a4)
1378 // However, don't subtract d until after the while loop so r doesn't become negative
1379 r.leftShift(32 * n);
1380 r.addLower(this, n);
1381
1382 // step 6: add b until r>=d
1383 while (r.compare(d) < 0) {
1384 r.add(b);
1385 quotient.subtract(MutableBigInteger.ONE);
1386 }
1387 r.subtract(d);
1388
1389 return r;
1390 }
1391
1392 /**
1393 * Returns a {@code MutableBigInteger} containing {@code blockLength} ints from
1394 * {@code this} number, starting at {@code index*blockLength}.<br/>
1395 * Used by Burnikel-Ziegler division.
1396 * @param index the block index
1397 * @param numBlocks the total number of blocks in {@code this} number
1398 * @param blockLength length of one block in units of 32 bits
1399 * @return
1400 */
1401 private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) {
1402 int blockStart = index * blockLength;
1403 if (blockStart >= intLen) {
1404 return new MutableBigInteger();
1405 }
1406
1407 int blockEnd;
1408 if (index == numBlocks-1) {
1409 blockEnd = intLen;
1410 } else {
1411 blockEnd = (index+1) * blockLength;
1412 }
1413 if (blockEnd > intLen) {
1414 return new MutableBigInteger();
1415 }
1416
1417 int[] newVal = Arrays.copyOfRange(value, offset+intLen-blockEnd, offset+intLen-blockStart);
1418 return new MutableBigInteger(newVal);
1419 }
1420
1421 /** @see BigInteger#bitLength() */
1422 int bitLength() {
1423 if (intLen == 0)
1424 return 0;
1425 return intLen*32 - Integer.numberOfLeadingZeros(value[offset]);
1426 }
1427
1428 /**
1429 * Internally used to calculate the quotient of this div v and places the
1430 * quotient in the provided MutableBigInteger object and the remainder is
1431 * returned.
1432 *
1433 * @return the remainder of the division will be returned.
1434 */
1435 long divide(long v, MutableBigInteger quotient) {
1436 if (v == 0)
1437 throw new ArithmeticException("BigInteger divide by zero");
1438
1439 // Dividend is zero
1440 if (intLen == 0) {
1441 quotient.intLen = quotient.offset = 0;
1442 return 0;
1443 }
1444 if (v < 0)
1445 v = -v;
1446
1447 int d = (int)(v >>> 32);
1448 quotient.clear();
1449 // Special case on word divisor
1450 if (d == 0)
1451 return divideOneWord((int)v, quotient) & LONG_MASK;
1452 else {
1453 return divideLongMagnitude(v, quotient).toLong();
1454 }
1455 }
1456
1457 private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) {
1458 int n2 = 32 - shift;
1459 int c=src[srcFrom];
1460 for (int i=0; i < srcLen-1; i++) {
1461 int b = c;
1462 c = src[++srcFrom];
1463 dst[dstFrom+i] = (b << shift) | (c >>> n2);
1464 }
1465 dst[dstFrom+srcLen-1] = c << shift;
1466 }
1467
1468 /**
1469 * Divide this MutableBigInteger by the divisor.
1470 * The quotient will be placed into the provided quotient object &
1471 * the remainder object is returned.
1472 */
1473 private MutableBigInteger divideMagnitude(MutableBigInteger div,
1474 MutableBigInteger quotient,
1475 boolean needRemainder ) {
1476 // assert div.intLen > 1
1477 // D1 normalize the divisor
1478 int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
1479 // Copy divisor value to protect divisor
1480 final int dlen = div.intLen;
1481 int[] divisor;
1482 MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
1483 if (shift > 0) {
1484 divisor = new int[dlen];
1485 copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
1486 if (Integer.numberOfLeadingZeros(value[offset]) >= shift) {
1487 int[] remarr = new int[intLen + 1];
1488 rem = new MutableBigInteger(remarr);
1489 rem.intLen = intLen;
1490 rem.offset = 1;
1491 copyAndShift(value,offset,intLen,remarr,1,shift);
1492 } else {
1493 int[] remarr = new int[intLen + 2];
1494 rem = new MutableBigInteger(remarr);
1495 rem.intLen = intLen+1;
1496 rem.offset = 1;
1497 int rFrom = offset;
1498 int c=0;
1499 int n2 = 32 - shift;
1500 for (int i=1; i < intLen+1; i++,rFrom++) {
1501 int b = c;
1502 c = value[rFrom];
1503 remarr[i] = (b << shift) | (c >>> n2);
1504 }
1505 remarr[intLen+1] = c << shift;
1506 }
1507 } else {
1508 divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
1509 rem = new MutableBigInteger(new int[intLen + 1]);
1510 System.arraycopy(value, offset, rem.value, 1, intLen);
1511 rem.intLen = intLen;
1512 rem.offset = 1;
1513 }
1514
1515 int nlen = rem.intLen;
1516
1517 // Set the quotient size
1518 final int limit = nlen - dlen + 1;
1519 if (quotient.value.length < limit) {
1520 quotient.value = new int[limit];
1521 quotient.offset = 0;
1522 }
1523 quotient.intLen = limit;
1524 int[] q = quotient.value;
1525
1526
1527 // Must insert leading 0 in rem if its length did not change
1528 if (rem.intLen == nlen) {
1529 rem.offset = 0;
1530 rem.value[0] = 0;
1531 rem.intLen++;
1532 }
1533
1534 int dh = divisor[0];
1535 long dhLong = dh & LONG_MASK;
1536 int dl = divisor[1];
1537
1538 // D2 Initialize j
1539 for (int j=0; j < limit-1; j++) {
1540 // D3 Calculate qhat
1541 // estimate qhat
1542 int qhat = 0;
1543 int qrem = 0;
1544 boolean skipCorrection = false;
1545 int nh = rem.value[j+rem.offset];
1546 int nh2 = nh + 0x80000000;
1547 int nm = rem.value[j+1+rem.offset];
1548
1549 if (nh == dh) {
1550 qhat = ~0;
1551 qrem = nh + nm;
1552 skipCorrection = qrem + 0x80000000 < nh2;
1553 } else {
1554 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1555 if (nChunk >= 0) {
1556 qhat = (int) (nChunk / dhLong);
1557 qrem = (int) (nChunk - (qhat * dhLong));
1558 } else {
1559 long tmp = divWord(nChunk, dh);
1560 qhat = (int) (tmp & LONG_MASK);
1561 qrem = (int) (tmp >>> 32);
1562 }
1563 }
1564
1565 if (qhat == 0)
1566 continue;
1567
1568 if (!skipCorrection) { // Correct qhat
1569 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1570 long rs = ((qrem & LONG_MASK) << 32) | nl;
1571 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1572
1573 if (unsignedLongCompare(estProduct, rs)) {
1574 qhat--;
1575 qrem = (int)((qrem & LONG_MASK) + dhLong);
1576 if ((qrem & LONG_MASK) >= dhLong) {
1577 estProduct -= (dl & LONG_MASK);
1578 rs = ((qrem & LONG_MASK) << 32) | nl;
1579 if (unsignedLongCompare(estProduct, rs))
1580 qhat--;
1581 }
1582 }
1583 }
1584
1585 // D4 Multiply and subtract
1586 rem.value[j+rem.offset] = 0;
1587 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1588
1589 // D5 Test remainder
1590 if (borrow + 0x80000000 > nh2) {
1591 // D6 Add back
1592 divadd(divisor, rem.value, j+1+rem.offset);
1593 qhat--;
1594 }
1595
1596 // Store the quotient digit
1597 q[j] = qhat;
1598 } // D7 loop on j
1599 // D3 Calculate qhat
1600 // estimate qhat
1601 int qhat = 0;
1602 int qrem = 0;
1603 boolean skipCorrection = false;
1604 int nh = rem.value[limit - 1 + rem.offset];
1605 int nh2 = nh + 0x80000000;
1606 int nm = rem.value[limit + rem.offset];
1607
1608 if (nh == dh) {
1609 qhat = ~0;
1610 qrem = nh + nm;
1611 skipCorrection = qrem + 0x80000000 < nh2;
1612 } else {
1613 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1614 if (nChunk >= 0) {
1615 qhat = (int) (nChunk / dhLong);
1616 qrem = (int) (nChunk - (qhat * dhLong));
1617 } else {
1618 long tmp = divWord(nChunk, dh);
1619 qhat = (int) (tmp & LONG_MASK);
1620 qrem = (int) (tmp >>> 32);
1621 }
1622 }
1623 if (qhat != 0) {
1624 if (!skipCorrection) { // Correct qhat
1625 long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
1626 long rs = ((qrem & LONG_MASK) << 32) | nl;
1627 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1628
1629 if (unsignedLongCompare(estProduct, rs)) {
1630 qhat--;
1631 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1632 if ((qrem & LONG_MASK) >= dhLong) {
1633 estProduct -= (dl & LONG_MASK);
1634 rs = ((qrem & LONG_MASK) << 32) | nl;
1635 if (unsignedLongCompare(estProduct, rs))
1636 qhat--;
1637 }
1638 }
1639 }
1640
1641
1642 // D4 Multiply and subtract
1643 int borrow;
1644 rem.value[limit - 1 + rem.offset] = 0;
1645 if(needRemainder)
1646 borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1647 else
1648 borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1649
1650 // D5 Test remainder
1651 if (borrow + 0x80000000 > nh2) {
1652 // D6 Add back
1653 if(needRemainder)
1654 divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
1655 qhat--;
1656 }
1657
1658 // Store the quotient digit
1659 q[(limit - 1)] = qhat;
1660 }
1661
1662
1663 if (needRemainder) {
1664 // D8 Unnormalize
1665 if (shift > 0)
1666 rem.rightShift(shift);
1667 rem.normalize();
1668 }
1669 quotient.normalize();
1670 return needRemainder ? rem : null;
1671 }
1672
1673 /**
1674 * Divide this MutableBigInteger by the divisor represented by positive long
1675 * value. The quotient will be placed into the provided quotient object &
1676 * the remainder object is returned.
1677 */
1678 private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
1679 // Remainder starts as dividend with space for a leading zero
1680 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
1681 System.arraycopy(value, offset, rem.value, 1, intLen);
1682 rem.intLen = intLen;
1683 rem.offset = 1;
1684
1685 int nlen = rem.intLen;
1686
1687 int limit = nlen - 2 + 1;
1688 if (quotient.value.length < limit) {
1689 quotient.value = new int[limit];
1690 quotient.offset = 0;
1691 }
1692 quotient.intLen = limit;
1693 int[] q = quotient.value;
1694
1695 // D1 normalize the divisor
1696 int shift = Long.numberOfLeadingZeros(ldivisor);
1697 if (shift > 0) {
1698 ldivisor<<=shift;
1699 rem.leftShift(shift);
1700 }
1701
1702 // Must insert leading 0 in rem if its length did not change
1703 if (rem.intLen == nlen) {
1704 rem.offset = 0;
1705 rem.value[0] = 0;
1706 rem.intLen++;
1707 }
1708
1709 int dh = (int)(ldivisor >>> 32);
1710 long dhLong = dh & LONG_MASK;
1711 int dl = (int)(ldivisor & LONG_MASK);
1712
1713 // D2 Initialize j
1714 for (int j = 0; j < limit; j++) {
1715 // D3 Calculate qhat
1716 // estimate qhat
1717 int qhat = 0;
1718 int qrem = 0;
1719 boolean skipCorrection = false;
1720 int nh = rem.value[j + rem.offset];
1721 int nh2 = nh + 0x80000000;
1722 int nm = rem.value[j + 1 + rem.offset];
1723
1724 if (nh == dh) {
1725 qhat = ~0;
1726 qrem = nh + nm;
1727 skipCorrection = qrem + 0x80000000 < nh2;
1728 } else {
1729 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1730 if (nChunk >= 0) {
1731 qhat = (int) (nChunk / dhLong);
1732 qrem = (int) (nChunk - (qhat * dhLong));
1733 } else {
1734 long tmp = divWord(nChunk, dh);
1735 qhat =(int)(tmp & LONG_MASK);
1736 qrem = (int)(tmp>>>32);
1737 }
1738 }
1739
1740 if (qhat == 0)
1741 continue;
1742
1743 if (!skipCorrection) { // Correct qhat
1744 long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
1745 long rs = ((qrem & LONG_MASK) << 32) | nl;
1746 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1747
1748 if (unsignedLongCompare(estProduct, rs)) {
1749 qhat--;
1750 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1751 if ((qrem & LONG_MASK) >= dhLong) {
1752 estProduct -= (dl & LONG_MASK);
1753 rs = ((qrem & LONG_MASK) << 32) | nl;
1754 if (unsignedLongCompare(estProduct, rs))
1755 qhat--;
1756 }
1757 }
1758 }
1759
1760 // D4 Multiply and subtract
1761 rem.value[j + rem.offset] = 0;
1762 int borrow = mulsubLong(rem.value, dh, dl, qhat, j + rem.offset);
1763
1764 // D5 Test remainder
1765 if (borrow + 0x80000000 > nh2) {
1766 // D6 Add back
1767 divaddLong(dh,dl, rem.value, j + 1 + rem.offset);
1768 qhat--;
1769 }
1770
1771 // Store the quotient digit
1772 q[j] = qhat;
1773 } // D7 loop on j
1774
1775 // D8 Unnormalize
1776 if (shift > 0)
1777 rem.rightShift(shift);
1778
1779 quotient.normalize();
1780 rem.normalize();
1781 return rem;
1782 }
1783
1784 /**
1785 * A primitive used for division by long.
1786 * Specialized version of the method divadd.
1787 * dh is a high part of the divisor, dl is a low part
1788 */
1789 private int divaddLong(int dh, int dl, int[] result, int offset) {
1790 long carry = 0;
1791
1792 long sum = (dl & LONG_MASK) + (result[1+offset] & LONG_MASK);
1793 result[1+offset] = (int)sum;
1794
1795 sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry;
1796 result[offset] = (int)sum;
1797 carry = sum >>> 32;
1798 return (int)carry;
1799 }
1800
1801 /**
1802 * This method is used for division by long.
1803 * Specialized version of the method sulsub.
1804 * dh is a high part of the divisor, dl is a low part
1805 */
1806 private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
1807 long xLong = x & LONG_MASK;
1808 offset += 2;
1809 long product = (dl & LONG_MASK) * xLong;
1810 long difference = q[offset] - product;
1811 q[offset--] = (int)difference;
1812 long carry = (product >>> 32)
1813 + (((difference & LONG_MASK) >
1814 (((~(int)product) & LONG_MASK))) ? 1:0);
1815 product = (dh & LONG_MASK) * xLong + carry;
1816 difference = q[offset] - product;
1817 q[offset--] = (int)difference;
1818 carry = (product >>> 32)
1819 + (((difference & LONG_MASK) >
1820 (((~(int)product) & LONG_MASK))) ? 1:0);
1821 return (int)carry;
1822 }
1823
1824 /**
1825 * Compare two longs as if they were unsigned.
1826 * Returns true iff one is bigger than two.
1827 */
1828 private boolean unsignedLongCompare(long one, long two) {
1829 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1830 }
1831
1832 /**
1833 * This method divides a long quantity by an int to estimate
1834 * qhat for two multi precision numbers. It is used when
1835 * the signed value of n is less than zero.
1836 * Returns long value where high 32 bits contain remainder value and
1837 * low 32 bits contain quotient value.
1838 */
1839 static long divWord(long n, int d) {
1840 long dLong = d & LONG_MASK;
1841 long r;
1842 long q;
1843 if (dLong == 1) {
1844 q = (int)n;
1845 r = 0;
1846 return (r << 32) | (q & LONG_MASK);
1847 }
1848
1849 // Approximate the quotient and remainder
1850 q = (n >>> 1) / (dLong >>> 1);
1851 r = n - q*dLong;
1852
1853 // Correct the approximation
1854 while (r < 0) {
1855 r += dLong;
1856 q--;
1857 }
1858 while (r >= dLong) {
1859 r -= dLong;
1860 q++;
1861 }
1862 // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1863 return (r << 32) | (q & LONG_MASK);
1864 }
1865
1866 /**
1867 * Calculate GCD of this and b. This and b are changed by the computation.
1868 */
1869 MutableBigInteger hybridGCD(MutableBigInteger b) {
1870 // Use Euclid's algorithm until the numbers are approximately the
1871 // same length, then use the binary GCD algorithm to find the GCD.
1872 MutableBigInteger a = this;
1873 MutableBigInteger q = new MutableBigInteger();
1874
1875 while (b.intLen != 0) {
1876 if (Math.abs(a.intLen - b.intLen) < 2)
1877 return a.binaryGCD(b);
1878
1879 MutableBigInteger r = a.divide(b, q);
1880 a = b;
1881 b = r;
1882 }
1883 return a;
1884 }
1885
1886 /**
1887 * Calculate GCD of this and v.
1888 * Assumes that this and v are not zero.
1889 */
1890 private MutableBigInteger binaryGCD(MutableBigInteger v) {
1891 // Algorithm B from Knuth section 4.5.2
1892 MutableBigInteger u = this;
1893 MutableBigInteger r = new MutableBigInteger();
1894
1895 // step B1
1896 int s1 = u.getLowestSetBit();
1897 int s2 = v.getLowestSetBit();
1898 int k = (s1 < s2) ? s1 : s2;
1899 if (k != 0) {
1900 u.rightShift(k);
1901 v.rightShift(k);
1902 }
1903
1904 // step B2
1905 boolean uOdd = (k == s1);
1906 MutableBigInteger t = uOdd ? v: u;
1907 int tsign = uOdd ? -1 : 1;
1908
1909 int lb;
1910 while ((lb = t.getLowestSetBit()) >= 0) {
1911 // steps B3 and B4
1912 t.rightShift(lb);
1913 // step B5
1914 if (tsign > 0)
1915 u = t;
1916 else
1917 v = t;
1918
1919 // Special case one word numbers
1920 if (u.intLen < 2 && v.intLen < 2) {
1921 int x = u.value[u.offset];
1922 int y = v.value[v.offset];
1923 x = binaryGcd(x, y);
1924 r.value[0] = x;
1925 r.intLen = 1;
1926 r.offset = 0;
1927 if (k > 0)
1928 r.leftShift(k);
1929 return r;
1930 }
1931
1932 // step B6
1933 if ((tsign = u.difference(v)) == 0)
1934 break;
1935 t = (tsign >= 0) ? u : v;
1936 }
1937
1938 if (k > 0)
1939 u.leftShift(k);
1940 return u;
1941 }
1942
1943 /**
1944 * Calculate GCD of a and b interpreted as unsigned integers.
1945 */
1946 static int binaryGcd(int a, int b) {
1947 if (b == 0)
1948 return a;
1949 if (a == 0)
1950 return b;
1951
1952 // Right shift a & b till their last bits equal to 1.
1953 int aZeros = Integer.numberOfTrailingZeros(a);
1954 int bZeros = Integer.numberOfTrailingZeros(b);
1955 a >>>= aZeros;
1956 b >>>= bZeros;
1957
1958 int t = (aZeros < bZeros ? aZeros : bZeros);
1959
1960 while (a != b) {
1961 if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned
1962 a -= b;
1963 a >>>= Integer.numberOfTrailingZeros(a);
1964 } else {
1965 b -= a;
1966 b >>>= Integer.numberOfTrailingZeros(b);
1967 }
1968 }
1969 return a<<t;
1970 }
1971
1972 /**
1973 * Returns the modInverse of this mod p. This and p are not affected by
1974 * the operation.
1975 */
1976 MutableBigInteger mutableModInverse(MutableBigInteger p) {
1977 // Modulus is odd, use Schroeppel's algorithm
1978 if (p.isOdd())
1979 return modInverse(p);
1980
1981 // Base and modulus are even, throw exception
1982 if (isEven())
1983 throw new ArithmeticException("BigInteger not invertible.");
1984
1985 // Get even part of modulus expressed as a power of 2
1986 int powersOf2 = p.getLowestSetBit();
1987
1988 // Construct odd part of modulus
1989 MutableBigInteger oddMod = new MutableBigInteger(p);
1990 oddMod.rightShift(powersOf2);
1991
1992 if (oddMod.isOne())
1993 return modInverseMP2(powersOf2);
1994
1995 // Calculate 1/a mod oddMod
1996 MutableBigInteger oddPart = modInverse(oddMod);
1997
1998 // Calculate 1/a mod evenMod
1999 MutableBigInteger evenPart = modInverseMP2(powersOf2);
2000
2001 // Combine the results using Chinese Remainder Theorem
2002 MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
2003 MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
2004
2005 MutableBigInteger temp1 = new MutableBigInteger();
2006 MutableBigInteger temp2 = new MutableBigInteger();
2007 MutableBigInteger result = new MutableBigInteger();
2008
2009 oddPart.leftShift(powersOf2);
2010 oddPart.multiply(y1, result);
2011
2012 evenPart.multiply(oddMod, temp1);
2013 temp1.multiply(y2, temp2);
2014
2015 result.add(temp2);
2016 return result.divide(p, temp1);
2017 }
2018
2019 /*
2020 * Calculate the multiplicative inverse of this mod 2^k.
2021 */
2022 MutableBigInteger modInverseMP2(int k) {
2023 if (isEven())
2024 throw new ArithmeticException("Non-invertible. (GCD != 1)");
2025
2026 if (k > 64)
2027 return euclidModInverse(k);
2028
2029 int t = inverseMod32(value[offset+intLen-1]);
2030
2031 if (k < 33) {
2032 t = (k == 32 ? t : t & ((1 << k) - 1));
2033 return new MutableBigInteger(t);
2034 }
2035
2036 long pLong = (value[offset+intLen-1] & LONG_MASK);
2037 if (intLen > 1)
2038 pLong |= ((long)value[offset+intLen-2] << 32);
2039 long tLong = t & LONG_MASK;
2040 tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
2041 tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
2042
2043 MutableBigInteger result = new MutableBigInteger(new int[2]);
2044 result.value[0] = (int)(tLong >>> 32);
2045 result.value[1] = (int)tLong;
2046 result.intLen = 2;
2047 result.normalize();
2048 return result;
2049 }
2050
2051 /**
2052 * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
2053 */
2054 static int inverseMod32(int val) {
2055 // Newton's iteration!
2056 int t = val;
2057 t *= 2 - val*t;
2058 t *= 2 - val*t;
2059 t *= 2 - val*t;
2060 t *= 2 - val*t;
2061 return t;
2062 }
2063
2064 /**
2065 * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
2066 */
2067 static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
2068 // Copy the mod to protect original
2069 return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
2070 }
2071
2072 /**
2073 * Calculate the multiplicative inverse of this mod mod, where mod is odd.
2074 * This and mod are not changed by the calculation.
2075 *
2076 * This method implements an algorithm due to Richard Schroeppel, that uses
2077 * the same intermediate representation as Montgomery Reduction
2078 * ("Montgomery Form"). The algorithm is described in an unpublished
2079 * manuscript entitled "Fast Modular Reciprocals."
2080 */
2081 private MutableBigInteger modInverse(MutableBigInteger mod) {
2082 MutableBigInteger p = new MutableBigInteger(mod);
2083 MutableBigInteger f = new MutableBigInteger(this);
2084 MutableBigInteger g = new MutableBigInteger(p);
2085 SignedMutableBigInteger c = new SignedMutableBigInteger(1);
2086 SignedMutableBigInteger d = new SignedMutableBigInteger();
2087 MutableBigInteger temp = null;
2088 SignedMutableBigInteger sTemp = null;
2089
2090 int k = 0;
2091 // Right shift f k times until odd, left shift d k times
2092 if (f.isEven()) {
2093 int trailingZeros = f.getLowestSetBit();
2094 f.rightShift(trailingZeros);
2095 d.leftShift(trailingZeros);
2096 k = trailingZeros;
2097 }
2098
2099 // The Almost Inverse Algorithm
2100 while (!f.isOne()) {
2101 // If gcd(f, g) != 1, number is not invertible modulo mod
2102 if (f.isZero())
2103 throw new ArithmeticException("BigInteger not invertible.");
2104
2105 // If f < g exchange f, g and c, d
2106 if (f.compare(g) < 0) {
2107 temp = f; f = g; g = temp;
2108 sTemp = d; d = c; c = sTemp;
2109 }
2110
2111 // If f == g (mod 4)
2112 if (((f.value[f.offset + f.intLen - 1] ^
2113 g.value[g.offset + g.intLen - 1]) & 3) == 0) {
2114 f.subtract(g);
2115 c.signedSubtract(d);
2116 } else { // If f != g (mod 4)
2117 f.add(g);
2118 c.signedAdd(d);
2119 }
2120
2121 // Right shift f k times until odd, left shift d k times
2122 int trailingZeros = f.getLowestSetBit();
2123 f.rightShift(trailingZeros);
2124 d.leftShift(trailingZeros);
2125 k += trailingZeros;
2126 }
2127
2128 while (c.sign < 0)
2129 c.signedAdd(p);
2130
2131 return fixup(c, p, k);
2132 }
2133
2134 /**
2135 * The Fixup Algorithm
2136 * Calculates X such that X = C * 2^(-k) (mod P)
2137 * Assumes C<P and P is odd.
2138 */
2139 static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
2140 int k) {
2141 MutableBigInteger temp = new MutableBigInteger();
2142 // Set r to the multiplicative inverse of p mod 2^32
2143 int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
2144
2145 for (int i=0, numWords = k >> 5; i < numWords; i++) {
2146 // V = R * c (mod 2^j)
2147 int v = r * c.value[c.offset + c.intLen-1];
2148 // c = c + (v * p)
2149 p.mul(v, temp);
2150 c.add(temp);
2151 // c = c / 2^j
2152 c.intLen--;
2153 }
2154 int numBits = k & 0x1f;
2155 if (numBits != 0) {
2156 // V = R * c (mod 2^j)
2157 int v = r * c.value[c.offset + c.intLen-1];
2158 v &= ((1<<numBits) - 1);
2159 // c = c + (v * p)
2160 p.mul(v, temp);
2161 c.add(temp);
2162 // c = c / 2^j
2163 c.rightShift(numBits);
2164 }
2165
2166 // In theory, c may be greater than p at this point (Very rare!)
2167 while (c.compare(p) >= 0)
2168 c.subtract(p);
2169
2170 return c;
2171 }
2172
2173 /**
2174 * Uses the extended Euclidean algorithm to compute the modInverse of base
2175 * mod a modulus that is a power of 2. The modulus is 2^k.
2176 */
2177 MutableBigInteger euclidModInverse(int k) {
2178 MutableBigInteger b = new MutableBigInteger(1);
2179 b.leftShift(k);
2180 MutableBigInteger mod = new MutableBigInteger(b);
2181
2182 MutableBigInteger a = new MutableBigInteger(this);
2183 MutableBigInteger q = new MutableBigInteger();
2184 MutableBigInteger r = b.divide(a, q);
2185
2186 MutableBigInteger swapper = b;
2187 // swap b & r
2188 b = r;
2189 r = swapper;
2190
2191 MutableBigInteger t1 = new MutableBigInteger(q);
2192 MutableBigInteger t0 = new MutableBigInteger(1);
2193 MutableBigInteger temp = new MutableBigInteger();
2194
2195 while (!b.isOne()) {
2196 r = a.divide(b, q);
2197
2198 if (r.intLen == 0)
2199 throw new ArithmeticException("BigInteger not invertible.");
2200
2201 swapper = r;
2202 a = swapper;
2203
2204 if (q.intLen == 1)
2205 t1.mul(q.value[q.offset], temp);
2206 else
2207 q.multiply(t1, temp);
2208 swapper = q;
2209 q = temp;
2210 temp = swapper;
2211 t0.add(q);
2212
2213 if (a.isOne())
2214 return t0;
2215
2216 r = b.divide(a, q);
2217
2218 if (r.intLen == 0)
2219 throw new ArithmeticException("BigInteger not invertible.");
2220
2221 swapper = b;
2222 b = r;
2223
2224 if (q.intLen == 1)
2225 t0.mul(q.value[q.offset], temp);
2226 else
2227 q.multiply(t0, temp);
2228 swapper = q; q = temp; temp = swapper;
2229
2230 t1.add(q);
2231 }
2232 mod.subtract(t1);
2233 return mod;
2234 }
2235 }