8014319: Faster division of large integers
Summary: Implement Burnickel-Ziegler division algorithm in BigInteger
Reviewed-by: bpb, martin
Contributed-by: Tim Buktu <tbuktu@hotmail.com>
--- a/jdk/src/share/classes/java/math/BigInteger.java Fri Jul 26 17:23:20 2013 -0700
+++ b/jdk/src/share/classes/java/math/BigInteger.java Fri Jul 26 17:03:19 2013 -0700
@@ -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.length<BURNIKEL_ZIEGLER_THRESHOLD || val.mag.length<BURNIKEL_ZIEGLER_THRESHOLD)
+ return divideKnuth(val);
+ else
+ return divideBurnikelZiegler(val);
+ }
+
+ /**
+ * Returns a BigInteger whose value is {@code (this / val)} using an O(n^2) algorithm from Knuth.
+ *
+ * @param val value by which this BigInteger is to be divided.
+ * @return {@code this / val}
+ * @throws ArithmeticException if {@code val} is zero.
+ * @see MutableBigInteger#divideKnuth(MutableBigInteger, MutableBigInteger, boolean)
+ */
+ private BigInteger divideKnuth(BigInteger val) {
MutableBigInteger q = new MutableBigInteger(),
a = new MutableBigInteger(this.mag),
b = new MutableBigInteger(val.mag);
- a.divide(b, q, false);
+ a.divideKnuth(b, q, false);
return q.toBigInteger(this.signum * val.signum);
}
@@ -1956,11 +1979,19 @@
* @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 divideAndRemainderKnuth(val);
+ else
+ return divideAndRemainderBurnikelZiegler(val);
+ }
+
+ /** Long division */
+ private BigInteger[] divideAndRemainderKnuth(BigInteger val) {
BigInteger[] result = new BigInteger[2];
MutableBigInteger q = new MutableBigInteger(),
a = new MutableBigInteger(this.mag),
b = new MutableBigInteger(val.mag);
- MutableBigInteger r = a.divide(b, q);
+ MutableBigInteger r = a.divideKnuth(b, q);
result[0] = q.toBigInteger(this.signum == val.signum ? 1 : -1);
result[1] = r.toBigInteger(this.signum);
return result;
@@ -1975,11 +2006,51 @@
* @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 remainderKnuth(val);
+ else
+ return remainderBurnikelZiegler(val);
+ }
+
+ /** Long division */
+ private BigInteger remainderKnuth(BigInteger val) {
MutableBigInteger q = new MutableBigInteger(),
a = new MutableBigInteger(this.mag),
b = new MutableBigInteger(val.mag);
- return a.divide(b, q).toBigInteger(this.signum);
+ return a.divideKnuth(b, q).toBigInteger(this.signum);
+ }
+
+ /**
+ * Calculates {@code this / val} using the Burnikel-Ziegler algorithm.
+ * @param val the divisor
+ * @return {@code this / val}
+ */
+ private BigInteger divideBurnikelZiegler(BigInteger val) {
+ return divideAndRemainderBurnikelZiegler(val)[0];
+ }
+
+ /**
+ * Calculates {@code this % val} using the Burnikel-Ziegler algorithm.
+ * @param val the divisor
+ * @return {@code this % val}
+ */
+ private BigInteger remainderBurnikelZiegler(BigInteger val) {
+ return divideAndRemainderBurnikelZiegler(val)[1];
+ }
+
+ /**
+ * Computes {@code this / val} and {@code this % val} using the
+ * Burnikel-Ziegler algorithm.
+ * @param val the divisor
+ * @return an array containing the quotient and remainder
+ */
+ private BigInteger[] divideAndRemainderBurnikelZiegler(BigInteger val) {
+ MutableBigInteger q = new MutableBigInteger();
+ MutableBigInteger r = new MutableBigInteger(this).divideAndRemainderBurnikelZiegler(new MutableBigInteger(val), q);
+ BigInteger qBigInt = q.isZero() ? ZERO : q.toBigInteger(signum*val.signum);
+ BigInteger rBigInt = r.isZero() ? ZERO : r.toBigInteger(signum);
+ return new BigInteger[] {qBigInt, rBigInt};
}
/**
@@ -3399,7 +3470,7 @@
/**
* Converts the specified BigInteger to a string and appends to
- * <code>sb</code>. This implements the recursive Schoenhage algorithm
+ * {@code sb}. This implements the recursive Schoenhage algorithm
* for base conversions.
* <p/>
* 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.
* <p/>
* This could be changed to a more complicated caching method using
- * <code>Future</code>.
+ * {@code Future}.
*/
private static BigInteger getRadixConversionCache(int radix, int exponent) {
BigInteger[] cacheLine = powerCache[radix]; // volatile read
--- a/jdk/src/share/classes/java/math/MutableBigInteger.java Fri Jul 26 17:23:20 2013 -0700
+++ b/jdk/src/share/classes/java/math/MutableBigInteger.java Fri Jul 26 17:03:19 2013 -0700
@@ -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<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
@@ -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<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD || b.intLen<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD)
+ 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.
@@ -917,38 +1168,34 @@
* changed.
*
*/
- MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
- return divide(b,quotient,true);
- }
-
- MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needReminder) {
+ 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;
- return needReminder ? new MutableBigInteger() : null;
+ 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 needReminder ? new MutableBigInteger(this) : null;
+ return needRemainder ? new MutableBigInteger(this) : null;
}
// Dividend equal to divisor
if (cmp == 0) {
quotient.value[0] = quotient.intLen = 1;
quotient.offset = 0;
- return needReminder ? new MutableBigInteger() : null;
+ 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(needReminder) {
+ if(needRemainder) {
if (r == 0)
return new MutableBigInteger();
return new MutableBigInteger(r);
@@ -956,7 +1203,216 @@
return null;
}
}
- return divideMagnitude(b, quotient, needReminder);
+
+ // 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;
+
+ 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.<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<BigInteger.BURNIKEL_ZIEGLER_THRESHOLD)
+ return divideKnuth(b, quotient);
+
+ // step 2: view this as [a1,a2,a3,a4] where each ai is n/2 ints or less
+ MutableBigInteger aUpper = new MutableBigInteger(this);
+ aUpper.safeRightShift(32*(n/2)); // aUpper = [a1,a2,a3]
+ keepLower(n/2); // this = a4
+
+ // step 3: q1=aUpper/b, r1=aUpper%b
+ MutableBigInteger q1 = new MutableBigInteger();
+ MutableBigInteger r1 = aUpper.divide3n2n(b, q1);
+
+ // step 4: quotient=[r1,this]/b, r2=[r1,this]%b
+ addDisjoint(r1, n/2); // this = [r1,this]
+ MutableBigInteger r2 = divide3n2n(b, quotient);
+
+ // step 5: let quotient=[q1,quotient] and return r2
+ quotient.addDisjoint(q1, n/2);
+ return r2;
+ }
+
+ /**
+ * This method implements algorithm 2 from pg. 5 of the Burnikel-Ziegler paper.
+ * It divides a 3n-digit number by a 2n-digit number.<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);
+ 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() */
+ 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<P and P is odd.
--- a/jdk/test/java/math/BigInteger/BigIntegerTest.java Fri Jul 26 17:23:20 2013 -0700
+++ b/jdk/test/java/math/BigInteger/BigIntegerTest.java Fri Jul 26 17:03:19 2013 -0700
@@ -43,12 +43,12 @@
* this test is a strong assurance that the BigInteger operations
* are working correctly.
*
- * Three arguments may be specified which give the number of
- * decimal digits you desire in the three batches of test numbers.
+ * Four arguments may be specified which give the number of
+ * decimal digits you desire in the four batches of test numbers.
*
* The tests are performed on arrays of random numbers which are
* generated by a Random class as well as special cases which
- * throw in boundary numbers such as 0, 1, maximum sized, etc.
+ * throw in boundary numbers such as 0, 1, maximum SIZEd, etc.
*
*/
public class BigIntegerTest {
@@ -63,11 +63,14 @@
//
// SCHOENHAGE_BASE_CONVERSION_THRESHOLD = 8 ints = 256 bits
//
+ // BURNIKEL_ZIEGLER_THRESHOLD = 50 ints = 1600 bits
+ //
static final int BITS_KARATSUBA = 1600;
static final int BITS_TOOM_COOK = 2400;
static final int BITS_KARATSUBA_SQUARE = 2880;
static final int BITS_TOOM_COOK_SQUARE = 4480;
static final int BITS_SCHOENHAGE_BASE = 256;
+ static final int BITS_BURNIKEL_ZIEGLER = 1600;
static final int ORDER_SMALL = 60;
static final int ORDER_MEDIUM = 100;
@@ -80,14 +83,15 @@
// #bits for testing Toom-Cook squaring
static final int ORDER_TOOM_COOK_SQUARE = 4600;
+ static final int SIZE = 1000; // numbers per batch
+
static Random rnd = new Random();
- static int size = 1000; // numbers per batch
static boolean failure = false;
public static void pow(int order) {
int failCount1 = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
// Test identity x^power == x*x*x ... *x
int power = rnd.nextInt(6) + 2;
BigInteger x = fetchNumber(order);
@@ -106,7 +110,7 @@
public static void square(int order) {
int failCount1 = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
// Test identity x^2 == x*x
BigInteger x = fetchNumber(order);
BigInteger xx = x.multiply(x);
@@ -121,7 +125,7 @@
public static void arithmetic(int order) {
int failCount = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(order);
while(x.compareTo(BigInteger.ZERO) != 1)
x = fetchNumber(order);
@@ -187,7 +191,7 @@
int failCount = 0;
BigInteger base = BigInteger.ONE.shiftLeft(BITS_KARATSUBA - 32 - 1);
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(BITS_KARATSUBA - 32 - 1);
BigInteger u = base.add(x);
int a = 1 + rnd.nextInt(31);
@@ -210,7 +214,7 @@
failCount = 0;
base = base.shiftLeft(BITS_TOOM_COOK - BITS_KARATSUBA);
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(BITS_TOOM_COOK - 32 - 1);
BigInteger u = base.add(x);
BigInteger u2 = u.shiftLeft(1);
@@ -241,7 +245,7 @@
int failCount = 0;
BigInteger base = BigInteger.ONE.shiftLeft(BITS_KARATSUBA_SQUARE - 32 - 1);
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(BITS_KARATSUBA_SQUARE - 32 - 1);
BigInteger u = base.add(x);
int a = 1 + rnd.nextInt(31);
@@ -259,7 +263,7 @@
failCount = 0;
base = base.shiftLeft(BITS_TOOM_COOK_SQUARE - BITS_KARATSUBA_SQUARE);
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(BITS_TOOM_COOK_SQUARE - 32 - 1);
BigInteger u = base.add(x);
int a = 1 + rnd.nextInt(31);
@@ -276,10 +280,61 @@
report("squareLarge Toom-Cook", failCount);
}
+ /**
+ * Sanity test for Burnikel-Ziegler division. The Burnikel-Ziegler division
+ * algorithm is used when each of the dividend and the divisor has at least
+ * a specified number of ints in its representation. This test is based on
+ * the observation that if {@code w = u*pow(2,a)} and {@code z = v*pow(2,b)}
+ * where {@code abs(u) > 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<SIZE; i++) {
+ BigInteger addend = new BigInteger(BITS_BURNIKEL_ZIEGLER - 34, rnd);
+ BigInteger v = base.add(addend);
+
+ BigInteger u = v.multiply(BigInteger.valueOf(2 + rnd.nextInt(Short.MAX_VALUE - 1)));
+
+ if(rnd.nextBoolean()) {
+ u = u.negate();
+ }
+ if(rnd.nextBoolean()) {
+ v = v.negate();
+ }
+
+ int a = 17 + rnd.nextInt(16);
+ int b = 1 + rnd.nextInt(16);
+ BigInteger w = u.multiply(BigInteger.valueOf(1L << a));
+ BigInteger z = v.multiply(BigInteger.valueOf(1L << b));
+
+ BigInteger[] divideResult = u.divideAndRemainder(v);
+ divideResult[0] = divideResult[0].multiply(BigInteger.valueOf(1L << (a - b)));
+ divideResult[1] = divideResult[1].multiply(BigInteger.valueOf(1L << b));
+ BigInteger[] bzResult = w.divideAndRemainder(z);
+
+ if (divideResult[0].compareTo(bzResult[0]) != 0 ||
+ divideResult[1].compareTo(bzResult[1]) != 0) {
+ failCount++;
+ }
+ }
+
+ report("divideLarge", failCount);
+ }
+
public static void bitCount() {
int failCount = 0;
- for (int i=0; i<size*10; i++) {
+ for (int i=0; i<SIZE*10; i++) {
int x = rnd.nextInt();
BigInteger bigX = BigInteger.valueOf((long)x);
int bit = (x < 0 ? 0 : 1);
@@ -300,7 +355,7 @@
public static void bitLength() {
int failCount = 0;
- for (int i=0; i<size*10; i++) {
+ for (int i=0; i<SIZE*10; i++) {
int x = rnd.nextInt();
BigInteger bigX = BigInteger.valueOf((long)x);
int signBit = (x < 0 ? 0x80000000 : 0);
@@ -321,7 +376,7 @@
public static void bitOps(int order) {
int failCount1 = 0, failCount2 = 0, failCount3 = 0;
- for (int i=0; i<size*5; i++) {
+ for (int i=0; i<SIZE*5; i++) {
BigInteger x = fetchNumber(order);
BigInteger y;
@@ -351,7 +406,7 @@
report("clearBit/testBit for " + order + " bits", failCount1);
report("flipBit/testBit for " + order + " bits", failCount2);
- for (int i=0; i<size*5; i++) {
+ for (int i=0; i<SIZE*5; i++) {
BigInteger x = fetchNumber(order);
// Test getLowestSetBit()
@@ -375,7 +430,7 @@
// Test identity x^y == x|y &~ x&y
int failCount = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(order);
BigInteger y = fetchNumber(order);
BigInteger z = x.xor(y);
@@ -387,7 +442,7 @@
// Test identity x &~ y == ~(~x | y)
failCount = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(order);
BigInteger y = fetchNumber(order);
BigInteger z = x.andNot(y);
@@ -440,7 +495,7 @@
public static void divideAndRemainder(int order) {
int failCount1 = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(order).abs();
while(x.compareTo(BigInteger.valueOf(3L)) != 1)
x = fetchNumber(order).abs();
@@ -519,7 +574,7 @@
public static void byteArrayConv(int order) {
int failCount = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(order);
while (x.equals(BigInteger.ZERO))
x = fetchNumber(order);
@@ -536,7 +591,7 @@
public static void modInv(int order) {
int failCount = 0, successCount = 0, nonInvCount = 0;
- for (int i=0; i<size; i++) {
+ for (int i=0; i<SIZE; i++) {
BigInteger x = fetchNumber(order);
while(x.equals(BigInteger.ZERO))
x = fetchNumber(order);
@@ -565,7 +620,7 @@
public static void modExp(int order1, int order2) {
int failCount = 0;
- for (int i=0; i<size/10; i++) {
+ for (int i=0; i<SIZE/10; i++) {
BigInteger m = fetchNumber(order1).abs();
while(m.compareTo(BigInteger.ONE) != 1)
m = fetchNumber(order1).abs();
@@ -883,8 +938,8 @@
/**
* Main to interpret arguments and run several tests.
*
- * Up to three arguments may be given to specify the size of BigIntegers
- * used for call parameters 1, 2, and 3. The size is interpreted as
+ * Up to three arguments may be given to specify the SIZE of BigIntegers
+ * used for call parameters 1, 2, and 3. The SIZE is interpreted as
* the maximum number of decimal digits that the parameters will have.
*
*/
@@ -945,6 +1000,7 @@
multiplyLarge();
squareLarge();
+ divideLarge();
if (failure)
throw new RuntimeException("Failure in BigIntegerTest.");
@@ -952,7 +1008,7 @@
/*
* Get a random or boundary-case number. This is designed to provide
- * a lot of numbers that will find failure points, such as max sized
+ * a lot of numbers that will find failure points, such as max SIZEd
* numbers, empty BigIntegers, etc.
*
* If order is less than 2, order is changed to 2.
@@ -987,13 +1043,13 @@
break;
case 4: // Random bit density
- int iterations = rnd.nextInt(order-1);
- result = BigInteger.ONE.shiftLeft(rnd.nextInt(order));
- for(int i=0; i<iterations; i++) {
- BigInteger temp = BigInteger.ONE.shiftLeft(
- rnd.nextInt(order));
- result = result.or(temp);
+ byte[] val = new byte[(order+7)/8];
+ int iterations = rnd.nextInt(order);
+ for (int i=0; i<iterations; i++) {
+ int bitIdx = rnd.nextInt(order);
+ val[bitIdx/8] |= 1 << (bitIdx%8);
}
+ result = new BigInteger(1, val);
break;
case 5: // Runs of consecutive ones and zeros
result = ZERO;