Skip to content

Instantly share code, notes, and snippets.

@juanmf

juanmf/MutableBigInteger.java

Last active March 24, 2016 12:58
Show Gist options
  • Save juanmf/5384e81c4852f234b300 to your computer and use it in GitHub Desktop.
Save juanmf/5384e81c4852f234b300 to your computer and use it in GitHub Desktop.
Download ZIP
MutableBigInteger
Raw
MutableBigInteger.java
/*
* The MIT License
*
* Copyright 2016 juan.fernandez.
*
* Permission is hereby granted, free of charge, to any person obtaining a copy
* of this software and associated documentation files (the "Software"), to deal
* in the Software without restriction, including without limitation the rights
* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
* copies of the Software, and to permit persons to whom the Software is
* furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included in
* all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
* THE SOFTWARE.
*/
package ar.com.docdigital.util;
/**
* 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 java.math.BigDecimal;
import java.math.BigInteger;
import java.nio.ByteBuffer;
import java.nio.IntBuffer;
import java.util.Arrays;
public class MutableBigInteger {
/**
* The MutableBigInteger constant zero.
*
* @since 1.2
*/
public static final MutableBigInteger ZERO = new MutableBigInteger(0);
public static final MutableBigInteger UNSAFE_AUX_VALUE = new MutableBigInteger(0);
public static final MutableBigInteger UNUSED_QUOTIENT_HOLDER = new MutableBigInteger(0);
/**
* This mask is used to obtain the value of an int as if it were unsigned.
*/
final static long LONG_MASK = 0xffffffffL;
/**
* Sentinel value for {@link #intCompact} indicating the
* significand information is only available from {@code intVal}.
*/
static final long INFLATED = Long.MIN_VALUE;
/**
* The threshold value for using Burnikel-Ziegler division. If the number
* of ints in the divisor are larger than this value, Burnikel-Ziegler
* division may be used. This value is found experimentally to work well.
*/
static final int BURNIKEL_ZIEGLER_THRESHOLD = 80;
/**
* The offset value for using Burnikel-Ziegler division. If the number
* of ints in the divisor exceeds the Burnikel-Ziegler threshold, and the
* number of ints in the dividend is greater than the number of ints in the
* divisor plus this value, Burnikel-Ziegler division will be used. This
* value is found experimentally to work well.
*/
static final int BURNIKEL_ZIEGLER_OFFSET = 40;
/**
* 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;
int[] backupValue;
/**
* 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;
int backupIntLen;
/**
* 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;
// Cache of common small BigDecimal values.
private static final BigDecimal zeroThroughTen[] = {
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ONE),
new BigDecimal(BigInteger.valueOf(2)),
new BigDecimal(BigInteger.valueOf(3)),
new BigDecimal(BigInteger.valueOf(4)),
new BigDecimal(BigInteger.valueOf(5)),
new BigDecimal(BigInteger.valueOf(6)),
new BigDecimal(BigInteger.valueOf(7)),
new BigDecimal(BigInteger.valueOf(8)),
new BigDecimal(BigInteger.valueOf(9)),
new BigDecimal(BigInteger.TEN),
};
// Cache of zero scaled by 0 - 15
private static final BigDecimal[] ZERO_SCALED_BY = {
zeroThroughTen[0],
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
new BigDecimal(BigInteger.ZERO),
};
/**
* Package private method to return bit length for an integer.
*/
static int bitLengthForInt(int n) {
return 32 - Integer.numberOfLeadingZeros(n);
}
static BigDecimal zeroValueOf(int scale) {
if (scale >= 0 && scale < ZERO_SCALED_BY.length)
return ZERO_SCALED_BY[scale];
else
return new BigDecimal(BigInteger.ZERO);
}
/**
* Compares this BigInteger with the specified Object for equality.
*
* @param x Object to which this BigInteger is to be compared.
* @return {@code true} if and only if the specified Object is a
* BigInteger whose value is numerically equal to this BigInteger.
*/
public boolean equals(Object x) {
// This test is just an optimization, which may or may not help
if (x == this)
return true;
if (!(x instanceof MutableBigInteger))
return false;
MutableBigInteger xInt = (MutableBigInteger) x;
// if (xInt.signum != signum) // siempre es positivo.
// return false;
int[] m = value;
int len = m.length;
int[] xm = xInt.value;
if (len != xm.length)
return false;
for (int i = 0; i < len; i++)
if (xm[i] != m[i])
return false;
return true;
}
/**
* Returns the hash code for this BigInteger.
*
* @return hash code for this BigInteger.
*/
public int hashCode() {
int hashCode = 0;
for (int i=0; i < value.length; i++)
hashCode = (int)(31*hashCode + (value[i] & LONG_MASK));
return hashCode;
}
// Constructors
/**
* The default constructor. An empty MutableBigInteger is created with
* a one word capacity.
*/
public MutableBigInteger() {
value = new int[1];
intLen = 0;
}
/**
* Construct a new MutableBigInteger with a magnitude specified by
* the int val.
*/
public MutableBigInteger(long val) {
if (val < 0) {
throw new UnsupportedOperationException("No negatives.");
// val = -val;
// signum = -1;
}
init(val);
}
public void init(long val) {
int highWord = (int)(val >>> 32);
if (highWord == 0) {
if (null == value || value.length != 1) {
value = new int[1];
}
value[0] = (int)val;
intLen = 1;
} else {
if (null == value || value.length != 2) {
value = new int[2];
}
value[0] = highWord;
value[1] = (int)val;
intLen = 2;
}
}
/**
* 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.
*/
public MutableBigInteger(java.math.BigInteger b) {
byte[] array = getBytesForIntArray(b);
ByteBuffer bbuf = ByteBuffer.wrap(array);
IntBuffer ibuf = bbuf.asIntBuffer();
int[] intArray = new int[ibuf.limit()];
ibuf.get(intArray);
intLen = intArray.length;
value = Arrays.copyOf(intArray, intLen);
}
/**
* BigInteger.toByteArray() yields the minimum amount of bytes that represents
* the value, that's not necessary multiple of 4, which makes ByteBuffer.asIntBuffer()
* loose info for values > Integer.MAX_VALUE.
*
* @param b
* @return a byte[] possibly with leading 0 bytes.
*/
private byte[] getBytesForIntArray(java.math.BigInteger b) {
byte[] array = b.toByteArray();
if (0 != array.length % 4) {
byte[] newArray = new byte[4 * (int)Math.ceil((float)array.length / 4)];
for (int i = array.length - 1, j = newArray.length - 1; i >= 0; i--, j--) {
newArray[j] = array[i];
}
array = newArray;
}
return array;
}
/**
* Construct a new MutableBigInteger with a magnitude equal to the
* specified MutableBigInteger.
*/
public 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.
*/
java.math.BigInteger toBigInteger(int sign) {
// throw new UnsupportedOperationException();
if (intLen == 0 || sign == 0)
return java.math.BigInteger.ZERO;
int[] mag = getMagnitudeArray();
ByteBuffer newBb = ByteBuffer.allocate(mag.length*4);
newBb.asIntBuffer().put(mag);
byte[] byteArray = newBb.array();
return new java.math.BigInteger(byteArray);
}
/**
* Converts this number to a nonnegative {@code BigInteger}.
*/
public java.math.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) {
throw new UnsupportedOperationException();
// if (intLen == 0 || sign == 0)
// return 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 java.math.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 <tt>b</tt>.
*/
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() {
java.math.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 = 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 = 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 java.math.BigInteger getLower(int n) {
throw new UnsupportedOperationException();
// if (isZero()) {
// return java.math.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 java.math.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.
*/
public 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;
}
public void addAndBackup(MutableBigInteger addend) {
if (null == backupValue || backupValue.length != value.length) {
backupValue = new int[value.length];
}
for (int w = 0; w < value.length; w++) {
backupValue[w] = value[w];
}
backupIntLen = intLen;
add(addend);
}
public void restoreBackup() {
if (null == backupValue || backupValue.length != value.length) {
value = new int[backupValue.length];
}
for (int w = 0; w < value.length; w++) {
value[w] = backupValue[w];
}
// value = backupValue;
intLen = backupIntLen;
}
/**
* 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.
*/
public 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.
*
*/
public MutableBigInteger remainder(MutableBigInteger b) {
return divide(b, UNUSED_QUOTIENT_HOLDER,true);
}
/**
* Calculates the quotient of this div b and places the quotient in the
* provided MutableBigInteger objects and the remainder object is returned.
*
*/
public MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
return divide(b,quotient,true);
}
MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
if (b.intLen < BURNIKEL_ZIEGLER_THRESHOLD ||
intLen - b.intLen < 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
* <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
* This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.
* The parameter beta was chosen to b 2<sup>32</sup> so almost all shifts are
* multiples of 32 bits.<br/>
* {@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/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.<br/>
* The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.
* <br/>
* {@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 < 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.<br/>
* The parameter beta is 2<sup>32</sup> so all shifts are multiples of 32 bits.<br/>
* <br/>
* {@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);
java.math.BigInteger b2 = b.getLower(n);
MutableBigInteger r;
MutableBigInteger d;
if (compareShifted(b, n) < 0) {
// step 3a: if a1<b1, let quotient=a12/b1 and r=a12%b1
r = a12.divide2n1n(b1, quotient);
// step 4: d=quotient*b2
d = new MutableBigInteger(quotient.toBigInteger().multiply(b2));
} else {
// step 3b: 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}.<br/>
* 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 <dLong, hence we're done.
return (r << 32) | (q & LONG_MASK);
}
/**
* 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<<t;
}
/**
* Returns the modInverse of this mod p. This and p are not affected by
* the operation.
*/
MutableBigInteger mutableModInverse(MutableBigInteger p) {
// Modulus is odd, use Schroeppel's algorithm
if (p.isOdd())
return modInverse(p);
// Base and modulus are even, throw exception
if (isEven())
throw new ArithmeticException("BigInteger not invertible.");
// Get even part of modulus expressed as a power of 2
int powersOf2 = p.getLowestSetBit();
// Construct odd part of modulus
MutableBigInteger oddMod = new MutableBigInteger(p);
oddMod.rightShift(powersOf2);
if (oddMod.isOne())
return modInverseMP2(powersOf2);
// Calculate 1/a mod oddMod
MutableBigInteger oddPart = modInverse(oddMod);
// Calculate 1/a mod evenMod
MutableBigInteger evenPart = modInverseMP2(powersOf2);
// Combine the results using Chinese Remainder Theorem
MutableBigInteger y1 = modInverseBP2(oddMod, powersOf2);
MutableBigInteger y2 = oddMod.modInverseMP2(powersOf2);
MutableBigInteger temp1 = new MutableBigInteger();
MutableBigInteger temp2 = new MutableBigInteger();
MutableBigInteger result = new MutableBigInteger();
oddPart.leftShift(powersOf2);
oddPart.multiply(y1, result);
evenPart.multiply(oddMod, temp1);
temp1.multiply(y2, temp2);
result.add(temp2);
return result.divide(p, temp1);
}
/*
* Calculate the multiplicative inverse of this mod 2^k.
*/
MutableBigInteger modInverseMP2(int k) {
if (isEven())
throw new ArithmeticException("Non-invertible. (GCD != 1)");
if (k > 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;
}
/**
* 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) {
throw new UnsupportedOperationException();
// 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<P and P is odd.
*/
static MutableBigInteger fixup(MutableBigInteger c, MutableBigInteger p,
int k) {
MutableBigInteger temp = new MutableBigInteger();
// Set r to the multiplicative inverse of p mod 2^32
int r = -inverseMod32(p.value[p.offset+p.intLen-1]);
for (int i=0, numWords = k >> 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<<numBits) - 1);
// c = c + (v * p)
p.mul(v, temp);
c.add(temp);
// c = c / 2^j
c.rightShift(numBits);
}
// In theory, c may be greater than p at this point (Very rare!)
while (c.compare(p) >= 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;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment