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