jaroslav@1258: /*
jaroslav@1258:  * Copyright (c) 1999, 2007, Oracle and/or its affiliates. All rights reserved.
jaroslav@1258:  * DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
jaroslav@1258:  *
jaroslav@1258:  * This code is free software; you can redistribute it and/or modify it
jaroslav@1258:  * under the terms of the GNU General Public License version 2 only, as
jaroslav@1258:  * published by the Free Software Foundation.  Oracle designates this
jaroslav@1258:  * particular file as subject to the "Classpath" exception as provided
jaroslav@1258:  * by Oracle in the LICENSE file that accompanied this code.
jaroslav@1258:  *
jaroslav@1258:  * This code is distributed in the hope that it will be useful, but WITHOUT
jaroslav@1258:  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
jaroslav@1258:  * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
jaroslav@1258:  * version 2 for more details (a copy is included in the LICENSE file that
jaroslav@1258:  * accompanied this code).
jaroslav@1258:  *
jaroslav@1258:  * You should have received a copy of the GNU General Public License version
jaroslav@1258:  * 2 along with this work; if not, write to the Free Software Foundation,
jaroslav@1258:  * Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
jaroslav@1258:  *
jaroslav@1258:  * Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
jaroslav@1258:  * or visit www.oracle.com if you need additional information or have any
jaroslav@1258:  * questions.
jaroslav@1258:  */
jaroslav@1258: 
jaroslav@1258: package java.math;
jaroslav@1258: 
jaroslav@1258: /**
jaroslav@1258:  * A class used to represent multiprecision integers that makes efficient
jaroslav@1258:  * use of allocated space by allowing a number to occupy only part of
jaroslav@1258:  * an array so that the arrays do not have to be reallocated as often.
jaroslav@1258:  * When performing an operation with many iterations the array used to
jaroslav@1258:  * hold a number is only reallocated when necessary and does not have to
jaroslav@1258:  * be the same size as the number it represents. A mutable number allows
jaroslav@1258:  * calculations to occur on the same number without having to create
jaroslav@1258:  * a new number for every step of the calculation as occurs with
jaroslav@1258:  * BigIntegers.
jaroslav@1258:  *
jaroslav@1258:  * @see     BigInteger
jaroslav@1258:  * @author  Michael McCloskey
jaroslav@1258:  * @since   1.3
jaroslav@1258:  */
jaroslav@1258: 
jaroslav@1258: import java.util.Arrays;
jaroslav@1258: 
jaroslav@1258: import static java.math.BigInteger.LONG_MASK;
jaroslav@1258: import static java.math.BigDecimal.INFLATED;
jaroslav@1258: 
jaroslav@1258: class MutableBigInteger {
jaroslav@1258:     /**
jaroslav@1258:      * Holds the magnitude of this MutableBigInteger in big endian order.
jaroslav@1258:      * The magnitude may start at an offset into the value array, and it may
jaroslav@1258:      * end before the length of the value array.
jaroslav@1258:      */
jaroslav@1258:     int[] value;
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * The number of ints of the value array that are currently used
jaroslav@1258:      * to hold the magnitude of this MutableBigInteger. The magnitude starts
jaroslav@1258:      * at an offset and offset + intLen may be less than value.length.
jaroslav@1258:      */
jaroslav@1258:     int intLen;
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * The offset into the value array where the magnitude of this
jaroslav@1258:      * MutableBigInteger begins.
jaroslav@1258:      */
jaroslav@1258:     int offset = 0;
jaroslav@1258: 
jaroslav@1258:     // Constants
jaroslav@1258:     /**
jaroslav@1258:      * MutableBigInteger with one element value array with the value 1. Used by
jaroslav@1258:      * BigDecimal divideAndRound to increment the quotient. Use this constant
jaroslav@1258:      * only when the method is not going to modify this object.
jaroslav@1258:      */
jaroslav@1258:     static final MutableBigInteger ONE = new MutableBigInteger(1);
jaroslav@1258: 
jaroslav@1258:     // Constructors
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * The default constructor. An empty MutableBigInteger is created with
jaroslav@1258:      * a one word capacity.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger() {
jaroslav@1258:         value = new int[1];
jaroslav@1258:         intLen = 0;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Construct a new MutableBigInteger with a magnitude specified by
jaroslav@1258:      * the int val.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger(int val) {
jaroslav@1258:         value = new int[1];
jaroslav@1258:         intLen = 1;
jaroslav@1258:         value[0] = val;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Construct a new MutableBigInteger with the specified value array
jaroslav@1258:      * up to the length of the array supplied.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger(int[] val) {
jaroslav@1258:         value = val;
jaroslav@1258:         intLen = val.length;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Construct a new MutableBigInteger with a magnitude equal to the
jaroslav@1258:      * specified BigInteger.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger(BigInteger b) {
jaroslav@1258:         intLen = b.mag.length;
jaroslav@1258:         value = Arrays.copyOf(b.mag, intLen);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Construct a new MutableBigInteger with a magnitude equal to the
jaroslav@1258:      * specified MutableBigInteger.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger(MutableBigInteger val) {
jaroslav@1258:         intLen = val.intLen;
jaroslav@1258:         value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Internal helper method to return the magnitude array. The caller is not
jaroslav@1258:      * supposed to modify the returned array.
jaroslav@1258:      */
jaroslav@1258:     private int[] getMagnitudeArray() {
jaroslav@1258:         if (offset > 0 || value.length != intLen)
jaroslav@1258:             return Arrays.copyOfRange(value, offset, offset + intLen);
jaroslav@1258:         return value;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Convert this MutableBigInteger to a long value. The caller has to make
jaroslav@1258:      * sure this MutableBigInteger can be fit into long.
jaroslav@1258:      */
jaroslav@1258:     private long toLong() {
jaroslav@1258:         assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
jaroslav@1258:         if (intLen == 0)
jaroslav@1258:             return 0;
jaroslav@1258:         long d = value[offset] & LONG_MASK;
jaroslav@1258:         return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Convert this MutableBigInteger to a BigInteger object.
jaroslav@1258:      */
jaroslav@1258:     BigInteger toBigInteger(int sign) {
jaroslav@1258:         if (intLen == 0 || sign == 0)
jaroslav@1258:             return BigInteger.ZERO;
jaroslav@1258:         return new BigInteger(getMagnitudeArray(), sign);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Convert this MutableBigInteger to BigDecimal object with the specified sign
jaroslav@1258:      * and scale.
jaroslav@1258:      */
jaroslav@1258:     BigDecimal toBigDecimal(int sign, int scale) {
jaroslav@1258:         if (intLen == 0 || sign == 0)
jaroslav@1258:             return BigDecimal.valueOf(0, scale);
jaroslav@1258:         int[] mag = getMagnitudeArray();
jaroslav@1258:         int len = mag.length;
jaroslav@1258:         int d = mag[0];
jaroslav@1258:         // If this MutableBigInteger can't be fit into long, we need to
jaroslav@1258:         // make a BigInteger object for the resultant BigDecimal object.
jaroslav@1258:         if (len > 2 || (d < 0 && len == 2))
jaroslav@1258:             return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
jaroslav@1258:         long v = (len == 2) ?
jaroslav@1258:             ((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
jaroslav@1258:             d & LONG_MASK;
jaroslav@1258:         return new BigDecimal(null, sign == -1 ? -v : v, scale, 0);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Clear out a MutableBigInteger for reuse.
jaroslav@1258:      */
jaroslav@1258:     void clear() {
jaroslav@1258:         offset = intLen = 0;
jaroslav@1258:         for (int index=0, n=value.length; index < n; index++)
jaroslav@1258:             value[index] = 0;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Set a MutableBigInteger to zero, removing its offset.
jaroslav@1258:      */
jaroslav@1258:     void reset() {
jaroslav@1258:         offset = intLen = 0;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
jaroslav@1258:      * as this MutableBigInteger is numerically less than, equal to, or
jaroslav@1258:      * greater than <tt>b</tt>.
jaroslav@1258:      */
jaroslav@1258:     final int compare(MutableBigInteger b) {
jaroslav@1258:         int blen = b.intLen;
jaroslav@1258:         if (intLen < blen)
jaroslav@1258:             return -1;
jaroslav@1258:         if (intLen > blen)
jaroslav@1258:            return 1;
jaroslav@1258: 
jaroslav@1258:         // Add Integer.MIN_VALUE to make the comparison act as unsigned integer
jaroslav@1258:         // comparison.
jaroslav@1258:         int[] bval = b.value;
jaroslav@1258:         for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
jaroslav@1258:             int b1 = value[i] + 0x80000000;
jaroslav@1258:             int b2 = bval[j]  + 0x80000000;
jaroslav@1258:             if (b1 < b2)
jaroslav@1258:                 return -1;
jaroslav@1258:             if (b1 > b2)
jaroslav@1258:                 return 1;
jaroslav@1258:         }
jaroslav@1258:         return 0;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Compare this against half of a MutableBigInteger object (Needed for
jaroslav@1258:      * remainder tests).
jaroslav@1258:      * Assumes no leading unnecessary zeros, which holds for results
jaroslav@1258:      * from divide().
jaroslav@1258:      */
jaroslav@1258:     final int compareHalf(MutableBigInteger b) {
jaroslav@1258:         int blen = b.intLen;
jaroslav@1258:         int len = intLen;
jaroslav@1258:         if (len <= 0)
jaroslav@1258:             return blen <=0 ? 0 : -1;
jaroslav@1258:         if (len > blen)
jaroslav@1258:             return 1;
jaroslav@1258:         if (len < blen - 1)
jaroslav@1258:             return -1;
jaroslav@1258:         int[] bval = b.value;
jaroslav@1258:         int bstart = 0;
jaroslav@1258:         int carry = 0;
jaroslav@1258:         // Only 2 cases left:len == blen or len == blen - 1
jaroslav@1258:         if (len != blen) { // len == blen - 1
jaroslav@1258:             if (bval[bstart] == 1) {
jaroslav@1258:                 ++bstart;
jaroslav@1258:                 carry = 0x80000000;
jaroslav@1258:             } else
jaroslav@1258:                 return -1;
jaroslav@1258:         }
jaroslav@1258:         // compare values with right-shifted values of b,
jaroslav@1258:         // carrying shifted-out bits across words
jaroslav@1258:         int[] val = value;
jaroslav@1258:         for (int i = offset, j = bstart; i < len + offset;) {
jaroslav@1258:             int bv = bval[j++];
jaroslav@1258:             long hb = ((bv >>> 1) + carry) & LONG_MASK;
jaroslav@1258:             long v = val[i++] & LONG_MASK;
jaroslav@1258:             if (v != hb)
jaroslav@1258:                 return v < hb ? -1 : 1;
jaroslav@1258:             carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
jaroslav@1258:         }
jaroslav@1258:         return carry == 0? 0 : -1;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Return the index of the lowest set bit in this MutableBigInteger. If the
jaroslav@1258:      * magnitude of this MutableBigInteger is zero, -1 is returned.
jaroslav@1258:      */
jaroslav@1258:     private final int getLowestSetBit() {
jaroslav@1258:         if (intLen == 0)
jaroslav@1258:             return -1;
jaroslav@1258:         int j, b;
jaroslav@1258:         for (j=intLen-1; (j>0) && (value[j+offset]==0); j--)
jaroslav@1258:             ;
jaroslav@1258:         b = value[j+offset];
jaroslav@1258:         if (b==0)
jaroslav@1258:             return -1;
jaroslav@1258:         return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Return the int in use in this MutableBigInteger at the specified
jaroslav@1258:      * index. This method is not used because it is not inlined on all
jaroslav@1258:      * platforms.
jaroslav@1258:      */
jaroslav@1258:     private final int getInt(int index) {
jaroslav@1258:         return value[offset+index];
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Return a long which is equal to the unsigned value of the int in
jaroslav@1258:      * use in this MutableBigInteger at the specified index. This method is
jaroslav@1258:      * not used because it is not inlined on all platforms.
jaroslav@1258:      */
jaroslav@1258:     private final long getLong(int index) {
jaroslav@1258:         return value[offset+index] & LONG_MASK;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Ensure that the MutableBigInteger is in normal form, specifically
jaroslav@1258:      * making sure that there are no leading zeros, and that if the
jaroslav@1258:      * magnitude is zero, then intLen is zero.
jaroslav@1258:      */
jaroslav@1258:     final void normalize() {
jaroslav@1258:         if (intLen == 0) {
jaroslav@1258:             offset = 0;
jaroslav@1258:             return;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         int index = offset;
jaroslav@1258:         if (value[index] != 0)
jaroslav@1258:             return;
jaroslav@1258: 
jaroslav@1258:         int indexBound = index+intLen;
jaroslav@1258:         do {
jaroslav@1258:             index++;
jaroslav@1258:         } while(index < indexBound && value[index]==0);
jaroslav@1258: 
jaroslav@1258:         int numZeros = index - offset;
jaroslav@1258:         intLen -= numZeros;
jaroslav@1258:         offset = (intLen==0 ?  0 : offset+numZeros);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * If this MutableBigInteger cannot hold len words, increase the size
jaroslav@1258:      * of the value array to len words.
jaroslav@1258:      */
jaroslav@1258:     private final void ensureCapacity(int len) {
jaroslav@1258:         if (value.length < len) {
jaroslav@1258:             value = new int[len];
jaroslav@1258:             offset = 0;
jaroslav@1258:             intLen = len;
jaroslav@1258:         }
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Convert this MutableBigInteger into an int array with no leading
jaroslav@1258:      * zeros, of a length that is equal to this MutableBigInteger's intLen.
jaroslav@1258:      */
jaroslav@1258:     int[] toIntArray() {
jaroslav@1258:         int[] result = new int[intLen];
jaroslav@1258:         for(int i=0; i<intLen; i++)
jaroslav@1258:             result[i] = value[offset+i];
jaroslav@1258:         return result;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Sets the int at index+offset in this MutableBigInteger to val.
jaroslav@1258:      * This does not get inlined on all platforms so it is not used
jaroslav@1258:      * as often as originally intended.
jaroslav@1258:      */
jaroslav@1258:     void setInt(int index, int val) {
jaroslav@1258:         value[offset + index] = val;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Sets this MutableBigInteger's value array to the specified array.
jaroslav@1258:      * The intLen is set to the specified length.
jaroslav@1258:      */
jaroslav@1258:     void setValue(int[] val, int length) {
jaroslav@1258:         value = val;
jaroslav@1258:         intLen = length;
jaroslav@1258:         offset = 0;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Sets this MutableBigInteger's value array to a copy of the specified
jaroslav@1258:      * array. The intLen is set to the length of the new array.
jaroslav@1258:      */
jaroslav@1258:     void copyValue(MutableBigInteger src) {
jaroslav@1258:         int len = src.intLen;
jaroslav@1258:         if (value.length < len)
jaroslav@1258:             value = new int[len];
jaroslav@1258:         System.arraycopy(src.value, src.offset, value, 0, len);
jaroslav@1258:         intLen = len;
jaroslav@1258:         offset = 0;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Sets this MutableBigInteger's value array to a copy of the specified
jaroslav@1258:      * array. The intLen is set to the length of the specified array.
jaroslav@1258:      */
jaroslav@1258:     void copyValue(int[] val) {
jaroslav@1258:         int len = val.length;
jaroslav@1258:         if (value.length < len)
jaroslav@1258:             value = new int[len];
jaroslav@1258:         System.arraycopy(val, 0, value, 0, len);
jaroslav@1258:         intLen = len;
jaroslav@1258:         offset = 0;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Returns true iff this MutableBigInteger has a value of one.
jaroslav@1258:      */
jaroslav@1258:     boolean isOne() {
jaroslav@1258:         return (intLen == 1) && (value[offset] == 1);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Returns true iff this MutableBigInteger has a value of zero.
jaroslav@1258:      */
jaroslav@1258:     boolean isZero() {
jaroslav@1258:         return (intLen == 0);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Returns true iff this MutableBigInteger is even.
jaroslav@1258:      */
jaroslav@1258:     boolean isEven() {
jaroslav@1258:         return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Returns true iff this MutableBigInteger is odd.
jaroslav@1258:      */
jaroslav@1258:     boolean isOdd() {
jaroslav@1258:         return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Returns true iff this MutableBigInteger is in normal form. A
jaroslav@1258:      * MutableBigInteger is in normal form if it has no leading zeros
jaroslav@1258:      * after the offset, and intLen + offset <= value.length.
jaroslav@1258:      */
jaroslav@1258:     boolean isNormal() {
jaroslav@1258:         if (intLen + offset > value.length)
jaroslav@1258:             return false;
jaroslav@1258:         if (intLen ==0)
jaroslav@1258:             return true;
jaroslav@1258:         return (value[offset] != 0);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Returns a String representation of this MutableBigInteger in radix 10.
jaroslav@1258:      */
jaroslav@1258:     public String toString() {
jaroslav@1258:         BigInteger b = toBigInteger(1);
jaroslav@1258:         return b.toString();
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Right shift this MutableBigInteger n bits. The MutableBigInteger is left
jaroslav@1258:      * in normal form.
jaroslav@1258:      */
jaroslav@1258:     void rightShift(int n) {
jaroslav@1258:         if (intLen == 0)
jaroslav@1258:             return;
jaroslav@1258:         int nInts = n >>> 5;
jaroslav@1258:         int nBits = n & 0x1F;
jaroslav@1258:         this.intLen -= nInts;
jaroslav@1258:         if (nBits == 0)
jaroslav@1258:             return;
jaroslav@1258:         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
jaroslav@1258:         if (nBits >= bitsInHighWord) {
jaroslav@1258:             this.primitiveLeftShift(32 - nBits);
jaroslav@1258:             this.intLen--;
jaroslav@1258:         } else {
jaroslav@1258:             primitiveRightShift(nBits);
jaroslav@1258:         }
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Left shift this MutableBigInteger n bits.
jaroslav@1258:      */
jaroslav@1258:     void leftShift(int n) {
jaroslav@1258:         /*
jaroslav@1258:          * If there is enough storage space in this MutableBigInteger already
jaroslav@1258:          * the available space will be used. Space to the right of the used
jaroslav@1258:          * ints in the value array is faster to utilize, so the extra space
jaroslav@1258:          * will be taken from the right if possible.
jaroslav@1258:          */
jaroslav@1258:         if (intLen == 0)
jaroslav@1258:            return;
jaroslav@1258:         int nInts = n >>> 5;
jaroslav@1258:         int nBits = n&0x1F;
jaroslav@1258:         int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
jaroslav@1258: 
jaroslav@1258:         // If shift can be done without moving words, do so
jaroslav@1258:         if (n <= (32-bitsInHighWord)) {
jaroslav@1258:             primitiveLeftShift(nBits);
jaroslav@1258:             return;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         int newLen = intLen + nInts +1;
jaroslav@1258:         if (nBits <= (32-bitsInHighWord))
jaroslav@1258:             newLen--;
jaroslav@1258:         if (value.length < newLen) {
jaroslav@1258:             // The array must grow
jaroslav@1258:             int[] result = new int[newLen];
jaroslav@1258:             for (int i=0; i<intLen; i++)
jaroslav@1258:                 result[i] = value[offset+i];
jaroslav@1258:             setValue(result, newLen);
jaroslav@1258:         } else if (value.length - offset >= newLen) {
jaroslav@1258:             // Use space on right
jaroslav@1258:             for(int i=0; i<newLen - intLen; i++)
jaroslav@1258:                 value[offset+intLen+i] = 0;
jaroslav@1258:         } else {
jaroslav@1258:             // Must use space on left
jaroslav@1258:             for (int i=0; i<intLen; i++)
jaroslav@1258:                 value[i] = value[offset+i];
jaroslav@1258:             for (int i=intLen; i<newLen; i++)
jaroslav@1258:                 value[i] = 0;
jaroslav@1258:             offset = 0;
jaroslav@1258:         }
jaroslav@1258:         intLen = newLen;
jaroslav@1258:         if (nBits == 0)
jaroslav@1258:             return;
jaroslav@1258:         if (nBits <= (32-bitsInHighWord))
jaroslav@1258:             primitiveLeftShift(nBits);
jaroslav@1258:         else
jaroslav@1258:             primitiveRightShift(32 -nBits);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * A primitive used for division. This method adds in one multiple of the
jaroslav@1258:      * divisor a back to the dividend result at a specified offset. It is used
jaroslav@1258:      * when qhat was estimated too large, and must be adjusted.
jaroslav@1258:      */
jaroslav@1258:     private int divadd(int[] a, int[] result, int offset) {
jaroslav@1258:         long carry = 0;
jaroslav@1258: 
jaroslav@1258:         for (int j=a.length-1; j >= 0; j--) {
jaroslav@1258:             long sum = (a[j] & LONG_MASK) +
jaroslav@1258:                        (result[j+offset] & LONG_MASK) + carry;
jaroslav@1258:             result[j+offset] = (int)sum;
jaroslav@1258:             carry = sum >>> 32;
jaroslav@1258:         }
jaroslav@1258:         return (int)carry;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * This method is used for division. It multiplies an n word input a by one
jaroslav@1258:      * word input x, and subtracts the n word product from q. This is needed
jaroslav@1258:      * when subtracting qhat*divisor from dividend.
jaroslav@1258:      */
jaroslav@1258:     private int mulsub(int[] q, int[] a, int x, int len, int offset) {
jaroslav@1258:         long xLong = x & LONG_MASK;
jaroslav@1258:         long carry = 0;
jaroslav@1258:         offset += len;
jaroslav@1258: 
jaroslav@1258:         for (int j=len-1; j >= 0; j--) {
jaroslav@1258:             long product = (a[j] & LONG_MASK) * xLong + carry;
jaroslav@1258:             long difference = q[offset] - product;
jaroslav@1258:             q[offset--] = (int)difference;
jaroslav@1258:             carry = (product >>> 32)
jaroslav@1258:                      + (((difference & LONG_MASK) >
jaroslav@1258:                          (((~(int)product) & LONG_MASK))) ? 1:0);
jaroslav@1258:         }
jaroslav@1258:         return (int)carry;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Right shift this MutableBigInteger n bits, where n is
jaroslav@1258:      * less than 32.
jaroslav@1258:      * Assumes that intLen > 0, n > 0 for speed
jaroslav@1258:      */
jaroslav@1258:     private final void primitiveRightShift(int n) {
jaroslav@1258:         int[] val = value;
jaroslav@1258:         int n2 = 32 - n;
jaroslav@1258:         for (int i=offset+intLen-1, c=val[i]; i>offset; i--) {
jaroslav@1258:             int b = c;
jaroslav@1258:             c = val[i-1];
jaroslav@1258:             val[i] = (c << n2) | (b >>> n);
jaroslav@1258:         }
jaroslav@1258:         val[offset] >>>= n;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Left shift this MutableBigInteger n bits, where n is
jaroslav@1258:      * less than 32.
jaroslav@1258:      * Assumes that intLen > 0, n > 0 for speed
jaroslav@1258:      */
jaroslav@1258:     private final void primitiveLeftShift(int n) {
jaroslav@1258:         int[] val = value;
jaroslav@1258:         int n2 = 32 - n;
jaroslav@1258:         for (int i=offset, c=val[i], m=i+intLen-1; i<m; i++) {
jaroslav@1258:             int b = c;
jaroslav@1258:             c = val[i+1];
jaroslav@1258:             val[i] = (b << n) | (c >>> n2);
jaroslav@1258:         }
jaroslav@1258:         val[offset+intLen-1] <<= n;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Adds the contents of two MutableBigInteger objects.The result
jaroslav@1258:      * is placed within this MutableBigInteger.
jaroslav@1258:      * The contents of the addend are not changed.
jaroslav@1258:      */
jaroslav@1258:     void add(MutableBigInteger addend) {
jaroslav@1258:         int x = intLen;
jaroslav@1258:         int y = addend.intLen;
jaroslav@1258:         int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
jaroslav@1258:         int[] result = (value.length < resultLen ? new int[resultLen] : value);
jaroslav@1258: 
jaroslav@1258:         int rstart = result.length-1;
jaroslav@1258:         long sum;
jaroslav@1258:         long carry = 0;
jaroslav@1258: 
jaroslav@1258:         // Add common parts of both numbers
jaroslav@1258:         while(x>0 && y>0) {
jaroslav@1258:             x--; y--;
jaroslav@1258:             sum = (value[x+offset] & LONG_MASK) +
jaroslav@1258:                 (addend.value[y+addend.offset] & LONG_MASK) + carry;
jaroslav@1258:             result[rstart--] = (int)sum;
jaroslav@1258:             carry = sum >>> 32;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // Add remainder of the longer number
jaroslav@1258:         while(x>0) {
jaroslav@1258:             x--;
jaroslav@1258:             if (carry == 0 && result == value && rstart == (x + offset))
jaroslav@1258:                 return;
jaroslav@1258:             sum = (value[x+offset] & LONG_MASK) + carry;
jaroslav@1258:             result[rstart--] = (int)sum;
jaroslav@1258:             carry = sum >>> 32;
jaroslav@1258:         }
jaroslav@1258:         while(y>0) {
jaroslav@1258:             y--;
jaroslav@1258:             sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
jaroslav@1258:             result[rstart--] = (int)sum;
jaroslav@1258:             carry = sum >>> 32;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         if (carry > 0) { // Result must grow in length
jaroslav@1258:             resultLen++;
jaroslav@1258:             if (result.length < resultLen) {
jaroslav@1258:                 int temp[] = new int[resultLen];
jaroslav@1258:                 // Result one word longer from carry-out; copy low-order
jaroslav@1258:                 // bits into new result.
jaroslav@1258:                 System.arraycopy(result, 0, temp, 1, result.length);
jaroslav@1258:                 temp[0] = 1;
jaroslav@1258:                 result = temp;
jaroslav@1258:             } else {
jaroslav@1258:                 result[rstart--] = 1;
jaroslav@1258:             }
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         value = result;
jaroslav@1258:         intLen = resultLen;
jaroslav@1258:         offset = result.length - resultLen;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Subtracts the smaller of this and b from the larger and places the
jaroslav@1258:      * result into this MutableBigInteger.
jaroslav@1258:      */
jaroslav@1258:     int subtract(MutableBigInteger b) {
jaroslav@1258:         MutableBigInteger a = this;
jaroslav@1258: 
jaroslav@1258:         int[] result = value;
jaroslav@1258:         int sign = a.compare(b);
jaroslav@1258: 
jaroslav@1258:         if (sign == 0) {
jaroslav@1258:             reset();
jaroslav@1258:             return 0;
jaroslav@1258:         }
jaroslav@1258:         if (sign < 0) {
jaroslav@1258:             MutableBigInteger tmp = a;
jaroslav@1258:             a = b;
jaroslav@1258:             b = tmp;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         int resultLen = a.intLen;
jaroslav@1258:         if (result.length < resultLen)
jaroslav@1258:             result = new int[resultLen];
jaroslav@1258: 
jaroslav@1258:         long diff = 0;
jaroslav@1258:         int x = a.intLen;
jaroslav@1258:         int y = b.intLen;
jaroslav@1258:         int rstart = result.length - 1;
jaroslav@1258: 
jaroslav@1258:         // Subtract common parts of both numbers
jaroslav@1258:         while (y>0) {
jaroslav@1258:             x--; y--;
jaroslav@1258: 
jaroslav@1258:             diff = (a.value[x+a.offset] & LONG_MASK) -
jaroslav@1258:                    (b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
jaroslav@1258:             result[rstart--] = (int)diff;
jaroslav@1258:         }
jaroslav@1258:         // Subtract remainder of longer number
jaroslav@1258:         while (x>0) {
jaroslav@1258:             x--;
jaroslav@1258:             diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
jaroslav@1258:             result[rstart--] = (int)diff;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         value = result;
jaroslav@1258:         intLen = resultLen;
jaroslav@1258:         offset = value.length - resultLen;
jaroslav@1258:         normalize();
jaroslav@1258:         return sign;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Subtracts the smaller of a and b from the larger and places the result
jaroslav@1258:      * into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
jaroslav@1258:      * operation was performed.
jaroslav@1258:      */
jaroslav@1258:     private int difference(MutableBigInteger b) {
jaroslav@1258:         MutableBigInteger a = this;
jaroslav@1258:         int sign = a.compare(b);
jaroslav@1258:         if (sign ==0)
jaroslav@1258:             return 0;
jaroslav@1258:         if (sign < 0) {
jaroslav@1258:             MutableBigInteger tmp = a;
jaroslav@1258:             a = b;
jaroslav@1258:             b = tmp;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         long diff = 0;
jaroslav@1258:         int x = a.intLen;
jaroslav@1258:         int y = b.intLen;
jaroslav@1258: 
jaroslav@1258:         // Subtract common parts of both numbers
jaroslav@1258:         while (y>0) {
jaroslav@1258:             x--; y--;
jaroslav@1258:             diff = (a.value[a.offset+ x] & LONG_MASK) -
jaroslav@1258:                 (b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
jaroslav@1258:             a.value[a.offset+x] = (int)diff;
jaroslav@1258:         }
jaroslav@1258:         // Subtract remainder of longer number
jaroslav@1258:         while (x>0) {
jaroslav@1258:             x--;
jaroslav@1258:             diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
jaroslav@1258:             a.value[a.offset+x] = (int)diff;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         a.normalize();
jaroslav@1258:         return sign;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Multiply the contents of two MutableBigInteger objects. The result is
jaroslav@1258:      * placed into MutableBigInteger z. The contents of y are not changed.
jaroslav@1258:      */
jaroslav@1258:     void multiply(MutableBigInteger y, MutableBigInteger z) {
jaroslav@1258:         int xLen = intLen;
jaroslav@1258:         int yLen = y.intLen;
jaroslav@1258:         int newLen = xLen + yLen;
jaroslav@1258: 
jaroslav@1258:         // Put z into an appropriate state to receive product
jaroslav@1258:         if (z.value.length < newLen)
jaroslav@1258:             z.value = new int[newLen];
jaroslav@1258:         z.offset = 0;
jaroslav@1258:         z.intLen = newLen;
jaroslav@1258: 
jaroslav@1258:         // The first iteration is hoisted out of the loop to avoid extra add
jaroslav@1258:         long carry = 0;
jaroslav@1258:         for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
jaroslav@1258:                 long product = (y.value[j+y.offset] & LONG_MASK) *
jaroslav@1258:                                (value[xLen-1+offset] & LONG_MASK) + carry;
jaroslav@1258:                 z.value[k] = (int)product;
jaroslav@1258:                 carry = product >>> 32;
jaroslav@1258:         }
jaroslav@1258:         z.value[xLen-1] = (int)carry;
jaroslav@1258: 
jaroslav@1258:         // Perform the multiplication word by word
jaroslav@1258:         for (int i = xLen-2; i >= 0; i--) {
jaroslav@1258:             carry = 0;
jaroslav@1258:             for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
jaroslav@1258:                 long product = (y.value[j+y.offset] & LONG_MASK) *
jaroslav@1258:                                (value[i+offset] & LONG_MASK) +
jaroslav@1258:                                (z.value[k] & LONG_MASK) + carry;
jaroslav@1258:                 z.value[k] = (int)product;
jaroslav@1258:                 carry = product >>> 32;
jaroslav@1258:             }
jaroslav@1258:             z.value[i] = (int)carry;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // Remove leading zeros from product
jaroslav@1258:         z.normalize();
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Multiply the contents of this MutableBigInteger by the word y. The
jaroslav@1258:      * result is placed into z.
jaroslav@1258:      */
jaroslav@1258:     void mul(int y, MutableBigInteger z) {
jaroslav@1258:         if (y == 1) {
jaroslav@1258:             z.copyValue(this);
jaroslav@1258:             return;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         if (y == 0) {
jaroslav@1258:             z.clear();
jaroslav@1258:             return;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // Perform the multiplication word by word
jaroslav@1258:         long ylong = y & LONG_MASK;
jaroslav@1258:         int[] zval = (z.value.length<intLen+1 ? new int[intLen + 1]
jaroslav@1258:                                               : z.value);
jaroslav@1258:         long carry = 0;
jaroslav@1258:         for (int i = intLen-1; i >= 0; i--) {
jaroslav@1258:             long product = ylong * (value[i+offset] & LONG_MASK) + carry;
jaroslav@1258:             zval[i+1] = (int)product;
jaroslav@1258:             carry = product >>> 32;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         if (carry == 0) {
jaroslav@1258:             z.offset = 1;
jaroslav@1258:             z.intLen = intLen;
jaroslav@1258:         } else {
jaroslav@1258:             z.offset = 0;
jaroslav@1258:             z.intLen = intLen + 1;
jaroslav@1258:             zval[0] = (int)carry;
jaroslav@1258:         }
jaroslav@1258:         z.value = zval;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:      /**
jaroslav@1258:      * This method is used for division of an n word dividend by a one word
jaroslav@1258:      * divisor. The quotient is placed into quotient. The one word divisor is
jaroslav@1258:      * specified by divisor.
jaroslav@1258:      *
jaroslav@1258:      * @return the remainder of the division is returned.
jaroslav@1258:      *
jaroslav@1258:      */
jaroslav@1258:     int divideOneWord(int divisor, MutableBigInteger quotient) {
jaroslav@1258:         long divisorLong = divisor & LONG_MASK;
jaroslav@1258: 
jaroslav@1258:         // Special case of one word dividend
jaroslav@1258:         if (intLen == 1) {
jaroslav@1258:             long dividendValue = value[offset] & LONG_MASK;
jaroslav@1258:             int q = (int) (dividendValue / divisorLong);
jaroslav@1258:             int r = (int) (dividendValue - q * divisorLong);
jaroslav@1258:             quotient.value[0] = q;
jaroslav@1258:             quotient.intLen = (q == 0) ? 0 : 1;
jaroslav@1258:             quotient.offset = 0;
jaroslav@1258:             return r;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         if (quotient.value.length < intLen)
jaroslav@1258:             quotient.value = new int[intLen];
jaroslav@1258:         quotient.offset = 0;
jaroslav@1258:         quotient.intLen = intLen;
jaroslav@1258: 
jaroslav@1258:         // Normalize the divisor
jaroslav@1258:         int shift = Integer.numberOfLeadingZeros(divisor);
jaroslav@1258: 
jaroslav@1258:         int rem = value[offset];
jaroslav@1258:         long remLong = rem & LONG_MASK;
jaroslav@1258:         if (remLong < divisorLong) {
jaroslav@1258:             quotient.value[0] = 0;
jaroslav@1258:         } else {
jaroslav@1258:             quotient.value[0] = (int)(remLong / divisorLong);
jaroslav@1258:             rem = (int) (remLong - (quotient.value[0] * divisorLong));
jaroslav@1258:             remLong = rem & LONG_MASK;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         int xlen = intLen;
jaroslav@1258:         int[] qWord = new int[2];
jaroslav@1258:         while (--xlen > 0) {
jaroslav@1258:             long dividendEstimate = (remLong<<32) |
jaroslav@1258:                 (value[offset + intLen - xlen] & LONG_MASK);
jaroslav@1258:             if (dividendEstimate >= 0) {
jaroslav@1258:                 qWord[0] = (int) (dividendEstimate / divisorLong);
jaroslav@1258:                 qWord[1] = (int) (dividendEstimate - qWord[0] * divisorLong);
jaroslav@1258:             } else {
jaroslav@1258:                 divWord(qWord, dividendEstimate, divisor);
jaroslav@1258:             }
jaroslav@1258:             quotient.value[intLen - xlen] = qWord[0];
jaroslav@1258:             rem = qWord[1];
jaroslav@1258:             remLong = rem & LONG_MASK;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         quotient.normalize();
jaroslav@1258:         // Unnormalize
jaroslav@1258:         if (shift > 0)
jaroslav@1258:             return rem % divisor;
jaroslav@1258:         else
jaroslav@1258:             return rem;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Calculates the quotient of this div b and places the quotient in the
jaroslav@1258:      * provided MutableBigInteger objects and the remainder object is returned.
jaroslav@1258:      *
jaroslav@1258:      * Uses Algorithm D in Knuth section 4.3.1.
jaroslav@1258:      * Many optimizations to that algorithm have been adapted from the Colin
jaroslav@1258:      * Plumb C library.
jaroslav@1258:      * It special cases one word divisors for speed. The content of b is not
jaroslav@1258:      * changed.
jaroslav@1258:      *
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
jaroslav@1258:         if (b.intLen == 0)
jaroslav@1258:             throw new ArithmeticException("BigInteger divide by zero");
jaroslav@1258: 
jaroslav@1258:         // Dividend is zero
jaroslav@1258:         if (intLen == 0) {
jaroslav@1258:             quotient.intLen = quotient.offset;
jaroslav@1258:             return new MutableBigInteger();
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         int cmp = compare(b);
jaroslav@1258:         // Dividend less than divisor
jaroslav@1258:         if (cmp < 0) {
jaroslav@1258:             quotient.intLen = quotient.offset = 0;
jaroslav@1258:             return new MutableBigInteger(this);
jaroslav@1258:         }
jaroslav@1258:         // Dividend equal to divisor
jaroslav@1258:         if (cmp == 0) {
jaroslav@1258:             quotient.value[0] = quotient.intLen = 1;
jaroslav@1258:             quotient.offset = 0;
jaroslav@1258:             return new MutableBigInteger();
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         quotient.clear();
jaroslav@1258:         // Special case one word divisor
jaroslav@1258:         if (b.intLen == 1) {
jaroslav@1258:             int r = divideOneWord(b.value[b.offset], quotient);
jaroslav@1258:             if (r == 0)
jaroslav@1258:                 return new MutableBigInteger();
jaroslav@1258:             return new MutableBigInteger(r);
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // Copy divisor value to protect divisor
jaroslav@1258:         int[] div = Arrays.copyOfRange(b.value, b.offset, b.offset + b.intLen);
jaroslav@1258:         return divideMagnitude(div, quotient);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Internally used  to calculate the quotient of this div v and places the
jaroslav@1258:      * quotient in the provided MutableBigInteger object and the remainder is
jaroslav@1258:      * returned.
jaroslav@1258:      *
jaroslav@1258:      * @return the remainder of the division will be returned.
jaroslav@1258:      */
jaroslav@1258:     long divide(long v, MutableBigInteger quotient) {
jaroslav@1258:         if (v == 0)
jaroslav@1258:             throw new ArithmeticException("BigInteger divide by zero");
jaroslav@1258: 
jaroslav@1258:         // Dividend is zero
jaroslav@1258:         if (intLen == 0) {
jaroslav@1258:             quotient.intLen = quotient.offset = 0;
jaroslav@1258:             return 0;
jaroslav@1258:         }
jaroslav@1258:         if (v < 0)
jaroslav@1258:             v = -v;
jaroslav@1258: 
jaroslav@1258:         int d = (int)(v >>> 32);
jaroslav@1258:         quotient.clear();
jaroslav@1258:         // Special case on word divisor
jaroslav@1258:         if (d == 0)
jaroslav@1258:             return divideOneWord((int)v, quotient) & LONG_MASK;
jaroslav@1258:         else {
jaroslav@1258:             int[] div = new int[]{ d, (int)(v & LONG_MASK) };
jaroslav@1258:             return divideMagnitude(div, quotient).toLong();
jaroslav@1258:         }
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Divide this MutableBigInteger by the divisor represented by its magnitude
jaroslav@1258:      * array. The quotient will be placed into the provided quotient object &
jaroslav@1258:      * the remainder object is returned.
jaroslav@1258:      */
jaroslav@1258:     private MutableBigInteger divideMagnitude(int[] divisor,
jaroslav@1258:                                               MutableBigInteger quotient) {
jaroslav@1258: 
jaroslav@1258:         // Remainder starts as dividend with space for a leading zero
jaroslav@1258:         MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
jaroslav@1258:         System.arraycopy(value, offset, rem.value, 1, intLen);
jaroslav@1258:         rem.intLen = intLen;
jaroslav@1258:         rem.offset = 1;
jaroslav@1258: 
jaroslav@1258:         int nlen = rem.intLen;
jaroslav@1258: 
jaroslav@1258:         // Set the quotient size
jaroslav@1258:         int dlen = divisor.length;
jaroslav@1258:         int limit = nlen - dlen + 1;
jaroslav@1258:         if (quotient.value.length < limit) {
jaroslav@1258:             quotient.value = new int[limit];
jaroslav@1258:             quotient.offset = 0;
jaroslav@1258:         }
jaroslav@1258:         quotient.intLen = limit;
jaroslav@1258:         int[] q = quotient.value;
jaroslav@1258: 
jaroslav@1258:         // D1 normalize the divisor
jaroslav@1258:         int shift = Integer.numberOfLeadingZeros(divisor[0]);
jaroslav@1258:         if (shift > 0) {
jaroslav@1258:             // First shift will not grow array
jaroslav@1258:             BigInteger.primitiveLeftShift(divisor, dlen, shift);
jaroslav@1258:             // But this one might
jaroslav@1258:             rem.leftShift(shift);
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // Must insert leading 0 in rem if its length did not change
jaroslav@1258:         if (rem.intLen == nlen) {
jaroslav@1258:             rem.offset = 0;
jaroslav@1258:             rem.value[0] = 0;
jaroslav@1258:             rem.intLen++;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         int dh = divisor[0];
jaroslav@1258:         long dhLong = dh & LONG_MASK;
jaroslav@1258:         int dl = divisor[1];
jaroslav@1258:         int[] qWord = new int[2];
jaroslav@1258: 
jaroslav@1258:         // D2 Initialize j
jaroslav@1258:         for(int j=0; j<limit; j++) {
jaroslav@1258:             // D3 Calculate qhat
jaroslav@1258:             // estimate qhat
jaroslav@1258:             int qhat = 0;
jaroslav@1258:             int qrem = 0;
jaroslav@1258:             boolean skipCorrection = false;
jaroslav@1258:             int nh = rem.value[j+rem.offset];
jaroslav@1258:             int nh2 = nh + 0x80000000;
jaroslav@1258:             int nm = rem.value[j+1+rem.offset];
jaroslav@1258: 
jaroslav@1258:             if (nh == dh) {
jaroslav@1258:                 qhat = ~0;
jaroslav@1258:                 qrem = nh + nm;
jaroslav@1258:                 skipCorrection = qrem + 0x80000000 < nh2;
jaroslav@1258:             } else {
jaroslav@1258:                 long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
jaroslav@1258:                 if (nChunk >= 0) {
jaroslav@1258:                     qhat = (int) (nChunk / dhLong);
jaroslav@1258:                     qrem = (int) (nChunk - (qhat * dhLong));
jaroslav@1258:                 } else {
jaroslav@1258:                     divWord(qWord, nChunk, dh);
jaroslav@1258:                     qhat = qWord[0];
jaroslav@1258:                     qrem = qWord[1];
jaroslav@1258:                 }
jaroslav@1258:             }
jaroslav@1258: 
jaroslav@1258:             if (qhat == 0)
jaroslav@1258:                 continue;
jaroslav@1258: 
jaroslav@1258:             if (!skipCorrection) { // Correct qhat
jaroslav@1258:                 long nl = rem.value[j+2+rem.offset] & LONG_MASK;
jaroslav@1258:                 long rs = ((qrem & LONG_MASK) << 32) | nl;
jaroslav@1258:                 long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
jaroslav@1258: 
jaroslav@1258:                 if (unsignedLongCompare(estProduct, rs)) {
jaroslav@1258:                     qhat--;
jaroslav@1258:                     qrem = (int)((qrem & LONG_MASK) + dhLong);
jaroslav@1258:                     if ((qrem & LONG_MASK) >=  dhLong) {
jaroslav@1258:                         estProduct -= (dl & LONG_MASK);
jaroslav@1258:                         rs = ((qrem & LONG_MASK) << 32) | nl;
jaroslav@1258:                         if (unsignedLongCompare(estProduct, rs))
jaroslav@1258:                             qhat--;
jaroslav@1258:                     }
jaroslav@1258:                 }
jaroslav@1258:             }
jaroslav@1258: 
jaroslav@1258:             // D4 Multiply and subtract
jaroslav@1258:             rem.value[j+rem.offset] = 0;
jaroslav@1258:             int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
jaroslav@1258: 
jaroslav@1258:             // D5 Test remainder
jaroslav@1258:             if (borrow + 0x80000000 > nh2) {
jaroslav@1258:                 // D6 Add back
jaroslav@1258:                 divadd(divisor, rem.value, j+1+rem.offset);
jaroslav@1258:                 qhat--;
jaroslav@1258:             }
jaroslav@1258: 
jaroslav@1258:             // Store the quotient digit
jaroslav@1258:             q[j] = qhat;
jaroslav@1258:         } // D7 loop on j
jaroslav@1258: 
jaroslav@1258:         // D8 Unnormalize
jaroslav@1258:         if (shift > 0)
jaroslav@1258:             rem.rightShift(shift);
jaroslav@1258: 
jaroslav@1258:         quotient.normalize();
jaroslav@1258:         rem.normalize();
jaroslav@1258:         return rem;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Compare two longs as if they were unsigned.
jaroslav@1258:      * Returns true iff one is bigger than two.
jaroslav@1258:      */
jaroslav@1258:     private boolean unsignedLongCompare(long one, long two) {
jaroslav@1258:         return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * This method divides a long quantity by an int to estimate
jaroslav@1258:      * qhat for two multi precision numbers. It is used when
jaroslav@1258:      * the signed value of n is less than zero.
jaroslav@1258:      */
jaroslav@1258:     private void divWord(int[] result, long n, int d) {
jaroslav@1258:         long dLong = d & LONG_MASK;
jaroslav@1258: 
jaroslav@1258:         if (dLong == 1) {
jaroslav@1258:             result[0] = (int)n;
jaroslav@1258:             result[1] = 0;
jaroslav@1258:             return;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // Approximate the quotient and remainder
jaroslav@1258:         long q = (n >>> 1) / (dLong >>> 1);
jaroslav@1258:         long r = n - q*dLong;
jaroslav@1258: 
jaroslav@1258:         // Correct the approximation
jaroslav@1258:         while (r < 0) {
jaroslav@1258:             r += dLong;
jaroslav@1258:             q--;
jaroslav@1258:         }
jaroslav@1258:         while (r >= dLong) {
jaroslav@1258:             r -= dLong;
jaroslav@1258:             q++;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // n - q*dlong == r && 0 <= r <dLong, hence we're done.
jaroslav@1258:         result[0] = (int)q;
jaroslav@1258:         result[1] = (int)r;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Calculate GCD of this and b. This and b are changed by the computation.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger hybridGCD(MutableBigInteger b) {
jaroslav@1258:         // Use Euclid's algorithm until the numbers are approximately the
jaroslav@1258:         // same length, then use the binary GCD algorithm to find the GCD.
jaroslav@1258:         MutableBigInteger a = this;
jaroslav@1258:         MutableBigInteger q = new MutableBigInteger();
jaroslav@1258: 
jaroslav@1258:         while (b.intLen != 0) {
jaroslav@1258:             if (Math.abs(a.intLen - b.intLen) < 2)
jaroslav@1258:                 return a.binaryGCD(b);
jaroslav@1258: 
jaroslav@1258:             MutableBigInteger r = a.divide(b, q);
jaroslav@1258:             a = b;
jaroslav@1258:             b = r;
jaroslav@1258:         }
jaroslav@1258:         return a;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Calculate GCD of this and v.
jaroslav@1258:      * Assumes that this and v are not zero.
jaroslav@1258:      */
jaroslav@1258:     private MutableBigInteger binaryGCD(MutableBigInteger v) {
jaroslav@1258:         // Algorithm B from Knuth section 4.5.2
jaroslav@1258:         MutableBigInteger u = this;
jaroslav@1258:         MutableBigInteger r = new MutableBigInteger();
jaroslav@1258: 
jaroslav@1258:         // step B1
jaroslav@1258:         int s1 = u.getLowestSetBit();
jaroslav@1258:         int s2 = v.getLowestSetBit();
jaroslav@1258:         int k = (s1 < s2) ? s1 : s2;
jaroslav@1258:         if (k != 0) {
jaroslav@1258:             u.rightShift(k);
jaroslav@1258:             v.rightShift(k);
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // step B2
jaroslav@1258:         boolean uOdd = (k==s1);
jaroslav@1258:         MutableBigInteger t = uOdd ? v: u;
jaroslav@1258:         int tsign = uOdd ? -1 : 1;
jaroslav@1258: 
jaroslav@1258:         int lb;
jaroslav@1258:         while ((lb = t.getLowestSetBit()) >= 0) {
jaroslav@1258:             // steps B3 and B4
jaroslav@1258:             t.rightShift(lb);
jaroslav@1258:             // step B5
jaroslav@1258:             if (tsign > 0)
jaroslav@1258:                 u = t;
jaroslav@1258:             else
jaroslav@1258:                 v = t;
jaroslav@1258: 
jaroslav@1258:             // Special case one word numbers
jaroslav@1258:             if (u.intLen < 2 && v.intLen < 2) {
jaroslav@1258:                 int x = u.value[u.offset];
jaroslav@1258:                 int y = v.value[v.offset];
jaroslav@1258:                 x  = binaryGcd(x, y);
jaroslav@1258:                 r.value[0] = x;
jaroslav@1258:                 r.intLen = 1;
jaroslav@1258:                 r.offset = 0;
jaroslav@1258:                 if (k > 0)
jaroslav@1258:                     r.leftShift(k);
jaroslav@1258:                 return r;
jaroslav@1258:             }
jaroslav@1258: 
jaroslav@1258:             // step B6
jaroslav@1258:             if ((tsign = u.difference(v)) == 0)
jaroslav@1258:                 break;
jaroslav@1258:             t = (tsign >= 0) ? u : v;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         if (k > 0)
jaroslav@1258:             u.leftShift(k);
jaroslav@1258:         return u;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Calculate GCD of a and b interpreted as unsigned integers.
jaroslav@1258:      */
jaroslav@1258:     static int binaryGcd(int a, int b) {
jaroslav@1258:         if (b==0)
jaroslav@1258:             return a;
jaroslav@1258:         if (a==0)
jaroslav@1258:             return b;
jaroslav@1258: 
jaroslav@1258:         // Right shift a & b till their last bits equal to 1.
jaroslav@1258:         int aZeros = Integer.numberOfTrailingZeros(a);
jaroslav@1258:         int bZeros = Integer.numberOfTrailingZeros(b);
jaroslav@1258:         a >>>= aZeros;
jaroslav@1258:         b >>>= bZeros;
jaroslav@1258: 
jaroslav@1258:         int t = (aZeros < bZeros ? aZeros : bZeros);
jaroslav@1258: 
jaroslav@1258:         while (a != b) {
jaroslav@1258:             if ((a+0x80000000) > (b+0x80000000)) {  // a > b as unsigned
jaroslav@1258:                 a -= b;
jaroslav@1258:                 a >>>= Integer.numberOfTrailingZeros(a);
jaroslav@1258:             } else {
jaroslav@1258:                 b -= a;
jaroslav@1258:                 b >>>= Integer.numberOfTrailingZeros(b);
jaroslav@1258:             }
jaroslav@1258:         }
jaroslav@1258:         return a<<t;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Returns the modInverse of this mod p. This and p are not affected by
jaroslav@1258:      * the operation.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger mutableModInverse(MutableBigInteger p) {
jaroslav@1258:         // Modulus is odd, use Schroeppel's algorithm
jaroslav@1258:         if (p.isOdd())
jaroslav@1258:             return modInverse(p);
jaroslav@1258: 
jaroslav@1258:         // Base and modulus are even, throw exception
jaroslav@1258:         if (isEven())
jaroslav@1258:             throw new ArithmeticException("BigInteger not invertible.");
jaroslav@1258: 
jaroslav@1258:         // Get even part of modulus expressed as a power of 2
jaroslav@1258:         int powersOf2 = p.getLowestSetBit();
jaroslav@1258: 
jaroslav@1258:         // Construct odd part of modulus
jaroslav@1258:         MutableBigInteger oddMod = new MutableBigInteger(p);
jaroslav@1258:         oddMod.rightShift(powersOf2);
jaroslav@1258: 
jaroslav@1258:         if (oddMod.isOne())
jaroslav@1258:             return modInverseMP2(powersOf2);
jaroslav@1258: 
jaroslav@1258:         // Calculate 1/a mod oddMod
jaroslav@1258:         MutableBigInteger oddPart = modInverse(oddMod);
jaroslav@1258: 
jaroslav@1258:         // Calculate 1/a mod evenMod
jaroslav@1258:         MutableBigInteger evenPart = modInverseMP2(powersOf2);
jaroslav@1258: 
jaroslav@1258:         // Combine the results using Chinese Remainder Theorem
jaroslav@1258:         MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
jaroslav@1258:         MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
jaroslav@1258: 
jaroslav@1258:         MutableBigInteger temp1 = new MutableBigInteger();
jaroslav@1258:         MutableBigInteger temp2 = new MutableBigInteger();
jaroslav@1258:         MutableBigInteger result = new MutableBigInteger();
jaroslav@1258: 
jaroslav@1258:         oddPart.leftShift(powersOf2);
jaroslav@1258:         oddPart.multiply(y1, result);
jaroslav@1258: 
jaroslav@1258:         evenPart.multiply(oddMod, temp1);
jaroslav@1258:         temp1.multiply(y2, temp2);
jaroslav@1258: 
jaroslav@1258:         result.add(temp2);
jaroslav@1258:         return result.divide(p, temp1);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /*
jaroslav@1258:      * Calculate the multiplicative inverse of this mod 2^k.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger modInverseMP2(int k) {
jaroslav@1258:         if (isEven())
jaroslav@1258:             throw new ArithmeticException("Non-invertible. (GCD != 1)");
jaroslav@1258: 
jaroslav@1258:         if (k > 64)
jaroslav@1258:             return euclidModInverse(k);
jaroslav@1258: 
jaroslav@1258:         int t = inverseMod32(value[offset+intLen-1]);
jaroslav@1258: 
jaroslav@1258:         if (k < 33) {
jaroslav@1258:             t = (k == 32 ? t : t & ((1 << k) - 1));
jaroslav@1258:             return new MutableBigInteger(t);
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         long pLong = (value[offset+intLen-1] & LONG_MASK);
jaroslav@1258:         if (intLen > 1)
jaroslav@1258:             pLong |=  ((long)value[offset+intLen-2] << 32);
jaroslav@1258:         long tLong = t & LONG_MASK;
jaroslav@1258:         tLong = tLong * (2 - pLong * tLong);  // 1 more Newton iter step
jaroslav@1258:         tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
jaroslav@1258: 
jaroslav@1258:         MutableBigInteger result = new MutableBigInteger(new int[2]);
jaroslav@1258:         result.value[0] = (int)(tLong >>> 32);
jaroslav@1258:         result.value[1] = (int)tLong;
jaroslav@1258:         result.intLen = 2;
jaroslav@1258:         result.normalize();
jaroslav@1258:         return result;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /*
jaroslav@1258:      * Returns the multiplicative inverse of val mod 2^32.  Assumes val is odd.
jaroslav@1258:      */
jaroslav@1258:     static int inverseMod32(int val) {
jaroslav@1258:         // Newton's iteration!
jaroslav@1258:         int t = val;
jaroslav@1258:         t *= 2 - val*t;
jaroslav@1258:         t *= 2 - val*t;
jaroslav@1258:         t *= 2 - val*t;
jaroslav@1258:         t *= 2 - val*t;
jaroslav@1258:         return t;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /*
jaroslav@1258:      * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
jaroslav@1258:      */
jaroslav@1258:     static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
jaroslav@1258:         // Copy the mod to protect original
jaroslav@1258:         return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Calculate the multiplicative inverse of this mod mod, where mod is odd.
jaroslav@1258:      * This and mod are not changed by the calculation.
jaroslav@1258:      *
jaroslav@1258:      * This method implements an algorithm due to Richard Schroeppel, that uses
jaroslav@1258:      * the same intermediate representation as Montgomery Reduction
jaroslav@1258:      * ("Montgomery Form").  The algorithm is described in an unpublished
jaroslav@1258:      * manuscript entitled "Fast Modular Reciprocals."
jaroslav@1258:      */
jaroslav@1258:     private MutableBigInteger modInverse(MutableBigInteger mod) {
jaroslav@1258:         MutableBigInteger p = new MutableBigInteger(mod);
jaroslav@1258:         MutableBigInteger f = new MutableBigInteger(this);
jaroslav@1258:         MutableBigInteger g = new MutableBigInteger(p);
jaroslav@1258:         SignedMutableBigInteger c = new SignedMutableBigInteger(1);
jaroslav@1258:         SignedMutableBigInteger d = new SignedMutableBigInteger();
jaroslav@1258:         MutableBigInteger temp = null;
jaroslav@1258:         SignedMutableBigInteger sTemp = null;
jaroslav@1258: 
jaroslav@1258:         int k = 0;
jaroslav@1258:         // Right shift f k times until odd, left shift d k times
jaroslav@1258:         if (f.isEven()) {
jaroslav@1258:             int trailingZeros = f.getLowestSetBit();
jaroslav@1258:             f.rightShift(trailingZeros);
jaroslav@1258:             d.leftShift(trailingZeros);
jaroslav@1258:             k = trailingZeros;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // The Almost Inverse Algorithm
jaroslav@1258:         while(!f.isOne()) {
jaroslav@1258:             // If gcd(f, g) != 1, number is not invertible modulo mod
jaroslav@1258:             if (f.isZero())
jaroslav@1258:                 throw new ArithmeticException("BigInteger not invertible.");
jaroslav@1258: 
jaroslav@1258:             // If f < g exchange f, g and c, d
jaroslav@1258:             if (f.compare(g) < 0) {
jaroslav@1258:                 temp = f; f = g; g = temp;
jaroslav@1258:                 sTemp = d; d = c; c = sTemp;
jaroslav@1258:             }
jaroslav@1258: 
jaroslav@1258:             // If f == g (mod 4)
jaroslav@1258:             if (((f.value[f.offset + f.intLen - 1] ^
jaroslav@1258:                  g.value[g.offset + g.intLen - 1]) & 3) == 0) {
jaroslav@1258:                 f.subtract(g);
jaroslav@1258:                 c.signedSubtract(d);
jaroslav@1258:             } else { // If f != g (mod 4)
jaroslav@1258:                 f.add(g);
jaroslav@1258:                 c.signedAdd(d);
jaroslav@1258:             }
jaroslav@1258: 
jaroslav@1258:             // Right shift f k times until odd, left shift d k times
jaroslav@1258:             int trailingZeros = f.getLowestSetBit();
jaroslav@1258:             f.rightShift(trailingZeros);
jaroslav@1258:             d.leftShift(trailingZeros);
jaroslav@1258:             k += trailingZeros;
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         while (c.sign < 0)
jaroslav@1258:            c.signedAdd(p);
jaroslav@1258: 
jaroslav@1258:         return fixup(c, p, k);
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /*
jaroslav@1258:      * The Fixup Algorithm
jaroslav@1258:      * Calculates X such that X = C * 2^(-k) (mod P)
jaroslav@1258:      * Assumes C<P and P is odd.
jaroslav@1258:      */
jaroslav@1258:     static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
jaroslav@1258:                                                                       int k) {
jaroslav@1258:         MutableBigInteger temp = new MutableBigInteger();
jaroslav@1258:         // Set r to the multiplicative inverse of p mod 2^32
jaroslav@1258:         int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
jaroslav@1258: 
jaroslav@1258:         for(int i=0, numWords = k >> 5; i<numWords; i++) {
jaroslav@1258:             // V = R * c (mod 2^j)
jaroslav@1258:             int  v = r * c.value[c.offset + c.intLen-1];
jaroslav@1258:             // c = c + (v * p)
jaroslav@1258:             p.mul(v, temp);
jaroslav@1258:             c.add(temp);
jaroslav@1258:             // c = c / 2^j
jaroslav@1258:             c.intLen--;
jaroslav@1258:         }
jaroslav@1258:         int numBits = k & 0x1f;
jaroslav@1258:         if (numBits != 0) {
jaroslav@1258:             // V = R * c (mod 2^j)
jaroslav@1258:             int v = r * c.value[c.offset + c.intLen-1];
jaroslav@1258:             v &= ((1<<numBits) - 1);
jaroslav@1258:             // c = c + (v * p)
jaroslav@1258:             p.mul(v, temp);
jaroslav@1258:             c.add(temp);
jaroslav@1258:             // c = c / 2^j
jaroslav@1258:             c.rightShift(numBits);
jaroslav@1258:         }
jaroslav@1258: 
jaroslav@1258:         // In theory, c may be greater than p at this point (Very rare!)
jaroslav@1258:         while (c.compare(p) >= 0)
jaroslav@1258:             c.subtract(p);
jaroslav@1258: 
jaroslav@1258:         return c;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258:     /**
jaroslav@1258:      * Uses the extended Euclidean algorithm to compute the modInverse of base
jaroslav@1258:      * mod a modulus that is a power of 2. The modulus is 2^k.
jaroslav@1258:      */
jaroslav@1258:     MutableBigInteger euclidModInverse(int k) {
jaroslav@1258:         MutableBigInteger b = new MutableBigInteger(1);
jaroslav@1258:         b.leftShift(k);
jaroslav@1258:         MutableBigInteger mod = new MutableBigInteger(b);
jaroslav@1258: 
jaroslav@1258:         MutableBigInteger a = new MutableBigInteger(this);
jaroslav@1258:         MutableBigInteger q = new MutableBigInteger();
jaroslav@1258:         MutableBigInteger r = b.divide(a, q);
jaroslav@1258: 
jaroslav@1258:         MutableBigInteger swapper = b;
jaroslav@1258:         // swap b & r
jaroslav@1258:         b = r;
jaroslav@1258:         r = swapper;
jaroslav@1258: 
jaroslav@1258:         MutableBigInteger t1 = new MutableBigInteger(q);
jaroslav@1258:         MutableBigInteger t0 = new MutableBigInteger(1);
jaroslav@1258:         MutableBigInteger temp = new MutableBigInteger();
jaroslav@1258: 
jaroslav@1258:         while (!b.isOne()) {
jaroslav@1258:             r = a.divide(b, q);
jaroslav@1258: 
jaroslav@1258:             if (r.intLen == 0)
jaroslav@1258:                 throw new ArithmeticException("BigInteger not invertible.");
jaroslav@1258: 
jaroslav@1258:             swapper = r;
jaroslav@1258:             a = swapper;
jaroslav@1258: 
jaroslav@1258:             if (q.intLen == 1)
jaroslav@1258:                 t1.mul(q.value[q.offset], temp);
jaroslav@1258:             else
jaroslav@1258:                 q.multiply(t1, temp);
jaroslav@1258:             swapper = q;
jaroslav@1258:             q = temp;
jaroslav@1258:             temp = swapper;
jaroslav@1258:             t0.add(q);
jaroslav@1258: 
jaroslav@1258:             if (a.isOne())
jaroslav@1258:                 return t0;
jaroslav@1258: 
jaroslav@1258:             r = b.divide(a, q);
jaroslav@1258: 
jaroslav@1258:             if (r.intLen == 0)
jaroslav@1258:                 throw new ArithmeticException("BigInteger not invertible.");
jaroslav@1258: 
jaroslav@1258:             swapper = b;
jaroslav@1258:             b =  r;
jaroslav@1258: 
jaroslav@1258:             if (q.intLen == 1)
jaroslav@1258:                 t0.mul(q.value[q.offset], temp);
jaroslav@1258:             else
jaroslav@1258:                 q.multiply(t0, temp);
jaroslav@1258:             swapper = q; q = temp; temp = swapper;
jaroslav@1258: 
jaroslav@1258:             t1.add(q);
jaroslav@1258:         }
jaroslav@1258:         mod.subtract(t1);
jaroslav@1258:         return mod;
jaroslav@1258:     }
jaroslav@1258: 
jaroslav@1258: }