META-INF.modules.java.base.classes.java.math.MutableBigInteger Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of java.base Show documentation
Show all versions of java.base Show documentation
Bytecoder java.base Module
/*
* Copyright (c) 1999, 2015, Oracle and/or its affiliates. All rights reserved.
* DO NOT ALTER OR REMOVE COPYRIGHT NOTICES OR THIS FILE HEADER.
*
* This code is free software; you can redistribute it and/or modify it
* under the terms of the GNU General Public License version 2 only, as
* published by the Free Software Foundation. Oracle designates this
* particular file as subject to the "Classpath" exception as provided
* by Oracle in the LICENSE file that accompanied this code.
*
* This code is distributed in the hope that it will be useful, but WITHOUT
* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
* FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
* version 2 for more details (a copy is included in the LICENSE file that
* accompanied this code).
*
* You should have received a copy of the GNU General Public License version
* 2 along with this work; if not, write to the Free Software Foundation,
* Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA.
*
* Please contact Oracle, 500 Oracle Parkway, Redwood Shores, CA 94065 USA
* or visit www.oracle.com if you need additional information or have any
* questions.
*/
package java.math;
/**
* A class used to represent multiprecision integers that makes efficient
* use of allocated space by allowing a number to occupy only part of
* an array so that the arrays do not have to be reallocated as often.
* When performing an operation with many iterations the array used to
* hold a number is only reallocated when necessary and does not have to
* be the same size as the number it represents. A mutable number allows
* calculations to occur on the same number without having to create
* a new number for every step of the calculation as occurs with
* BigIntegers.
*
* @see BigInteger
* @author Michael McCloskey
* @author Timothy Buktu
* @since 1.3
*/
import static java.math.BigDecimal.INFLATED;
import static java.math.BigInteger.LONG_MASK;
import java.util.Arrays;
class MutableBigInteger {
/**
* Holds the magnitude of this MutableBigInteger in big endian order.
* The magnitude may start at an offset into the value array, and it may
* end before the length of the value array.
*/
int[] value;
/**
* The number of ints of the value array that are currently used
* to hold the magnitude of this MutableBigInteger. The magnitude starts
* at an offset and offset + intLen may be less than value.length.
*/
int intLen;
/**
* The offset into the value array where the magnitude of this
* MutableBigInteger begins.
*/
int offset = 0;
// Constants
/**
* MutableBigInteger with one element value array with the value 1. Used by
* BigDecimal divideAndRound to increment the quotient. Use this constant
* only when the method is not going to modify this object.
*/
static final MutableBigInteger ONE = new MutableBigInteger(1);
/**
* The minimum {@code intLen} for cancelling powers of two before
* dividing.
* If the number of ints is less than this threshold,
* {@code divideKnuth} does not eliminate common powers of two from
* the dividend and divisor.
*/
static final int KNUTH_POW2_THRESH_LEN = 6;
/**
* The minimum number of trailing zero ints for cancelling powers of two
* before dividing.
* If the dividend and divisor don't share at least this many zero ints
* at the end, {@code divideKnuth} does not eliminate common powers
* of two from the dividend and divisor.
*/
static final int KNUTH_POW2_THRESH_ZEROS = 3;
// Constructors
/**
* The default constructor. An empty MutableBigInteger is created with
* a one word capacity.
*/
MutableBigInteger() {
value = new int[1];
intLen = 0;
}
/**
* Construct a new MutableBigInteger with a magnitude specified by
* the int val.
*/
MutableBigInteger(int val) {
value = new int[1];
intLen = 1;
value[0] = val;
}
/**
* Construct a new MutableBigInteger with the specified value array
* up to the length of the array supplied.
*/
MutableBigInteger(int[] val) {
value = val;
intLen = val.length;
}
/**
* Construct a new MutableBigInteger with a magnitude equal to the
* specified BigInteger.
*/
MutableBigInteger(BigInteger b) {
intLen = b.mag.length;
value = Arrays.copyOf(b.mag, intLen);
}
/**
* Construct a new MutableBigInteger with a magnitude equal to the
* specified MutableBigInteger.
*/
MutableBigInteger(MutableBigInteger val) {
intLen = val.intLen;
value = Arrays.copyOfRange(val.value, val.offset, val.offset + intLen);
}
/**
* Makes this number an {@code n}-int number all of whose bits are ones.
* Used by Burnikel-Ziegler division.
* @param n number of ints in the {@code value} array
* @return a number equal to {@code ((1<<(32*n)))-1}
*/
private void ones(int n) {
if (n > value.length)
value = new int[n];
Arrays.fill(value, -1);
offset = 0;
intLen = n;
}
/**
* Internal helper method to return the magnitude array. The caller is not
* supposed to modify the returned array.
*/
private int[] getMagnitudeArray() {
if (offset > 0 || value.length != intLen)
return Arrays.copyOfRange(value, offset, offset + intLen);
return value;
}
/**
* Convert this MutableBigInteger to a long value. The caller has to make
* sure this MutableBigInteger can be fit into long.
*/
private long toLong() {
assert (intLen <= 2) : "this MutableBigInteger exceeds the range of long";
if (intLen == 0)
return 0;
long d = value[offset] & LONG_MASK;
return (intLen == 2) ? d << 32 | (value[offset + 1] & LONG_MASK) : d;
}
/**
* Convert this MutableBigInteger to a BigInteger object.
*/
BigInteger toBigInteger(int sign) {
if (intLen == 0 || sign == 0)
return BigInteger.ZERO;
return new BigInteger(getMagnitudeArray(), sign);
}
/**
* Converts this number to a nonnegative {@code BigInteger}.
*/
BigInteger toBigInteger() {
normalize();
return toBigInteger(isZero() ? 0 : 1);
}
/**
* Convert this MutableBigInteger to BigDecimal object with the specified sign
* and scale.
*/
BigDecimal toBigDecimal(int sign, int scale) {
if (intLen == 0 || sign == 0)
return BigDecimal.zeroValueOf(scale);
int[] mag = getMagnitudeArray();
int len = mag.length;
int d = mag[0];
// If this MutableBigInteger can't be fit into long, we need to
// make a BigInteger object for the resultant BigDecimal object.
if (len > 2 || (d < 0 && len == 2))
return new BigDecimal(new BigInteger(mag, sign), INFLATED, scale, 0);
long v = (len == 2) ?
((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
d & LONG_MASK;
return BigDecimal.valueOf(sign == -1 ? -v : v, scale);
}
/**
* This is for internal use in converting from a MutableBigInteger
* object into a long value given a specified sign.
* returns INFLATED if value is not fit into long
*/
long toCompactValue(int sign) {
if (intLen == 0 || sign == 0)
return 0L;
int[] mag = getMagnitudeArray();
int len = mag.length;
int d = mag[0];
// If this MutableBigInteger can not be fitted into long, we need to
// make a BigInteger object for the resultant BigDecimal object.
if (len > 2 || (d < 0 && len == 2))
return INFLATED;
long v = (len == 2) ?
((mag[1] & LONG_MASK) | (d & LONG_MASK) << 32) :
d & LONG_MASK;
return sign == -1 ? -v : v;
}
/**
* Clear out a MutableBigInteger for reuse.
*/
void clear() {
offset = intLen = 0;
for (int index=0, n=value.length; index < n; index++)
value[index] = 0;
}
/**
* Set a MutableBigInteger to zero, removing its offset.
*/
void reset() {
offset = intLen = 0;
}
/**
* Compare the magnitude of two MutableBigIntegers. Returns -1, 0 or 1
* as this MutableBigInteger is numerically less than, equal to, or
* greater than {@code b}.
*/
final int compare(MutableBigInteger b) {
int blen = b.intLen;
if (intLen < blen)
return -1;
if (intLen > blen)
return 1;
// Add Integer.MIN_VALUE to make the comparison act as unsigned integer
// comparison.
int[] bval = b.value;
for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
int b1 = value[i] + 0x80000000;
int b2 = bval[j] + 0x80000000;
if (b1 < b2)
return -1;
if (b1 > b2)
return 1;
}
return 0;
}
/**
* Returns a value equal to what {@code b.leftShift(32*ints); return compare(b);}
* would return, but doesn't change the value of {@code b}.
*/
private int compareShifted(MutableBigInteger b, int ints) {
int blen = b.intLen;
int alen = intLen - ints;
if (alen < blen)
return -1;
if (alen > blen)
return 1;
// Add Integer.MIN_VALUE to make the comparison act as unsigned integer
// comparison.
int[] bval = b.value;
for (int i = offset, j = b.offset; i < alen + offset; i++, j++) {
int b1 = value[i] + 0x80000000;
int b2 = bval[j] + 0x80000000;
if (b1 < b2)
return -1;
if (b1 > b2)
return 1;
}
return 0;
}
/**
* Compare this against half of a MutableBigInteger object (Needed for
* remainder tests).
* Assumes no leading unnecessary zeros, which holds for results
* from divide().
*/
final int compareHalf(MutableBigInteger b) {
int blen = b.intLen;
int len = intLen;
if (len <= 0)
return blen <= 0 ? 0 : -1;
if (len > blen)
return 1;
if (len < blen - 1)
return -1;
int[] bval = b.value;
int bstart = 0;
int carry = 0;
// Only 2 cases left:len == blen or len == blen - 1
if (len != blen) { // len == blen - 1
if (bval[bstart] == 1) {
++bstart;
carry = 0x80000000;
} else
return -1;
}
// compare values with right-shifted values of b,
// carrying shifted-out bits across words
int[] val = value;
for (int i = offset, j = bstart; i < len + offset;) {
int bv = bval[j++];
long hb = ((bv >>> 1) + carry) & LONG_MASK;
long v = val[i++] & LONG_MASK;
if (v != hb)
return v < hb ? -1 : 1;
carry = (bv & 1) << 31; // carray will be either 0x80000000 or 0
}
return carry == 0 ? 0 : -1;
}
/**
* Return the index of the lowest set bit in this MutableBigInteger. If the
* magnitude of this MutableBigInteger is zero, -1 is returned.
*/
private final int getLowestSetBit() {
if (intLen == 0)
return -1;
int j, b;
for (j=intLen-1; (j > 0) && (value[j+offset] == 0); j--)
;
b = value[j+offset];
if (b == 0)
return -1;
return ((intLen-1-j)<<5) + Integer.numberOfTrailingZeros(b);
}
/**
* Return the int in use in this MutableBigInteger at the specified
* index. This method is not used because it is not inlined on all
* platforms.
*/
private final int getInt(int index) {
return value[offset+index];
}
/**
* Return a long which is equal to the unsigned value of the int in
* use in this MutableBigInteger at the specified index. This method is
* not used because it is not inlined on all platforms.
*/
private final long getLong(int index) {
return value[offset+index] & LONG_MASK;
}
/**
* Ensure that the MutableBigInteger is in normal form, specifically
* making sure that there are no leading zeros, and that if the
* magnitude is zero, then intLen is zero.
*/
final void normalize() {
if (intLen == 0) {
offset = 0;
return;
}
int index = offset;
if (value[index] != 0)
return;
int indexBound = index+intLen;
do {
index++;
} while(index < indexBound && value[index] == 0);
int numZeros = index - offset;
intLen -= numZeros;
offset = (intLen == 0 ? 0 : offset+numZeros);
}
/**
* If this MutableBigInteger cannot hold len words, increase the size
* of the value array to len words.
*/
private final void ensureCapacity(int len) {
if (value.length < len) {
value = new int[len];
offset = 0;
intLen = len;
}
}
/**
* Convert this MutableBigInteger into an int array with no leading
* zeros, of a length that is equal to this MutableBigInteger's intLen.
*/
int[] toIntArray() {
int[] result = new int[intLen];
for(int i=0; i < intLen; i++)
result[i] = value[offset+i];
return result;
}
/**
* Sets the int at index+offset in this MutableBigInteger to val.
* This does not get inlined on all platforms so it is not used
* as often as originally intended.
*/
void setInt(int index, int val) {
value[offset + index] = val;
}
/**
* Sets this MutableBigInteger's value array to the specified array.
* The intLen is set to the specified length.
*/
void setValue(int[] val, int length) {
value = val;
intLen = length;
offset = 0;
}
/**
* Sets this MutableBigInteger's value array to a copy of the specified
* array. The intLen is set to the length of the new array.
*/
void copyValue(MutableBigInteger src) {
int len = src.intLen;
if (value.length < len)
value = new int[len];
System.arraycopy(src.value, src.offset, value, 0, len);
intLen = len;
offset = 0;
}
/**
* Sets this MutableBigInteger's value array to a copy of the specified
* array. The intLen is set to the length of the specified array.
*/
void copyValue(int[] val) {
int len = val.length;
if (value.length < len)
value = new int[len];
System.arraycopy(val, 0, value, 0, len);
intLen = len;
offset = 0;
}
/**
* Returns true iff this MutableBigInteger has a value of one.
*/
boolean isOne() {
return (intLen == 1) && (value[offset] == 1);
}
/**
* Returns true iff this MutableBigInteger has a value of zero.
*/
boolean isZero() {
return (intLen == 0);
}
/**
* Returns true iff this MutableBigInteger is even.
*/
boolean isEven() {
return (intLen == 0) || ((value[offset + intLen - 1] & 1) == 0);
}
/**
* Returns true iff this MutableBigInteger is odd.
*/
boolean isOdd() {
return isZero() ? false : ((value[offset + intLen - 1] & 1) == 1);
}
/**
* Returns true iff this MutableBigInteger is in normal form. A
* MutableBigInteger is in normal form if it has no leading zeros
* after the offset, and intLen + offset <= value.length.
*/
boolean isNormal() {
if (intLen + offset > value.length)
return false;
if (intLen == 0)
return true;
return (value[offset] != 0);
}
/**
* Returns a String representation of this MutableBigInteger in radix 10.
*/
public String toString() {
BigInteger b = toBigInteger(1);
return b.toString();
}
/**
* Like {@link #rightShift(int)} but {@code n} can be greater than the length of the number.
*/
void safeRightShift(int n) {
if (n/32 >= intLen) {
reset();
} else {
rightShift(n);
}
}
/**
* Right shift this MutableBigInteger n bits. The MutableBigInteger is left
* in normal form.
*/
void rightShift(int n) {
if (intLen == 0)
return;
int nInts = n >>> 5;
int nBits = n & 0x1F;
this.intLen -= nInts;
if (nBits == 0)
return;
int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
if (nBits >= bitsInHighWord) {
this.primitiveLeftShift(32 - nBits);
this.intLen--;
} else {
primitiveRightShift(nBits);
}
}
/**
* Like {@link #leftShift(int)} but {@code n} can be zero.
*/
void safeLeftShift(int n) {
if (n > 0) {
leftShift(n);
}
}
/**
* Left shift this MutableBigInteger n bits.
*/
void leftShift(int n) {
/*
* If there is enough storage space in this MutableBigInteger already
* the available space will be used. Space to the right of the used
* ints in the value array is faster to utilize, so the extra space
* will be taken from the right if possible.
*/
if (intLen == 0)
return;
int nInts = n >>> 5;
int nBits = n&0x1F;
int bitsInHighWord = BigInteger.bitLengthForInt(value[offset]);
// If shift can be done without moving words, do so
if (n <= (32-bitsInHighWord)) {
primitiveLeftShift(nBits);
return;
}
int newLen = intLen + nInts +1;
if (nBits <= (32-bitsInHighWord))
newLen--;
if (value.length < newLen) {
// The array must grow
int[] result = new int[newLen];
for (int i=0; i < intLen; i++)
result[i] = value[offset+i];
setValue(result, newLen);
} else if (value.length - offset >= newLen) {
// Use space on right
for(int i=0; i < newLen - intLen; i++)
value[offset+intLen+i] = 0;
} else {
// Must use space on left
for (int i=0; i < intLen; i++)
value[i] = value[offset+i];
for (int i=intLen; i < newLen; i++)
value[i] = 0;
offset = 0;
}
intLen = newLen;
if (nBits == 0)
return;
if (nBits <= (32-bitsInHighWord))
primitiveLeftShift(nBits);
else
primitiveRightShift(32 -nBits);
}
/**
* A primitive used for division. This method adds in one multiple of the
* divisor a back to the dividend result at a specified offset. It is used
* when qhat was estimated too large, and must be adjusted.
*/
private int divadd(int[] a, int[] result, int offset) {
long carry = 0;
for (int j=a.length-1; j >= 0; j--) {
long sum = (a[j] & LONG_MASK) +
(result[j+offset] & LONG_MASK) + carry;
result[j+offset] = (int)sum;
carry = sum >>> 32;
}
return (int)carry;
}
/**
* This method is used for division. It multiplies an n word input a by one
* word input x, and subtracts the n word product from q. This is needed
* when subtracting qhat*divisor from dividend.
*/
private int mulsub(int[] q, int[] a, int x, int len, int offset) {
long xLong = x & LONG_MASK;
long carry = 0;
offset += len;
for (int j=len-1; j >= 0; j--) {
long product = (a[j] & LONG_MASK) * xLong + carry;
long difference = q[offset] - product;
q[offset--] = (int)difference;
carry = (product >>> 32)
+ (((difference & LONG_MASK) >
(((~(int)product) & LONG_MASK))) ? 1:0);
}
return (int)carry;
}
/**
* The method is the same as mulsun, except the fact that q array is not
* updated, the only result of the method is borrow flag.
*/
private int mulsubBorrow(int[] q, int[] a, int x, int len, int offset) {
long xLong = x & LONG_MASK;
long carry = 0;
offset += len;
for (int j=len-1; j >= 0; j--) {
long product = (a[j] & LONG_MASK) * xLong + carry;
long difference = q[offset--] - product;
carry = (product >>> 32)
+ (((difference & LONG_MASK) >
(((~(int)product) & LONG_MASK))) ? 1:0);
}
return (int)carry;
}
/**
* Right shift this MutableBigInteger n bits, where n is
* less than 32.
* Assumes that intLen > 0, n > 0 for speed
*/
private final void primitiveRightShift(int n) {
int[] val = value;
int n2 = 32 - n;
for (int i=offset+intLen-1, c=val[i]; i > offset; i--) {
int b = c;
c = val[i-1];
val[i] = (c << n2) | (b >>> n);
}
val[offset] >>>= n;
}
/**
* Left shift this MutableBigInteger n bits, where n is
* less than 32.
* Assumes that intLen > 0, n > 0 for speed
*/
private final void primitiveLeftShift(int n) {
int[] val = value;
int n2 = 32 - n;
for (int i=offset, c=val[i], m=i+intLen-1; i < m; i++) {
int b = c;
c = val[i+1];
val[i] = (b << n) | (c >>> n2);
}
val[offset+intLen-1] <<= n;
}
/**
* Returns a {@code BigInteger} equal to the {@code n}
* low ints of this number.
*/
private BigInteger getLower(int n) {
if (isZero()) {
return BigInteger.ZERO;
} else if (intLen < n) {
return toBigInteger(1);
} else {
// strip zeros
int len = n;
while (len > 0 && value[offset+intLen-len] == 0)
len--;
int sign = len > 0 ? 1 : 0;
return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
}
}
/**
* Discards all ints whose index is greater than {@code n}.
*/
private void keepLower(int n) {
if (intLen >= n) {
offset += intLen - n;
intLen = n;
}
}
/**
* Adds the contents of two MutableBigInteger objects.The result
* is placed within this MutableBigInteger.
* The contents of the addend are not changed.
*/
void add(MutableBigInteger addend) {
int x = intLen;
int y = addend.intLen;
int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
int[] result = (value.length < resultLen ? new int[resultLen] : value);
int rstart = result.length-1;
long sum;
long carry = 0;
// Add common parts of both numbers
while(x > 0 && y > 0) {
x--; y--;
sum = (value[x+offset] & LONG_MASK) +
(addend.value[y+addend.offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
// Add remainder of the longer number
while(x > 0) {
x--;
if (carry == 0 && result == value && rstart == (x + offset))
return;
sum = (value[x+offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
while(y > 0) {
y--;
sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
if (carry > 0) { // Result must grow in length
resultLen++;
if (result.length < resultLen) {
int temp[] = new int[resultLen];
// Result one word longer from carry-out; copy low-order
// bits into new result.
System.arraycopy(result, 0, temp, 1, result.length);
temp[0] = 1;
result = temp;
} else {
result[rstart--] = 1;
}
}
value = result;
intLen = resultLen;
offset = result.length - resultLen;
}
/**
* Adds the value of {@code addend} shifted {@code n} ints to the left.
* Has the same effect as {@code addend.leftShift(32*ints); add(addend);}
* but doesn't change the value of {@code addend}.
*/
void addShifted(MutableBigInteger addend, int n) {
if (addend.isZero()) {
return;
}
int x = intLen;
int y = addend.intLen + n;
int resultLen = (intLen > y ? intLen : y);
int[] result = (value.length < resultLen ? new int[resultLen] : value);
int rstart = result.length-1;
long sum;
long carry = 0;
// Add common parts of both numbers
while (x > 0 && y > 0) {
x--; y--;
int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
sum = (value[x+offset] & LONG_MASK) +
(bval & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
// Add remainder of the longer number
while (x > 0) {
x--;
if (carry == 0 && result == value && rstart == (x + offset)) {
return;
}
sum = (value[x+offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
while (y > 0) {
y--;
int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
sum = (bval & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
if (carry > 0) { // Result must grow in length
resultLen++;
if (result.length < resultLen) {
int temp[] = new int[resultLen];
// Result one word longer from carry-out; copy low-order
// bits into new result.
System.arraycopy(result, 0, temp, 1, result.length);
temp[0] = 1;
result = temp;
} else {
result[rstart--] = 1;
}
}
value = result;
intLen = resultLen;
offset = result.length - resultLen;
}
/**
* Like {@link #addShifted(MutableBigInteger, int)} but {@code this.intLen} must
* not be greater than {@code n}. In other words, concatenates {@code this}
* and {@code addend}.
*/
void addDisjoint(MutableBigInteger addend, int n) {
if (addend.isZero())
return;
int x = intLen;
int y = addend.intLen + n;
int resultLen = (intLen > y ? intLen : y);
int[] result;
if (value.length < resultLen)
result = new int[resultLen];
else {
result = value;
Arrays.fill(value, offset+intLen, value.length, 0);
}
int rstart = result.length-1;
// copy from this if needed
System.arraycopy(value, offset, result, rstart+1-x, x);
y -= x;
rstart -= x;
int len = Math.min(y, addend.value.length-addend.offset);
System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len);
// zero the gap
for (int i=rstart+1-y+len; i < rstart+1; i++)
result[i] = 0;
value = result;
intLen = resultLen;
offset = result.length - resultLen;
}
/**
* Adds the low {@code n} ints of {@code addend}.
*/
void addLower(MutableBigInteger addend, int n) {
MutableBigInteger a = new MutableBigInteger(addend);
if (a.offset + a.intLen >= n) {
a.offset = a.offset + a.intLen - n;
a.intLen = n;
}
a.normalize();
add(a);
}
/**
* Subtracts the smaller of this and b from the larger and places the
* result into this MutableBigInteger.
*/
int subtract(MutableBigInteger b) {
MutableBigInteger a = this;
int[] result = value;
int sign = a.compare(b);
if (sign == 0) {
reset();
return 0;
}
if (sign < 0) {
MutableBigInteger tmp = a;
a = b;
b = tmp;
}
int resultLen = a.intLen;
if (result.length < resultLen)
result = new int[resultLen];
long diff = 0;
int x = a.intLen;
int y = b.intLen;
int rstart = result.length - 1;
// Subtract common parts of both numbers
while (y > 0) {
x--; y--;
diff = (a.value[x+a.offset] & LONG_MASK) -
(b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
result[rstart--] = (int)diff;
}
// Subtract remainder of longer number
while (x > 0) {
x--;
diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
result[rstart--] = (int)diff;
}
value = result;
intLen = resultLen;
offset = value.length - resultLen;
normalize();
return sign;
}
/**
* Subtracts the smaller of a and b from the larger and places the result
* into the larger. Returns 1 if the answer is in a, -1 if in b, 0 if no
* operation was performed.
*/
private int difference(MutableBigInteger b) {
MutableBigInteger a = this;
int sign = a.compare(b);
if (sign == 0)
return 0;
if (sign < 0) {
MutableBigInteger tmp = a;
a = b;
b = tmp;
}
long diff = 0;
int x = a.intLen;
int y = b.intLen;
// Subtract common parts of both numbers
while (y > 0) {
x--; y--;
diff = (a.value[a.offset+ x] & LONG_MASK) -
(b.value[b.offset+ y] & LONG_MASK) - ((int)-(diff>>32));
a.value[a.offset+x] = (int)diff;
}
// Subtract remainder of longer number
while (x > 0) {
x--;
diff = (a.value[a.offset+ x] & LONG_MASK) - ((int)-(diff>>32));
a.value[a.offset+x] = (int)diff;
}
a.normalize();
return sign;
}
/**
* Multiply the contents of two MutableBigInteger objects. The result is
* placed into MutableBigInteger z. The contents of y are not changed.
*/
void multiply(MutableBigInteger y, MutableBigInteger z) {
int xLen = intLen;
int yLen = y.intLen;
int newLen = xLen + yLen;
// Put z into an appropriate state to receive product
if (z.value.length < newLen)
z.value = new int[newLen];
z.offset = 0;
z.intLen = newLen;
// The first iteration is hoisted out of the loop to avoid extra add
long carry = 0;
for (int j=yLen-1, k=yLen+xLen-1; j >= 0; j--, k--) {
long product = (y.value[j+y.offset] & LONG_MASK) *
(value[xLen-1+offset] & LONG_MASK) + carry;
z.value[k] = (int)product;
carry = product >>> 32;
}
z.value[xLen-1] = (int)carry;
// Perform the multiplication word by word
for (int i = xLen-2; i >= 0; i--) {
carry = 0;
for (int j=yLen-1, k=yLen+i; j >= 0; j--, k--) {
long product = (y.value[j+y.offset] & LONG_MASK) *
(value[i+offset] & LONG_MASK) +
(z.value[k] & LONG_MASK) + carry;
z.value[k] = (int)product;
carry = product >>> 32;
}
z.value[i] = (int)carry;
}
// Remove leading zeros from product
z.normalize();
}
/**
* Multiply the contents of this MutableBigInteger by the word y. The
* result is placed into z.
*/
void mul(int y, MutableBigInteger z) {
if (y == 1) {
z.copyValue(this);
return;
}
if (y == 0) {
z.clear();
return;
}
// Perform the multiplication word by word
long ylong = y & LONG_MASK;
int[] zval = (z.value.length < intLen+1 ? new int[intLen + 1]
: z.value);
long carry = 0;
for (int i = intLen-1; i >= 0; i--) {
long product = ylong * (value[i+offset] & LONG_MASK) + carry;
zval[i+1] = (int)product;
carry = product >>> 32;
}
if (carry == 0) {
z.offset = 1;
z.intLen = intLen;
} else {
z.offset = 0;
z.intLen = intLen + 1;
zval[0] = (int)carry;
}
z.value = zval;
}
/**
* This method is used for division of an n word dividend by a one word
* divisor. The quotient is placed into quotient. The one word divisor is
* specified by divisor.
*
* @return the remainder of the division is returned.
*
*/
int divideOneWord(int divisor, MutableBigInteger quotient) {
long divisorLong = divisor & LONG_MASK;
// Special case of one word dividend
if (intLen == 1) {
long dividendValue = value[offset] & LONG_MASK;
int q = (int) (dividendValue / divisorLong);
int r = (int) (dividendValue - q * divisorLong);
quotient.value[0] = q;
quotient.intLen = (q == 0) ? 0 : 1;
quotient.offset = 0;
return r;
}
if (quotient.value.length < intLen)
quotient.value = new int[intLen];
quotient.offset = 0;
quotient.intLen = intLen;
// Normalize the divisor
int shift = Integer.numberOfLeadingZeros(divisor);
int rem = value[offset];
long remLong = rem & LONG_MASK;
if (remLong < divisorLong) {
quotient.value[0] = 0;
} else {
quotient.value[0] = (int)(remLong / divisorLong);
rem = (int) (remLong - (quotient.value[0] * divisorLong));
remLong = rem & LONG_MASK;
}
int xlen = intLen;
while (--xlen > 0) {
long dividendEstimate = (remLong << 32) |
(value[offset + intLen - xlen] & LONG_MASK);
int q;
if (dividendEstimate >= 0) {
q = (int) (dividendEstimate / divisorLong);
rem = (int) (dividendEstimate - q * divisorLong);
} else {
long tmp = divWord(dividendEstimate, divisor);
q = (int) (tmp & LONG_MASK);
rem = (int) (tmp >>> 32);
}
quotient.value[intLen - xlen] = q;
remLong = rem & LONG_MASK;
}
quotient.normalize();
// Unnormalize
if (shift > 0)
return rem % divisor;
else
return rem;
}
/**
* Calculates the quotient of this div b and places the quotient in the
* provided MutableBigInteger objects and the remainder object is returned.
*
*/
MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
return divide(b,quotient,true);
}
MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
if (b.intLen < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD ||
intLen - b.intLen < BigInteger.BURNIKEL_ZIEGLER_OFFSET) {
return divideKnuth(b, quotient, needRemainder);
} else {
return divideAndRemainderBurnikelZiegler(b, quotient);
}
}
/**
* @see #divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
*/
MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient) {
return divideKnuth(b,quotient,true);
}
/**
* Calculates the quotient of this div b and places the quotient in the
* provided MutableBigInteger objects and the remainder object is returned.
*
* Uses Algorithm D in Knuth section 4.3.1.
* Many optimizations to that algorithm have been adapted from the Colin
* Plumb C library.
* It special cases one word divisors for speed. The content of b is not
* changed.
*
*/
MutableBigInteger divideKnuth(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
if (b.intLen == 0)
throw new ArithmeticException("BigInteger divide by zero");
// Dividend is zero
if (intLen == 0) {
quotient.intLen = quotient.offset = 0;
return needRemainder ? new MutableBigInteger() : null;
}
int cmp = compare(b);
// Dividend less than divisor
if (cmp < 0) {
quotient.intLen = quotient.offset = 0;
return needRemainder ? new MutableBigInteger(this) : null;
}
// Dividend equal to divisor
if (cmp == 0) {
quotient.value[0] = quotient.intLen = 1;
quotient.offset = 0;
return needRemainder ? new MutableBigInteger() : null;
}
quotient.clear();
// Special case one word divisor
if (b.intLen == 1) {
int r = divideOneWord(b.value[b.offset], quotient);
if(needRemainder) {
if (r == 0)
return new MutableBigInteger();
return new MutableBigInteger(r);
} else {
return null;
}
}
// Cancel common powers of two if we're above the KNUTH_POW2_* thresholds
if (intLen >= KNUTH_POW2_THRESH_LEN) {
int trailingZeroBits = Math.min(getLowestSetBit(), b.getLowestSetBit());
if (trailingZeroBits >= KNUTH_POW2_THRESH_ZEROS*32) {
MutableBigInteger a = new MutableBigInteger(this);
b = new MutableBigInteger(b);
a.rightShift(trailingZeroBits);
b.rightShift(trailingZeroBits);
MutableBigInteger r = a.divideKnuth(b, quotient);
r.leftShift(trailingZeroBits);
return r;
}
}
return divideMagnitude(b, quotient, needRemainder);
}
/**
* Computes {@code this/b} and {@code this%b} using the
* Burnikel-Ziegler algorithm.
* This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
* The parameter beta was chosen to b 232 so almost all shifts are
* multiples of 32 bits.
* {@code this} and {@code b} must be nonnegative.
* @param b the divisor
* @param quotient output parameter for {@code this/b}
* @return the remainder
*/
MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) {
int r = intLen;
int s = b.intLen;
// Clear the quotient
quotient.offset = quotient.intLen = 0;
if (r < s) {
return this;
} else {
// Unlike Knuth division, we don't check for common powers of two here because
// BZ already runs faster if both numbers contain powers of two and cancelling them has no
// additional benefit.
// step 1: let m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));
int j = (s+m-1) / m; // step 2a: j = ceil(s/m)
int n = j * m; // step 2b: block length in 32-bit units
long n32 = 32L * n; // block length in bits
int sigma = (int) Math.max(0, n32 - b.bitLength()); // step 3: sigma = max{T | (2^T)*B < beta^n}
MutableBigInteger bShifted = new MutableBigInteger(b);
bShifted.safeLeftShift(sigma); // step 4a: shift b so its length is a multiple of n
MutableBigInteger aShifted = new MutableBigInteger (this);
aShifted.safeLeftShift(sigma); // step 4b: shift a by the same amount
// step 5: t is the number of blocks needed to accommodate a plus one additional bit
int t = (int) ((aShifted.bitLength()+n32) / n32);
if (t < 2) {
t = 2;
}
// step 6: conceptually split a into blocks a[t-1], ..., a[0]
MutableBigInteger a1 = aShifted.getBlock(t-1, t, n); // the most significant block of a
// step 7: z[t-2] = [a[t-1], a[t-2]]
MutableBigInteger z = aShifted.getBlock(t-2, t, n); // the second to most significant block
z.addDisjoint(a1, n); // z[t-2]
// do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers
MutableBigInteger qi = new MutableBigInteger();
MutableBigInteger ri;
for (int i=t-2; i > 0; i--) {
// step 8a: compute (qi,ri) such that z=b*qi+ri
ri = z.divide2n1n(bShifted, qi);
// step 8b: z = [ri, a[i-1]]
z = aShifted.getBlock(i-1, t, n); // a[i-1]
z.addDisjoint(ri, n);
quotient.addShifted(qi, i*n); // update q (part of step 9)
}
// final iteration of step 8: do the loop one more time for i=0 but leave z unchanged
ri = z.divide2n1n(bShifted, qi);
quotient.add(qi);
ri.rightShift(sigma); // step 9: a and b were shifted, so shift back
return ri;
}
}
/**
* This method implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.
* It divides a 2n-digit number by a n-digit number.
* The parameter beta is 232 so all shifts are multiples of 32 bits.
*
* {@code this} must be a nonnegative number such that {@code this.bitLength() <= 2*b.bitLength()}
* @param b a positive number such that {@code b.bitLength()} is even
* @param quotient output parameter for {@code this/b}
* @return {@code this%b}
*/
private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
int n = b.intLen;
// step 1: base case
if (n%2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
return divideKnuth(b, quotient);
}
// step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less
MutableBigInteger aUpper = new MutableBigInteger(this);
aUpper.safeRightShift(32*(n/2)); // aUpper = [a1,a2,a3]
keepLower(n/2); // this = a4
// step 3: q1=aUpper/b, r1=aUpper%b
MutableBigInteger q1 = new MutableBigInteger();
MutableBigInteger r1 = aUpper.divide3n2n(b, q1);
// step 4: quotient=[r1,this]/b, r2=[r1,this]%b
addDisjoint(r1, n/2); // this = [r1,this]
MutableBigInteger r2 = divide3n2n(b, quotient);
// step 5: let quotient=[q1,quotient] and return r2
quotient.addDisjoint(q1, n/2);
return r2;
}
/**
* This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.
* It divides a 3n-digit number by a 2n-digit number.
* The parameter beta is 232 so all shifts are multiples of 32 bits.
*
* {@code this} must be a nonnegative number such that {@code 2*this.bitLength() <= 3*b.bitLength()}
* @param quotient output parameter for {@code this/b}
* @return {@code this%b}
*/
private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) {
int n = b.intLen / 2; // half the length of b in ints
// step 1: view this as [a1,a2,a3] where each ai is n ints or less; let a12=[a1,a2]
MutableBigInteger a12 = new MutableBigInteger(this);
a12.safeRightShift(32*n);
// step 2: view b as [b1,b2] where each bi is n ints or less
MutableBigInteger b1 = new MutableBigInteger(b);
b1.safeRightShift(n * 32);
BigInteger b2 = b.getLower(n);
MutableBigInteger r;
MutableBigInteger d;
if (compareShifted(b, n) < 0) {
// step 3a: if a1=b1, let quotient=beta^n-1 and r=a12-b1*2^n+b1
quotient.ones(n);
a12.add(b1);
b1.leftShift(32*n);
a12.subtract(b1);
r = a12;
// step 4: d=quotient*b2=(b2 << 32*n) - b2
d = new MutableBigInteger(b2);
d.leftShift(32 * n);
d.subtract(new MutableBigInteger(b2));
}
// step 5: r = r*beta^n + a3 - d (paper says a4)
// However, don't subtract d until after the while loop so r doesn't become negative
r.leftShift(32 * n);
r.addLower(this, n);
// step 6: add b until r>=d
while (r.compare(d) < 0) {
r.add(b);
quotient.subtract(MutableBigInteger.ONE);
}
r.subtract(d);
return r;
}
/**
* Returns a {@code MutableBigInteger} containing {@code blockLength} ints from
* {@code this} number, starting at {@code index*blockLength}.
* Used by Burnikel-Ziegler division.
* @param index the block index
* @param numBlocks the total number of blocks in {@code this} number
* @param blockLength length of one block in units of 32 bits
* @return
*/
private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) {
int blockStart = index * blockLength;
if (blockStart >= intLen) {
return new MutableBigInteger();
}
int blockEnd;
if (index == numBlocks-1) {
blockEnd = intLen;
} else {
blockEnd = (index+1) * blockLength;
}
if (blockEnd > intLen) {
return new MutableBigInteger();
}
int[] newVal = Arrays.copyOfRange(value, offset+intLen-blockEnd, offset+intLen-blockStart);
return new MutableBigInteger(newVal);
}
/** @see BigInteger#bitLength() */
long bitLength() {
if (intLen == 0)
return 0;
return intLen*32L - Integer.numberOfLeadingZeros(value[offset]);
}
/**
* Internally used to calculate the quotient of this div v and places the
* quotient in the provided MutableBigInteger object and the remainder is
* returned.
*
* @return the remainder of the division will be returned.
*/
long divide(long v, MutableBigInteger quotient) {
if (v == 0)
throw new ArithmeticException("BigInteger divide by zero");
// Dividend is zero
if (intLen == 0) {
quotient.intLen = quotient.offset = 0;
return 0;
}
if (v < 0)
v = -v;
int d = (int)(v >>> 32);
quotient.clear();
// Special case on word divisor
if (d == 0)
return divideOneWord((int)v, quotient) & LONG_MASK;
else {
return divideLongMagnitude(v, quotient).toLong();
}
}
private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) {
int n2 = 32 - shift;
int c=src[srcFrom];
for (int i=0; i < srcLen-1; i++) {
int b = c;
c = src[++srcFrom];
dst[dstFrom+i] = (b << shift) | (c >>> n2);
}
dst[dstFrom+srcLen-1] = c << shift;
}
/**
* Divide this MutableBigInteger by the divisor.
* The quotient will be placed into the provided quotient object &
* the remainder object is returned.
*/
private MutableBigInteger divideMagnitude(MutableBigInteger div,
MutableBigInteger quotient,
boolean needRemainder ) {
// assert div.intLen > 1
// D1 normalize the divisor
int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
// Copy divisor value to protect divisor
final int dlen = div.intLen;
int[] divisor;
MutableBigInteger rem; // Remainder starts as dividend with space for a leading zero
if (shift > 0) {
divisor = new int[dlen];
copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
if (Integer.numberOfLeadingZeros(value[offset]) >= shift) {
int[] remarr = new int[intLen + 1];
rem = new MutableBigInteger(remarr);
rem.intLen = intLen;
rem.offset = 1;
copyAndShift(value,offset,intLen,remarr,1,shift);
} else {
int[] remarr = new int[intLen + 2];
rem = new MutableBigInteger(remarr);
rem.intLen = intLen+1;
rem.offset = 1;
int rFrom = offset;
int c=0;
int n2 = 32 - shift;
for (int i=1; i < intLen+1; i++,rFrom++) {
int b = c;
c = value[rFrom];
remarr[i] = (b << shift) | (c >>> n2);
}
remarr[intLen+1] = c << shift;
}
} else {
divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
rem = new MutableBigInteger(new int[intLen + 1]);
System.arraycopy(value, offset, rem.value, 1, intLen);
rem.intLen = intLen;
rem.offset = 1;
}
int nlen = rem.intLen;
// Set the quotient size
final int limit = nlen - dlen + 1;
if (quotient.value.length < limit) {
quotient.value = new int[limit];
quotient.offset = 0;
}
quotient.intLen = limit;
int[] q = quotient.value;
// Must insert leading 0 in rem if its length did not change
if (rem.intLen == nlen) {
rem.offset = 0;
rem.value[0] = 0;
rem.intLen++;
}
int dh = divisor[0];
long dhLong = dh & LONG_MASK;
int dl = divisor[1];
// D2 Initialize j
for (int j=0; j < limit-1; j++) {
// D3 Calculate qhat
// estimate qhat
int qhat = 0;
int qrem = 0;
boolean skipCorrection = false;
int nh = rem.value[j+rem.offset];
int nh2 = nh + 0x80000000;
int nm = rem.value[j+1+rem.offset];
if (nh == dh) {
qhat = ~0;
qrem = nh + nm;
skipCorrection = qrem + 0x80000000 < nh2;
} else {
long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
if (nChunk >= 0) {
qhat = (int) (nChunk / dhLong);
qrem = (int) (nChunk - (qhat * dhLong));
} else {
long tmp = divWord(nChunk, dh);
qhat = (int) (tmp & LONG_MASK);
qrem = (int) (tmp >>> 32);
}
}
if (qhat == 0)
continue;
if (!skipCorrection) { // Correct qhat
long nl = rem.value[j+2+rem.offset] & LONG_MASK;
long rs = ((qrem & LONG_MASK) << 32) | nl;
long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
if (unsignedLongCompare(estProduct, rs)) {
qhat--;
qrem = (int)((qrem & LONG_MASK) + dhLong);
if ((qrem & LONG_MASK) >= dhLong) {
estProduct -= (dl & LONG_MASK);
rs = ((qrem & LONG_MASK) << 32) | nl;
if (unsignedLongCompare(estProduct, rs))
qhat--;
}
}
}
// D4 Multiply and subtract
rem.value[j+rem.offset] = 0;
int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
// D5 Test remainder
if (borrow + 0x80000000 > nh2) {
// D6 Add back
divadd(divisor, rem.value, j+1+rem.offset);
qhat--;
}
// Store the quotient digit
q[j] = qhat;
} // D7 loop on j
// D3 Calculate qhat
// estimate qhat
int qhat = 0;
int qrem = 0;
boolean skipCorrection = false;
int nh = rem.value[limit - 1 + rem.offset];
int nh2 = nh + 0x80000000;
int nm = rem.value[limit + rem.offset];
if (nh == dh) {
qhat = ~0;
qrem = nh + nm;
skipCorrection = qrem + 0x80000000 < nh2;
} else {
long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
if (nChunk >= 0) {
qhat = (int) (nChunk / dhLong);
qrem = (int) (nChunk - (qhat * dhLong));
} else {
long tmp = divWord(nChunk, dh);
qhat = (int) (tmp & LONG_MASK);
qrem = (int) (tmp >>> 32);
}
}
if (qhat != 0) {
if (!skipCorrection) { // Correct qhat
long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
long rs = ((qrem & LONG_MASK) << 32) | nl;
long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
if (unsignedLongCompare(estProduct, rs)) {
qhat--;
qrem = (int) ((qrem & LONG_MASK) + dhLong);
if ((qrem & LONG_MASK) >= dhLong) {
estProduct -= (dl & LONG_MASK);
rs = ((qrem & LONG_MASK) << 32) | nl;
if (unsignedLongCompare(estProduct, rs))
qhat--;
}
}
}
// D4 Multiply and subtract
int borrow;
rem.value[limit - 1 + rem.offset] = 0;
if(needRemainder)
borrow = mulsub(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
else
borrow = mulsubBorrow(rem.value, divisor, qhat, dlen, limit - 1 + rem.offset);
// D5 Test remainder
if (borrow + 0x80000000 > nh2) {
// D6 Add back
if(needRemainder)
divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
qhat--;
}
// Store the quotient digit
q[(limit - 1)] = qhat;
}
if (needRemainder) {
// D8 Unnormalize
if (shift > 0)
rem.rightShift(shift);
rem.normalize();
}
quotient.normalize();
return needRemainder ? rem : null;
}
/**
* Divide this MutableBigInteger by the divisor represented by positive long
* value. The quotient will be placed into the provided quotient object &
* the remainder object is returned.
*/
private MutableBigInteger divideLongMagnitude(long ldivisor, MutableBigInteger quotient) {
// Remainder starts as dividend with space for a leading zero
MutableBigInteger rem = new MutableBigInteger(new int[intLen + 1]);
System.arraycopy(value, offset, rem.value, 1, intLen);
rem.intLen = intLen;
rem.offset = 1;
int nlen = rem.intLen;
int limit = nlen - 2 + 1;
if (quotient.value.length < limit) {
quotient.value = new int[limit];
quotient.offset = 0;
}
quotient.intLen = limit;
int[] q = quotient.value;
// D1 normalize the divisor
int shift = Long.numberOfLeadingZeros(ldivisor);
if (shift > 0) {
ldivisor<<=shift;
rem.leftShift(shift);
}
// Must insert leading 0 in rem if its length did not change
if (rem.intLen == nlen) {
rem.offset = 0;
rem.value[0] = 0;
rem.intLen++;
}
int dh = (int)(ldivisor >>> 32);
long dhLong = dh & LONG_MASK;
int dl = (int)(ldivisor & LONG_MASK);
// D2 Initialize j
for (int j = 0; j < limit; j++) {
// D3 Calculate qhat
// estimate qhat
int qhat = 0;
int qrem = 0;
boolean skipCorrection = false;
int nh = rem.value[j + rem.offset];
int nh2 = nh + 0x80000000;
int nm = rem.value[j + 1 + rem.offset];
if (nh == dh) {
qhat = ~0;
qrem = nh + nm;
skipCorrection = qrem + 0x80000000 < nh2;
} else {
long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
if (nChunk >= 0) {
qhat = (int) (nChunk / dhLong);
qrem = (int) (nChunk - (qhat * dhLong));
} else {
long tmp = divWord(nChunk, dh);
qhat =(int)(tmp & LONG_MASK);
qrem = (int)(tmp>>>32);
}
}
if (qhat == 0)
continue;
if (!skipCorrection) { // Correct qhat
long nl = rem.value[j + 2 + rem.offset] & LONG_MASK;
long rs = ((qrem & LONG_MASK) << 32) | nl;
long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
if (unsignedLongCompare(estProduct, rs)) {
qhat--;
qrem = (int) ((qrem & LONG_MASK) + dhLong);
if ((qrem & LONG_MASK) >= dhLong) {
estProduct -= (dl & LONG_MASK);
rs = ((qrem & LONG_MASK) << 32) | nl;
if (unsignedLongCompare(estProduct, rs))
qhat--;
}
}
}
// D4 Multiply and subtract
rem.value[j + rem.offset] = 0;
int borrow = mulsubLong(rem.value, dh, dl, qhat, j + rem.offset);
// D5 Test remainder
if (borrow + 0x80000000 > nh2) {
// D6 Add back
divaddLong(dh,dl, rem.value, j + 1 + rem.offset);
qhat--;
}
// Store the quotient digit
q[j] = qhat;
} // D7 loop on j
// D8 Unnormalize
if (shift > 0)
rem.rightShift(shift);
quotient.normalize();
rem.normalize();
return rem;
}
/**
* A primitive used for division by long.
* Specialized version of the method divadd.
* dh is a high part of the divisor, dl is a low part
*/
private int divaddLong(int dh, int dl, int[] result, int offset) {
long carry = 0;
long sum = (dl & LONG_MASK) + (result[1+offset] & LONG_MASK);
result[1+offset] = (int)sum;
sum = (dh & LONG_MASK) + (result[offset] & LONG_MASK) + carry;
result[offset] = (int)sum;
carry = sum >>> 32;
return (int)carry;
}
/**
* This method is used for division by long.
* Specialized version of the method sulsub.
* dh is a high part of the divisor, dl is a low part
*/
private int mulsubLong(int[] q, int dh, int dl, int x, int offset) {
long xLong = x & LONG_MASK;
offset += 2;
long product = (dl & LONG_MASK) * xLong;
long difference = q[offset] - product;
q[offset--] = (int)difference;
long carry = (product >>> 32)
+ (((difference & LONG_MASK) >
(((~(int)product) & LONG_MASK))) ? 1:0);
product = (dh & LONG_MASK) * xLong + carry;
difference = q[offset] - product;
q[offset--] = (int)difference;
carry = (product >>> 32)
+ (((difference & LONG_MASK) >
(((~(int)product) & LONG_MASK))) ? 1:0);
return (int)carry;
}
/**
* Compare two longs as if they were unsigned.
* Returns true iff one is bigger than two.
*/
private boolean unsignedLongCompare(long one, long two) {
return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
}
/**
* This method divides a long quantity by an int to estimate
* qhat for two multi precision numbers. It is used when
* the signed value of n is less than zero.
* Returns long value where high 32 bits contain remainder value and
* low 32 bits contain quotient value.
*/
static long divWord(long n, int d) {
long dLong = d & LONG_MASK;
long r;
long q;
if (dLong == 1) {
q = (int)n;
r = 0;
return (r << 32) | (q & LONG_MASK);
}
// Approximate the quotient and remainder
q = (n >>> 1) / (dLong >>> 1);
r = n - q*dLong;
// Correct the approximation
while (r < 0) {
r += dLong;
q--;
}
while (r >= dLong) {
r -= dLong;
q++;
}
// n - q*dlong == r && 0 <= r not changed. The value of {@code this} is assumed
* to be non-negative.
*
* @implNote The implementation is based on the material in Henry S. Warren,
* Jr., Hacker's Delight (2nd ed.) (Addison Wesley, 2013), 279-282.
*
* @throws ArithmeticException if the value returned by {@code bitLength()}
* overflows the range of {@code int}.
* @return the integer square root of {@code this}
* @since 9
*/
MutableBigInteger sqrt() {
// Special cases.
if (this.isZero()) {
return new MutableBigInteger(0);
} else if (this.value.length == 1
&& (this.value[0] & LONG_MASK) < 4) { // result is unity
return ONE;
}
if (bitLength() <= 63) {
// Initial estimate is the square root of the positive long value.
long v = new BigInteger(this.value, 1).longValueExact();
long xk = (long)Math.floor(Math.sqrt(v));
// Refine the estimate.
do {
long xk1 = (xk + v/xk)/2;
// Terminate when non-decreasing.
if (xk1 >= xk) {
return new MutableBigInteger(new int[] {
(int)(xk >>> 32), (int)(xk & LONG_MASK)
});
}
xk = xk1;
} while (true);
} else {
// Set up the initial estimate of the iteration.
// Obtain the bitLength > 63.
int bitLength = (int) this.bitLength();
if (bitLength != this.bitLength()) {
throw new ArithmeticException("bitLength() integer overflow");
}
// Determine an even valued right shift into positive long range.
int shift = bitLength - 63;
if (shift % 2 == 1) {
shift++;
}
// Shift the value into positive long range.
MutableBigInteger xk = new MutableBigInteger(this);
xk.rightShift(shift);
xk.normalize();
// Use the square root of the shifted value as an approximation.
double d = new BigInteger(xk.value, 1).doubleValue();
BigInteger bi = BigInteger.valueOf((long)Math.ceil(Math.sqrt(d)));
xk = new MutableBigInteger(bi.mag);
// Shift the approximate square root back into the original range.
xk.leftShift(shift / 2);
// Refine the estimate.
MutableBigInteger xk1 = new MutableBigInteger();
do {
// xk1 = (xk + n/xk)/2
this.divide(xk, xk1, false);
xk1.add(xk);
xk1.rightShift(1);
// Terminate when non-decreasing.
if (xk1.compare(xk) >= 0) {
return xk;
}
// xk = xk1
xk.copyValue(xk1);
xk1.reset();
} while (true);
}
}
/**
* Calculate GCD of this and b. This and b are changed by the computation.
*/
MutableBigInteger hybridGCD(MutableBigInteger b) {
// Use Euclid's algorithm until the numbers are approximately the
// same length, then use the binary GCD algorithm to find the GCD.
MutableBigInteger a = this;
MutableBigInteger q = new MutableBigInteger();
while (b.intLen != 0) {
if (Math.abs(a.intLen - b.intLen) < 2)
return a.binaryGCD(b);
MutableBigInteger r = a.divide(b, q);
a = b;
b = r;
}
return a;
}
/**
* Calculate GCD of this and v.
* Assumes that this and v are not zero.
*/
private MutableBigInteger binaryGCD(MutableBigInteger v) {
// Algorithm B from Knuth section 4.5.2
MutableBigInteger u = this;
MutableBigInteger r = new MutableBigInteger();
// step B1
int s1 = u.getLowestSetBit();
int s2 = v.getLowestSetBit();
int k = (s1 < s2) ? s1 : s2;
if (k != 0) {
u.rightShift(k);
v.rightShift(k);
}
// step B2
boolean uOdd = (k == s1);
MutableBigInteger t = uOdd ? v: u;
int tsign = uOdd ? -1 : 1;
int lb;
while ((lb = t.getLowestSetBit()) >= 0) {
// steps B3 and B4
t.rightShift(lb);
// step B5
if (tsign > 0)
u = t;
else
v = t;
// Special case one word numbers
if (u.intLen < 2 && v.intLen < 2) {
int x = u.value[u.offset];
int y = v.value[v.offset];
x = binaryGcd(x, y);
r.value[0] = x;
r.intLen = 1;
r.offset = 0;
if (k > 0)
r.leftShift(k);
return r;
}
// step B6
if ((tsign = u.difference(v)) == 0)
break;
t = (tsign >= 0) ? u : v;
}
if (k > 0)
u.leftShift(k);
return u;
}
/**
* Calculate GCD of a and b interpreted as unsigned integers.
*/
static int binaryGcd(int a, int b) {
if (b == 0)
return a;
if (a == 0)
return b;
// Right shift a & b till their last bits equal to 1.
int aZeros = Integer.numberOfTrailingZeros(a);
int bZeros = Integer.numberOfTrailingZeros(b);
a >>>= aZeros;
b >>>= bZeros;
int t = (aZeros < bZeros ? aZeros : bZeros);
while (a != b) {
if ((a+0x80000000) > (b+0x80000000)) { // a > b as unsigned
a -= b;
a >>>= Integer.numberOfTrailingZeros(a);
} else {
b -= a;
b >>>= Integer.numberOfTrailingZeros(b);
}
}
return a< 64)
return euclidModInverse(k);
int t = inverseMod32(value[offset+intLen-1]);
if (k < 33) {
t = (k == 32 ? t : t & ((1 << k) - 1));
return new MutableBigInteger(t);
}
long pLong = (value[offset+intLen-1] & LONG_MASK);
if (intLen > 1)
pLong |= ((long)value[offset+intLen-2] << 32);
long tLong = t & LONG_MASK;
tLong = tLong * (2 - pLong * tLong); // 1 more Newton iter step
tLong = (k == 64 ? tLong : tLong & ((1L << k) - 1));
MutableBigInteger result = new MutableBigInteger(new int[2]);
result.value[0] = (int)(tLong >>> 32);
result.value[1] = (int)tLong;
result.intLen = 2;
result.normalize();
return result;
}
/**
* Returns the multiplicative inverse of val mod 2^32. Assumes val is odd.
*/
static int inverseMod32(int val) {
// Newton's iteration!
int t = val;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
return t;
}
/**
* Returns the multiplicative inverse of val mod 2^64. Assumes val is odd.
*/
static long inverseMod64(long val) {
// Newton's iteration!
long t = val;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
t *= 2 - val*t;
assert(t * val == 1);
return t;
}
/**
* Calculate the multiplicative inverse of 2^k mod mod, where mod is odd.
*/
static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) {
// Copy the mod to protect original
return fixup(new MutableBigInteger(1), new MutableBigInteger(mod), k);
}
/**
* Calculate the multiplicative inverse of this mod mod, where mod is odd.
* This and mod are not changed by the calculation.
*
* This method implements an algorithm due to Richard Schroeppel, that uses
* the same intermediate representation as Montgomery Reduction
* ("Montgomery Form"). The algorithm is described in an unpublished
* manuscript entitled "Fast Modular Reciprocals."
*/
private MutableBigInteger modInverse(MutableBigInteger mod) {
MutableBigInteger p = new MutableBigInteger(mod);
MutableBigInteger f = new MutableBigInteger(this);
MutableBigInteger g = new MutableBigInteger(p);
SignedMutableBigInteger c = new SignedMutableBigInteger(1);
SignedMutableBigInteger d = new SignedMutableBigInteger();
MutableBigInteger temp = null;
SignedMutableBigInteger sTemp = null;
int k = 0;
// Right shift f k times until odd, left shift d k times
if (f.isEven()) {
int trailingZeros = f.getLowestSetBit();
f.rightShift(trailingZeros);
d.leftShift(trailingZeros);
k = trailingZeros;
}
// The Almost Inverse Algorithm
while (!f.isOne()) {
// If gcd(f, g) != 1, number is not invertible modulo mod
if (f.isZero())
throw new ArithmeticException("BigInteger not invertible.");
// If f < g exchange f, g and c, d
if (f.compare(g) < 0) {
temp = f; f = g; g = temp;
sTemp = d; d = c; c = sTemp;
}
// If f == g (mod 4)
if (((f.value[f.offset + f.intLen - 1] ^
g.value[g.offset + g.intLen - 1]) & 3) == 0) {
f.subtract(g);
c.signedSubtract(d);
} else { // If f != g (mod 4)
f.add(g);
c.signedAdd(d);
}
// Right shift f k times until odd, left shift d k times
int trailingZeros = f.getLowestSetBit();
f.rightShift(trailingZeros);
d.leftShift(trailingZeros);
k += trailingZeros;
}
while (c.sign < 0)
c.signedAdd(p);
return fixup(c, p, k);
}
/**
* The Fixup Algorithm
* Calculates X such that X = C * 2^(-k) (mod P)
* Assumes C> 5; i < numWords; i++) {
// V = R * c (mod 2^j)
int v = r * c.value[c.offset + c.intLen-1];
// c = c + (v * p)
p.mul(v, temp);
c.add(temp);
// c = c / 2^j
c.intLen--;
}
int numBits = k & 0x1f;
if (numBits != 0) {
// V = R * c (mod 2^j)
int v = r * c.value[c.offset + c.intLen-1];
v &= ((1<= 0)
c.subtract(p);
return c;
}
/**
* Uses the extended Euclidean algorithm to compute the modInverse of base
* mod a modulus that is a power of 2. The modulus is 2^k.
*/
MutableBigInteger euclidModInverse(int k) {
MutableBigInteger b = new MutableBigInteger(1);
b.leftShift(k);
MutableBigInteger mod = new MutableBigInteger(b);
MutableBigInteger a = new MutableBigInteger(this);
MutableBigInteger q = new MutableBigInteger();
MutableBigInteger r = b.divide(a, q);
MutableBigInteger swapper = b;
// swap b & r
b = r;
r = swapper;
MutableBigInteger t1 = new MutableBigInteger(q);
MutableBigInteger t0 = new MutableBigInteger(1);
MutableBigInteger temp = new MutableBigInteger();
while (!b.isOne()) {
r = a.divide(b, q);
if (r.intLen == 0)
throw new ArithmeticException("BigInteger not invertible.");
swapper = r;
a = swapper;
if (q.intLen == 1)
t1.mul(q.value[q.offset], temp);
else
q.multiply(t1, temp);
swapper = q;
q = temp;
temp = swapper;
t0.add(q);
if (a.isOne())
return t0;
r = b.divide(a, q);
if (r.intLen == 0)
throw new ArithmeticException("BigInteger not invertible.");
swapper = b;
b = r;
if (q.intLen == 1)
t0.mul(q.value[q.offset], temp);
else
q.multiply(t0, temp);
swapper = q; q = temp; temp = swapper;
t1.add(q);
}
mod.subtract(t1);
return mod;
}
}