BigInteger patch for efficient multiplication and division (bug #4837946) (original) (raw)

--- jdk8/jdk/src/share/classes/java/math/BigInteger.java 2012-12-28 20:43:25.314218154 +0100

+++ src/main/java/java/math/BigInteger.java 2012-12-28 14:59:56.049558654 +0100

@@ -94,6 +94,9 @@

* @see BigDecimal

* @author Josh Bloch

* @author Michael McCloskey

• * @author Alan Eliasen

• * @author Timothy Buktu

* @since JDK1.1

*/

@@ -174,6 +177,30 @@

*/

final static long LONG_MASK = 0xffffffffL;

• /**

• * The threshold value for using Karatsuba multiplication. If the number

• * of ints in both mag arrays are greater than this number, then

• * Karatsuba multiplication will be used. This value is found

• * experimentally to work well.

• */

• private static final int KARATSUBA_THRESHOLD = 50;

• /**

• * The threshold value for using 3-way Toom-Cook multiplication.

• * If the number of ints in both mag arrays are greater than this number,

• * then Toom-Cook multiplication will be used. This value is found

• * experimentally to work well.

• */

• private static final int TOOM_COOK_THRESHOLD = 75;

• /**

• * The threshold value for using Burnikel-Ziegler division. If the number

• * of ints in the number are larger than this value,

• * Burnikel-Ziegler division will be used. This value is found

• * experimentally to work well.

• */

• private static final int BURNIKEL_ZIEGLER_THRESHOLD = 50;

//Constructors

/**

@@ -978,6 +1005,19 @@

return (val[0]>0 ? new BigInteger(val, 1) : new BigInteger(val));

}

• /**

• * Returns an n-bit number all of whose bits are ones

• * @param n a non-negative number

• * @return ONE.shiftLeft(32*n).subtract(ONE)

• */

• private static BigInteger ones(int n) {

• int[] mag = new int[(n+31)/32];

• Arrays.fill(mag, -1);

• if (n%32 != 0)

• mag[0] >>>= 32 - (n%32);

• return new BigInteger(mag, 1);

• }

// Constants

/**

@@ -1290,17 +1330,28 @@

public BigInteger multiply(BigInteger val) {

if (val.signum == 0 || signum == 0)

return ZERO;

- int resultSign = signum == val.signum ? 1 : -1;

- if (val.mag.length == 1) {

- return multiplyByInt(mag,val.mag[0], resultSign);

- }

- if(mag.length == 1) {

- return multiplyByInt(val.mag,mag[0], resultSign);

- }

- int[] result = multiplyToLen(mag, mag.length,

- val.mag, val.mag.length, null);

- result = trustedStripLeadingZeroInts(result);

- return new BigInteger(result, resultSign);

• int xlen = mag.length;

• int ylen = val.mag.length;

• if ((xlen < KARATSUBA_THRESHOLD) || (ylen < KARATSUBA_THRESHOLD))

• {

• int resultSign = signum == val.signum ? 1 : -1;

• if (val.mag.length == 1) {

• return multiplyByInt(mag,val.mag[0], resultSign);

• }

• if(mag.length == 1) {

• return multiplyByInt(val.mag,mag[0], resultSign);

• }

• int[] result = multiplyToLen(mag, xlen,

• val.mag, ylen, null);

• result = trustedStripLeadingZeroInts(result);

• return new BigInteger(result, resultSign);

• }

• else

• if ((xlen < TOOM_COOK_THRESHOLD) && (ylen < TOOM_COOK_THRESHOLD))

• return multiplyKaratsuba(this, val);

• else

• return multiplyToomCook3(this, val);

}

private static BigInteger multiplyByInt(int[] x, int y, int sign) {

@@ -1402,6 +1453,265 @@

}

/**

• * Multiplies two BigIntegers using the Karatsuba multiplication

• * algorithm. This is a recursive divide-and-conquer algorithm which is

• * more efficient for large numbers than what is commonly called the

• * "grade-school" algorithm used in multiplyToLen. If the numbers to be

• * multiplied have length n, the "grade-school" algorithm has an

• * asymptotic complexity of O(n^2). In contrast, the Karatsuba algorithm

• * has complexity of O(n^(log2(3))), or O(n^1.585). It achieves this

• * increased performance by doing 3 multiplies instead of 4 when

• * evaluating the product. As it has some overhead, should be used when

• * both numbers are larger than a certain threshold (found

• * experimentally).

• *

• * See: http://en.wikipedia.org/wiki/Karatsuba\_algorithm

• */

• private static BigInteger multiplyKaratsuba(BigInteger x, BigInteger y)

• {

• int xlen = x.mag.length;

• int ylen = y.mag.length;

• // The number of ints in each half of the number.

• int half = (Math.max(xlen, ylen)+1) / 2;

• // xl and yl are the lower halves of x and y respectively,

• // xh and yh are the upper halves.

• BigInteger xl = x.getLower(half);

• BigInteger xh = x.getUpper(half);

• BigInteger yl = y.getLower(half);

• BigInteger yh = y.getUpper(half);

• BigInteger p1 = xh.multiply(yh); // p1 = xh*yh

• BigInteger p2 = xl.multiply(yl); // p2 = xl*yl

• // p3=(xh+xl)*(yh+yl)

• BigInteger p3 = xh.add(xl).multiply(yh.add(yl));

• // result = p1 * 2^(32*2*half) + (p3 - p1 - p2) * 2^(32*half) + p2

• BigInteger result = p1.shiftLeft(32*half).add(p3.subtract(p1).subtract(p2)).shiftLeft(32*half).add(p2);

• if (x.signum != y.signum)

• return result.negate();

• else

• return result;

• }

• /**

• * Multiplies two BigIntegers using a 3-way Toom-Cook multiplication

• * algorithm. This is a recursive divide-and-conquer algorithm which is

• * more efficient for large numbers than what is commonly called the

• * "grade-school" algorithm used in multiplyToLen. If the numbers to be

• * multiplied have length n, the "grade-school" algorithm has an

• * asymptotic complexity of O(n^2). In contrast, 3-way Toom-Cook has a

• * complexity of about O(n^1.465). It achieves this increased asymptotic

• * performance by breaking each number into three parts and by doing 5

• * multiplies instead of 9 when evaluating the product. Due to overhead

• * (additions, shifts, and one division) in the Toom-Cook algorithm, it

• * should only be used when both numbers are larger than a certain

• * threshold (found experimentally). This threshold is generally larger

• * than that for Karatsuba multiplication, so this algorithm is generally

• * only used when numbers become significantly larger.

• *

• * The algorithm used is the "optimal" 3-way Toom-Cook algorithm outlined

• * by Marco Bodrato.

• *

• * See: http://bodrato.it/toom-cook/

• * http://bodrato.it/papers/#WAIFI2007

• *

• * "Towards Optimal Toom-Cook Multiplication for Univariate and

• * Multivariate Polynomials in Characteristic 2 and 0." by Marco BODRATO;

• * In C.Carlet and B.Sunar, Eds., "WAIFI'07 proceedings", p. 116-133,

• * LNCS #4547. Springer, Madrid, Spain, June 21-22, 2007.

• *

• */

• private static BigInteger multiplyToomCook3(BigInteger a, BigInteger b)

• {

• int alen = a.mag.length;

• int blen = b.mag.length;

• int largest = Math.max(alen, blen);

• // k is the size (in ints) of the lower-order slices.

• int k = (largest+2)/3; // Equal to ceil(largest/3)

• // r is the size (in ints) of the highest-order slice.

• int r = largest - 2*k;

• // Obtain slices of the numbers. a2 and b2 are the most significant

• // bits of the numbers a and b, and a0 and b0 the least significant.

• BigInteger a0, a1, a2, b0, b1, b2;

• a2 = a.getToomSlice(k, r, 0, largest);

• a1 = a.getToomSlice(k, r, 1, largest);

• a0 = a.getToomSlice(k, r, 2, largest);

• b2 = b.getToomSlice(k, r, 0, largest);

• b1 = b.getToomSlice(k, r, 1, largest);

• b0 = b.getToomSlice(k, r, 2, largest);

• BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1, db1;

• v0 = a0.multiply(b0);

• da1 = a2.add(a0);

• db1 = b2.add(b0);

• vm1 = da1.subtract(a1).multiply(db1.subtract(b1));

• da1 = da1.add(a1);

• db1 = db1.add(b1);

• v1 = da1.multiply(db1);

• v2 = da1.add(a2).shiftLeft(1).subtract(a0).multiply(

• db1.add(b2).shiftLeft(1).subtract(b0));

• vinf = a2.multiply(b2);

• /* The algorithm requires two divisions by 2 and one by 3.

• All divisions are known to be exact, that is, they do not produce

• remainders, and all results are positive. The divisions by 2 are

• implemented as right shifts which are relatively efficient, leaving

• only an exact division by 3, which is done by a specialized

• linear-time algorithm. */

• t2 = v2.subtract(vm1).exactDivideBy3();

• tm1 = v1.subtract(vm1).shiftRight(1);

• t1 = v1.subtract(v0);

• t2 = t2.subtract(t1).shiftRight(1);

• t1 = t1.subtract(tm1).subtract(vinf);

• t2 = t2.subtract(vinf.shiftLeft(1));

• tm1 = tm1.subtract(t2);

• // Number of bits to shift left.

• int ss = k*32;

• BigInteger result = vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);

• if (a.signum != b.signum)

• return result.negate();

• else

• return result;

• }

• /** Returns a slice of a BigInteger for use in Toom-Cook multiplication.

• @param lowerSize The size of the lower-order bit slices.

• @param upperSize The size of the higher-order bit slices.

• @param slice The index of which slice is requested, which must be a

• number from 0 to size-1. Slice 0 is the highest-order bits,

• and slice size-1 are the lowest-order bits.

• Slice 0 may be of different size than the other slices.

• @param fullsize The size of the larger integer array, used to align

• slices to the appropriate position when multiplying different-sized

• numbers.

• */

• private BigInteger getToomSlice(int lowerSize, int upperSize, int slice,

• int fullsize)

• {

• int start, end, sliceSize, len, offset;

• len = mag.length;

• offset = fullsize - len;

• if (slice == 0)

• {

• start = 0 - offset;

• end = upperSize - 1 - offset;

• }

• else

• {

• start = upperSize + (slice-1)*lowerSize - offset;

• end = start + lowerSize - 1;

• }

• if (start < 0)

• start = 0;

• if (end < 0)

• return ZERO;

• sliceSize = (end-start) + 1;

• if (sliceSize <= 0)

• return ZERO;

• // While performing Toom-Cook, all slices are positive and

• // the sign is adjusted when the final number is composed.

• if (start==0 && sliceSize >= len)

• return this.abs();

• int intSlice[] = new int[sliceSize];

• System.arraycopy(mag, start, intSlice, 0, sliceSize);

• return new BigInteger(trustedStripLeadingZeroInts(intSlice), 1);

• }

• /** Does an exact division (that is, the remainder is known to be zero)

• of the specified number by 3. This is used in Toom-Cook

• multiplication. This is an efficient algorithm that runs in linear

• time. If the argument is not exactly divisible by 3, results are

• undefined. Note that this is expected to be called with positive

• arguments only. */

• private BigInteger exactDivideBy3()

• {

• int len = mag.length;

• int[] result = new int[len];

• long x, w, q, borrow;

• borrow = 0L;

• for (int i=len-1; i>=0; i--)

• {

• x = (mag[i] & LONG_MASK);

• w = x - borrow;

• if (borrow > x) // Did we make the number go negative?

• borrow = 1L;

• else

• borrow = 0L;

• // 0xAAAAAAAB is the modular inverse of 3 (mod 2^32). Thus,

• // the effect of this is to divide by 3 (mod 2^32).

• // This is much faster than division on most architectures.

• q = (w * 0xAAAAAAABL) & LONG_MASK;

• result[i] = (int) q;

• // Now check the borrow. The second check can of course be

• // eliminated if the first fails.

• if (q >= 0x55555556L)

• {

• borrow++;

• if (q >= 0xAAAAAAABL)

• borrow++;

• }

• }

• result = trustedStripLeadingZeroInts(result);

• return new BigInteger(result, signum);

• }

• /**

• * Returns a new BigInteger representing n lower ints of the number.

• * This is used by Karatsuba multiplication and Burnikel-Ziegler division.

• */

• private BigInteger getLower(int n) {

• int len = mag.length;

• if (len <= n)

• return this;

• int lowerInts[] = new int[n];

• System.arraycopy(mag, len-n, lowerInts, 0, n);

• return new BigInteger(trustedStripLeadingZeroInts(lowerInts), 1);

• }

• /**

• * Returns a new BigInteger representing mag.length-n upper

• * ints of the number. This is used by Karatsuba multiplication.

• */

• private BigInteger getUpper(int n) {

• int len = mag.length;

• if (len <= n)

• return ZERO;

• int upperLen = len - n;

• int upperInts[] = new int[upperLen];

• System.arraycopy(mag, 0, upperInts, 0, upperLen);

• return new BigInteger(trustedStripLeadingZeroInts(upperInts), 1);

• }

• /**

* Returns a BigInteger whose value is {@code (this²)}.

*

* @return {@code this²}

@@ -1488,6 +1798,14 @@

* @throws ArithmeticException if {@code val} is zero.

*/

public BigInteger divide(BigInteger val) {

• if (mag.length<BURNIKEL_ZIEGLER_THRESHOLD || val.mag.length<BURNIKEL_ZIEGLER_THRESHOLD)

• return divideLong(val);

• else

• return divideBurnikelZiegler(val);

• }

• /** Long division */

• private BigInteger divideLong(BigInteger val) {

MutableBigInteger q = new MutableBigInteger(),

a = new MutableBigInteger(this.mag),

b = new MutableBigInteger(val.mag);

@@ -1508,6 +1826,14 @@

* @throws ArithmeticException if {@code val} is zero.

*/

public BigInteger[] divideAndRemainder(BigInteger val) {

• if (mag.length<BURNIKEL_ZIEGLER_THRESHOLD || val.mag.length<BURNIKEL_ZIEGLER_THRESHOLD)

• return divideAndRemainderLong(val);

• else

• return divideAndRemainderBurnikelZiegler(val);

• }

• /** Long division */

• private BigInteger[] divideAndRemainderLong(BigInteger val) {

BigInteger[] result = new BigInteger[2];

MutableBigInteger q = new MutableBigInteger(),

a = new MutableBigInteger(this.mag),

@@ -1527,6 +1853,14 @@

* @throws ArithmeticException if {@code val} is zero.

*/

public BigInteger remainder(BigInteger val) {

• if (mag.length<BURNIKEL_ZIEGLER_THRESHOLD || val.mag.length<BURNIKEL_ZIEGLER_THRESHOLD)

• return remainderLong(val);

• else

• return remainderBurnikelZiegler(val);

• }

• /** Long division */

• private BigInteger remainderLong(BigInteger val) {

MutableBigInteger q = new MutableBigInteger(),

a = new MutableBigInteger(this.mag),

b = new MutableBigInteger(val.mag);

@@ -1535,6 +1869,207 @@

}

/**

• * Calculates this / val using the Burnikel-Ziegler algorithm.

• * @param val the divisor

• * @return this / val

• */

• private BigInteger divideBurnikelZiegler(BigInteger val) {

• return divideAndRemainderBurnikelZiegler(val)[0];

• }

• /**

• * Calculates this % val using the Burnikel-Ziegler algorithm.

• * @param val the divisor

• * @return this % val

• */

• private BigInteger remainderBurnikelZiegler(BigInteger val) {

• return divideAndRemainderBurnikelZiegler(val)[1];

• }

• /**

• * Computes this / val and this % val using the

• * Burnikel-Ziegler algorithm.

• * @param val the divisor

• * @return an array containing the quotient and remainder

• */

• private BigInteger[] divideAndRemainderBurnikelZiegler(BigInteger val) {

• BigInteger[] c = divideAndRemainderBurnikelZieglerPositive(abs(), val.abs());

• // fix signs

• if (signum*val.signum < 0)

• c[0] = c[0].negate();

• if (signum < 0)

• c[1] = c[1].negate();

• return c;

• }

• /**

• * Computes a/b and a%b using the

• * Burnikel-Ziegler algorithm.

• * This method implements algorithm 3 from pg. 9 of the Burnikel-Ziegler paper.

• * The parameter β is 2^32 so all shifts are multiples of 32 bits.
  

• * a and b must be nonnegative.

• * @param a the dividend

• * @param b the divisor

• * @return an array containing the quotient and remainder

• */

• private BigInteger[] divideAndRemainderBurnikelZieglerPositive(BigInteger a, BigInteger b) {

• int r = (a.bitLength()+31) / 32;

• int s = (b.bitLength()+31) / 32;

• if (r < s)

• return new BigInteger[] {ZERO, a};

• else {

• // 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; // j = ceil(s/m)

• int n = j * m; // block length in 32-bit units

• int n32 = 32 * n; // block length in bits

• int sigma = Math.max(0, n32 - b.bitLength());

• b = b.shiftLeft(sigma); // shift b so its length is a multiple of n

• a = a.shiftLeft(sigma); // shift a by the same amount

• // t is the number of blocks needed to accommodate 'a' plus one additional bit

• int t = (a.bitLength()+n32) / n32;

• if (t < 2)

• t = 2;

• BigInteger a1 = a.getBlock(t-1, t, n); // the most significant block of a

• BigInteger a2 = a.getBlock(t-2, t, n); // the second to most significant block

• // do schoolbook division on blocks, dividing 2-block numbers by 1-block numbers

• BigInteger z = a1.shiftLeftInts(n).add(a2); // Z[t-2]

• BigInteger quotient = ZERO;

• BigInteger[] c;

• for (int i=t-2; i>0; i--) {

• c = burnikel21(z, b);

• z = a.getBlock(i-1, t, n);

• z = z.add(c[1].shiftLeftInts(n));

• quotient = quotient.add(c[0]).shiftLeftInts(n);

• }

• // do the loop one more time for i=0 but leave z unchanged

• c = burnikel21(z, b);

• quotient = quotient.add(c[0]);

• BigInteger remainder = c[1].shiftRight(sigma); // a and b were shifted, so shift back

• return new BigInteger[] {quotient, remainder};

• }

• }

• /**

• * Returns a BigInteger containing blockLength ints from

• * this number, starting at index*blockLength.
  

• * Used by Burnikel-Ziegler division.

• * @param index the block index

• * @param numBlocks the total number of blocks in this number

• * @param blockLength length of one block in units of 32 bits

• * @return

• */

• private BigInteger getBlock(int index, int numBlocks, int blockLength) {

• int blockStart = index * blockLength;

• if (blockStart >= mag.length)

• return ZERO;

• int blockEnd;

• if (index == numBlocks-1)

• blockEnd = (bitLength()+31) / 32;

• else

• blockEnd = (index+1) * blockLength;

• if (blockEnd > mag.length)

• return ZERO;

• int[] newMag = trustedStripLeadingZeroInts(Arrays.copyOfRange(mag, mag.length-blockEnd, mag.length-blockStart));

• return new BigInteger(newMag, signum);

• }

• /**

• * Implements algorithm 1 from pg. 4 of the Burnikel-Ziegler paper.

• * The parameter β is 2^32 so all shifts are multiples of 32 bits.

• * @param a a nonnegative number such that a.bitLength() <= 2*b.bitLength()

• * @param b a positive number such that b.bitLength() is even

• * @return a/b and a%b

• */

• private BigInteger[] burnikel21(BigInteger a, BigInteger b) {

• int n = (b.bitLength()+31) / 32;

• if (n%2!=0 || n<BURNIKEL_ZIEGLER_THRESHOLD)

• return a.divideAndRemainderLong(b);

• // view a as [a1,a2,a3,a4] and divide [a1,a2,a3] by b

• BigInteger[] c1 = burnikel32(a.shiftRightInts(n/2), b);

• // divide the concatenation of c1[1] and a4 by b

• BigInteger a4 = a.getLower(n/2);

• BigInteger[] c2 = burnikel32(c1[1].shiftLeftInts(n/2).add(a4), b);

• // quotient = the concatentation of the two above quotients

• return new BigInteger[] {c1[0].shiftLeftInts(n/2).add(c2[0]), c2[1]};

• }

• /**

• * Implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.

• * The parameter β is 2^32 so all shifts are multiples of 32 bits.
  

• * @param a a nonnegative number such that 2*a.bitLength() <= 3*b.bitLength()

• * @param b a positive number such that b.bitLength() is even

• * @return a/b and a%b

• */

• private BigInteger[] burnikel32(BigInteger a, BigInteger b) {

• int n = (b.bitLength()+63) / 64; // half the length of b in ints

• // split a in 3 parts of length n or less

• BigInteger a1 = a.shiftRightInts(2*n);

• BigInteger a2 = a.shiftAndTruncate(n);

• BigInteger a3 = a.getLower(n);

• // split a in 2 parts of length n or less

• BigInteger b1 = b.shiftRightInts(n);

• BigInteger b2 = b.getLower(n);

• BigInteger q, r1;

• BigInteger a12 = a1.shiftLeftInts(n).add(a2); // concatenation of a1 and a2

• if (a1.compareTo(b1) < 0) {

• // q=a12/b1, r=a12%b1

• BigInteger[] c = burnikel21(a12, b1);

• q = c[0];

• r1 = c[1];

• }

• else {

• // q=β^n-1, r=a12-b1*2^n+b1

• q = ones(32*n);

• r1 = a12.subtract(b1.shiftLeftInts(n)).add(b1);

• }

• BigInteger d = q.multiply(b2);

• BigInteger r = r1.shiftLeftInts(n).add(a3).subtract(d); // r = r1*β^n + a3 - d (paper says a4)

• // add b until r>=0

• while (r.signum() < 0) {

• r = r.add(b);

• q = q.subtract(ONE);

• }

• return new BigInteger[] {q, r};

• }

• /**

• * Returns a number equal to this.shiftRightInts(n).getLower(n).
  

• * Used by Burnikel-Ziegler division.

• * @param n a non-negative number

• * @return n bits of this starting at bit n

• */

• private BigInteger shiftAndTruncate(int n) {

• if (mag.length <= n)

• return ZERO;

• if (mag.length <= 2*n) {

• int[] newMag = trustedStripLeadingZeroInts(Arrays.copyOfRange(mag, 0, mag.length-n));

• return new BigInteger(newMag, signum);

• }

• else {

• int[] newMag = trustedStripLeadingZeroInts(Arrays.copyOfRange(mag, mag.length-2*n, mag.length-n));

• return new BigInteger(newMag, signum);

• }

• }

• /**

* Returns a BigInteger whose value is (this^exponent).

* Note that {@code exponent} is an integer rather than a BigInteger.

*

@@ -2323,6 +2858,28 @@

return val;

}

• /**

• * Shifts a number to the left by a multiple of 32.

• * @param n a non-negative number

• * @return this.shiftLeft(32*n)

• */

• private BigInteger shiftLeftInts(int n) {

• int[] newMag = trustedStripLeadingZeroInts(Arrays.copyOf(mag, mag.length+n));

• return new BigInteger(newMag, signum);

• }

• /**

• * Shifts a number to the right by a multiple of 32.

• * @param n a non-negative number

• * @return this.shiftRight(32*n)

• */

• private BigInteger shiftRightInts(int n) {

• if (n >= mag.length)

• return ZERO;

• else

• return new BigInteger(Arrays.copyOf(mag, mag.length-n), signum);

• }

// Bitwise Operations

/**