1 /*
2 * Copyright (c) 1999, 2011, 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 * @since 1.3
42 */
43
44 import java.util.Arrays;
45
46 import static java.math.BigInteger.LONG_MASK;
47 import static java.math.BigDecimal.INFLATED;
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 // Constructors
79
80 /**
81 * The default constructor. An empty MutableBigInteger is created with
82 * a one word capacity.
83 */
84 MutableBigInteger() {
85 value = new int[1];
86 intLen = 0;
87 }
88
89 /**
90 * Construct a new MutableBigInteger with a magnitude specified by
91 * the int val.
92 */
93 MutableBigInteger(int val) {
94 value = new int[1];
95 intLen = 1;
96 value[0] = val;
97 }
98
99 /**
100 * Construct a new MutableBigInteger with the specified value array
101 * up to the length of the array supplied.
102 */
103 MutableBigInteger(int[] val) {
104 value = val;
105 intLen = val.length;
106 }
107
108 /**
109 * Construct a new MutableBigInteger with a magnitude equal to the
110 * specified BigInteger.
111 */
112 MutableBigInteger(BigInteger b) {
113 intLen = b.mag.length;
114 value = Arrays.copyOf(b.mag, intLen);
115 }
116
117 /**
118 * Construct a new MutableBigInteger with a magnitude equal to the
119 * specified MutableBigInteger.
120 */
121 MutableBigInteger(MutableBigInteger val) {
122 intLen = val.intLen;
123 value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
124 }
125
126 /**
127 * Internal helper method to return the magnitude array. The caller is not
128 * supposed to modify the returned array.
129 */
130 private int[] getMagnitudeArray() {
131 if (offset > 0 || value.length != intLen)
132 return Arrays.copyOfRange(value, offset, offset + intLen);
133 return value;
134 }
135
136 /**
137 * Convert this MutableBigInteger to a long value. The caller has to make
138 * sure this MutableBigInteger can be fit into long.
139 */
140 private long toLong() {
141 assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
142 if (intLen == 0)
143 return 0;
144 long d = value[offset] & LONG_MASK;
145 return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
146 }
147
148 /**
149 * Convert this MutableBigInteger to a BigInteger object.
150 */
151 BigInteger toBigInteger(int sign) {
152 if (intLen == 0 || sign == 0)
153 return BigInteger.ZERO;
154 return new BigInteger(getMagnitudeArray(), sign);
155 }
156
157 /**
158 * Convert this MutableBigInteger to BigDecimal object with the specified sign
159 * and scale.
160 */
161 BigDecimal toBigDecimal(int sign, int scale) {
162 if (intLen == 0 || sign == 0)
163 return BigDecimal.zeroValueOf(scale);
164 int[] mag = getMagnitudeArray();
165 int len = mag.length;
166 int d = mag[0];
167 // If this MutableBigInteger can't be fit into long, we need to
168 // make a BigInteger object for the resultant BigDecimal object.
169 if (len > 2 || (d < 0 && len == 2))
170 return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
171 long v = (len == 2) ?
172 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
173 d & LONG_MASK;
174 return BigDecimal.valueOf(sign == -1 ? -v : v, scale);
175 }
176
177 /**
178 * This is for internal use in converting from a MutableBigInteger
179 * object into a long value given a specified sign.
180 * returns INFLATED if value is not fit into long
181 */
182 long toCompactValue(int sign) {
183 if (intLen == 0 || sign == 0)
184 return 0L;
185 int[] mag = getMagnitudeArray();
186 int len = mag.length;
187 int d = mag[0];
188 // If this MutableBigInteger can not be fitted into long, we need to
189 // make a BigInteger object for the resultant BigDecimal object.
190 if (len > 2 || (d < 0 && len == 2))
191 return INFLATED;
192 long v = (len == 2) ?
193 ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
194 d & LONG_MASK;
195 return sign == -1 ? -v : v;
196 }
197
198 /**
199 * Clear out a MutableBigInteger for reuse.
200 */
201 void clear() {
202 offset = intLen = 0;
203 for (int index=0, n=value.length; index < n; index++)
204 value[index] = 0;
205 }
206
207 /**
208 * Set a MutableBigInteger to zero, removing its offset.
209 */
210 void reset() {
211 offset = intLen = 0;
212 }
213
214 /**
215 * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
216 * as this MutableBigInteger is numerically less than, equal to, or
217 * greater than <tt>b</tt>.
218 */
219 final int compare(MutableBigInteger b) {
220 int blen = b.intLen;
221 if (intLen < blen)
222 return -1;
223 if (intLen > blen)
224 return 1;
225
226 // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
227 // comparison.
228 int[] bval = b.value;
229 for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
230 int b1 = value[i] + 0x80000000;
231 int b2 = bval[j] + 0x80000000;
232 if (b1 < b2)
233 return -1;
234 if (b1 > b2)
235 return 1;
236 }
237 return 0;
238 }
239
240 /**
241 * Compare this against half of a MutableBigInteger object (Needed for
242 * remainder tests).
243 * Assumes no leading unnecessary zeros, which holds for results
244 * from divide().
245 */
246 final int compareHalf(MutableBigInteger b) {
247 int blen = b.intLen;
248 int len = intLen;
249 if (len <= 0)
250 return blen <=0 ? 0 : -1;
251 if (len > blen)
252 return 1;
253 if (len < blen - 1)
254 return -1;
255 int[] bval = b.value;
256 int bstart = 0;
257 int carry = 0;
258 // Only 2 cases left:len == blen or len == blen - 1
259 if (len != blen) { // len == blen - 1
260 if (bval[bstart] == 1) {
261 ++bstart;
262 carry = 0x80000000;
263 } else
264 return -1;
265 }
266 // compare values with right-shifted values of b,
267 // carrying shifted-out bits across words
268 int[] val = value;
269 for (int i = offset, j = bstart; i < len + offset;) {
270 int bv = bval[j++];
271 long hb = ((bv >>> 1) + carry) & LONG_MASK;
272 long v = val[i++] & LONG_MASK;
273 if (v != hb)
274 return v < hb ? -1 : 1;
275 carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
276 }
277 return carry == 0? 0 : -1;
278 }
279
280 /**
281 * Return the index of the lowest set bit in this MutableBigInteger. If the
282 * magnitude of this MutableBigInteger is zero, -1 is returned.
283 */
284 private final int getLowestSetBit() {
285 if (intLen == 0)
286 return -1;
287 int j, b;
288 for (j=intLen-1; (j>0) && (value[j+offset]==0); j--)
289 ;
290 b = value[j+offset];
291 if (b==0)
292 return -1;
293 return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
294 }
295
296 /**
297 * Return the int in use in this MutableBigInteger at the specified
298 * index. This method is not used because it is not inlined on all
299 * platforms.
300 */
301 private final int getInt(int index) {
302 return value[offset+index];
303 }
304
305 /**
306 * Return a long which is equal to the unsigned value of the int in
307 * use in this MutableBigInteger at the specified index. This method is
308 * not used because it is not inlined on all platforms.
309 */
310 private final long getLong(int index) {
311 return value[offset+index] & LONG_MASK;
312 }
313
314 /**
315 * Ensure that the MutableBigInteger is in normal form, specifically
316 * making sure that there are no leading zeros, and that if the
317 * magnitude is zero, then intLen is zero.
318 */
319 final void normalize() {
320 if (intLen == 0) {
321 offset = 0;
322 return;
323 }
324
325 int index = offset;
326 if (value[index] != 0)
327 return;
328
329 int indexBound = index+intLen;
330 do {
331 index++;
332 } while(index < indexBound && value[index]==0);
333
334 int numZeros = index - offset;
335 intLen -= numZeros;
336 offset = (intLen==0 ? 0 : offset+numZeros);
337 }
338
339 /**
340 * If this MutableBigInteger cannot hold len words, increase the size
341 * of the value array to len words.
342 */
343 private final void ensureCapacity(int len) {
344 if (value.length < len) {
345 value = new int[len];
346 offset = 0;
347 intLen = len;
348 }
349 }
350
351 /**
352 * Convert this MutableBigInteger into an int array with no leading
353 * zeros, of a length that is equal to this MutableBigInteger's intLen.
354 */
355 int[] toIntArray() {
356 int[] result = new int[intLen];
357 for(int i=0; i<intLen; i++)
358 result[i] = value[offset+i];
359 return result;
360 }
361
362 /**
363 * Sets the int at index+offset in this MutableBigInteger to val.
364 * This does not get inlined on all platforms so it is not used
365 * as often as originally intended.
366 */
367 void setInt(int index, int val) {
368 value[offset + index] = val;
369 }
370
371 /**
372 * Sets this MutableBigInteger's value array to the specified array.
373 * The intLen is set to the specified length.
374 */
375 void setValue(int[] val, int length) {
376 value = val;
377 intLen = length;
378 offset = 0;
379 }
380
381 /**
382 * Sets this MutableBigInteger's value array to a copy of the specified
383 * array. The intLen is set to the length of the new array.
384 */
385 void copyValue(MutableBigInteger src) {
386 int len = src.intLen;
387 if (value.length < len)
388 value = new int[len];
389 System.arraycopy(src.value, src.offset, value, 0, len);
390 intLen = len;
391 offset = 0;
392 }
393
394 /**
395 * Sets this MutableBigInteger's value array to a copy of the specified
396 * array. The intLen is set to the length of the specified array.
397 */
398 void copyValue(int[] val) {
399 int len = val.length;
400 if (value.length < len)
401 value = new int[len];
402 System.arraycopy(val, 0, value, 0, len);
403 intLen = len;
404 offset = 0;
405 }
406
407 /**
408 * Returns true iff this MutableBigInteger has a value of one.
409 */
410 boolean isOne() {
411 return (intLen == 1) && (value[offset] == 1);
412 }
413
414 /**
415 * Returns true iff this MutableBigInteger has a value of zero.
416 */
417 boolean isZero() {
418 return (intLen == 0);
419 }
420
421 /**
422 * Returns true iff this MutableBigInteger is even.
423 */
424 boolean isEven() {
425 return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
426 }
427
428 /**
429 * Returns true iff this MutableBigInteger is odd.
430 */
431 boolean isOdd() {
432 return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
433 }
434
435 /**
436 * Returns true iff this MutableBigInteger is in normal form. A
437 * MutableBigInteger is in normal form if it has no leading zeros
438 * after the offset, and intLen + offset <= value.length.
439 */
440 boolean isNormal() {
441 if (intLen + offset > value.length)
442 return false;
443 if (intLen ==0)
444 return true;
445 return (value[offset] != 0);
446 }
447
448 /**
449 * Returns a String representation of this MutableBigInteger in radix 10.
450 */
451 public String toString() {
452 BigInteger b = toBigInteger(1);
453 return b.toString();
454 }
455
456 /**
457 * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
458 * in normal form.
459 */
460 void rightShift(int n) {
461 if (intLen == 0)
462 return;
463 int nInts = n >>> 5;
464 int nBits = n & 0x1F;
465 this.intLen -= nInts;
466 if (nBits == 0)
467 return;
468 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
469 if (nBits >= bitsInHighWord) {
470 this.primitiveLeftShift(32 - nBits);
471 this.intLen--;
472 } else {
473 primitiveRightShift(nBits);
474 }
475 }
476
477 /**
478 * Left shift this MutableBigInteger n bits.
479 */
480 void leftShift(int n) {
481 /*
482 * If there is enough storage space in this MutableBigInteger already
483 * the available space will be used. Space to the right of the used
484 * ints in the value array is faster to utilize, so the extra space
485 * will be taken from the right if possible.
486 */
487 if (intLen == 0)
488 return;
489 int nInts = n >>> 5;
490 int nBits = n&0x1F;
491 int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
492
493 // If shift can be done without moving words, do so
494 if (n <= (32-bitsInHighWord)) {
495 primitiveLeftShift(nBits);
496 return;
497 }
498
499 int newLen = intLen + nInts +1;
500 if (nBits <= (32-bitsInHighWord))
501 newLen--;
502 if (value.length < newLen) {
503 // The array must grow
504 int[] result = new int[newLen];
505 for (int i=0; i<intLen; i++)
506 result[i] = value[offset+i];
507 setValue(result, newLen);
508 } else if (value.length - offset >= newLen) {
509 // Use space on right
510 for(int i=0; i<newLen - intLen; i++)
511 value[offset+intLen+i] = 0;
512 } else {
513 // Must use space on left
514 for (int i=0; i<intLen; i++)
515 value[i] = value[offset+i];
516 for (int i=intLen; i<newLen; i++)
517 value[i] = 0;
518 offset = 0;
519 }
520 intLen = newLen;
521 if (nBits == 0)
522 return;
523 if (nBits <= (32-bitsInHighWord))
524 primitiveLeftShift(nBits);
525 else
526 primitiveRightShift(32 -nBits);
527 }
528
529 /**
530 * A primitive used for division. This method adds in one multiple of the
531 * divisor a back to the dividend result at a specified offset. It is used
532 * when qhat was estimated too large, and must be adjusted.
533 */
534 private int divadd(int[] a, int[] result, int offset) {
535 long carry = 0;
536
537 for (int j=a.length-1; j >= 0; j--) {
538 long sum = (a[j] & LONG_MASK) +
539 (result[j+offset] & LONG_MASK) + carry;
540 result[j+offset] = (int)sum;
541 carry = sum >>> 32;
542 }
543 return (int)carry;
544 }
545
546 /**
547 * This method is used for division. It multiplies an n word input a by one
548 * word input x, and subtracts the n word product from q. This is needed
549 * when subtracting qhat*divisor from dividend.
550 */
551 private int mulsub(int[] q, int[] a, int x, int len, int offset) {
552 long xLong = x & LONG_MASK;
553 long carry = 0;
554 offset += len;
555
556 for (int j=len-1; j >= 0; j--) {
557 long product = (a[j] & LONG_MASK) * xLong + carry;
558 long difference = q[offset] - product;
559 q[offset--] = (int)difference;
560 carry = (product >>> 32)
561 + (((difference & LONG_MASK) >
562 (((~(int)product) & LONG_MASK))) ? 1:0);
563 }
564 return (int)carry;
565 }
566
567 /**
568 * The method is the same as mulsun, except the fact that q array is not
569 * updated, the only result of the method is borrow flag.
570 */
571 private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
572 long xLong = x & LONG_MASK;
573 long carry = 0;
574 offset += len;
575 for (int j=len-1; j >= 0; j--) {
576 long product = (a[j] & LONG_MASK) * xLong + carry;
577 long difference = q[offset--] - product;
578 carry = (product >>> 32)
579 + (((difference & LONG_MASK) >
580 (((~(int)product) & LONG_MASK))) ? 1:0);
581 }
582 return (int)carry;
583 }
584
585 /**
586 * Right shift this MutableBigInteger n bits, where n is
587 * less than 32.
588 * Assumes that intLen > 0, n > 0 for speed
589 */
590 private final void primitiveRightShift(int n) {
591 int[] val = value;
592 int n2 = 32 - n;
593 for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
594 int b = c;
595 c = val[i-1];
596 val[i] = (c << n2) | (b >>> n);
597 }
598 val[offset] >>>= n;
599 }
600
601 /**
602 * Left shift this MutableBigInteger n bits, where n is
603 * less than 32.
604 * Assumes that intLen > 0, n > 0 for speed
605 */
606 private final void primitiveLeftShift(int n) {
607 int[] val = value;
608 int n2 = 32 - n;
609 for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) {
610 int b = c;
611 c = val[i+1];
612 val[i] = (b << n) | (c >>> n2);
613 }
614 val[offset+intLen-1] <<= n;
615 }
616
617 /**
618 * Adds the contents of two MutableBigInteger objects.The result
619 * is placed within this MutableBigInteger.
620 * The contents of the addend are not changed.
621 */
622 void add(MutableBigInteger addend) {
623 int x = intLen;
624 int y = addend.intLen;
625 int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
626 int[] result = (value.length < resultLen ? new int[resultLen] : value);
627
628 int rstart = result.length-1;
629 long sum;
630 long carry = 0;
631
632 // Add common parts of both numbers
633 while(x>0 && y>0) {
634 x--; y--;
635 sum = (value[x+offset] & LONG_MASK) +
636 (addend.value[y+addend.offset] & LONG_MASK) + carry;
637 result[rstart--] = (int)sum;
638 carry = sum >>> 32;
639 }
640
641 // Add remainder of the longer number
642 while(x>0) {
643 x--;
644 if (carry == 0 && result == value && rstart == (x + offset))
645 return;
646 sum = (value[x+offset] & LONG_MASK) + carry;
647 result[rstart--] = (int)sum;
648 carry = sum >>> 32;
649 }
650 while(y>0) {
651 y--;
652 sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
653 result[rstart--] = (int)sum;
654 carry = sum >>> 32;
655 }
656
657 if (carry > 0) { // Result must grow in length
658 resultLen++;
659 if (result.length < resultLen) {
660 int temp[] = new int[resultLen];
661 // Result one word longer from carry-out; copy low-order
662 // bits into new result.
663 System.arraycopy(result, 0, temp, 1, result.length);
664 temp[0] = 1;
665 result = temp;
666 } else {
667 result[rstart--] = 1;
668 }
669 }
670
671 value = result;
672 intLen = resultLen;
673 offset = result.length - resultLen;
674 }
675
676
677 /**
678 * Subtracts the smaller of this and b from the larger and places the
679 * result into this MutableBigInteger.
680 */
681 int subtract(MutableBigInteger b) {
682 MutableBigInteger a = this;
683
684 int[] result = value;
685 int sign = a.compare(b);
686
687 if (sign == 0) {
688 reset();
689 return 0;
690 }
691 if (sign < 0) {
692 MutableBigInteger tmp = a;
693 a = b;
694 b = tmp;
695 }
696
697 int resultLen = a.intLen;
698 if (result.length < resultLen)
699 result = new int[resultLen];
700
701 long diff = 0;
702 int x = a.intLen;
703 int y = b.intLen;
704 int rstart = result.length - 1;
705
706 // Subtract common parts of both numbers
707 while (y>0) {
708 x--; y--;
709
710 diff = (a.value[x+a.offset] & LONG_MASK) -
711 (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
712 result[rstart--] = (int)diff;
713 }
714 // Subtract remainder of longer number
715 while (x>0) {
716 x--;
717 diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
718 result[rstart--] = (int)diff;
719 }
720
721 value = result;
722 intLen = resultLen;
723 offset = value.length - resultLen;
724 normalize();
725 return sign;
726 }
727
728 /**
729 * Subtracts the smaller of a and b from the larger and places the result
730 * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
731 * operation was performed.
732 */
733 private int difference(MutableBigInteger b) {
734 MutableBigInteger a = this;
735 int sign = a.compare(b);
736 if (sign ==0)
737 return 0;
738 if (sign < 0) {
739 MutableBigInteger tmp = a;
740 a = b;
741 b = tmp;
742 }
743
744 long diff = 0;
745 int x = a.intLen;
746 int y = b.intLen;
747
748 // Subtract common parts of both numbers
749 while (y>0) {
750 x--; y--;
751 diff = (a.value[a.offset+ x] & LONG_MASK) -
752 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
753 a.value[a.offset+x] = (int)diff;
754 }
755 // Subtract remainder of longer number
756 while (x>0) {
757 x--;
758 diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
759 a.value[a.offset+x] = (int)diff;
760 }
761
762 a.normalize();
763 return sign;
764 }
765
766 /**
767 * Multiply the contents of two MutableBigInteger objects. The result is
768 * placed into MutableBigInteger z. The contents of y are not changed.
769 */
770 void multiply(MutableBigInteger y, MutableBigInteger z) {
771 int xLen = intLen;
772 int yLen = y.intLen;
773 int newLen = xLen + yLen;
774
775 // Put z into an appropriate state to receive product
776 if (z.value.length < newLen)
777 z.value = new int[newLen];
778 z.offset = 0;
779 z.intLen = newLen;
780
781 // The first iteration is hoisted out of the loop to avoid extra add
782 long carry = 0;
783 for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
784 long product = (y.value[j+y.offset] & LONG_MASK) *
785 (value[xLen-1+offset] & LONG_MASK) + carry;
786 z.value[k] = (int)product;
787 carry = product >>> 32;
788 }
789 z.value[xLen-1] = (int)carry;
790
791 // Perform the multiplication word by word
792 for (int i = xLen-2; i >= 0; i--) {
793 carry = 0;
794 for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
795 long product = (y.value[j+y.offset] & LONG_MASK) *
796 (value[i+offset] & LONG_MASK) +
797 (z.value[k] & LONG_MASK) + carry;
798 z.value[k] = (int)product;
799 carry = product >>> 32;
800 }
801 z.value[i] = (int)carry;
802 }
803
804 // Remove leading zeros from product
805 z.normalize();
806 }
807
808 /**
809 * Multiply the contents of this MutableBigInteger by the word y. The
810 * result is placed into z.
811 */
812 void mul(int y, MutableBigInteger z) {
813 if (y == 1) {
814 z.copyValue(this);
815 return;
816 }
817
818 if (y == 0) {
819 z.clear();
820 return;
821 }
822
823 // Perform the multiplication word by word
824 long ylong = y & LONG_MASK;
825 int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1]
826 : z.value);
827 long carry = 0;
828 for (int i = intLen-1; i >= 0; i--) {
829 long product = ylong * (value[i+offset] & LONG_MASK) + carry;
830 zval[i+1] = (int)product;
831 carry = product >>> 32;
832 }
833
834 if (carry == 0) {
835 z.offset = 1;
836 z.intLen = intLen;
837 } else {
838 z.offset = 0;
839 z.intLen = intLen + 1;
840 zval[0] = (int)carry;
841 }
842 z.value = zval;
843 }
844
845 /**
846 * This method is used for division of an n word dividend by a one word
847 * divisor. The quotient is placed into quotient. The one word divisor is
848 * specified by divisor.
849 *
850 * @return the remainder of the division is returned.
851 *
852 */
853 int divideOneWord(int divisor, MutableBigInteger quotient) {
854 long divisorLong = divisor & LONG_MASK;
855
856 // Special case of one word dividend
857 if (intLen == 1) {
858 long dividendValue = value[offset] & LONG_MASK;
859 int q = (int) (dividendValue / divisorLong);
860 int r = (int) (dividendValue - q * divisorLong);
861 quotient.value[0] = q;
862 quotient.intLen = (q == 0) ? 0 : 1;
863 quotient.offset = 0;
864 return r;
865 }
866
867 if (quotient.value.length < intLen)
868 quotient.value = new int[intLen];
869 quotient.offset = 0;
870 quotient.intLen = intLen;
871
872 // Normalize the divisor
873 int shift = Integer.numberOfLeadingZeros(divisor);
874
875 int rem = value[offset];
876 long remLong = rem & LONG_MASK;
877 if (remLong < divisorLong) {
878 quotient.value[0] = 0;
879 } else {
880 quotient.value[0] = (int)(remLong / divisorLong);
881 rem = (int) (remLong - (quotient.value[0] * divisorLong));
882 remLong = rem & LONG_MASK;
883 }
884 int xlen = intLen;
885 while (--xlen > 0) {
886 long dividendEstimate = (remLong << 32) |
887 (value[offset + intLen - xlen] & LONG_MASK);
888 int q;
889 if (dividendEstimate >= 0) {
890 q = (int) (dividendEstimate / divisorLong);
891 rem = (int) (dividendEstimate - q * divisorLong);
892 } else {
893 long tmp = divWord(dividendEstimate, divisor);
894 q = (int) (tmp & LONG_MASK);
895 rem = (int) (tmp >>> 32);
896 }
897 quotient.value[intLen - xlen] = q;
898 remLong = rem & LONG_MASK;
899 }
900
901 quotient.normalize();
902 // Unnormalize
903 if (shift > 0)
904 return rem % divisor;
905 else
906 return rem;
907 }
908
909 /**
910 * Calculates the quotient of this div b and places the quotient in the
911 * provided MutableBigInteger objects and the remainder object is returned.
912 *
913 * Uses Algorithm D in Knuth section 4.3.1.
914 * Many optimizations to that algorithm have been adapted from the Colin
915 * Plumb C library.
916 * It special cases one word divisors for speed. The content of b is not
917 * changed.
918 *
919 */
920 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
921 return divide(b,quotient,true);
922 }
923
924 MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needReminder) {
925 if (b.intLen == 0)
926 throw new ArithmeticException("BigInteger divide by zero");
927
928 // Dividend is zero
929 if (intLen == 0) {
930 quotient.intLen = quotient.offset;
931 return needReminder ? new MutableBigInteger() : null;
932 }
933
934 int cmp = compare(b);
935 // Dividend less than divisor
936 if (cmp < 0) {
937 quotient.intLen = quotient.offset = 0;
938 return needReminder ? new MutableBigInteger(this) : null;
939 }
940 // Dividend equal to divisor
941 if (cmp == 0) {
942 quotient.value[0] = quotient.intLen = 1;
943 quotient.offset = 0;
944 return needReminder ? new MutableBigInteger() : null;
945 }
946
947 quotient.clear();
948 // Special case one word divisor
949 if (b.intLen == 1) {
950 int r = divideOneWord(b.value[b.offset], quotient);
951 if(needReminder) {
952 if (r == 0)
953 return new MutableBigInteger();
954 return new MutableBigInteger(r);
955 } else {
956 return null;
957 }
958 }
959 return divideMagnitude(b, quotient, needReminder);
960 }
961
962 /**
963 * Internally used to calculate the quotient of this div v and places the
964 * quotient in the provided MutableBigInteger object and the remainder is
965 * returned.
966 *
967 * @return the remainder of the division will be returned.
968 */
969 long divide(long v, MutableBigInteger quotient) {
970 if (v == 0)
971 throw new ArithmeticException("BigInteger divide by zero");
972
973 // Dividend is zero
974 if (intLen == 0) {
975 quotient.intLen = quotient.offset = 0;
976 return 0;
977 }
978 if (v < 0)
979 v = -v;
980
981 int d = (int)(v >>> 32);
982 quotient.clear();
983 // Special case on word divisor
984 if (d == 0)
985 return divideOneWord((int)v, quotient) & LONG_MASK;
986 else {
987 return divideLongMagnitude(v, quotient).toLong();
988 }
989 }
990
991 private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) {
992 int n2 = 32 - shift;
993 int c=src[srcFrom];
994 for (int i=0; i < srcLen-1; i++) {
995 int b = c;
996 c = src[++srcFrom];
997 dst[dstFrom+i] = (b << shift) | (c >>> n2);
998 }
999 dst[dstFrom+srcLen-1] = c << shift;
1000 }
1001
1002 /**
1003 * Divide this MutableBigInteger by the divisor.
1004 * The quotient will be placed into the provided quotient object &
1005 * the remainder object is returned.
1006 */
1007 private MutableBigInteger divideMagnitude(MutableBigInteger div,
1008 MutableBigInteger quotient,
1009 boolean needReminder ) {
1010 // assert div.intLen > 1
1011 // D1 normalize the divisor
1012 int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
1013 // Copy divisor value to protect divisor
1014 final int dlen = div.intLen;
1015 int[] divisor;
1016 MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
1017 if (shift > 0) {
1018 divisor = new int[dlen];
1019 copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
1020 if(Integer.numberOfLeadingZeros(value[offset])>=shift) {
1021 int[] remarr = new int[intLen + 1];
1022 rem = new MutableBigInteger(remarr);
1023 rem.intLen = intLen;
1024 rem.offset = 1;
1025 copyAndShift(value,offset,intLen,remarr,1,shift);
1026 } else {
1027 int[] remarr = new int[intLen + 2];
1028 rem = new MutableBigInteger(remarr);
1029 rem.intLen = intLen+1;
1030 rem.offset = 1;
1031 int rFrom = offset;
1032 int c=0;
1033 int n2 = 32 - shift;
1034 for (int i=1; i < intLen+1; i++,rFrom++) {
1035 int b = c;
1036 c = value[rFrom];
1037 remarr[i] = (b << shift) | (c >>> n2);
1038 }
1039 remarr[intLen+1] = c << shift;
1040 }
1041 } else {
1042 divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
1043 rem = new MutableBigInteger(new int[intLen + 1]);
1044 System.arraycopy(value, offset, rem.value, 1, intLen);
1045 rem.intLen = intLen;
1046 rem.offset = 1;
1047 }
1048
1049 int nlen = rem.intLen;
1050
1051 // Set the quotient size
1052 final int limit = nlen - dlen + 1;
1053 if (quotient.value.length < limit) {
1054 quotient.value = new int[limit];
1055 quotient.offset = 0;
1056 }
1057 quotient.intLen = limit;
1058 int[] q = quotient.value;
1059
1060
1061 // Must insert leading 0 in rem if its length did not change
1062 if (rem.intLen == nlen) {
1063 rem.offset = 0;
1064 rem.value[0] = 0;
1065 rem.intLen++;
1066 }
1067
1068 int dh = divisor[0];
1069 long dhLong = dh & LONG_MASK;
1070 int dl = divisor[1];
1071
1072 // D2 Initialize j
1073 for(int j=0; j<limit-1; j++) {
1074 // D3 Calculate qhat
1075 // estimate qhat
1076 int qhat = 0;
1077 int qrem = 0;
1078 boolean skipCorrection = false;
1079 int nh = rem.value[j+rem.offset];
1080 int nh2 = nh + 0x80000000;
1081 int nm = rem.value[j+1+rem.offset];
1082
1083 if (nh == dh) {
1084 qhat = ~0;
1085 qrem = nh + nm;
1086 skipCorrection = qrem + 0x80000000 < nh2;
1087 } else {
1088 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
1089 if (nChunk >= 0) {
1090 qhat = (int) (nChunk / dhLong);
1091 qrem = (int) (nChunk - (qhat * dhLong));
1092 } else {
1093 long tmp = divWord(nChunk, dh);
1094 qhat = (int) (tmp & LONG_MASK);
1095 qrem = (int) (tmp >>> 32);
1096 }
1097 }
1098
1099 if (qhat == 0)
1100 continue;
1101
1102 if (!skipCorrection) { // Correct qhat
1103 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
1104 long rs = ((qrem & LONG_MASK) << 32) | nl;
1105 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1106
1107 if (unsignedLongCompare(estProduct, rs)) {
1108 qhat--;
1109 qrem = (int)((qrem & LONG_MASK) + dhLong);
1110 if ((qrem & LONG_MASK) >= dhLong) {
1111 estProduct -= (dl & LONG_MASK);
1112 rs = ((qrem & LONG_MASK) << 32) | nl;
1113 if (unsignedLongCompare(estProduct, rs))
1114 qhat--;
1115 }
1116 }
1117 }
1118
1119 // D4 Multiply and subtract
1120 rem.value[j+rem.offset] = 0;
1121 int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
1122
1123 // D5 Test remainder
1124 if (borrow + 0x80000000 > nh2) {
1125 // D6 Add back
1126 divadd(divisor, rem.value, j+1+rem.offset);
1127 qhat--;
1128 }
1129
1130 // Store the quotient digit
1131 q[j] = qhat;
1132 } // D7 loop on j
1133 // D3 Calculate qhat
1134 // estimate qhat
1135 int qhat = 0;
1136 int qrem = 0;
1137 boolean skipCorrection = false;
1138 int nh = rem.value[limit - 1 + rem.offset];
1139 int nh2 = nh + 0x80000000;
1140 int nm = rem.value[limit + rem.offset];
1141
1142 if (nh == dh) {
1143 qhat = ~0;
1144 qrem = nh + nm;
1145 skipCorrection = qrem + 0x80000000 < nh2;
1146 } else {
1147 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1148 if (nChunk >= 0) {
1149 qhat = (int) (nChunk / dhLong);
1150 qrem = (int) (nChunk - (qhat * dhLong));
1151 } else {
1152 long tmp = divWord(nChunk, dh);
1153 qhat = (int) (tmp & LONG_MASK);
1154 qrem = (int) (tmp >>> 32);
1155 }
1156 }
1157 if (qhat != 0) {
1158 if (!skipCorrection) { // Correct qhat
1159 long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
1160 long rs = ((qrem & LONG_MASK) << 32) | nl;
1161 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1162
1163 if (unsignedLongCompare(estProduct, rs)) {
1164 qhat--;
1165 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1166 if ((qrem & LONG_MASK) >= dhLong) {
1167 estProduct -= (dl & LONG_MASK);
1168 rs = ((qrem & LONG_MASK) << 32) | nl;
1169 if (unsignedLongCompare(estProduct, rs))
1170 qhat--;
1171 }
1172 }
1173 }
1174
1175
1176 // D4 Multiply and subtract
1177 int borrow;
1178 rem.value[limit - 1 + rem.offset] = 0;
1179 if(needReminder)
1180 borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1181 else
1182 borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
1183
1184 // D5 Test remainder
1185 if (borrow + 0x80000000 > nh2) {
1186 // D6 Add back
1187 if(needReminder)
1188 divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
1189 qhat--;
1190 }
1191
1192 // Store the quotient digit
1193 q[(limit - 1)] = qhat;
1194 }
1195
1196
1197 if(needReminder) {
1198 // D8 Unnormalize
1199 if (shift > 0)
1200 rem.rightShift(shift);
1201 rem.normalize();
1202 }
1203 quotient.normalize();
1204 return needReminder ? rem : null;
1205 }
1206
1207 /**
1208 * Divide this MutableBigInteger by the divisor represented by positive long
1209 * value. The quotient will be placed into the provided quotient object &
1210 * the remainder object is returned.
1211 */
1212 private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
1213 // Remainder starts as dividend with space for a leading zero
1214 MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
1215 System.arraycopy(value, offset, rem.value, 1, intLen);
1216 rem.intLen = intLen;
1217 rem.offset = 1;
1218
1219 int nlen = rem.intLen;
1220
1221 int limit = nlen - 2 + 1;
1222 if (quotient.value.length < limit) {
1223 quotient.value = new int[limit];
1224 quotient.offset = 0;
1225 }
1226 quotient.intLen = limit;
1227 int[] q = quotient.value;
1228
1229 // D1 normalize the divisor
1230 int shift = Long.numberOfLeadingZeros(ldivisor);
1231 if (shift > 0) {
1232 ldivisor<<=shift;
1233 rem.leftShift(shift);
1234 }
1235
1236 // Must insert leading 0 in rem if its length did not change
1237 if (rem.intLen == nlen) {
1238 rem.offset = 0;
1239 rem.value[0] = 0;
1240 rem.intLen++;
1241 }
1242
1243 int dh = (int)(ldivisor >>> 32);
1244 long dhLong = dh & LONG_MASK;
1245 int dl = (int)(ldivisor & LONG_MASK);
1246
1247 // D2 Initialize j
1248 for (int j = 0; j < limit; j++) {
1249 // D3 Calculate qhat
1250 // estimate qhat
1251 int qhat = 0;
1252 int qrem = 0;
1253 boolean skipCorrection = false;
1254 int nh = rem.value[j + rem.offset];
1255 int nh2 = nh + 0x80000000;
1256 int nm = rem.value[j + 1 + rem.offset];
1257
1258 if (nh == dh) {
1259 qhat = ~0;
1260 qrem = nh + nm;
1261 skipCorrection = qrem + 0x80000000 < nh2;
1262 } else {
1263 long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
1264 if (nChunk >= 0) {
1265 qhat = (int) (nChunk / dhLong);
1266 qrem = (int) (nChunk - (qhat * dhLong));
1267 } else {
1268 long tmp = divWord(nChunk, dh);
1269 qhat =(int)(tmp & LONG_MASK);
1270 qrem = (int)(tmp>>>32);
1271 }
1272 }
1273
1274 if (qhat == 0)
1275 continue;
1276
1277 if (!skipCorrection) { // Correct qhat
1278 long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
1279 long rs = ((qrem & LONG_MASK) << 32) | nl;
1280 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
1281
1282 if (unsignedLongCompare(estProduct, rs)) {
1283 qhat--;
1284 qrem = (int) ((qrem & LONG_MASK) + dhLong);
1285 if ((qrem & LONG_MASK) >= dhLong) {
1286 estProduct -= (dl & LONG_MASK);
1287 rs = ((qrem & LONG_MASK) << 32) | nl;
1288 if (unsignedLongCompare(estProduct, rs))
1289 qhat--;
1290 }
1291 }
1292 }
1293
1294 // D4 Multiply and subtract
1295 rem.value[j + rem.offset] = 0;
1296 int borrow = mulsubLong(rem.value, dh, dl, qhat, j + rem.offset);
1297
1298 // D5 Test remainder
1299 if (borrow + 0x80000000 > nh2) {
1300 // D6 Add back
1301 divaddLong(dh,dl, rem.value, j + 1 + rem.offset);
1302 qhat--;
1303 }
1304
1305 // Store the quotient digit
1306 q[j] = qhat;
1307 } // D7 loop on j
1308
1309 // D8 Unnormalize
1310 if (shift > 0)
1311 rem.rightShift(shift);
1312
1313 quotient.normalize();
1314 rem.normalize();
1315 return rem;
1316 }
1317
1318 /**
1319 * A primitive used for division by long.
1320 * Specialized version of the method divadd.
1321 * dh is a high part of the divisor, dl is a low part
1322 */
1323 private int divaddLong(int dh, int dl, int[] result, int offset) {
1324 long carry = 0;
1325
1326 long sum = (dl & LONG_MASK) + (result[1+offset] & LONG_MASK);
1327 result[1+offset] = (int)sum;
1328
1329 sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry;
1330 result[offset] = (int)sum;
1331 carry = sum >>> 32;
1332 return (int)carry;
1333 }
1334
1335 /**
1336 * This method is used for division by long.
1337 * Specialized version of the method sulsub.
1338 * dh is a high part of the divisor, dl is a low part
1339 */
1340 private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
1341 long xLong = x & LONG_MASK;
1342 offset += 2;
1343 long product = (dl & LONG_MASK) * xLong;
1344 long difference = q[offset] - product;
1345 q[offset--] = (int)difference;
1346 long carry = (product >>> 32)
1347 + (((difference & LONG_MASK) >
1348 (((~(int)product) & LONG_MASK))) ? 1:0);
1349 product = (dh & LONG_MASK) * xLong + carry;
1350 difference = q[offset] - product;
1351 q[offset--] = (int)difference;
1352 carry = (product >>> 32)
1353 + (((difference & LONG_MASK) >
1354 (((~(int)product) & LONG_MASK))) ? 1:0);
1355 return (int)carry;
1356 }
1357
1358 /**
1359 * Compare two longs as if they were unsigned.
1360 * Returns true iff one is bigger than two.
1361 */
1362 private boolean unsignedLongCompare(long one, long two) {
1363 return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
1364 }
1365
1366 /**
1367 * This method divides a long quantity by an int to estimate
1368 * qhat for two multi precision numbers. It is used when
1369 * the signed value of n is less than zero.
1370 * Returns long value where high 32 bits contain reminder value and
1371 * low 32 bits contain quotient value.
1372 */
1373 static long divWord(long n, int d) {
1374 long dLong = d & LONG_MASK;
1375 long r;
1376 long q;
1377 if (dLong == 1) {
1378 q = (int)n;
1379 r = 0;
1380 return (r << 32) | (q & LONG_MASK);
1381 }
1382
1383 // Approximate the quotient and remainder
1384 q = (n >>> 1) / (dLong >>> 1);
1385 r = n - q*dLong;
1386
1387 // Correct the approximation
1388 while (r < 0) {
1389 r += dLong;
1390 q--;
1391 }
1392 while (r >= dLong) {
1393 r -= dLong;
1394 q++;
1395 }
1396 // n - q*dlong == r && 0 <= r <dLong, hence we're done.
1397 return (r << 32) | (q & LONG_MASK);
1398 }
1399
1400 /**
1401 * Calculate GCD of this and b. This and b are changed by the computation.
1402 */
1403 MutableBigInteger hybridGCD(MutableBigInteger b) {
1404 // Use Euclid's algorithm until the numbers are approximately the
1405 // same length, then use the binary GCD algorithm to find the GCD.
1406 MutableBigInteger a = this;
1407 MutableBigInteger q = new MutableBigInteger();
1408
1409 while (b.intLen != 0) {
1410 if (Math.abs(a.intLen - b.intLen) < 2)
1411 return a.binaryGCD(b);
1412
1413 MutableBigInteger r = a.divide(b, q);
1414 a = b;
1415 b = r;
1416 }
1417 return a;
1418 }
1419
1420 /**
1421 * Calculate GCD of this and v.
1422 * Assumes that this and v are not zero.
1423 */
1424 private MutableBigInteger binaryGCD(MutableBigInteger v) {
1425 // Algorithm B from Knuth section 4.5.2
1426 MutableBigInteger u = this;
1427 MutableBigInteger r = new MutableBigInteger();
1428
1429 // step B1
1430 int s1 = u.getLowestSetBit();
1431 int s2 = v.getLowestSetBit();
1432 int k = (s1 < s2) ? s1 : s2;
1433 if (k != 0) {
1434 u.rightShift(k);
1435 v.rightShift(k);
1436 }
1437
1438 // step B2
1439 boolean uOdd = (k==s1);
1440 MutableBigInteger t = uOdd ? v: u;
1441 int tsign = uOdd ? -1 : 1;
1442
1443 int lb;
1444 while ((lb = t.getLowestSetBit()) >= 0) {
1445 // steps B3 and B4
1446 t.rightShift(lb);
1447 // step B5
1448 if (tsign > 0)
1449 u = t;
1450 else
1451 v = t;
1452
1453 // Special case one word numbers
1454 if (u.intLen < 2 && v.intLen < 2) {
1455 int x = u.value[u.offset];
1456 int y = v.value[v.offset];
1457 x = binaryGcd(x, y);
1458 r.value[0] = x;
1459 r.intLen = 1;
1460 r.offset = 0;
1461 if (k > 0)
1462 r.leftShift(k);
1463 return r;
1464 }
1465
1466 // step B6
1467 if ((tsign = u.difference(v)) == 0)
1468 break;
1469 t = (tsign >= 0) ? u : v;
1470 }
1471
1472 if (k > 0)
1473 u.leftShift(k);
1474 return u;
1475 }
1476
1477 /**
1478 * Calculate GCD of a and b interpreted as unsigned integers.
1479 */
1480 static int binaryGcd(int a, int b) {
1481 if (b==0)
1482 return a;
1483 if (a==0)
1484 return b;
1485
1486 // Right shift a & b till their last bits equal to 1.
1487 int aZeros = Integer.numberOfTrailingZeros(a);
1488 int bZeros = Integer.numberOfTrailingZeros(b);
1489 a >>>= aZeros;
1490 b >>>= bZeros;
1491
1492 int t = (aZeros < bZeros ? aZeros : bZeros);
1493
1494 while (a != b) {
1495 if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned
1496 a -= b;
1497 a >>>= Integer.numberOfTrailingZeros(a);
1498 } else {
1499 b -= a;
1500 b >>>= Integer.numberOfTrailingZeros(b);
1501 }
1502 }
1503 return a<<t;
1504 }
1505
1506 /**
1507 * Returns the modInverse of this mod p. This and p are not affected by
1508 * the operation.
1509 */
1510 MutableBigInteger mutableModInverse(MutableBigInteger p) {
1511 // Modulus is odd, use Schroeppel's algorithm
1512 if (p.isOdd())
1513 return modInverse(p);
1514
1515 // Base and modulus are even, throw exception
1516 if (isEven())
1517 throw new ArithmeticException("BigInteger not invertible.");
1518
1519 // Get even part of modulus expressed as a power of 2
1520 int powersOf2 = p.getLowestSetBit();
1521
1522 // Construct odd part of modulus
1523 MutableBigInteger oddMod = new MutableBigInteger(p);
1524 oddMod.rightShift(powersOf2);
1525
1526 if (oddMod.isOne())
1527 return modInverseMP2(powersOf2);
1528
1529 // Calculate 1/a mod oddMod
1530 MutableBigInteger oddPart = modInverse(oddMod);
1531
1532 // Calculate 1/a mod evenMod
1533 MutableBigInteger evenPart = modInverseMP2(powersOf2);
1534
1535 // Combine the results using Chinese Remainder Theorem
1536 MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
1537 MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
1538
1539 MutableBigInteger temp1 = new MutableBigInteger();
1540 MutableBigInteger temp2 = new MutableBigInteger();
1541 MutableBigInteger result = new MutableBigInteger();
1542
1543 oddPart.leftShift(powersOf2);
1544 oddPart.multiply(y1, result);
1545
1546 evenPart.multiply(oddMod, temp1);
1547 temp1.multiply(y2, temp2);
1548
1549 result.add(temp2);
1550 return result.divide(p, temp1);
1551 }
1552
1553 /*
1554 * Calculate the multiplicative inverse of this mod 2^k.
1555 */
1556 MutableBigInteger modInverseMP2(int k) {
1557 if (isEven())
1558 throw new ArithmeticException("Non-invertible. (GCD != 1)");
1559
1560 if (k > 64)
1561 return euclidModInverse(k);
1562
1563 int t = inverseMod32(value[offset+intLen-1]);
1564
1565 if (k < 33) {
1566 t = (k == 32 ? t : t & ((1 << k) - 1));
1567 return new MutableBigInteger(t);
1568 }
1569
1570 long pLong = (value[offset+intLen-1] & LONG_MASK);
1571 if (intLen > 1)
1572 pLong |= ((long)value[offset+intLen-2] << 32);
1573 long tLong = t & LONG_MASK;
1574 tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
1575 tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
1576
1577 MutableBigInteger result = new MutableBigInteger(new int[2]);
1578 result.value[0] = (int)(tLong >>> 32);
1579 result.value[1] = (int)tLong;
1580 result.intLen = 2;
1581 result.normalize();
1582 return result;
1583 }
1584
1585 /*
1586 * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
1587 */
1588 static int inverseMod32(int val) {
1589 // Newton's iteration!
1590 int t = val;
1591 t *= 2 - val*t;
1592 t *= 2 - val*t;
1593 t *= 2 - val*t;
1594 t *= 2 - val*t;
1595 return t;
1596 }
1597
1598 /*
1599 * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
1600 */
1601 static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
1602 // Copy the mod to protect original
1603 return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
1604 }
1605
1606 /**
1607 * Calculate the multiplicative inverse of this mod mod, where mod is odd.
1608 * This and mod are not changed by the calculation.
1609 *
1610 * This method implements an algorithm due to Richard Schroeppel, that uses
1611 * the same intermediate representation as Montgomery Reduction
1612 * ("Montgomery Form"). The algorithm is described in an unpublished
1613 * manuscript entitled "Fast Modular Reciprocals."
1614 */
1615 private MutableBigInteger modInverse(MutableBigInteger mod) {
1616 MutableBigInteger p = new MutableBigInteger(mod);
1617 MutableBigInteger f = new MutableBigInteger(this);
1618 MutableBigInteger g = new MutableBigInteger(p);
1619 SignedMutableBigInteger c = new SignedMutableBigInteger(1);
1620 SignedMutableBigInteger d = new SignedMutableBigInteger();
1621 MutableBigInteger temp = null;
1622 SignedMutableBigInteger sTemp = null;
1623
1624 int k = 0;
1625 // Right shift f k times until odd, left shift d k times
1626 if (f.isEven()) {
1627 int trailingZeros = f.getLowestSetBit();
1628 f.rightShift(trailingZeros);
1629 d.leftShift(trailingZeros);
1630 k = trailingZeros;
1631 }
1632
1633 // The Almost Inverse Algorithm
1634 while(!f.isOne()) {
1635 // If gcd(f, g) != 1, number is not invertible modulo mod
1636 if (f.isZero())
1637 throw new ArithmeticException("BigInteger not invertible.");
1638
1639 // If f < g exchange f, g and c, d
1640 if (f.compare(g) < 0) {
1641 temp = f; f = g; g = temp;
1642 sTemp = d; d = c; c = sTemp;
1643 }
1644
1645 // If f == g (mod 4)
1646 if (((f.value[f.offset + f.intLen - 1] ^
1647 g.value[g.offset + g.intLen - 1]) & 3) == 0) {
1648 f.subtract(g);
1649 c.signedSubtract(d);
1650 } else { // If f != g (mod 4)
1651 f.add(g);
1652 c.signedAdd(d);
1653 }
1654
1655 // Right shift f k times until odd, left shift d k times
1656 int trailingZeros = f.getLowestSetBit();
1657 f.rightShift(trailingZeros);
1658 d.leftShift(trailingZeros);
1659 k += trailingZeros;
1660 }
1661
1662 while (c.sign < 0)
1663 c.signedAdd(p);
1664
1665 return fixup(c, p, k);
1666 }
1667
1668 /*
1669 * The Fixup Algorithm
1670 * Calculates X such that X = C * 2^(-k) (mod P)
1671 * Assumes C<P and P is odd.
1672 */
1673 static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
1674 int k) {
1675 MutableBigInteger temp = new MutableBigInteger();
1676 // Set r to the multiplicative inverse of p mod 2^32
1677 int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
1678
1679 for(int i=0, numWords = k >> 5; i<numWords; i++) {
1680 // V = R * c (mod 2^j)
1681 int v = r * c.value[c.offset + c.intLen-1];
1682 // c = c + (v * p)
1683 p.mul(v, temp);
1684 c.add(temp);
1685 // c = c / 2^j
1686 c.intLen--;
1687 }
1688 int numBits = k & 0x1f;
1689 if (numBits != 0) {
1690 // V = R * c (mod 2^j)
1691 int v = r * c.value[c.offset + c.intLen-1];
1692 v &= ((1<<numBits) - 1);
1693 // c = c + (v * p)
1694 p.mul(v, temp);
1695 c.add(temp);
1696 // c = c / 2^j
1697 c.rightShift(numBits);
1698 }
1699
1700 // In theory, c may be greater than p at this point (Very rare!)
1701 while (c.compare(p) >= 0)
1702 c.subtract(p);
1703
1704 return c;
1705 }
1706
1707 /**
1708 * Uses the extended Euclidean algorithm to compute the modInverse of base
1709 * mod a modulus that is a power of 2. The modulus is 2^k.
1710 */
1711 MutableBigInteger euclidModInverse(int k) {
1712 MutableBigInteger b = new MutableBigInteger(1);
1713 b.leftShift(k);
1714 MutableBigInteger mod = new MutableBigInteger(b);
1715
1716 MutableBigInteger a = new MutableBigInteger(this);
1717 MutableBigInteger q = new MutableBigInteger();
1718 MutableBigInteger r = b.divide(a, q);
1719
1720 MutableBigInteger swapper = b;
1721 // swap b & r
1722 b = r;
1723 r = swapper;
1724
1725 MutableBigInteger t1 = new MutableBigInteger(q);
1726 MutableBigInteger t0 = new MutableBigInteger(1);
1727 MutableBigInteger temp = new MutableBigInteger();
1728
1729 while (!b.isOne()) {
1730 r = a.divide(b, q);
1731
1732 if (r.intLen == 0)
1733 throw new ArithmeticException("BigInteger not invertible.");
1734
1735 swapper = r;
1736 a = swapper;
1737
1738 if (q.intLen == 1)
1739 t1.mul(q.value[q.offset], temp);
1740 else
1741 q.multiply(t1, temp);
1742 swapper = q;
1743 q = temp;
1744 temp = swapper;
1745 t0.add(q);
1746
1747 if (a.isOne())
1748 return t0;
1749
1750 r = b.divide(a, q);
1751
1752 if (r.intLen == 0)
1753 throw new ArithmeticException("BigInteger not invertible.");
1754
1755 swapper = b;
1756 b = r;
1757
1758 if (q.intLen == 1)
1759 t0.mul(q.value[q.offset], temp);
1760 else
1761 q.multiply(t0, temp);
1762 swapper = q; q = temp; temp = swapper;
1763
1764 t1.add(q);
1765 }
1766 mod.subtract(t1);
1767 return mod;
1768 }
1769 }