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 * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
534 * in normal form.
535 */
536 void rightShift(int n) {
537 if (intLen == 0)
538 return;
539 int nInts = n >>> 5;
540 int nBits = n & 0x1F;
541 this.intLen -= nInts;
542 if (nBits == 0)
543 return;
544 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
545 if (nBits >= bitsInHighWord) {
546 this.primitiveLeftShift(32 - nBits);
547 this.intLen--;
548 } else {
549 primitiveRightShift(nBits);
550 }
551 }
552
553 /**
554 * Like {@link #leftShift(int)} but {@code n} can be zero.
555 */
556 void safeLeftShift(int n) {
557 if (n > 0)
558 leftShift(n);
559 }
560
561 /**
562 * Left shift this MutableBigInteger n bits.
563 */
564 void leftShift(int n) {
565 /*
566 * If there is enough storage space in this MutableBigInteger already
567 * the available space will be used. Space to the right of the used
568 * ints in the value array is faster to utilize, so the extra space
569 * will be taken from the right if possible.
570 */
571 if (intLen == 0)
572 return;
573 int nInts = n >>> 5;
574 int nBits = n&0x1F;
575 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
576
577 // If shift can be done without moving words, do so
578 if (n <= (32-bitsInHighWord)) {
579 primitiveLeftShift(nBits);
580 return;
581 }
582
583 int newLen = intLen + nInts +1;
584 if (nBits <= (32-bitsInHighWord))
585 newLen--;
586 if (value.length < newLen) {
587 // The array must grow
588 int[] result = new int[newLen];
589 for (int i=0; i<intLen; i++)
590 result[i] = value[offset+i];
591 setValue(result, newLen);
592 } else if (value.length - offset >= newLen) {
593 // Use space on right
594 for(int i=0; i<newLen - intLen; i++)
595 value[offset+intLen+i] = 0;
596 } else {
597 // Must use space on left
598 for (int i=0; i<intLen; i++)
599 value[i] = value[offset+i];
600 for (int i=intLen; i<newLen; i++)
601 value[i] = 0;
602 offset = 0;
603 }
604 intLen = newLen;
605 if (nBits == 0)
606 return;
607 if (nBits <= (32-bitsInHighWord))
608 primitiveLeftShift(nBits);
609 else
610 primitiveRightShift(32 -nBits);
611 }
612
613 /**
614 * A primitive used for division. This method adds in one multiple of the
615 * divisor a back to the dividend result at a specified offset. It is used
616 * when qhat was estimated too large, and must be adjusted.
617 */
618 private int divadd(int[] a, int[] result, int offset) {
619 long carry = 0;
620
621 for (int j=a.length-1; j >= 0; j--) {
622 long sum = (a[j] & LONG_MASK) +
623 (result[j+offset] & LONG_MASK) + carry;
624 result[j+offset] = (int)sum;
625 carry = sum >>> 32;
626 }
627 return (int)carry;
628 }
629
630 /**
631 * This method is used for division. It multiplies an n word input a by one
632 * word input x, and subtracts the n word product from q. This is needed
633 * when subtracting qhat*divisor from dividend.
634 */
635 private int mulsub(int[] q, int[] a, int x, int len, int offset) {
636 long xLong = x & LONG_MASK;
637 long carry = 0;
638 offset += len;
639
640 for (int j=len-1; j >= 0; j--) {
641 long product = (a[j] & LONG_MASK) * xLong + carry;
642 long difference = q[offset] - product;
643 q[offset--] = (int)difference;
644 carry = (product >>> 32)
645 + (((difference & LONG_MASK) >
646 (((~(int)product) & LONG_MASK))) ? 1:0);
647 }
648 return (int)carry;
649 }
650
651 /**
652 * The method is the same as mulsun, except the fact that q array is not
653 * updated, the only result of the method is borrow flag.
654 */
655 private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
656 long xLong = x & LONG_MASK;
657 long carry = 0;
658 offset += len;
659 for (int j=len-1; j >= 0; j--) {
660 long product = (a[j] & LONG_MASK) * xLong + carry;
661 long difference = q[offset--] - product;
662 carry = (product >>> 32)
663 + (((difference & LONG_MASK) >
664 (((~(int)product) & LONG_MASK))) ? 1:0);
665 }
666 return (int)carry;
667 }
668
669 /**
670 * Right shift this MutableBigInteger n bits, where n is
671 * less than 32.
672 * Assumes that intLen > 0, n > 0 for speed
673 */
674 private final void primitiveRightShift(int n) {
675 int[] val = value;
676 int n2 = 32 - n;
677 for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
678 int b = c;
679 c = val[i-1];
680 val[i] = (c << n2) | (b >>> n);
681 }
682 val[offset] >>>= n;
683 }
684
685 /**
686 * Left shift this MutableBigInteger n bits, where n is
687 * less than 32.
688 * Assumes that intLen > 0, n > 0 for speed
689 */
690 private final void primitiveLeftShift(int n) {
691 int[] val = value;
692 int n2 = 32 - n;
693 for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) {
694 int b = c;
695 c = val[i+1];
696 val[i] = (b << n) | (c >>> n2);
697 }
698 val[offset+intLen-1] <<= n;
699 }
700
701 /**
702 * Returns a {@code BigInteger} equal to the {@code n}
703 * low ints of this number.
704 */
705 private BigInteger getLower(int n) {
706 if (isZero())
707 return BigInteger.ZERO;
708 else if (intLen < n)
709 return toBigInteger(1);
710 else {
711 // strip zeros
712 int len = n;
713 while (len>0 && value[offset+intLen-len]==0)
714 len--;
715 int sign = len>0 ? 1 : 0;
716 return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
717 }
718 }
719
720 /**
721 * Discards all ints whose index is greater than {@code n}.
722 */
723 private void keepLower(int n) {
724 if (intLen >= n) {
725 offset += intLen - n;
726 intLen = n;
727 }
728 }
729
730 /**
731 * Adds the contents of two MutableBigInteger objects.The result
732 * is placed within this MutableBigInteger.
733 * The contents of the addend are not changed.
734 */
735 void add(MutableBigInteger addend) {
736 int x = intLen;
737 int y = addend.intLen;
738 int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
739 int[] result = (value.length < resultLen ? new int[resultLen] : value);
740
741 int rstart = result.length-1;
742 long sum;
743 long carry = 0;
744
745 // Add common parts of both numbers
746 while(x>0 && y>0) {
747 x--; y--;
748 sum = (value[x+offset] & LONG_MASK) +
749 (addend.value[y+addend.offset] & LONG_MASK) + carry;
750 result[rstart--] = (int)sum;
751 carry = sum >>> 32;
752 }
753
754 // Add remainder of the longer number
755 while(x>0) {
756 x--;
757 if (carry == 0 && result == value && rstart == (x + offset))
758 return;
759 sum = (value[x+offset] & LONG_MASK) + carry;
760 result[rstart--] = (int)sum;
761 carry = sum >>> 32;
762 }
763 while(y>0) {
764 y--;
765 sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
766 result[rstart--] = (int)sum;
767 carry = sum >>> 32;
768 }
769
770 if (carry > 0) { // Result must grow in length
771 resultLen++;
772 if (result.length < resultLen) {
773 int temp[] = new int[resultLen];
774 // Result one word longer from carry-out; copy low-order
775 // bits into new result.
776 System.arraycopy(result, 0, temp, 1, result.length);
777 temp[0] = 1;
778 result = temp;
779 } else {
780 result[rstart--] = 1;
781 }
782 }
783
784 value = result;
785 intLen = resultLen;
786 offset = result.length - resultLen;
787 }
788
789 /**
790 * Adds the value of {@code addend} shifted {@code n} ints to the left.
791 * Has the same effect as {@code addend.leftShift(32*ints); add(b);}
792 * but doesn't change the value of {@code b}.
793 */
794 void addShifted(MutableBigInteger addend, int n) {
795 if (addend.isZero())
796 return;
797
798 int x = intLen;
799 int y = addend.intLen + n;
800 int resultLen = (intLen > y ? intLen : y);
801 int[] result = (value.length < resultLen ? new int[resultLen] : value);
802
803 int rstart = result.length-1;
804 long sum;
805 long carry = 0;
806
807 // Add common parts of both numbers
808 while(x>0 && y>0) {
809 x--; y--;
810 int bval = y+addend.offset<addend.value.length ? addend.value[y+addend.offset] : 0;
811 sum = (value[x+offset] & LONG_MASK) +
812 (bval & LONG_MASK) + carry;
813 result[rstart--] = (int)sum;
814 carry = sum >>> 32;
815 }
816
817 // Add remainder of the longer number
818 while(x>0) {
819 x--;
820 if (carry == 0 && result == value && rstart == (x + offset))
821 return;
822 sum = (value[x+offset] & LONG_MASK) + carry;
823 result[rstart--] = (int)sum;
824 carry = sum >>> 32;
825 }
826 while(y>0) {
827 y--;
828 int bval = y+addend.offset<addend.value.length ? addend.value[y+addend.offset] : 0;
829 sum = (bval & LONG_MASK) + carry;
830 result[rstart--] = (int)sum;
831 carry = sum >>> 32;
832 }
833
834 if (carry > 0) { // Result must grow in length
835 resultLen++;
836 if (result.length < resultLen) {
837 int temp[] = new int[resultLen];
838 // Result one word longer from carry-out; copy low-order
839 // bits into new result.
840 System.arraycopy(result, 0, temp, 1, result.length);
841 temp[0] = 1;
842 result = temp;
843 } else {
844 result[rstart--] = 1;
845 }
846 }
847
848 value = result;
849 intLen = resultLen;
850 offset = result.length - resultLen;
851 }
852
853 /**
854 * Like {@link #addShifted(MutableBigInteger, int)} but {@code this.intLen} must
855 * not be greater than {@code n}. In other words, concatenates {@code this}
856 * and {@code addend}.
857 */
858 void addDisjoint(MutableBigInteger addend, int n) {
859 if (addend.isZero())
860 return;
861
862 int x = intLen;
863 int y = addend.intLen + n;
864 int resultLen = (intLen > y ? intLen : y);
865 int[] result;
866 if (value.length < resultLen)
867 result = new int[resultLen];
868 else {
869 result = value;
870 Arrays.fill(value, offset+intLen, value.length, 0);
871 }
872
873 int rstart = result.length-1;
874
875 // copy from this if needed
876 System.arraycopy(value, offset, result, rstart+1-x, x);
877 y -= x;
878 rstart -= x;
879
880 int len = Math.min(y, addend.value.length-addend.offset);
881 System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len);
882
883 // zero the gap
884 for (int i=rstart+1-y+len; i<rstart+1; i++)
885 result[i] = 0;
886
887 value = result;
888 intLen = resultLen;
889 offset = result.length - resultLen;
890 }
891
892 /**
893 * Adds the low {@code n} ints of {@code addend}.
894 */
895 void addLower(MutableBigInteger addend, int n) {
896 MutableBigInteger a = new MutableBigInteger(addend);
897 if (a.offset + a.intLen >= n) {
898 a.offset = a.offset + a.intLen - n;
899 a.intLen = n;
900 }
901 a.normalize();
902 add(a);
903 }
904
905 /**
906 * Subtracts the smaller of this and b from the larger and places the
907 * result into this MutableBigInteger.
908 */
909 int subtract(MutableBigInteger b) {
910 MutableBigInteger a = this;
911
912 int[] result = value;
913 int sign = a.compare(b);
914
915 if (sign == 0) {
916 reset();
917 return 0;
918 }
919 if (sign < 0) {
920 MutableBigInteger tmp = a;
921 a = b;
922 b = tmp;
923 }
924
925 int resultLen = a.intLen;
926 if (result.length < resultLen)
927 result = new int[resultLen];
928
929 long diff = 0;
930 int x = a.intLen;
931 int y = b.intLen;
932 int rstart = result.length - 1;
933
934 // Subtract common parts of both numbers
935 while (y>0) {
936 x--; y--;
937
938 diff = (a.value[x+a.offset] & LONG_MASK) -
939 (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
940 result[rstart--] = (int)diff;
941 }
942 // Subtract remainder of longer number
943 while (x>0) {
944 x--;
945 diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
946 result[rstart--] = (int)diff;
947 }
948
949 value = result;
950 intLen = resultLen;
951 offset = value.length - resultLen;
952 normalize();
953 return sign;
954 }
955
956 /**
957 * Subtracts the smaller of a and b from the larger and places the result
958 * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
959 * operation was performed.
960 */
961 private int difference(MutableBigInteger b) {
962 MutableBigInteger a = this;
963 int sign = a.compare(b);
964 if (sign ==0)
965 return 0;
966 if (sign < 0) {
967 MutableBigInteger tmp = a;
968 a = b;
969 b = tmp;
970 }
971
972 long diff = 0;
973 int x = a.intLen;
974 int y = b.intLen;
975
976 // Subtract common parts of both numbers
977 while (y>0) {
978 x--; y--;
979 diff = (a.value[a.offset+ x] & LONG_MASK) -
980 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
981 a.value[a.offset+x] = (int)diff;
982 }
983 // Subtract remainder of longer number
984 while (x>0) {
985 x--;
986 diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
987 a.value[a.offset+x] = (int)diff;
988 }
989
990 a.normalize();
991 return sign;
992 }
993
994 /**
995 * Multiply the contents of two MutableBigInteger objects. The result is
996 * placed into MutableBigInteger z. The contents of y are not changed.
997 */
998 void multiply(MutableBigInteger y, MutableBigInteger z) {
999 int xLen = intLen;
1000 int yLen = y.intLen;
1001 int newLen = xLen + yLen;
1002
1003 // Put z into an appropriate state to receive product
1004 if (z.value.length < newLen)
1005 z.value = new int[newLen];
1006 z.offset = 0;
1007 z.intLen = newLen;
1008
1009 // The first iteration is hoisted out of the loop to avoid extra add
1010 long carry = 0;
1011 for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
1012 long product = (y.value[j+y.offset] & LONG_MASK) *
1013 (value[xLen-1+offset] & LONG_MASK) + carry;
1014 z.value[k] = (int)product;
1015 carry = product >>> 32;
1016 }
1017 z.value[xLen-1] = (int)carry;
1018
1019 // Perform the multiplication word by word
1020 for (int i = xLen-2; i >= 0; i--) {
1021 carry = 0;
1022 for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
1023 long product = (y.value[j+y.offset] & LONG_MASK) *
1024 (value[i+offset] & LONG_MASK) +
1025 (z.value[k] & LONG_MASK) + carry;
1026 z.value[k] = (int)product;
1027 carry = product >>> 32;
1028 }
1029 z.value[i] = (int)carry;
1030 }
1031
1032 // Remove leading zeros from product
1033 z.normalize();
1034 }
1035
1036 /**
1037 * Multiply the contents of this MutableBigInteger by the word y. The
1038 * result is placed into z.
1039 */
1040 void mul(int y, MutableBigInteger z) {
1041 if (y == 1) {
1042 z.copyValue(this);
1043 return;
1044 }
1045
1046 if (y == 0) {
1047 z.clear();
1048 return;
1049 }
1050
1051 // Perform the multiplication word by word
1052 long ylong = y & LONG_MASK;
1053 int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1]
1054 : z.value);
1055 long carry = 0;
1056 for (int i = intLen-1; i >= 0; i--) {
1057 long product = ylong * (value[i+offset] & LONG_MASK) + carry;
1058 zval[i+1] = (int)product;
1059 carry = product >>> 32;
1060 }
1061
1062 if (carry == 0) {
1063 z.offset = 1;
1064 z.intLen = intLen;
1065 } else {
1066 z.offset = 0;
1067 z.intLen = intLen + 1;
1068 zval[0] = (int)carry;
1069 }
1070 z.value = zval;
1071 }
1072
1073 /**
1074 * This method is used for division of an n word dividend by a one word
1075 * divisor. The quotient is placed into quotient. The one word divisor is
1076 * specified by divisor.
1077 *
1078 * @return the remainder of the division is returned.
1079 *
1080 */
1081 int divideOneWord(int divisor, MutableBigInteger quotient) {
1082 long divisorLong = divisor & LONG_MASK;
1083
1084 // Special case of one word dividend
1085 if (intLen == 1) {
1086 long dividendValue = value[offset] & LONG_MASK;
1087 int q = (int) (dividendValue / divisorLong);
1088 int r = (int) (dividendValue - q * divisorLong);
1089 quotient.value[0] = q;
1090 quotient.intLen = (q == 0) ? 0 : 1;
1091 quotient.offset = 0;
1092 return r;
1093 }
1094
1095 if (quotient.value.length < intLen)
1096 quotient.value = new int[intLen];
1097 quotient.offset = 0;
1098 quotient.intLen = intLen;
1099
1100 // Normalize the divisor
1101 int shift = Integer.numberOfLeadingZeros(divisor);
1102
1103 int rem = value[offset];
1104 long remLong = rem & LONG_MASK;
1105 if (remLong < divisorLong) {
1106 quotient.value[0] = 0;
1107 } else {
1108 quotient.value[0] = (int)(remLong / divisorLong);
1109 rem = (int) (remLong - (quotient.value[0] * divisorLong));
1110 remLong = rem & LONG_MASK;
1111 }
1112 int xlen = intLen;
1113 while (--xlen > 0) {
1114 long dividendEstimate = (remLong << 32) |
1115 (value[offset + intLen - xlen] & LONG_MASK);
1116 int q;
1117 if (dividendEstimate >= 0) {
1118 q = (int) (dividendEstimate / divisorLong);
1119 rem = (int) (dividendEstimate - q * divisorLong);
1120 } else {
1121 long tmp = divWord(dividendEstimate, divisor);
1122 q = (int) (tmp & LONG_MASK);
1123 rem = (int) (tmp >>> 32);
1124 }
1125 quotient.value[intLen - xlen] = q;
1126 remLong = rem & LONG_MASK;
1127 }
1128
1129 quotient.normalize();
1130 // Unnormalize
1131 if (shift > 0)
1132 return rem % divisor;
1133 else
1134 return rem;
1135 }
1136
1137 /**
1138 * Calculates the quotient of this div b and places the quotient in the
1139 * provided MutableBigInteger objects and the remainder object is returned.
1140 *
1141 */
1142 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
1143 return divide(b,quotient,true);
1144 }
1145
1146 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1147 if (intLen<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD || b.intLen<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD)
1148 return divideKnuth(b, quotient, needRemainder);
1149 else
1150 return divideAndRemainderBurnikelZiegler(b, quotient);
1151 }
1152
1153 /**
1154 * @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
1155 */
1156 MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) {
1157 return divideKnuth(b,quotient,true);
1158 }
1159
1160 /**
1161 * Calculates the quotient of this div b and places the quotient in the
1162 * provided MutableBigInteger objects and the remainder object is returned.
1163 *
1164 * Uses Algorithm D in Knuth section 4.3.1.
1165 * Many optimizations to that algorithm have been adapted from the Colin
1166 * Plumb C library.
1167 * It special cases one word divisors for speed. The content of b is not
1168 * changed.
1169 *
1170 */
1171 MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1172 if (b.intLen == 0)
1173 throw new ArithmeticException("BigInteger divide by zero");
1174
1175 // Dividend is zero
1176 if (intLen == 0) {
1177 quotient.intLen = quotient.offset = 0;
1178 return needRemainder ? new MutableBigInteger() : null;
1179 }
1180
1181 int cmp = compare(b);
1182 // Dividend less than divisor
1183 if (cmp < 0) {
1184 quotient.intLen = quotient.offset = 0;
1185 return needRemainder ? new MutableBigInteger(this) : null;
1186 }
1187 // Dividend equal to divisor
1188 if (cmp == 0) {
1189 quotient.value[0] = quotient.intLen = 1;
1190 quotient.offset = 0;
1191 return needRemainder ? new MutableBigInteger() : null;
1192 }
1193
1194 quotient.clear();
1195 // Special case one word divisor
1196 if (b.intLen == 1) {
1197 int r = divideOneWord(b.value[b.offset], quotient);
1198 if(needRemainder) {
1199 if (r == 0)
1200 return new MutableBigInteger();
1201 return new MutableBigInteger(r);
1202 } else {
1203 return null;
1204 }
1205 }
1206
1207 // Cancel common powers of two if we're above the KNUTH_POW2_* thresholds
1208 if (intLen >= KNUTH_POW2_THRESH_LEN) {
1209 int trailingZeroBits = Math.min(getLowestSetBit(), b.getLowestSetBit());
1210 if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS*32) {
1211 MutableBigInteger a = new MutableBigInteger(this);
1212 b = new MutableBigInteger(b);
1213 a.rightShift(trailingZeroBits);
1214 b.rightShift(trailingZeroBits);
1215 MutableBigInteger r = a.divideKnuth(b, quotient);
1216 r.leftShift(trailingZeroBits);
1217 return r;
1218 }
1219 }
1220
1221 return divideMagnitude(b, quotient, needRemainder);
1222 }
1223
1224 /**
1225 * Computes {@code this/b} and {@code this%b} using the
1226 * <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
1227 * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
1228 * The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are
1229 * multiples of 32 bits.<br/>
1230 * {@code this} and {@code b} must be nonnegative.
1231 * @param b the divisor
1232 * @param quotient output parameter for {@code this/b}
1233 * @return the remainder
1234 */
1235 MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) {
1236 int r = intLen;
1237 int s = b.intLen;
1238
1239 if (r < s)
1240 return this;
1241 else {
1242 // Unlike Knuth division, we don't check for common powers of two here because
1243 // BZ already runs faster if both numbers contain powers of two and cancelling them has no
1244 // additional benefit.
1245
1246 // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
1247 int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
1248
1249 int j = (s+m-1) / m; // step 2a: j = ceil(s/m)
1250 int n = j * m; // step 2b: block length in 32-bit units
1251 int n32 = 32 * n; // block length in bits
1252 int sigma = Math.max(0, n32 - b.bitLength()); // step 3: sigma = max{T | (2^T)*B < beta^n}
1253 MutableBigInteger bShifted = new MutableBigInteger(b);
1254 bShifted.safeLeftShift(sigma); // step 4a: shift b so its length is a multiple of n
1255 safeLeftShift(sigma); // step 4b: shift this by the same amount
1256
1257 // step 5: t is the number of blocks needed to accommodate this plus one additional bit
1258 int t = (bitLength()+n32) / n32;
1259 if (t < 2)
1260 t = 2;
1261
1262 // step 6: conceptually split this into blocks a[t-1], ..., a[0]
1263 MutableBigInteger a1 = getBlock(t-1, t, n); // the most significant block of this
1264
1265 // step 7: z[t-2] = [a[t-1], a[t-2]]
1266 MutableBigInteger z = getBlock(t-2, t, n); // the second to most significant block
1267 z.addDisjoint(a1, n); // z[t-2]
1268
1269 // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
1270 MutableBigInteger qi = new MutableBigInteger();
1271 MutableBigInteger ri;
1272 quotient.offset = quotient.intLen = 0;
1273 for (int i=t-2; i>0; i--) {
1274 // step 8a: compute (qi,ri) such that z=b*qi+ri
1275 ri = z.divide2n1n(bShifted, qi);
1276
1277 // step 8b: z = [ri, a[i-1]]
1278 z = getBlock(i-1, t, n); // a[i-1]
1279 z.addDisjoint(ri, n);
1280 quotient.addShifted(qi, i*n); // update q (part of step 9)
1281 }
1282 // final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
1283 ri = z.divide2n1n(bShifted, qi);
1284 quotient.add(qi);
1285
1286 ri.rightShift(sigma); // step 9: this and b were shifted, so shift back
1287 return ri;
1288 }
1289 }
1290
1291 /**
1292 * This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
1293 * It divides a 2n-digit number by a n-digit number.<br/>
1294 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
1295 * <br/>
1296 * {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
1297 * @param b a positive number such that {@code b.bitLength()} is even
1298 * @param quotient output parameter for {@code this/b}
1299 * @return {@code this%b}
1300 */
1301 private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
1302 int n = b.intLen;
1303
1304 // step 1: base case
1305 if (n%2!=0 || n<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD)
1306 return divideKnuth(b, quotient);
1307
1308 // step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less
1309 MutableBigInteger aUpper = new MutableBigInteger(this);
1310 aUpper.safeRightShift(32*(n/2)); // aUpper = [a1,a2,a3]
1311 keepLower(n/2); // this = a4
1312
1313 // step 3: q1=aUpper/b, r1=aUpper%b
1314 MutableBigInteger q1 = new MutableBigInteger();
1315 MutableBigInteger r1 = aUpper.divide3n2n(b, q1);
1316
1317 // step 4: quotient=[r1,this]/b, r2=[r1,this]%b
1318 addDisjoint(r1, n/2); // this = [r1,this]
1319 MutableBigInteger r2 = divide3n2n(b, quotient);
1320
1321 // step 5: let quotient=[q1,quotient] and return r2
1322 quotient.addDisjoint(q1, n/2);
1323 return r2;
1324 }
1325
1326 /**
1327 * This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.
1328 * It divides a 3n-digit number by a 2n-digit number.<br/>
1329 * The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.<br/>
1330 * <br/>
1331 * {@code this} must be a nonnegative number such that {@code 2*this.bitLength() <= 3*b.bitLength()}
1332 * @param quotient output parameter for {@code this/b}
1333 * @return {@code this%b}
1334 */
1335 private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) {
1336 int n = b.intLen / 2; // half the length of b in ints
1337
1338 // step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2]
1339 MutableBigInteger a12 = new MutableBigInteger(this);
1340 a12.safeRightShift(32*n);
1341
1342 // step 2: view b as [b1,b2] where each bi is n ints or less
1343 MutableBigInteger b1 = new MutableBigInteger(b);
1344 b1.safeRightShift(n * 32);
1345 BigInteger b2 = b.getLower(n);
1346
1347 MutableBigInteger r;
1348 MutableBigInteger d;
1349 if (compareShifted(b, n) < 0) {
1350 // step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1
1351 r = a12.divide2n1n(b1, quotient);
1352
1353 // step 4: d=quotient*b2
1354 d = new MutableBigInteger(quotient.toBigInteger().multiply(b2));
1355 }
1356 else {
1357 // step 3b: if a1>=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1
1358 quotient.ones(n);
1359 a12.add(b1);
1360 b1.leftShift(32*n);
1361 a12.subtract(b1);
1362 r = a12;
1363
1364 // step 4: d=quotient*b2=(b2 << 32*n) - b2
1365 d = new MutableBigInteger(b2);
1366 d.leftShift(32 * n);
1367 d.subtract(new MutableBigInteger(b2));
1368 }
1369
1370 // step 5: r = r*beta^n + a3 - d (paper says a4)
1371 // However, don't subtract d until after the while loop so r doesn't become negative
1372 r.leftShift(32 * n);
1373 r.addLower(this, n);
1374
1375 // step 6: add b until r>=d
1376 while (r.compare(d) < 0) {
1377 r.add(b);
1378 quotient.subtract(MutableBigInteger.ONE);
1379 }
1380 r.subtract(d);
1381
1382 return r;
1383 }
1384
1385 /**
1386 * Returns a {@code MutableBigInteger} containing {@code blockLength} ints from
1387 * {@code this} number, starting at {@code index*blockLength}.<br/>
1388 * Used by Burnikel-Ziegler division.
1389 * @param index the block index
1390 * @param numBlocks the total number of blocks in {@code this} number
1391 * @param blockLength length of one block in units of 32 bits
1392 * @return
1393 */
1394 private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) {
1395 int blockStart = index * blockLength;
1396 if (blockStart >= intLen)
1397 return new MutableBigInteger();
1398
1399 int blockEnd;
1400 if (index == numBlocks-1)
1401 blockEnd = intLen;
1402 else
1403 blockEnd = (index+1) * blockLength;
1404 if (blockEnd > intLen)
1405 return new MutableBigInteger();
1406
1407 int[] newVal = Arrays.copyOfRange(value, offset+intLen-blockEnd, offset+intLen-blockStart);
1408 return new MutableBigInteger(newVal);
1409 }
1410
1411 /** @see BigInteger#bitLength() */
1412 int bitLength() {
1413 if (intLen == 0)
1414 return 0;
1415 return intLen*32 - Integer.numberOfLeadingZeros(value[offset]);
1416 }
1417
1418 /**
1419 * Internally used to calculate the quotient of this div v and places the
1420 * quotient in the provided MutableBigInteger object and the remainder is
1421 * returned.
1422 *
1423 * @return the remainder of the division will be returned.
1424 */
1425 long divide(long v, MutableBigInteger quotient) {
1426 if (v == 0)
1427 throw new ArithmeticException("BigInteger divide by zero");
1428
1429 // Dividend is zero
1430 if (intLen == 0) {
1431 quotient.intLen = quotient.offset = 0;
1432 return 0;
1433 }
1434 if (v < 0)
1435 v = -v;
1436
1437 int d = (int)(v >>> 32);
1438 quotient.clear();
1439 // Special case on word divisor
1440 if (d == 0)
1441 return divideOneWord((int)v, quotient) & LONG_MASK;
1442 else {
1443 return divideLongMagnitude(v, quotient).toLong();
1444 }
1445 }
1446
1447 private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) {
1448 int n2 = 32 - shift;
1449 int c=src[srcFrom];
1450 for (int i=0; i < srcLen-1; i++) {
1451 int b = c;
1452 c = src[++srcFrom];
1453 dst[dstFrom+i] = (b << shift) | (c >>> n2);
1454 }
1455 dst[dstFrom+srcLen-1] = c << shift;
1456 }
1457
1458 /**
1459 * Divide this MutableBigInteger by the divisor.
1460 * The quotient will be placed into the provided quotient object &
1461 * the remainder object is returned.
1462 */
1463 private MutableBigInteger divideMagnitude(MutableBigInteger div,
1464 MutableBigInteger quotient,
1465 boolean needRemainder ) {
1466 // assert div.intLen > 1
1467 // D1 normalize the divisor
1468 int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
1469 // Copy divisor value to protect divisor
1470 final int dlen = div.intLen;
1471 int[] divisor;
1472 MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
1473 if (shift > 0) {
1474 divisor = new int[dlen];
1475 copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
1476 if(Integer.numberOfLeadingZeros(value[offset])>=shift) {
1477 int[] remarr = new int[intLen + 1];
1478 rem = new MutableBigInteger(remarr);
1479 rem.intLen = intLen;
1480 rem.offset = 1;
1481 copyAndShift(value,offset,intLen,remarr,1,shift);
1482 } else {
1483 int[] remarr = new int[intLen + 2];
1484 rem = new MutableBigInteger(remarr);
1485 rem.intLen = intLen+1;
1486 rem.offset = 1;
1487 int rFrom = offset;
1488 int c=0;
1489 int n2 = 32 - shift;
1490 for (int i=1; i < intLen+1; i++,rFrom++) {
1491 int b = c;
1492 c = value[rFrom];
1493 remarr[i] = (b << shift) | (c >>> n2);
1494 }
1495 remarr[intLen+1] = c << shift;
1496 }
1497 } else {
1498 divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
1499 rem = new MutableBigInteger(new int[intLen + 1]);
1500 System.arraycopy(value, offset, rem.value, 1, intLen);
1501 rem.intLen = intLen;
1502 rem.offset = 1;
1503 }
1504
1505 int nlen = rem.intLen;
1506
1507 // Set the quotient size
1508 final int limit = nlen - dlen + 1;
1509 if (quotient.value.length < limit) {
1510 quotient.value = new int[limit];
1511 quotient.offset = 0;
1512 }
1513 quotient.intLen = limit;
1514 int[] q = quotient.value;
1515
1516
1517 // Must insert leading 0 in rem if its length did not change
1518 if (rem.intLen == nlen) {
1519 rem.offset = 0;
1520 rem.value[0] = 0;
1521 rem.intLen++;
1522 }
1523
1524 int dh = divisor[0];
1525 long dhLong = dh & LONG_MASK;
1526 int dl = divisor[1];
1527
1528 // D2 Initialize j
1529 for(int j=0; j<limit-1; j++) {
1530 // D3 Calculate qhat
1531 // estimate qhat
1532 int qhat = 0;
1533 int qrem = 0;
1534 boolean skipCorrection = false;
1535 int nh = rem.value[j+rem.offset];
1536 int nh2 = nh + 0x80000000;
1537 int nm = rem.value[j+1+rem.offset];
1538
1539 if (nh == dh) {
1540 qhat = ~0;
1541 qrem = nh + nm;
1542 skipCorrection = qrem + 0x80000000 < nh2;
1543 } else {
1544 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1545 if (nChunk >= 0) {
1546 qhat = (int) (nChunk / dhLong);
1547 qrem = (int) (nChunk - (qhat * dhLong));
1548 } else {
1549 long tmp = divWord(nChunk, dh);
1550 qhat = (int) (tmp & LONG_MASK);
1551 qrem = (int) (tmp >>> 32);
1552 }
1553 }
1554
1555 if (qhat == 0)
1556 continue;
1557
1558 if (!skipCorrection) { // Correct qhat
1559 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1560 long rs = ((qrem & LONG_MASK) << 32) | nl;
1561 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1562
1563 if (unsignedLongCompare(estProduct, rs)) {
1564 qhat--;
1565 qrem = (int)((qrem & LONG_MASK) + dhLong);
1566 if ((qrem & LONG_MASK) >= dhLong) {
1567 estProduct -= (dl & LONG_MASK);
1568 rs = ((qrem & LONG_MASK) << 32) | nl;
1569 if (unsignedLongCompare(estProduct, rs))
1570 qhat--;
1571 }
1572 }
1573 }
1574
1575 // D4 Multiply and subtract
1576 rem.value[j+rem.offset] = 0;
1577 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1578
1579 // D5 Test remainder
1580 if (borrow + 0x80000000 > nh2) {
1581 // D6 Add back
1582 divadd(divisor, rem.value, j+1+rem.offset);
1583 qhat--;
1584 }
1585
1586 // Store the quotient digit
1587 q[j] = qhat;
1588 } // D7 loop on j
1589 // D3 Calculate qhat
1590 // estimate qhat
1591 int qhat = 0;
1592 int qrem = 0;
1593 boolean skipCorrection = false;
1594 int nh = rem.value[limit - 1 + rem.offset];
1595 int nh2 = nh + 0x80000000;
1596 int nm = rem.value[limit + rem.offset];
1597
1598 if (nh == dh) {
1599 qhat = ~0;
1600 qrem = nh + nm;
1601 skipCorrection = qrem + 0x80000000 < nh2;
1602 } else {
1603 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1604 if (nChunk >= 0) {
1605 qhat = (int) (nChunk / dhLong);
1606 qrem = (int) (nChunk - (qhat * dhLong));
1607 } else {
1608 long tmp = divWord(nChunk, dh);
1609 qhat = (int) (tmp & LONG_MASK);
1610 qrem = (int) (tmp >>> 32);
1611 }
1612 }
1613 if (qhat != 0) {
1614 if (!skipCorrection) { // Correct qhat
1615 long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
1616 long rs = ((qrem & LONG_MASK) << 32) | nl;
1617 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1618
1619 if (unsignedLongCompare(estProduct, rs)) {
1620 qhat--;
1621 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1622 if ((qrem & LONG_MASK) >= dhLong) {
1623 estProduct -= (dl & LONG_MASK);
1624 rs = ((qrem & LONG_MASK) << 32) | nl;
1625 if (unsignedLongCompare(estProduct, rs))
1626 qhat--;
1627 }
1628 }
1629 }
1630
1631
1632 // D4 Multiply and subtract
1633 int borrow;
1634 rem.value[limit - 1 + rem.offset] = 0;
1635 if(needRemainder)
1636 borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1637 else
1638 borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1639
1640 // D5 Test remainder
1641 if (borrow + 0x80000000 > nh2) {
1642 // D6 Add back
1643 if(needRemainder)
1644 divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
1645 qhat--;
1646 }
1647
1648 // Store the quotient digit
1649 q[(limit - 1)] = qhat;
1650 }
1651
1652
1653 if(needRemainder) {
1654 // D8 Unnormalize
1655 if (shift > 0)
1656 rem.rightShift(shift);
1657 rem.normalize();
1658 }
1659 quotient.normalize();
1660 return needRemainder ? rem : null;
1661 }
1662
1663 /**
1664 * Divide this MutableBigInteger by the divisor represented by positive long
1665 * value. The quotient will be placed into the provided quotient object &
1666 * the remainder object is returned.
1667 */
1668 private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
1669 // Remainder starts as dividend with space for a leading zero
1670 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
1671 System.arraycopy(value, offset, rem.value, 1, intLen);
1672 rem.intLen = intLen;
1673 rem.offset = 1;
1674
1675 int nlen = rem.intLen;
1676
1677 int limit = nlen - 2 + 1;
1678 if (quotient.value.length < limit) {
1679 quotient.value = new int[limit];
1680 quotient.offset = 0;
1681 }
1682 quotient.intLen = limit;
1683 int[] q = quotient.value;
1684
1685 // D1 normalize the divisor
1686 int shift = Long.numberOfLeadingZeros(ldivisor);
1687 if (shift > 0) {
1688 ldivisor<<=shift;
1689 rem.leftShift(shift);
1690 }
1691
1692 // Must insert leading 0 in rem if its length did not change
1693 if (rem.intLen == nlen) {
1694 rem.offset = 0;
1695 rem.value[0] = 0;
1696 rem.intLen++;
1697 }
1698
1699 int dh = (int)(ldivisor >>> 32);
1700 long dhLong = dh & LONG_MASK;
1701 int dl = (int)(ldivisor & LONG_MASK);
1702
1703 // D2 Initialize j
1704 for (int j = 0; j < limit; j++) {
1705 // D3 Calculate qhat
1706 // estimate qhat
1707 int qhat = 0;
1708 int qrem = 0;
1709 boolean skipCorrection = false;
1710 int nh = rem.value[j + rem.offset];
1711 int nh2 = nh + 0x80000000;
1712 int nm = rem.value[j + 1 + rem.offset];
1713
1714 if (nh == dh) {
1715 qhat = ~0;
1716 qrem = nh + nm;
1717 skipCorrection = qrem + 0x80000000 < nh2;
1718 } else {
1719 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1720 if (nChunk >= 0) {
1721 qhat = (int) (nChunk / dhLong);
1722 qrem = (int) (nChunk - (qhat * dhLong));
1723 } else {
1724 long tmp = divWord(nChunk, dh);
1725 qhat =(int)(tmp & LONG_MASK);
1726 qrem = (int)(tmp>>>32);
1727 }
1728 }
1729
1730 if (qhat == 0)
1731 continue;
1732
1733 if (!skipCorrection) { // Correct qhat
1734 long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
1735 long rs = ((qrem & LONG_MASK) << 32) | nl;
1736 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1737
1738 if (unsignedLongCompare(estProduct, rs)) {
1739 qhat--;
1740 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1741 if ((qrem & LONG_MASK) >= dhLong) {
1742 estProduct -= (dl & LONG_MASK);
1743 rs = ((qrem & LONG_MASK) << 32) | nl;
1744 if (unsignedLongCompare(estProduct, rs))
1745 qhat--;
1746 }
1747 }
1748 }
1749
1750 // D4 Multiply and subtract
1751 rem.value[j + rem.offset] = 0;
1752 int borrow = mulsubLong(rem.value, dh, dl, qhat, j + rem.offset);
1753
1754 // D5 Test remainder
1755 if (borrow + 0x80000000 > nh2) {
1756 // D6 Add back
1757 divaddLong(dh,dl, rem.value, j + 1 + rem.offset);
1758 qhat--;
1759 }
1760
1761 // Store the quotient digit
1762 q[j] = qhat;
1763 } // D7 loop on j
1764
1765 // D8 Unnormalize
1766 if (shift > 0)
1767 rem.rightShift(shift);
1768
1769 quotient.normalize();
1770 rem.normalize();
1771 return rem;
1772 }
1773
1774 /**
1775 * A primitive used for division by long.
1776 * Specialized version of the method divadd.
1777 * dh is a high part of the divisor, dl is a low part
1778 */
1779 private int divaddLong(int dh, int dl, int[] result, int offset) {
1780 long carry = 0;
1781
1782 long sum = (dl & LONG_MASK) + (result[1+offset] & LONG_MASK);
1783 result[1+offset] = (int)sum;
1784
1785 sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry;
1786 result[offset] = (int)sum;
1787 carry = sum >>> 32;
1788 return (int)carry;
1789 }
1790
1791 /**
1792 * This method is used for division by long.
1793 * Specialized version of the method sulsub.
1794 * dh is a high part of the divisor, dl is a low part
1795 */
1796 private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
1797 long xLong = x & LONG_MASK;
1798 offset += 2;
1799 long product = (dl & LONG_MASK) * xLong;
1800 long difference = q[offset] - product;
1801 q[offset--] = (int)difference;
1802 long carry = (product >>> 32)
1803 + (((difference & LONG_MASK) >
1804 (((~(int)product) & LONG_MASK))) ? 1:0);
1805 product = (dh & LONG_MASK) * xLong + carry;
1806 difference = q[offset] - product;
1807 q[offset--] = (int)difference;
1808 carry = (product >>> 32)
1809 + (((difference & LONG_MASK) >
1810 (((~(int)product) & LONG_MASK))) ? 1:0);
1811 return (int)carry;
1812 }
1813
1814 /**
1815 * Compare two longs as if they were unsigned.
1816 * Returns true iff one is bigger than two.
1817 */
1818 private boolean unsignedLongCompare(long one, long two) {
1819 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1820 }
1821
1822 /**
1823 * This method divides a long quantity by an int to estimate
1824 * qhat for two multi precision numbers. It is used when
1825 * the signed value of n is less than zero.
1826 * Returns long value where high 32 bits contain remainder value and
1827 * low 32 bits contain quotient value.
1828 */
1829 static long divWord(long n, int d) {
1830 long dLong = d & LONG_MASK;
1831 long r;
1832 long q;
1833 if (dLong == 1) {
1834 q = (int)n;
1835 r = 0;
1836 return (r << 32) | (q & LONG_MASK);
1837 }
1838
1839 // Approximate the quotient and remainder
1840 q = (n >>> 1) / (dLong >>> 1);
1841 r = n - q*dLong;
1842
1843 // Correct the approximation
1844 while (r < 0) {
1845 r += dLong;
1846 q--;
1847 }
1848 while (r >= dLong) {
1849 r -= dLong;
1850 q++;
1851 }
1852 // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1853 return (r << 32) | (q & LONG_MASK);
1854 }
1855
1856 /**
1857 * Calculate GCD of this and b. This and b are changed by the computation.
1858 */
1859 MutableBigInteger hybridGCD(MutableBigInteger b) {
1860 // Use Euclid's algorithm until the numbers are approximately the
1861 // same length, then use the binary GCD algorithm to find the GCD.
1862 MutableBigInteger a = this;
1863 MutableBigInteger q = new MutableBigInteger();
1864
1865 while (b.intLen != 0) {
1866 if (Math.abs(a.intLen - b.intLen) < 2)
1867 return a.binaryGCD(b);
1868
1869 MutableBigInteger r = a.divide(b, q);
1870 a = b;
1871 b = r;
1872 }
1873 return a;
1874 }
1875
1876 /**
1877 * Calculate GCD of this and v.
1878 * Assumes that this and v are not zero.
1879 */
1880 private MutableBigInteger binaryGCD(MutableBigInteger v) {
1881 // Algorithm B from Knuth section 4.5.2
1882 MutableBigInteger u = this;
1883 MutableBigInteger r = new MutableBigInteger();
1884
1885 // step B1
1886 int s1 = u.getLowestSetBit();
1887 int s2 = v.getLowestSetBit();
1888 int k = (s1 < s2) ? s1 : s2;
1889 if (k != 0) {
1890 u.rightShift(k);
1891 v.rightShift(k);
1892 }
1893
1894 // step B2
1895 boolean uOdd = (k==s1);
1896 MutableBigInteger t = uOdd ? v: u;
1897 int tsign = uOdd ? -1 : 1;
1898
1899 int lb;
1900 while ((lb = t.getLowestSetBit()) >= 0) {
1901 // steps B3 and B4
1902 t.rightShift(lb);
1903 // step B5
1904 if (tsign > 0)
1905 u = t;
1906 else
1907 v = t;
1908
1909 // Special case one word numbers
1910 if (u.intLen < 2 && v.intLen < 2) {
1911 int x = u.value[u.offset];
1912 int y = v.value[v.offset];
1913 x = binaryGcd(x, y);
1914 r.value[0] = x;
1915 r.intLen = 1;
1916 r.offset = 0;
1917 if (k > 0)
1918 r.leftShift(k);
1919 return r;
1920 }
1921
1922 // step B6
1923 if ((tsign = u.difference(v)) == 0)
1924 break;
1925 t = (tsign >= 0) ? u : v;
1926 }
1927
1928 if (k > 0)
1929 u.leftShift(k);
1930 return u;
1931 }
1932
1933 /**
1934 * Calculate GCD of a and b interpreted as unsigned integers.
1935 */
1936 static int binaryGcd(int a, int b) {
1937 if (b==0)
1938 return a;
1939 if (a==0)
1940 return b;
1941
1942 // Right shift a & b till their last bits equal to 1.
1943 int aZeros = Integer.numberOfTrailingZeros(a);
1944 int bZeros = Integer.numberOfTrailingZeros(b);
1945 a >>>= aZeros;
1946 b >>>= bZeros;
1947
1948 int t = (aZeros < bZeros ? aZeros : bZeros);
1949
1950 while (a != b) {
1951 if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned
1952 a -= b;
1953 a >>>= Integer.numberOfTrailingZeros(a);
1954 } else {
1955 b -= a;
1956 b >>>= Integer.numberOfTrailingZeros(b);
1957 }
1958 }
1959 return a<<t;
1960 }
1961
1962 /**
1963 * Returns the modInverse of this mod p. This and p are not affected by
1964 * the operation.
1965 */
1966 MutableBigInteger mutableModInverse(MutableBigInteger p) {
1967 // Modulus is odd, use Schroeppel's algorithm
1968 if (p.isOdd())
1969 return modInverse(p);
1970
1971 // Base and modulus are even, throw exception
1972 if (isEven())
1973 throw new ArithmeticException("BigInteger not invertible.");
1974
1975 // Get even part of modulus expressed as a power of 2
1976 int powersOf2 = p.getLowestSetBit();
1977
1978 // Construct odd part of modulus
1979 MutableBigInteger oddMod = new MutableBigInteger(p);
1980 oddMod.rightShift(powersOf2);
1981
1982 if (oddMod.isOne())
1983 return modInverseMP2(powersOf2);
1984
1985 // Calculate 1/a mod oddMod
1986 MutableBigInteger oddPart = modInverse(oddMod);
1987
1988 // Calculate 1/a mod evenMod
1989 MutableBigInteger evenPart = modInverseMP2(powersOf2);
1990
1991 // Combine the results using Chinese Remainder Theorem
1992 MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
1993 MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
1994
1995 MutableBigInteger temp1 = new MutableBigInteger();
1996 MutableBigInteger temp2 = new MutableBigInteger();
1997 MutableBigInteger result = new MutableBigInteger();
1998
1999 oddPart.leftShift(powersOf2);
2000 oddPart.multiply(y1, result);
2001
2002 evenPart.multiply(oddMod, temp1);
2003 temp1.multiply(y2, temp2);
2004
2005 result.add(temp2);
2006 return result.divide(p, temp1);
2007 }
2008
2009 /*
2010 * Calculate the multiplicative inverse of this mod 2^k.
2011 */
2012 MutableBigInteger modInverseMP2(int k) {
2013 if (isEven())
2014 throw new ArithmeticException("Non-invertible. (GCD != 1)");
2015
2016 if (k > 64)
2017 return euclidModInverse(k);
2018
2019 int t = inverseMod32(value[offset+intLen-1]);
2020
2021 if (k < 33) {
2022 t = (k == 32 ? t : t & ((1 << k) - 1));
2023 return new MutableBigInteger(t);
2024 }
2025
2026 long pLong = (value[offset+intLen-1] & LONG_MASK);
2027 if (intLen > 1)
2028 pLong |= ((long)value[offset+intLen-2] << 32);
2029 long tLong = t & LONG_MASK;
2030 tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
2031 tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
2032
2033 MutableBigInteger result = new MutableBigInteger(new int[2]);
2034 result.value[0] = (int)(tLong >>> 32);
2035 result.value[1] = (int)tLong;
2036 result.intLen = 2;
2037 result.normalize();
2038 return result;
2039 }
2040
2041 /**
2042 * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
2043 */
2044 static int inverseMod32(int val) {
2045 // Newton's iteration!
2046 int t = val;
2047 t *= 2 - val*t;
2048 t *= 2 - val*t;
2049 t *= 2 - val*t;
2050 t *= 2 - val*t;
2051 return t;
2052 }
2053
2054 /**
2055 * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
2056 */
2057 static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
2058 // Copy the mod to protect original
2059 return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
2060 }
2061
2062 /**
2063 * Calculate the multiplicative inverse of this mod mod, where mod is odd.
2064 * This and mod are not changed by the calculation.
2065 *
2066 * This method implements an algorithm due to Richard Schroeppel, that uses
2067 * the same intermediate representation as Montgomery Reduction
2068 * ("Montgomery Form"). The algorithm is described in an unpublished
2069 * manuscript entitled "Fast Modular Reciprocals."
2070 */
2071 private MutableBigInteger modInverse(MutableBigInteger mod) {
2072 MutableBigInteger p = new MutableBigInteger(mod);
2073 MutableBigInteger f = new MutableBigInteger(this);
2074 MutableBigInteger g = new MutableBigInteger(p);
2075 SignedMutableBigInteger c = new SignedMutableBigInteger(1);
2076 SignedMutableBigInteger d = new SignedMutableBigInteger();
2077 MutableBigInteger temp = null;
2078 SignedMutableBigInteger sTemp = null;
2079
2080 int k = 0;
2081 // Right shift f k times until odd, left shift d k times
2082 if (f.isEven()) {
2083 int trailingZeros = f.getLowestSetBit();
2084 f.rightShift(trailingZeros);
2085 d.leftShift(trailingZeros);
2086 k = trailingZeros;
2087 }
2088
2089 // The Almost Inverse Algorithm
2090 while(!f.isOne()) {
2091 // If gcd(f, g) != 1, number is not invertible modulo mod
2092 if (f.isZero())
2093 throw new ArithmeticException("BigInteger not invertible.");
2094
2095 // If f < g exchange f, g and c, d
2096 if (f.compare(g) < 0) {
2097 temp = f; f = g; g = temp;
2098 sTemp = d; d = c; c = sTemp;
2099 }
2100
2101 // If f == g (mod 4)
2102 if (((f.value[f.offset + f.intLen - 1] ^
2103 g.value[g.offset + g.intLen - 1]) & 3) == 0) {
2104 f.subtract(g);
2105 c.signedSubtract(d);
2106 } else { // If f != g (mod 4)
2107 f.add(g);
2108 c.signedAdd(d);
2109 }
2110
2111 // Right shift f k times until odd, left shift d k times
2112 int trailingZeros = f.getLowestSetBit();
2113 f.rightShift(trailingZeros);
2114 d.leftShift(trailingZeros);
2115 k += trailingZeros;
2116 }
2117
2118 while (c.sign < 0)
2119 c.signedAdd(p);
2120
2121 return fixup(c, p, k);
2122 }
2123
2124 /**
2125 * The Fixup Algorithm
2126 * Calculates X such that X = C * 2^(-k) (mod P)
2127 * Assumes C<P and P is odd.
2128 */
2129 static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
2130 int k) {
2131 MutableBigInteger temp = new MutableBigInteger();
2132 // Set r to the multiplicative inverse of p mod 2^32
2133 int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
2134
2135 for(int i=0, numWords = k >> 5; i<numWords; i++) {
2136 // V = R * c (mod 2^j)
2137 int v = r * c.value[c.offset + c.intLen-1];
2138 // c = c + (v * p)
2139 p.mul(v, temp);
2140 c.add(temp);
2141 // c = c / 2^j
2142 c.intLen--;
2143 }
2144 int numBits = k & 0x1f;
2145 if (numBits != 0) {
2146 // V = R * c (mod 2^j)
2147 int v = r * c.value[c.offset + c.intLen-1];
2148 v &= ((1<<numBits) - 1);
2149 // c = c + (v * p)
2150 p.mul(v, temp);
2151 c.add(temp);
2152 // c = c / 2^j
2153 c.rightShift(numBits);
2154 }
2155
2156 // In theory, c may be greater than p at this point (Very rare!)
2157 while (c.compare(p) >= 0)
2158 c.subtract(p);
2159
2160 return c;
2161 }
2162
2163 /**
2164 * Uses the extended Euclidean algorithm to compute the modInverse of base
2165 * mod a modulus that is a power of 2. The modulus is 2^k.
2166 */
2167 MutableBigInteger euclidModInverse(int k) {
2168 MutableBigInteger b = new MutableBigInteger(1);
2169 b.leftShift(k);
2170 MutableBigInteger mod = new MutableBigInteger(b);
2171
2172 MutableBigInteger a = new MutableBigInteger(this);
2173 MutableBigInteger q = new MutableBigInteger();
2174 MutableBigInteger r = b.divide(a, q);
2175
2176 MutableBigInteger swapper = b;
2177 // swap b & r
2178 b = r;
2179 r = swapper;
2180
2181 MutableBigInteger t1 = new MutableBigInteger(q);
2182 MutableBigInteger t0 = new MutableBigInteger(1);
2183 MutableBigInteger temp = new MutableBigInteger();
2184
2185 while (!b.isOne()) {
2186 r = a.divide(b, q);
2187
2188 if (r.intLen == 0)
2189 throw new ArithmeticException("BigInteger not invertible.");
2190
2191 swapper = r;
2192 a = swapper;
2193
2194 if (q.intLen == 1)
2195 t1.mul(q.value[q.offset], temp);
2196 else
2197 q.multiply(t1, temp);
2198 swapper = q;
2199 q = temp;
2200 temp = swapper;
2201 t0.add(q);
2202
2203 if (a.isOne())
2204 return t0;
2205
2206 r = b.divide(a, q);
2207
2208 if (r.intLen == 0)
2209 throw new ArithmeticException("BigInteger not invertible.");
2210
2211 swapper = b;
2212 b = r;
2213
2214 if (q.intLen == 1)
2215 t0.mul(q.value[q.offset], temp);
2216 else
2217 q.multiply(t0, temp);
2218 swapper = q; q = temp; temp = swapper;
2219
2220 t1.add(q);
2221 }
2222 mod.subtract(t1);
2223 return mod;
2224 }
2225 }