8014319: Faster division of large integers

Summary: Implement Burnickel-Ziegler division algorithm in BigInteger
Reviewed-by: bpb, martin
Contributed-by: NTim Buktu <tbuktu@hotmail.com>

src/share/classes/java/math/BigInteger.java

@@ -33,7 +33,6 @@ import java.io.IOException;
 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 @@ import sun.misc.FloatConsts;
  * @author  Josh Bloch
  * @author  Michael McCloskey
  * @author  Alan Eliasen
+ * @author  Timothy Buktu
  * @since JDK1.1
  */
 
@@ -214,6 +214,14 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
      */
     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,
@@ -1781,7 +1789,7 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
             if (len < TOOM_COOK_SQUARE_THRESHOLD)
                 return squareKaratsuba();
             else
-               return squareToomCook3();
+                return squareToomCook3();
     }
 
     /**
@@ -1936,11 +1944,26 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
      * @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 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
      * @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 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
      * @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 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
 
     /**
      * 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 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
      * 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

src/share/classes/java/math/MutableBigInteger.java

 /*
- * 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,13 +38,13 @@ package java.math;
  *
  * @see     BigInteger
  * @author  Michael McCloskey
+ * @author  Timothy Buktu
  * @since   1.3
  */
 
-import java.util.Arrays;
-
-import static java.math.BigInteger.LONG_MASK;
 import static java.math.BigDecimal.INFLATED;
+import static java.math.BigInteger.LONG_MASK;
+import java.util.Arrays;
 
 class MutableBigInteger {
     /**
@@ -75,6 +75,24 @@ class MutableBigInteger {
      */
     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
 
     /**
@@ -123,6 +141,20 @@ class MutableBigInteger {
         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.
@@ -154,6 +186,14 @@ class MutableBigInteger {
         return new BigInteger(getMagnitudeArray(), sign);
     }
 
+    /**
+     * Converts this number to a nonnegative {@code BigInteger}.
+     */
+    BigInteger toBigInteger() {
+        normalize();
+        return toBigInteger(isZero() ? 0 : 1);
+    }
+
     /**
      * Convert this MutableBigInteger to BigDecimal object with the specified sign
      * and scale.
@@ -237,6 +277,32 @@ class MutableBigInteger {
         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).
@@ -453,6 +519,16 @@ class MutableBigInteger {
         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.
@@ -474,6 +550,14 @@ class MutableBigInteger {
         }
     }
 
+    /**
+     * 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.
      */
@@ -614,6 +698,35 @@ class MutableBigInteger {
         val[offset+intLen-1] <<= n;
     }
 
+    /**
+     * Returns a {@code BigInteger} equal to the {@code n}
+     * low ints of this number.
+     */
+    private BigInteger getLower(int n) {
+        if (isZero())
+            return BigInteger.ZERO;
+        else if (intLen < n)
+            return toBigInteger(1);
+        else {
+            // strip zeros
+            int len = n;
+            while (len>0 && value[offset+intLen-len]==0)
+                len--;
+            int sign = len>0 ? 1 : 0;
+            return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
+        }
+    }
+
+    /**
+     * Discards all ints whose index is greater than {@code n}.
+     */
+    private void keepLower(int n) {
+        if (intLen >= n) {
+            offset += intLen - n;
+            intLen = n;
+        }
+    }
+
     /**
      * Adds the contents of two MutableBigInteger objects.The result
      * is placed within this MutableBigInteger.
@@ -673,6 +786,121 @@ class MutableBigInteger {
         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
@@ -906,6 +1134,29 @@ class MutableBigInteger {
             return rem;
     }
 
+    /**
+     * Calculates the quotient of this div b and places the quotient in the
+     * provided MutableBigInteger objects and the remainder object is returned.
+     *
+     */
+    MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient) {
+        return divide(b,quotient,true);
+    }
+
+    MutableBigInteger divide(MutableBigInteger b, MutableBigInteger quotient, boolean needRemainder) {
+        if (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.
@@ -917,38 +1168,34 @@ class MutableBigInteger {
      * 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 @@ class MutableBigInteger {
                 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 @@ class MutableBigInteger {
      */
     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 @@ class MutableBigInteger {
             // 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 @@ class MutableBigInteger {
             // 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 @@ class MutableBigInteger {
         }
 
 
-        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 @@ class MutableBigInteger {
      * 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 @@ class MutableBigInteger {
         return result;
     }
 
-    /*
+    /**
      * Returns the multiplicative inverse of val mod 2^32.  Assumes val is odd.
      */
     static int inverseMod32(int val) {
@@ -1595,7 +2051,7 @@ class MutableBigInteger {
         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 @@ class MutableBigInteger {
         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.

test/java/math/BigInteger/BigIntegerTest.java

@@ -43,12 +43,12 @@ import java.util.Random;
  * 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 @@ public class BigIntegerTest {
     //
     // 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 @@ public class BigIntegerTest {
     // #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 class BigIntegerTest {
     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 class BigIntegerTest {
     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 @@ public class BigIntegerTest {
         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 @@ public class BigIntegerTest {
 
         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 @@ public class BigIntegerTest {
         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 @@ public class BigIntegerTest {
 
         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 @@ public class BigIntegerTest {
         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 class BigIntegerTest {
     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 class BigIntegerTest {
     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 @@ public class BigIntegerTest {
         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 @@ public class BigIntegerTest {
 
         // 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 @@ public class BigIntegerTest {
 
         // 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 class BigIntegerTest {
     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 class BigIntegerTest {
     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 class BigIntegerTest {
     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 class BigIntegerTest {
     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 @@ public class BigIntegerTest {
     /**
      * 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 @@ public class BigIntegerTest {
 
         multiplyLarge();
         squareLarge();
+        divideLarge();
 
         if (failure)
             throw new RuntimeException("Failure in BigIntegerTest.");
@@ -952,7 +1008,7 @@ public class BigIntegerTest {
 
     /*
      * 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 @@ public class BigIntegerTest {
                 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;