# HG changeset patch # User bpb # Date 1374883399 25200 # Node ID 4699ce24107442e95f7d634fc4166e438d04066a # Parent 9f9ffe6be557392243dc98095c34cc5ae8c74618 8014319: Faster division of large integers Summary: Implement Burnickel-Ziegler division algorithm in BigInteger Reviewed-by: bpb, martin Contributed-by: Tim Buktu diff --git a/src/share/classes/java/math/BigInteger.java b/src/share/classes/java/math/BigInteger.java --- a/src/share/classes/java/math/BigInteger.java +++ b/src/share/classes/java/math/BigInteger.java @@ -33,7 +33,6 @@ import java.io.ObjectInputStream; import java.io.ObjectOutputStream; import java.io.ObjectStreamField; -import java.util.ArrayList; import java.util.Arrays; import java.util.Random; import sun.misc.DoubleConsts; @@ -101,6 +100,7 @@ * @author Josh Bloch * @author Michael McCloskey * @author Alan Eliasen + * @author Timothy Buktu * @since JDK1.1 */ @@ -215,6 +215,14 @@ private static final int TOOM_COOK_SQUARE_THRESHOLD = 140; /** + * 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. + */ + static final int BURNIKEL_ZIEGLER_THRESHOLD = 50; + + /** * The threshold value for using Schoenhage recursive base conversion. If * the number of ints in the number are larger than this value, * the Schoenhage algorithm will be used. In practice, it appears that the @@ -1781,7 +1789,7 @@ if (len < TOOM_COOK_SQUARE_THRESHOLD) return squareKaratsuba(); else - return squareToomCook3(); + return squareToomCook3(); } /** @@ -1936,11 +1944,26 @@ * @throws ArithmeticException if {@code val} is zero. */ public BigInteger divide(BigInteger val) { + if (mag.lengthsb. This implements the recursive Schoenhage algorithm + * {@code sb}. This implements the recursive Schoenhage algorithm * for base conversions. *

* See Knuth, Donald, _The Art of Computer Programming_, Vol. 2, @@ -3450,7 +3521,7 @@ * If this value doesn't already exist in the cache, it is added. *

* This could be changed to a more complicated caching method using - * Future. + * {@code Future}. */ private static BigInteger getRadixConversionCache(int radix, int exponent) { BigInteger[] cacheLine = powerCache[radix]; // volatile read diff --git a/src/share/classes/java/math/MutableBigInteger.java b/src/share/classes/java/math/MutableBigInteger.java --- a/src/share/classes/java/math/MutableBigInteger.java +++ b/src/share/classes/java/math/MutableBigInteger.java @@ -1,5 +1,5 @@ /* - * Copyright (c) 1999, 2011, Oracle and/or its affiliates. All rights reserved. + * Copyright (c) 1999, 2013, 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 @@ -38,14 +38,14 @@ * * @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; -import static java.math.BigInteger.LONG_MASK; -import static java.math.BigDecimal.INFLATED; - class MutableBigInteger { /** * Holds the magnitude of this MutableBigInteger in big endian order. @@ -75,6 +75,24 @@ */ 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 /** @@ -124,6 +142,20 @@ } /** + * 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. */ @@ -155,6 +187,14 @@ } /** + * 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. */ @@ -238,6 +278,32 @@ } /** + * 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 @@ -454,6 +520,16 @@ } /** + * 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. */ @@ -475,6 +551,14 @@ } /** + * 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) { @@ -615,6 +699,35 @@ } /** + * 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. @@ -673,6 +786,121 @@ 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(b);} + * but doesn't change the value of {@code b}. + */ + 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>> 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>> 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= 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 @@ -910,6 +1138,29 @@ * 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 (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; + + 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 + int n32 = 32 * n; // block length in bits + int sigma = 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 + safeLeftShift(sigma); // step 4b: shift this by the same amount + + // step 5: t is the number of blocks needed to accommodate this plus one additional bit + int t = (bitLength()+n32) / n32; + if (t < 2) + t = 2; + + // step 6: conceptually split this into blocks a[t-1], ..., a[0] + MutableBigInteger a1 = getBlock(t-1, t, n); // the most significant block of this + + // step 7: z[t-2] = [a[t-1], a[t-2]] + MutableBigInteger z = 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; + quotient.offset = quotient.intLen = 0; + 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 = 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: this 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 + * 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() */ + int bitLength() { + if (intLen == 0) + return 0; + return intLen*32 - Integer.numberOfLeadingZeros(value[offset]); } /** @@ -1006,7 +1462,7 @@ */ private MutableBigInteger divideMagnitude(MutableBigInteger div, MutableBigInteger quotient, - boolean needReminder ) { + boolean needRemainder ) { // assert div.intLen > 1 // D1 normalize the divisor int shift = Integer.numberOfLeadingZeros(div.value[div.offset]); @@ -1176,7 +1632,7 @@ // D4 Multiply and subtract int borrow; rem.value[limit - 1 + rem.offset] = 0; - if(needReminder) + 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); @@ -1184,7 +1640,7 @@ // D5 Test remainder if (borrow + 0x80000000 > nh2) { // D6 Add back - if(needReminder) + if(needRemainder) divadd(divisor, rem.value, limit - 1 + 1 + rem.offset); qhat--; } @@ -1194,14 +1650,14 @@ } - if(needReminder) { + if(needRemainder) { // D8 Unnormalize if (shift > 0) rem.rightShift(shift); rem.normalize(); } quotient.normalize(); - return needReminder ? rem : null; + return needRemainder ? rem : null; } /** @@ -1367,7 +1823,7 @@ * 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 reminder value and + * Returns long value where high 32 bits contain remainder value and * low 32 bits contain quotient value. */ static long divWord(long n, int d) { @@ -1582,7 +2038,7 @@ return result; } - /* + /** * Returns the multiplicative inverse of val mod 2^32. Assumes val is odd. */ static int inverseMod32(int val) { @@ -1595,7 +2051,7 @@ return t; } - /* + /** * Calculate the multiplicative inverse of 2^k mod mod, where mod is odd. */ static MutableBigInteger modInverseBP2(MutableBigInteger mod, int k) { @@ -1665,7 +2121,7 @@ return fixup(c, p, k); } - /* + /** * The Fixup Algorithm * Calculates X such that X = C * 2^(-k) (mod P) * Assumes C

abs(v)} and {@code a > b && b > 0}, then if + * {@code w/z = q1*z + r1} and {@code u/v = q2*v + r2}, then + * {@code q1 = q2*pow(2,a-b)} and {@code r1 = r2*pow(2,b)}. The test + * ensures that {@code v} is just under the B-Z threshold and that {@code w} + * and {@code z} are both over the threshold. This implies that {@code u/v} + * uses the standard division algorithm and {@code w/z} uses the B-Z + * algorithm. The results of the two algorithms are then compared using the + * observation described in the foregoing and if they are not equal a + * failure is logged. + */ + public static void divideLarge() { + int failCount = 0; + + BigInteger base = BigInteger.ONE.shiftLeft(BITS_BURNIKEL_ZIEGLER - 33); + for (int i=0; i