View Javadoc
1   /*
2    * Copyright (c) 1999, 2013, Oracle and/or its affiliates. All rights reserved.
3    * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
4    *
5    * This code is free software; you can redistribute it and/or modify it
6    * under the terms of the GNU General Public License version 2 only, as
7    * published by the Free Software Foundation.  Oracle designates this
8    * particular file as subject to the "Classpath" exception as provided
9    * by Oracle in the LICENSE file that accompanied this code.
10   *
11   * This code is distributed in the hope that it will be useful, but WITHOUT
12   * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13   * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
14   * version 2 for more details (a copy is included in the LICENSE file that
15   * accompanied this code).
16   *
17   * You should have received a copy of the GNU General Public License version
18   * 2 along with this work; if not, write to the Free Software Foundation,
19   * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
20   *
21   * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
22   * or visit www.oracle.com if you need additional information or have any
23   * questions.
24   */
25  
26  package java.math;
27  
28  /**
29   * A class used to represent multiprecision integers that makes efficient
30   * use of allocated space by allowing a number to occupy only part of
31   * an array so that the arrays do not have to be reallocated as often.
32   * When performing an operation with many iterations the array used to
33   * hold a number is only reallocated when necessary and does not have to
34   * be the same size as the number it represents. A mutable number allows
35   * calculations to occur on the same number without having to create
36   * a new number for every step of the calculation as occurs with
37   * BigIntegers.
38   *
39   * @see     BigInteger
40   * @author  Michael McCloskey
41   * @author  Timothy Buktu
42   * @since   1.3
43   */
44  
45  import static java.math.BigDecimal.INFLATED;
46  import static java.math.BigInteger.LONG_MASK;
47  import java.util.Arrays;
48  
49  class MutableBigInteger {
50      /**
51       * Holds the magnitude of this MutableBigInteger in big endian order.
52       * The magnitude may start at an offset into the value array, and it may
53       * end before the length of the value array.
54       */
55      int[] value;
56  
57      /**
58       * The number of ints of the value array that are currently used
59       * to hold the magnitude of this MutableBigInteger. The magnitude starts
60       * at an offset and offset + intLen may be less than value.length.
61       */
62      int intLen;
63  
64      /**
65       * The offset into the value array where the magnitude of this
66       * MutableBigInteger begins.
67       */
68      int offset = 0;
69  
70      // Constants
71      /**
72       * MutableBigInteger with one element value array with the value 1. Used by
73       * BigDecimal divideAndRound to increment the quotient. Use this constant
74       * only when the method is not going to modify this object.
75       */
76      static final MutableBigInteger ONE = new MutableBigInteger(1);
77  
78      /**
79       * The minimum {@code intLen} for cancelling powers of two before
80       * dividing.
81       * If the number of ints is less than this threshold,
82       * {@code divideKnuth} does not eliminate common powers of two from
83       * the dividend and divisor.
84       */
85      static final int KNUTH_POW2_THRESH_LEN = 6;
86  
87      /**
88       * The minimum number of trailing zero ints for cancelling powers of two
89       * before dividing.
90       * If the dividend and divisor don't share at least this many zero ints
91       * at the end, {@code divideKnuth} does not eliminate common powers
92       * of two from the dividend and divisor.
93       */
94      static final int KNUTH_POW2_THRESH_ZEROS = 3;
95  
96      // Constructors
97  
98      /**
99       * The default constructor. An empty MutableBigInteger is created with
100      * a one word capacity.
101      */
102     MutableBigInteger() {
103         value = new int[1];
104         intLen = 0;
105     }
106 
107     /**
108      * Construct a new MutableBigInteger with a magnitude specified by
109      * the int val.
110      */
111     MutableBigInteger(int val) {
112         value = new int[1];
113         intLen = 1;
114         value[0] = val;
115     }
116 
117     /**
118      * Construct a new MutableBigInteger with the specified value array
119      * up to the length of the array supplied.
120      */
121     MutableBigInteger(int[] val) {
122         value = val;
123         intLen = val.length;
124     }
125 
126     /**
127      * Construct a new MutableBigInteger with a magnitude equal to the
128      * specified BigInteger.
129      */
130     MutableBigInteger(BigInteger b) {
131         intLen = b.mag.length;
132         value = Arrays.copyOf(b.mag, intLen);
133     }
134 
135     /**
136      * Construct a new MutableBigInteger with a magnitude equal to the
137      * specified MutableBigInteger.
138      */
139     MutableBigInteger(MutableBigInteger val) {
140         intLen = val.intLen;
141         value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
142     }
143 
144     /**
145      * Makes this number an {@code n}-int number all of whose bits are ones.
146      * Used by Burnikel-Ziegler division.
147      * @param n number of ints in the {@code value} array
148      * @return a number equal to {@code ((1<<(32*n)))-1}
149      */
150     private void ones(int n) {
151         if (n > value.length)
152             value = new int[n];
153         Arrays.fill(value, -1);
154         offset = 0;
155         intLen = n;
156     }
157 
158     /**
159      * Internal helper method to return the magnitude array. The caller is not
160      * supposed to modify the returned array.
161      */
162     private int[] getMagnitudeArray() {
163         if (offset > 0 || value.length != intLen)
164             return Arrays.copyOfRange(value, offset, offset + intLen);
165         return value;
166     }
167 
168     /**
169      * Convert this MutableBigInteger to a long value. The caller has to make
170      * sure this MutableBigInteger can be fit into long.
171      */
172     private long toLong() {
173         assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
174         if (intLen == 0)
175             return 0;
176         long d = value[offset] & LONG_MASK;
177         return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
178     }
179 
180     /**
181      * Convert this MutableBigInteger to a BigInteger object.
182      */
183     BigInteger toBigInteger(int sign) {
184         if (intLen == 0 || sign == 0)
185             return BigInteger.ZERO;
186         return new BigInteger(getMagnitudeArray(), sign);
187     }
188 
189     /**
190      * Converts this number to a nonnegative {@code BigInteger}.
191      */
192     BigInteger toBigInteger() {
193         normalize();
194         return toBigInteger(isZero() ? 0 : 1);
195     }
196 
197     /**
198      * Convert this MutableBigInteger to BigDecimal object with the specified sign
199      * and scale.
200      */
201     BigDecimal toBigDecimal(int sign, int scale) {
202         if (intLen == 0 || sign == 0)
203             return BigDecimal.zeroValueOf(scale);
204         int[] mag = getMagnitudeArray();
205         int len = mag.length;
206         int d = mag[0];
207         // If this MutableBigInteger can't be fit into long, we need to
208         // make a BigInteger object for the resultant BigDecimal object.
209         if (len > 2 || (d < 0 && len == 2))
210             return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
211         long v = (len == 2) ?
212             ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
213             d & LONG_MASK;
214         return BigDecimal.valueOf(sign == -1 ? -v : v, scale);
215     }
216 
217     /**
218      * This is for internal use in converting from a MutableBigInteger
219      * object into a long value given a specified sign.
220      * returns INFLATED if value is not fit into long
221      */
222     long toCompactValue(int sign) {
223         if (intLen == 0 || sign == 0)
224             return 0L;
225         int[] mag = getMagnitudeArray();
226         int len = mag.length;
227         int d = mag[0];
228         // If this MutableBigInteger can not be fitted into long, we need to
229         // make a BigInteger object for the resultant BigDecimal object.
230         if (len > 2 || (d < 0 && len == 2))
231             return INFLATED;
232         long v = (len == 2) ?
233             ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
234             d & LONG_MASK;
235         return sign == -1 ? -v : v;
236     }
237 
238     /**
239      * Clear out a MutableBigInteger for reuse.
240      */
241     void clear() {
242         offset = intLen = 0;
243         for (int index=0, n=value.length; index < n; index++)
244             value[index] = 0;
245     }
246 
247     /**
248      * Set a MutableBigInteger to zero, removing its offset.
249      */
250     void reset() {
251         offset = intLen = 0;
252     }
253 
254     /**
255      * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
256      * as this MutableBigInteger is numerically less than, equal to, or
257      * greater than <tt>b</tt>.
258      */
259     final int compare(MutableBigInteger b) {
260         int blen = b.intLen;
261         if (intLen < blen)
262             return -1;
263         if (intLen > blen)
264            return 1;
265 
266         // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
267         // comparison.
268         int[] bval = b.value;
269         for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
270             int b1 = value[i] + 0x80000000;
271             int b2 = bval[j]  + 0x80000000;
272             if (b1 < b2)
273                 return -1;
274             if (b1 > b2)
275                 return 1;
276         }
277         return 0;
278     }
279 
280     /**
281      * Returns a value equal to what {@code b.leftShift(32*ints); return compare(b);}
282      * would return, but doesn't change the value of {@code b}.
283      */
284     private int compareShifted(MutableBigInteger b, int ints) {
285         int blen = b.intLen;
286         int alen = intLen - ints;
287         if (alen < blen)
288             return -1;
289         if (alen > blen)
290            return 1;
291 
292         // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
293         // comparison.
294         int[] bval = b.value;
295         for (int i = offset, j = b.offset; i < alen + offset; i++, j++) {
296             int b1 = value[i] + 0x80000000;
297             int b2 = bval[j]  + 0x80000000;
298             if (b1 < b2)
299                 return -1;
300             if (b1 > b2)
301                 return 1;
302         }
303         return 0;
304     }
305 
306     /**
307      * Compare this against half of a MutableBigInteger object (Needed for
308      * remainder tests).
309      * Assumes no leading unnecessary zeros, which holds for results
310      * from divide().
311      */
312     final int compareHalf(MutableBigInteger b) {
313         int blen = b.intLen;
314         int len = intLen;
315         if (len <= 0)
316             return blen <= 0 ? 0 : -1;
317         if (len > blen)
318             return 1;
319         if (len < blen - 1)
320             return -1;
321         int[] bval = b.value;
322         int bstart = 0;
323         int carry = 0;
324         // Only 2 cases left:len == blen or len == blen - 1
325         if (len != blen) { // len == blen - 1
326             if (bval[bstart] == 1) {
327                 ++bstart;
328                 carry = 0x80000000;
329             } else
330                 return -1;
331         }
332         // compare values with right-shifted values of b,
333         // carrying shifted-out bits across words
334         int[] val = value;
335         for (int i = offset, j = bstart; i < len + offset;) {
336             int bv = bval[j++];
337             long hb = ((bv >>> 1) + carry) & LONG_MASK;
338             long v = val[i++] & LONG_MASK;
339             if (v != hb)
340                 return v < hb ? -1 : 1;
341             carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
342         }
343         return carry == 0 ? 0 : -1;
344     }
345 
346     /**
347      * Return the index of the lowest set bit in this MutableBigInteger. If the
348      * magnitude of this MutableBigInteger is zero, -1 is returned.
349      */
350     private final int getLowestSetBit() {
351         if (intLen == 0)
352             return -1;
353         int j, b;
354         for (j=intLen-1; (j > 0) && (value[j+offset] == 0); j--)
355             ;
356         b = value[j+offset];
357         if (b == 0)
358             return -1;
359         return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
360     }
361 
362     /**
363      * Return the int in use in this MutableBigInteger at the specified
364      * index. This method is not used because it is not inlined on all
365      * platforms.
366      */
367     private final int getInt(int index) {
368         return value[offset+index];
369     }
370 
371     /**
372      * Return a long which is equal to the unsigned value of the int in
373      * use in this MutableBigInteger at the specified index. This method is
374      * not used because it is not inlined on all platforms.
375      */
376     private final long getLong(int index) {
377         return value[offset+index] & LONG_MASK;
378     }
379 
380     /**
381      * Ensure that the MutableBigInteger is in normal form, specifically
382      * making sure that there are no leading zeros, and that if the
383      * magnitude is zero, then intLen is zero.
384      */
385     final void normalize() {
386         if (intLen == 0) {
387             offset = 0;
388             return;
389         }
390 
391         int index = offset;
392         if (value[index] != 0)
393             return;
394 
395         int indexBound = index+intLen;
396         do {
397             index++;
398         } while(index < indexBound && value[index] == 0);
399 
400         int numZeros = index - offset;
401         intLen -= numZeros;
402         offset = (intLen == 0 ?  0 : offset+numZeros);
403     }
404 
405     /**
406      * If this MutableBigInteger cannot hold len words, increase the size
407      * of the value array to len words.
408      */
409     private final void ensureCapacity(int len) {
410         if (value.length < len) {
411             value = new int[len];
412             offset = 0;
413             intLen = len;
414         }
415     }
416 
417     /**
418      * Convert this MutableBigInteger into an int array with no leading
419      * zeros, of a length that is equal to this MutableBigInteger's intLen.
420      */
421     int[] toIntArray() {
422         int[] result = new int[intLen];
423         for(int i=0; i < intLen; i++)
424             result[i] = value[offset+i];
425         return result;
426     }
427 
428     /**
429      * Sets the int at index+offset in this MutableBigInteger to val.
430      * This does not get inlined on all platforms so it is not used
431      * as often as originally intended.
432      */
433     void setInt(int index, int val) {
434         value[offset + index] = val;
435     }
436 
437     /**
438      * Sets this MutableBigInteger's value array to the specified array.
439      * The intLen is set to the specified length.
440      */
441     void setValue(int[] val, int length) {
442         value = val;
443         intLen = length;
444         offset = 0;
445     }
446 
447     /**
448      * Sets this MutableBigInteger's value array to a copy of the specified
449      * array. The intLen is set to the length of the new array.
450      */
451     void copyValue(MutableBigInteger src) {
452         int len = src.intLen;
453         if (value.length < len)
454             value = new int[len];
455         System.arraycopy(src.value, src.offset, value, 0, len);
456         intLen = len;
457         offset = 0;
458     }
459 
460     /**
461      * Sets this MutableBigInteger's value array to a copy of the specified
462      * array. The intLen is set to the length of the specified array.
463      */
464     void copyValue(int[] val) {
465         int len = val.length;
466         if (value.length < len)
467             value = new int[len];
468         System.arraycopy(val, 0, value, 0, len);
469         intLen = len;
470         offset = 0;
471     }
472 
473     /**
474      * Returns true iff this MutableBigInteger has a value of one.
475      */
476     boolean isOne() {
477         return (intLen == 1) && (value[offset] == 1);
478     }
479 
480     /**
481      * Returns true iff this MutableBigInteger has a value of zero.
482      */
483     boolean isZero() {
484         return (intLen == 0);
485     }
486 
487     /**
488      * Returns true iff this MutableBigInteger is even.
489      */
490     boolean isEven() {
491         return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
492     }
493 
494     /**
495      * Returns true iff this MutableBigInteger is odd.
496      */
497     boolean isOdd() {
498         return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
499     }
500 
501     /**
502      * Returns true iff this MutableBigInteger is in normal form. A
503      * MutableBigInteger is in normal form if it has no leading zeros
504      * after the offset, and intLen + offset <= value.length.
505      */
506     boolean isNormal() {
507         if (intLen + offset > value.length)
508             return false;
509         if (intLen == 0)
510             return true;
511         return (value[offset] != 0);
512     }
513 
514     /**
515      * Returns a String representation of this MutableBigInteger in radix 10.
516      */
517     public String toString() {
518         BigInteger b = toBigInteger(1);
519         return b.toString();
520     }
521 
522     /**
523      * Like {@link #rightShift(int)} but {@code n} can be greater than the length of the number.
524      */
525     void safeRightShift(int n) {
526         if (n/32 >= intLen) {
527             reset();
528         } else {
529             rightShift(n);
530         }
531     }
532 
533     /**
534      * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
535      * in normal form.
536      */
537     void rightShift(int n) {
538         if (intLen == 0)
539             return;
540         int nInts = n >>> 5;
541         int nBits = n & 0x1F;
542         this.intLen -= nInts;
543         if (nBits == 0)
544             return;
545         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
546         if (nBits >= bitsInHighWord) {
547             this.primitiveLeftShift(32 - nBits);
548             this.intLen--;
549         } else {
550             primitiveRightShift(nBits);
551         }
552     }
553 
554     /**
555      * Like {@link #leftShift(int)} but {@code n} can be zero.
556      */
557     void safeLeftShift(int n) {
558         if (n > 0) {
559             leftShift(n);
560         }
561     }
562 
563     /**
564      * Left shift this MutableBigInteger n bits.
565      */
566     void leftShift(int n) {
567         /*
568          * If there is enough storage space in this MutableBigInteger already
569          * the available space will be used. Space to the right of the used
570          * ints in the value array is faster to utilize, so the extra space
571          * will be taken from the right if possible.
572          */
573         if (intLen == 0)
574            return;
575         int nInts = n >>> 5;
576         int nBits = n&0x1F;
577         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
578 
579         // If shift can be done without moving words, do so
580         if (n <= (32-bitsInHighWord)) {
581             primitiveLeftShift(nBits);
582             return;
583         }
584 
585         int newLen = intLen + nInts +1;
586         if (nBits <= (32-bitsInHighWord))
587             newLen--;
588         if (value.length < newLen) {
589             // The array must grow
590             int[] result = new int[newLen];
591             for (int i=0; i < intLen; i++)
592                 result[i] = value[offset+i];
593             setValue(result, newLen);
594         } else if (value.length - offset >= newLen) {
595             // Use space on right
596             for(int i=0; i < newLen - intLen; i++)
597                 value[offset+intLen+i] = 0;
598         } else {
599             // Must use space on left
600             for (int i=0; i < intLen; i++)
601                 value[i] = value[offset+i];
602             for (int i=intLen; i < newLen; i++)
603                 value[i] = 0;
604             offset = 0;
605         }
606         intLen = newLen;
607         if (nBits == 0)
608             return;
609         if (nBits <= (32-bitsInHighWord))
610             primitiveLeftShift(nBits);
611         else
612             primitiveRightShift(32 -nBits);
613     }
614 
615     /**
616      * A primitive used for division. This method adds in one multiple of the
617      * divisor a back to the dividend result at a specified offset. It is used
618      * when qhat was estimated too large, and must be adjusted.
619      */
620     private int divadd(int[] a, int[] result, int offset) {
621         long carry = 0;
622 
623         for (int j=a.length-1; j >= 0; j--) {
624             long sum = (a[j] & LONG_MASK) +
625                        (result[j+offset] & LONG_MASK) + carry;
626             result[j+offset] = (int)sum;
627             carry = sum >>> 32;
628         }
629         return (int)carry;
630     }
631 
632     /**
633      * This method is used for division. It multiplies an n word input a by one
634      * word input x, and subtracts the n word product from q. This is needed
635      * when subtracting qhat*divisor from dividend.
636      */
637     private int mulsub(int[] q, int[] a, int x, int len, int offset) {
638         long xLong = x & LONG_MASK;
639         long carry = 0;
640         offset += len;
641 
642         for (int j=len-1; j >= 0; j--) {
643             long product = (a[j] & LONG_MASK) * xLong + carry;
644             long difference = q[offset] - product;
645             q[offset--] = (int)difference;
646             carry = (product >>> 32)
647                      + (((difference & LONG_MASK) >
648                          (((~(int)product) & LONG_MASK))) ? 1:0);
649         }
650         return (int)carry;
651     }
652 
653     /**
654      * The method is the same as mulsun, except the fact that q array is not
655      * updated, the only result of the method is borrow flag.
656      */
657     private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
658         long xLong = x & LONG_MASK;
659         long carry = 0;
660         offset += len;
661         for (int j=len-1; j >= 0; j--) {
662             long product = (a[j] & LONG_MASK) * xLong + carry;
663             long difference = q[offset--] - product;
664             carry = (product >>> 32)
665                      + (((difference & LONG_MASK) >
666                          (((~(int)product) & LONG_MASK))) ? 1:0);
667         }
668         return (int)carry;
669     }
670 
671     /**
672      * Right shift this MutableBigInteger n bits, where n is
673      * less than 32.
674      * Assumes that intLen > 0, n > 0 for speed
675      */
676     private final void primitiveRightShift(int n) {
677         int[] val = value;
678         int n2 = 32 - n;
679         for (int i=offset+intLen-1, c=val[i]; i > offset; i--) {
680             int b = c;
681             c = val[i-1];
682             val[i] = (c << n2) | (b >>> n);
683         }
684         val[offset] >>>= n;
685     }
686 
687     /**
688      * Left shift this MutableBigInteger n bits, where n is
689      * less than 32.
690      * Assumes that intLen > 0, n > 0 for speed
691      */
692     private final void primitiveLeftShift(int n) {
693         int[] val = value;
694         int n2 = 32 - n;
695         for (int i=offset, c=val[i], m=i+intLen-1; i < m; i++) {
696             int b = c;
697             c = val[i+1];
698             val[i] = (b << n) | (c >>> n2);
699         }
700         val[offset+intLen-1] <<= n;
701     }
702 
703     /**
704      * Returns a {@code BigInteger} equal to the {@code n}
705      * low ints of this number.
706      */
707     private BigInteger getLower(int n) {
708         if (isZero()) {
709             return BigInteger.ZERO;
710         } else if (intLen < n) {
711             return toBigInteger(1);
712         } else {
713             // strip zeros
714             int len = n;
715             while (len > 0 && value[offset+intLen-len] == 0)
716                 len--;
717             int sign = len > 0 ? 1 : 0;
718             return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
719         }
720     }
721 
722     /**
723      * Discards all ints whose index is greater than {@code n}.
724      */
725     private void keepLower(int n) {
726         if (intLen >= n) {
727             offset += intLen - n;
728             intLen = n;
729         }
730     }
731 
732     /**
733      * Adds the contents of two MutableBigInteger objects.The result
734      * is placed within this MutableBigInteger.
735      * The contents of the addend are not changed.
736      */
737     void add(MutableBigInteger addend) {
738         int x = intLen;
739         int y = addend.intLen;
740         int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
741         int[] result = (value.length < resultLen ? new int[resultLen] : value);
742 
743         int rstart = result.length-1;
744         long sum;
745         long carry = 0;
746 
747         // Add common parts of both numbers
748         while(x > 0 && y > 0) {
749             x--; y--;
750             sum = (value[x+offset] & LONG_MASK) +
751                 (addend.value[y+addend.offset] & LONG_MASK) + carry;
752             result[rstart--] = (int)sum;
753             carry = sum >>> 32;
754         }
755 
756         // Add remainder of the longer number
757         while(x > 0) {
758             x--;
759             if (carry == 0 && result == value && rstart == (x + offset))
760                 return;
761             sum = (value[x+offset] & LONG_MASK) + carry;
762             result[rstart--] = (int)sum;
763             carry = sum >>> 32;
764         }
765         while(y > 0) {
766             y--;
767             sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
768             result[rstart--] = (int)sum;
769             carry = sum >>> 32;
770         }
771 
772         if (carry > 0) { // Result must grow in length
773             resultLen++;
774             if (result.length < resultLen) {
775                 int temp[] = new int[resultLen];
776                 // Result one word longer from carry-out; copy low-order
777                 // bits into new result.
778                 System.arraycopy(result, 0, temp, 1, result.length);
779                 temp[0] = 1;
780                 result = temp;
781             } else {
782                 result[rstart--] = 1;
783             }
784         }
785 
786         value = result;
787         intLen = resultLen;
788         offset = result.length - resultLen;
789     }
790 
791     /**
792      * Adds the value of {@code addend} shifted {@code n} ints to the left.
793      * Has the same effect as {@code addend.leftShift(32*ints); add(addend);}
794      * but doesn't change the value of {@code addend}.
795      */
796     void addShifted(MutableBigInteger addend, int n) {
797         if (addend.isZero()) {
798             return;
799         }
800 
801         int x = intLen;
802         int y = addend.intLen + n;
803         int resultLen = (intLen > y ? intLen : y);
804         int[] result = (value.length < resultLen ? new int[resultLen] : value);
805 
806         int rstart = result.length-1;
807         long sum;
808         long carry = 0;
809 
810         // Add common parts of both numbers
811         while (x > 0 && y > 0) {
812             x--; y--;
813             int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
814             sum = (value[x+offset] & LONG_MASK) +
815                 (bval & LONG_MASK) + carry;
816             result[rstart--] = (int)sum;
817             carry = sum >>> 32;
818         }
819 
820         // Add remainder of the longer number
821         while (x > 0) {
822             x--;
823             if (carry == 0 && result == value && rstart == (x + offset)) {
824                 return;
825             }
826             sum = (value[x+offset] & LONG_MASK) + carry;
827             result[rstart--] = (int)sum;
828             carry = sum >>> 32;
829         }
830         while (y > 0) {
831             y--;
832             int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
833             sum = (bval & LONG_MASK) + carry;
834             result[rstart--] = (int)sum;
835             carry = sum >>> 32;
836         }
837 
838         if (carry > 0) { // Result must grow in length
839             resultLen++;
840             if (result.length < resultLen) {
841                 int temp[] = new int[resultLen];
842                 // Result one word longer from carry-out; copy low-order
843                 // bits into new result.
844                 System.arraycopy(result, 0, temp, 1, result.length);
845                 temp[0] = 1;
846                 result = temp;
847             } else {
848                 result[rstart--] = 1;
849             }
850         }
851 
852         value = result;
853         intLen = resultLen;
854         offset = result.length - resultLen;
855     }
856 
857     /**
858      * Like {@link #addShifted(MutableBigInteger, int)} but {@code this.intLen} must
859      * not be greater than {@code n}. In other words, concatenates {@code this}
860      * and {@code addend}.
861      */
862     void addDisjoint(MutableBigInteger addend, int n) {
863         if (addend.isZero())
864             return;
865 
866         int x = intLen;
867         int y = addend.intLen + n;
868         int resultLen = (intLen > y ? intLen : y);
869         int[] result;
870         if (value.length < resultLen)
871             result = new int[resultLen];
872         else {
873             result = value;
874             Arrays.fill(value, offset+intLen, value.length, 0);
875         }
876 
877         int rstart = result.length-1;
878 
879         // copy from this if needed
880         System.arraycopy(value, offset, result, rstart+1-x, x);
881         y -= x;
882         rstart -= x;
883 
884         int len = Math.min(y, addend.value.length-addend.offset);
885         System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len);
886 
887         // zero the gap
888         for (int i=rstart+1-y+len; i < rstart+1; i++)
889             result[i] = 0;
890 
891         value = result;
892         intLen = resultLen;
893         offset = result.length - resultLen;
894     }
895 
896     /**
897      * Adds the low {@code n} ints of {@code addend}.
898      */
899     void addLower(MutableBigInteger addend, int n) {
900         MutableBigInteger a = new MutableBigInteger(addend);
901         if (a.offset + a.intLen >= n) {
902             a.offset = a.offset + a.intLen - n;
903             a.intLen = n;
904         }
905         a.normalize();
906         add(a);
907     }
908 
909     /**
910      * Subtracts the smaller of this and b from the larger and places the
911      * result into this MutableBigInteger.
912      */
913     int subtract(MutableBigInteger b) {
914         MutableBigInteger a = this;
915 
916         int[] result = value;
917         int sign = a.compare(b);
918 
919         if (sign == 0) {
920             reset();
921             return 0;
922         }
923         if (sign < 0) {
924             MutableBigInteger tmp = a;
925             a = b;
926             b = tmp;
927         }
928 
929         int resultLen = a.intLen;
930         if (result.length < resultLen)
931             result = new int[resultLen];
932 
933         long diff = 0;
934         int x = a.intLen;
935         int y = b.intLen;
936         int rstart = result.length - 1;
937 
938         // Subtract common parts of both numbers
939         while (y > 0) {
940             x--; y--;
941 
942             diff = (a.value[x+a.offset] & LONG_MASK) -
943                    (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
944             result[rstart--] = (int)diff;
945         }
946         // Subtract remainder of longer number
947         while (x > 0) {
948             x--;
949             diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
950             result[rstart--] = (int)diff;
951         }
952 
953         value = result;
954         intLen = resultLen;
955         offset = value.length - resultLen;
956         normalize();
957         return sign;
958     }
959 
960     /**
961      * Subtracts the smaller of a and b from the larger and places the result
962      * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
963      * operation was performed.
964      */
965     private int difference(MutableBigInteger b) {
966         MutableBigInteger a = this;
967         int sign = a.compare(b);
968         if (sign == 0)
969             return 0;
970         if (sign < 0) {
971             MutableBigInteger tmp = a;
972             a = b;
973             b = tmp;
974         }
975 
976         long diff = 0;
977         int x = a.intLen;
978         int y = b.intLen;
979 
980         // Subtract common parts of both numbers
981         while (y > 0) {
982             x--; y--;
983             diff = (a.value[a.offset+ x] & LONG_MASK) -
984                 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
985             a.value[a.offset+x] = (int)diff;
986         }
987         // Subtract remainder of longer number
988         while (x > 0) {
989             x--;
990             diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
991             a.value[a.offset+x] = (int)diff;
992         }
993 
994         a.normalize();
995         return sign;
996     }
997 
998     /**
999      * Multiply the contents of two MutableBigInteger objects. The result is
1000      * placed into MutableBigInteger z. The contents of y are not changed.
1001      */
1002     void multiply(MutableBigInteger y, MutableBigInteger z) {
1003         int xLen = intLen;
1004         int yLen = y.intLen;
1005         int newLen = xLen + yLen;
1006 
1007         // Put z into an appropriate state to receive product
1008         if (z.value.length < newLen)
1009             z.value = new int[newLen];
1010         z.offset = 0;
1011         z.intLen = newLen;
1012 
1013         // The first iteration is hoisted out of the loop to avoid extra add
1014         long carry = 0;
1015         for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
1016                 long product = (y.value[j+y.offset] & LONG_MASK) *
1017                                (value[xLen-1+offset] & LONG_MASK) + carry;
1018                 z.value[k] = (int)product;
1019                 carry = product >>> 32;
1020         }
1021         z.value[xLen-1] = (int)carry;
1022 
1023         // Perform the multiplication word by word
1024         for (int i = xLen-2; i >= 0; i--) {
1025             carry = 0;
1026             for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
1027                 long product = (y.value[j+y.offset] & LONG_MASK) *
1028                                (value[i+offset] & LONG_MASK) +
1029                                (z.value[k] & LONG_MASK) + carry;
1030                 z.value[k] = (int)product;
1031                 carry = product >>> 32;
1032             }
1033             z.value[i] = (int)carry;
1034         }
1035 
1036         // Remove leading zeros from product
1037         z.normalize();
1038     }
1039 
1040     /**
1041      * Multiply the contents of this MutableBigInteger by the word y. The
1042      * result is placed into z.
1043      */
1044     void mul(int y, MutableBigInteger z) {
1045         if (y == 1) {
1046             z.copyValue(this);
1047             return;
1048         }
1049 
1050         if (y == 0) {
1051             z.clear();
1052             return;
1053         }
1054 
1055         // Perform the multiplication word by word
1056         long ylong = y & LONG_MASK;
1057         int[] zval = (z.value.length < intLen+1 ? new int[intLen + 1]
1058                                               : z.value);
1059         long carry = 0;
1060         for (int i = intLen-1; i >= 0; i--) {
1061             long product = ylong * (value[i+offset] & LONG_MASK) + carry;
1062             zval[i+1] = (int)product;
1063             carry = product >>> 32;
1064         }
1065 
1066         if (carry == 0) {
1067             z.offset = 1;
1068             z.intLen = intLen;
1069         } else {
1070             z.offset = 0;
1071             z.intLen = intLen + 1;
1072             zval[0] = (int)carry;
1073         }
1074         z.value = zval;
1075     }
1076 
1077      /**
1078      * This method is used for division of an n word dividend by a one word
1079      * divisor. The quotient is placed into quotient. The one word divisor is
1080      * specified by divisor.
1081      *
1082      * @return the remainder of the division is returned.
1083      *
1084      */
1085     int divideOneWord(int divisor, MutableBigInteger quotient) {
1086         long divisorLong = divisor & LONG_MASK;
1087 
1088         // Special case of one word dividend
1089         if (intLen == 1) {
1090             long dividendValue = value[offset] & LONG_MASK;
1091             int q = (int) (dividendValue / divisorLong);
1092             int r = (int) (dividendValue - q * divisorLong);
1093             quotient.value[0] = q;
1094             quotient.intLen = (q == 0) ? 0 : 1;
1095             quotient.offset = 0;
1096             return r;
1097         }
1098 
1099         if (quotient.value.length < intLen)
1100             quotient.value = new int[intLen];
1101         quotient.offset = 0;
1102         quotient.intLen = intLen;
1103 
1104         // Normalize the divisor
1105         int shift = Integer.numberOfLeadingZeros(divisor);
1106 
1107         int rem = value[offset];
1108         long remLong = rem & LONG_MASK;
1109         if (remLong < divisorLong) {
1110             quotient.value[0] = 0;
1111         } else {
1112             quotient.value[0] = (int)(remLong / divisorLong);
1113             rem = (int) (remLong - (quotient.value[0] * divisorLong));
1114             remLong = rem & LONG_MASK;
1115         }
1116         int xlen = intLen;
1117         while (--xlen > 0) {
1118             long dividendEstimate = (remLong << 32) |
1119                     (value[offset + intLen - xlen] & LONG_MASK);
1120             int q;
1121             if (dividendEstimate >= 0) {
1122                 q = (int) (dividendEstimate / divisorLong);
1123                 rem = (int) (dividendEstimate - q * divisorLong);
1124             } else {
1125                 long tmp = divWord(dividendEstimate, divisor);
1126                 q = (int) (tmp & LONG_MASK);
1127                 rem = (int) (tmp >>> 32);
1128             }
1129             quotient.value[intLen - xlen] = q;
1130             remLong = rem & LONG_MASK;
1131         }
1132 
1133         quotient.normalize();
1134         // Unnormalize
1135         if (shift > 0)
1136             return rem % divisor;
1137         else
1138             return rem;
1139     }
1140 
1141     /**
1142      * Calculates the quotient of this div b and places the quotient in the
1143      * provided MutableBigInteger objects and the remainder object is returned.
1144      *
1145      */
1146     MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
1147         return divide(b,quotient,true);
1148     }
1149 
1150     MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1151         if (b.intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD ||
1152                 intLen - b.intLen < BigInteger.BURNIKEL_ZIEGLER_OFFSET) {
1153             return divideKnuth(b, quotient, needRemainder);
1154         } else {
1155             return divideAndRemainderBurnikelZiegler(b, quotient);
1156         }
1157     }
1158 
1159     /**
1160      * @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
1161      */
1162     MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) {
1163         return divideKnuth(b,quotient,true);
1164     }
1165 
1166     /**
1167      * Calculates the quotient of this div b and places the quotient in the
1168      * provided MutableBigInteger objects and the remainder object is returned.
1169      *
1170      * Uses Algorithm D in Knuth section 4.3.1.
1171      * Many optimizations to that algorithm have been adapted from the Colin
1172      * Plumb C library.
1173      * It special cases one word divisors for speed. The content of b is not
1174      * changed.
1175      *
1176      */
1177     MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
1178         if (b.intLen == 0)
1179             throw new ArithmeticException("BigInteger divide by zero");
1180 
1181         // Dividend is zero
1182         if (intLen == 0) {
1183             quotient.intLen = quotient.offset = 0;
1184             return needRemainder ? new MutableBigInteger() : null;
1185         }
1186 
1187         int cmp = compare(b);
1188         // Dividend less than divisor
1189         if (cmp < 0) {
1190             quotient.intLen = quotient.offset = 0;
1191             return needRemainder ? new MutableBigInteger(this) : null;
1192         }
1193         // Dividend equal to divisor
1194         if (cmp == 0) {
1195             quotient.value[0] = quotient.intLen = 1;
1196             quotient.offset = 0;
1197             return needRemainder ? new MutableBigInteger() : null;
1198         }
1199 
1200         quotient.clear();
1201         // Special case one word divisor
1202         if (b.intLen == 1) {
1203             int r = divideOneWord(b.value[b.offset], quotient);
1204             if(needRemainder) {
1205                 if (r == 0)
1206                     return new MutableBigInteger();
1207                 return new MutableBigInteger(r);
1208             } else {
1209                 return null;
1210             }
1211         }
1212 
1213         // Cancel common powers of two if we're above the KNUTH_POW2_* thresholds
1214         if (intLen >= KNUTH_POW2_THRESH_LEN) {
1215             int trailingZeroBits = Math.min(getLowestSetBit(), b.getLowestSetBit());
1216             if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS*32) {
1217                 MutableBigInteger a = new MutableBigInteger(this);
1218                 b = new MutableBigInteger(b);
1219                 a.rightShift(trailingZeroBits);
1220                 b.rightShift(trailingZeroBits);
1221                 MutableBigInteger r = a.divideKnuth(b, quotient);
1222                 r.leftShift(trailingZeroBits);
1223                 return r;
1224             }
1225         }
1226 
1227         return divideMagnitude(b, quotient, needRemainder);
1228     }
1229 
1230     /**
1231      * Computes {@code this/b} and {@code this%b} using the
1232      * <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
1233      * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
1234      * The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are
1235      * multiples of 32 bits.<br/>
1236      * {@code this} and {@code b} must be nonnegative.
1237      * @param b the divisor
1238      * @param quotient output parameter for {@code this/b}
1239      * @return the remainder
1240      */
1241     MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) {
1242         int r = intLen;
1243         int s = b.intLen;
1244 
1245         // Clear the quotient
1246         quotient.offset = quotient.intLen = 0;
1247 
1248         if (r < s) {
1249             return this;
1250         } else {
1251             // Unlike Knuth division, we don't check for common powers of two here because
1252             // BZ already runs faster if both numbers contain powers of two and cancelling them has no
1253             // additional benefit.
1254 
1255             // step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
1256             int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
1257 
1258             int j = (s+m-1) / m;      // step 2a: j = ceil(s/m)
1259             int n = j * m;            // step 2b: block length in 32-bit units
1260             long n32 = 32L * n;         // block length in bits
1261             int sigma = (int) Math.max(0, n32 - b.bitLength());   // step 3: sigma = max{T | (2^T)*B < beta^n}
1262             MutableBigInteger bShifted = new MutableBigInteger(b);
1263             bShifted.safeLeftShift(sigma);   // step 4a: shift b so its length is a multiple of n
1264             safeLeftShift(sigma);     // step 4b: shift this by the same amount
1265 
1266             // step 5: t is the number of blocks needed to accommodate this plus one additional bit
1267             int t = (int) ((bitLength()+n32) / n32);
1268             if (t < 2) {
1269                 t = 2;
1270             }
1271 
1272             // step 6: conceptually split this into blocks a[t-1], ..., a[0]
1273             MutableBigInteger a1 = getBlock(t-1, t, n);   // the most significant block of this
1274 
1275             // step 7: z[t-2] = [a[t-1], a[t-2]]
1276             MutableBigInteger z = 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 = 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: this 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 GCD of this and b. This and b are changed by the computation.
1870      */
1871     MutableBigInteger hybridGCD(MutableBigInteger b) {
1872         // Use Euclid's algorithm until the numbers are approximately the
1873         // same length, then use the binary GCD algorithm to find the GCD.
1874         MutableBigInteger a = this;
1875         MutableBigInteger q = new MutableBigInteger();
1876 
1877         while (b.intLen != 0) {
1878             if (Math.abs(a.intLen - b.intLen) < 2)
1879                 return a.binaryGCD(b);
1880 
1881             MutableBigInteger r = a.divide(b, q);
1882             a = b;
1883             b = r;
1884         }
1885         return a;
1886     }
1887 
1888     /**
1889      * Calculate GCD of this and v.
1890      * Assumes that this and v are not zero.
1891      */
1892     private MutableBigInteger binaryGCD(MutableBigInteger v) {
1893         // Algorithm B from Knuth section 4.5.2
1894         MutableBigInteger u = this;
1895         MutableBigInteger r = new MutableBigInteger();
1896 
1897         // step B1
1898         int s1 = u.getLowestSetBit();
1899         int s2 = v.getLowestSetBit();
1900         int k = (s1 < s2) ? s1 : s2;
1901         if (k != 0) {
1902             u.rightShift(k);
1903             v.rightShift(k);
1904         }
1905 
1906         // step B2
1907         boolean uOdd = (k == s1);
1908         MutableBigInteger t = uOdd ? v: u;
1909         int tsign = uOdd ? -1 : 1;
1910 
1911         int lb;
1912         while ((lb = t.getLowestSetBit()) >= 0) {
1913             // steps B3 and B4
1914             t.rightShift(lb);
1915             // step B5
1916             if (tsign > 0)
1917                 u = t;
1918             else
1919                 v = t;
1920 
1921             // Special case one word numbers
1922             if (u.intLen < 2 && v.intLen < 2) {
1923                 int x = u.value[u.offset];
1924                 int y = v.value[v.offset];
1925                 x  = binaryGcd(x, y);
1926                 r.value[0] = x;
1927                 r.intLen = 1;
1928                 r.offset = 0;
1929                 if (k > 0)
1930                     r.leftShift(k);
1931                 return r;
1932             }
1933 
1934             // step B6
1935             if ((tsign = u.difference(v)) == 0)
1936                 break;
1937             t = (tsign >= 0) ? u : v;
1938         }
1939 
1940         if (k > 0)
1941             u.leftShift(k);
1942         return u;
1943     }
1944 
1945     /**
1946      * Calculate GCD of a and b interpreted as unsigned integers.
1947      */
1948     static int binaryGcd(int a, int b) {
1949         if (b == 0)
1950             return a;
1951         if (a == 0)
1952             return b;
1953 
1954         // Right shift a & b till their last bits equal to 1.
1955         int aZeros = Integer.numberOfTrailingZeros(a);
1956         int bZeros = Integer.numberOfTrailingZeros(b);
1957         a >>>= aZeros;
1958         b >>>= bZeros;
1959 
1960         int t = (aZeros < bZeros ? aZeros : bZeros);
1961 
1962         while (a != b) {
1963             if ((a+0x80000000) > (b+0x80000000)) {  // a > b as unsigned
1964                 a -= b;
1965                 a >>>= Integer.numberOfTrailingZeros(a);
1966             } else {
1967                 b -= a;
1968                 b >>>= Integer.numberOfTrailingZeros(b);
1969             }
1970         }
1971         return a<<t;
1972     }
1973 
1974     /**
1975      * Returns the modInverse of this mod p. This and p are not affected by
1976      * the operation.
1977      */
1978     MutableBigInteger mutableModInverse(MutableBigInteger p) {
1979         // Modulus is odd, use Schroeppel's algorithm
1980         if (p.isOdd())
1981             return modInverse(p);
1982 
1983         // Base and modulus are even, throw exception
1984         if (isEven())
1985             throw new ArithmeticException("BigInteger not invertible.");
1986 
1987         // Get even part of modulus expressed as a power of 2
1988         int powersOf2 = p.getLowestSetBit();
1989 
1990         // Construct odd part of modulus
1991         MutableBigInteger oddMod = new MutableBigInteger(p);
1992         oddMod.rightShift(powersOf2);
1993 
1994         if (oddMod.isOne())
1995             return modInverseMP2(powersOf2);
1996 
1997         // Calculate 1/a mod oddMod
1998         MutableBigInteger oddPart = modInverse(oddMod);
1999 
2000         // Calculate 1/a mod evenMod
2001         MutableBigInteger evenPart = modInverseMP2(powersOf2);
2002 
2003         // Combine the results using Chinese Remainder Theorem
2004         MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
2005         MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
2006 
2007         MutableBigInteger temp1 = new MutableBigInteger();
2008         MutableBigInteger temp2 = new MutableBigInteger();
2009         MutableBigInteger result = new MutableBigInteger();
2010 
2011         oddPart.leftShift(powersOf2);
2012         oddPart.multiply(y1, result);
2013 
2014         evenPart.multiply(oddMod, temp1);
2015         temp1.multiply(y2, temp2);
2016 
2017         result.add(temp2);
2018         return result.divide(p, temp1);
2019     }
2020 
2021     /*
2022      * Calculate the multiplicative inverse of this mod 2^k.
2023      */
2024     MutableBigInteger modInverseMP2(int k) {
2025         if (isEven())
2026             throw new ArithmeticException("Non-invertible. (GCD != 1)");
2027 
2028         if (k > 64)
2029             return euclidModInverse(k);
2030 
2031         int t = inverseMod32(value[offset+intLen-1]);
2032 
2033         if (k < 33) {
2034             t = (k == 32 ? t : t & ((1 << k) - 1));
2035             return new MutableBigInteger(t);
2036         }
2037 
2038         long pLong = (value[offset+intLen-1] & LONG_MASK);
2039         if (intLen > 1)
2040             pLong |=  ((long)value[offset+intLen-2] << 32);
2041         long tLong = t & LONG_MASK;
2042         tLong = tLong * (2 - pLong * tLong);  // 1 more Newton iter step
2043         tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
2044 
2045         MutableBigInteger result = new MutableBigInteger(new int[2]);
2046         result.value[0] = (int)(tLong >>> 32);
2047         result.value[1] = (int)tLong;
2048         result.intLen = 2;
2049         result.normalize();
2050         return result;
2051     }
2052 
2053     /**
2054      * Returns the multiplicative inverse of val mod 2^32.  Assumes val is odd.
2055      */
2056     static int inverseMod32(int val) {
2057         // Newton's iteration!
2058         int t = val;
2059         t *= 2 - val*t;
2060         t *= 2 - val*t;
2061         t *= 2 - val*t;
2062         t *= 2 - val*t;
2063         return t;
2064     }
2065 
2066     /**
2067      * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
2068      */
2069     static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
2070         // Copy the mod to protect original
2071         return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
2072     }
2073 
2074     /**
2075      * Calculate the multiplicative inverse of this mod mod, where mod is odd.
2076      * This and mod are not changed by the calculation.
2077      *
2078      * This method implements an algorithm due to Richard Schroeppel, that uses
2079      * the same intermediate representation as Montgomery Reduction
2080      * ("Montgomery Form").  The algorithm is described in an unpublished
2081      * manuscript entitled "Fast Modular Reciprocals."
2082      */
2083     private MutableBigInteger modInverse(MutableBigInteger mod) {
2084         MutableBigInteger p = new MutableBigInteger(mod);
2085         MutableBigInteger f = new MutableBigInteger(this);
2086         MutableBigInteger g = new MutableBigInteger(p);
2087         SignedMutableBigInteger c = new SignedMutableBigInteger(1);
2088         SignedMutableBigInteger d = new SignedMutableBigInteger();
2089         MutableBigInteger temp = null;
2090         SignedMutableBigInteger sTemp = null;
2091 
2092         int k = 0;
2093         // Right shift f k times until odd, left shift d k times
2094         if (f.isEven()) {
2095             int trailingZeros = f.getLowestSetBit();
2096             f.rightShift(trailingZeros);
2097             d.leftShift(trailingZeros);
2098             k = trailingZeros;
2099         }
2100 
2101         // The Almost Inverse Algorithm
2102         while (!f.isOne()) {
2103             // If gcd(f, g) != 1, number is not invertible modulo mod
2104             if (f.isZero())
2105                 throw new ArithmeticException("BigInteger not invertible.");
2106 
2107             // If f < g exchange f, g and c, d
2108             if (f.compare(g) < 0) {
2109                 temp = f; f = g; g = temp;
2110                 sTemp = d; d = c; c = sTemp;
2111             }
2112 
2113             // If f == g (mod 4)
2114             if (((f.value[f.offset + f.intLen - 1] ^
2115                  g.value[g.offset + g.intLen - 1]) & 3) == 0) {
2116                 f.subtract(g);
2117                 c.signedSubtract(d);
2118             } else { // If f != g (mod 4)
2119                 f.add(g);
2120                 c.signedAdd(d);
2121             }
2122 
2123             // Right shift f k times until odd, left shift d k times
2124             int trailingZeros = f.getLowestSetBit();
2125             f.rightShift(trailingZeros);
2126             d.leftShift(trailingZeros);
2127             k += trailingZeros;
2128         }
2129 
2130         while (c.sign < 0)
2131            c.signedAdd(p);
2132 
2133         return fixup(c, p, k);
2134     }
2135 
2136     /**
2137      * The Fixup Algorithm
2138      * Calculates X such that X = C * 2^(-k) (mod P)
2139      * Assumes C<P and P is odd.
2140      */
2141     static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
2142                                                                       int k) {
2143         MutableBigInteger temp = new MutableBigInteger();
2144         // Set r to the multiplicative inverse of p mod 2^32
2145         int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
2146 
2147         for (int i=0, numWords = k >> 5; i < numWords; i++) {
2148             // V = R * c (mod 2^j)
2149             int  v = r * c.value[c.offset + c.intLen-1];
2150             // c = c + (v * p)
2151             p.mul(v, temp);
2152             c.add(temp);
2153             // c = c / 2^j
2154             c.intLen--;
2155         }
2156         int numBits = k & 0x1f;
2157         if (numBits != 0) {
2158             // V = R * c (mod 2^j)
2159             int v = r * c.value[c.offset + c.intLen-1];
2160             v &= ((1<<numBits) - 1);
2161             // c = c + (v * p)
2162             p.mul(v, temp);
2163             c.add(temp);
2164             // c = c / 2^j
2165             c.rightShift(numBits);
2166         }
2167 
2168         // In theory, c may be greater than p at this point (Very rare!)
2169         while (c.compare(p) >= 0)
2170             c.subtract(p);
2171 
2172         return c;
2173     }
2174 
2175     /**
2176      * Uses the extended Euclidean algorithm to compute the modInverse of base
2177      * mod a modulus that is a power of 2. The modulus is 2^k.
2178      */
2179     MutableBigInteger euclidModInverse(int k) {
2180         MutableBigInteger b = new MutableBigInteger(1);
2181         b.leftShift(k);
2182         MutableBigInteger mod = new MutableBigInteger(b);
2183 
2184         MutableBigInteger a = new MutableBigInteger(this);
2185         MutableBigInteger q = new MutableBigInteger();
2186         MutableBigInteger r = b.divide(a, q);
2187 
2188         MutableBigInteger swapper = b;
2189         // swap b & r
2190         b = r;
2191         r = swapper;
2192 
2193         MutableBigInteger t1 = new MutableBigInteger(q);
2194         MutableBigInteger t0 = new MutableBigInteger(1);
2195         MutableBigInteger temp = new MutableBigInteger();
2196 
2197         while (!b.isOne()) {
2198             r = a.divide(b, q);
2199 
2200             if (r.intLen == 0)
2201                 throw new ArithmeticException("BigInteger not invertible.");
2202 
2203             swapper = r;
2204             a = swapper;
2205 
2206             if (q.intLen == 1)
2207                 t1.mul(q.value[q.offset], temp);
2208             else
2209                 q.multiply(t1, temp);
2210             swapper = q;
2211             q = temp;
2212             temp = swapper;
2213             t0.add(q);
2214 
2215             if (a.isOne())
2216                 return t0;
2217 
2218             r = b.divide(a, q);
2219 
2220             if (r.intLen == 0)
2221                 throw new ArithmeticException("BigInteger not invertible.");
2222 
2223             swapper = b;
2224             b =  r;
2225 
2226             if (q.intLen == 1)
2227                 t0.mul(q.value[q.offset], temp);
2228             else
2229                 q.multiply(t0, temp);
2230             swapper = q; q = temp; temp = swapper;
2231 
2232             t1.add(q);
2233         }
2234         mod.subtract(t1);
2235         return mod;
2236     }
2237 }