4837946: Faster multiplication and exponentiation of large integers

4646474: BigInteger.pow() algorithm slow in 1.4.0
Summary: Implement Karatsuba and 3-way Toom-Cook multiplication as well as exponentiation using Karatsuba and Toom-Cook squaring.
Reviewed-by: alanb, bpb, martin
Contributed-by: Alan Eliasen <eliasen@mindspring.com>

2 files changed, 536 insertions, 67 deletions

diff --git a/src/share/classes/java/math/BigDecimal.java b/src/share/classes/java/math/BigDecimal.java
index b7f6380948..d2a4492237 100644
--- a/src/share/classes/java/math/BigDecimal.java
+++ b/src/share/classes/java/math/BigDecimal.java
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 1996, 2011, Oracle and/or its affiliates. All rights reserved.
+ * Copyright (c) 1996, 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
@@ -3538,24 +3538,7 @@ public class BigDecimal extends Number implements Comparable<BigDecimal> {
                 return expandBigIntegerTenPowers(n);
         }
 
-        if (n < 1024*524288) {
-            // BigInteger.pow is slow, so make 10**n by constructing a
-            // BigInteger from a character string (still not very fast)
-            // which occupies no more than 1GB (!) of memory.
-            char tenpow[] = new char[n + 1];
-            tenpow[0] = '1';
-            for (int i = 1; i <= n; i++) {
-                tenpow[i] = '0';
-            }
-            return new BigInteger(tenpow, 1, tenpow.length);
-        }
-
-        if ((n & 0x1) == 0x1) {
-            return BigInteger.TEN.multiply(bigTenToThe(n - 1));
-        } else {
-            BigInteger tmp = bigTenToThe(n/2);
-            return tmp.multiply(tmp);
-        }
+        return BigInteger.TEN.pow(n);
     }
 
     /**
diff --git a/src/share/classes/java/math/BigInteger.java b/src/share/classes/java/math/BigInteger.java
index a4988deaa4..4c8b7816d0 100644
--- a/src/share/classes/java/math/BigInteger.java
+++ b/src/share/classes/java/math/BigInteger.java
@@ -1,5 +1,5 @@
 /*
- * Copyright (c) 1996, 2011, Oracle and/or its affiliates. All rights reserved.
+ * Copyright (c) 1996, 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
@@ -29,9 +29,12 @@
 
 package java.math;
 
-import java.util.Random;
-import java.io.*;
+import java.io.IOException;
+import java.io.ObjectInputStream;
+import java.io.ObjectOutputStream;
+import java.io.ObjectStreamField;
 import java.util.Arrays;
+import java.util.Random;
 
 /**
  * Immutable arbitrary-precision integers.  All operations behave as if
@@ -94,6 +97,7 @@ import java.util.Arrays;
  * @see     BigDecimal
  * @author  Josh Bloch
  * @author  Michael McCloskey
+ * @author  Alan Eliasen
  * @since JDK1.1
  */
 
@@ -174,6 +178,39 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
      */
     final static long LONG_MASK = 0xffffffffL;
 
+    /**
+     * The threshold value for using Karatsuba multiplication.  If the number
+     * of ints in both mag arrays are greater than this number, then
+     * Karatsuba multiplication will be used.   This value is found
+     * experimentally to work well.
+     */
+    private static final int KARATSUBA_THRESHOLD = 50;
+
+    /**
+     * The threshold value for using 3-way Toom-Cook multiplication.
+     * If the number of ints in each mag array is greater than the
+     * Karatsuba threshold, and the number of ints in at least one of
+     * the mag arrays is greater than this threshold, then Toom-Cook
+     * multiplication will be used.
+     */
+    private static final int TOOM_COOK_THRESHOLD = 75;
+
+    /**
+     * The threshold value for using Karatsuba squaring.  If the number
+     * of ints in the number are larger than this value,
+     * Karatsuba squaring will be used.   This value is found
+     * experimentally to work well.
+     */
+    private static final int KARATSUBA_SQUARE_THRESHOLD = 90;
+
+    /**
+     * The threshold value for using Toom-Cook squaring.  If the number
+     * of ints in the number are larger than this value,
+     * Toom-Cook squaring will be used.   This value is found
+     * experimentally to work well.
+     */
+    private static final int TOOM_COOK_SQUARE_THRESHOLD = 140;
+
     //Constructors
 
     /**
@@ -522,15 +559,16 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
 
         if (bitLength < 2)
             throw new ArithmeticException("bitLength < 2");
-        // The cutoff of 95 was chosen empirically for best performance
-        prime = (bitLength < 95 ? smallPrime(bitLength, certainty, rnd)
+        prime = (bitLength < SMALL_PRIME_THRESHOLD
+                                ? smallPrime(bitLength, certainty, rnd)
                                 : largePrime(bitLength, certainty, rnd));
         signum = 1;
         mag = prime.mag;
     }
 
     // Minimum size in bits that the requested prime number has
-    // before we use the large prime number generating algorithms
+    // before we use the large prime number generating algorithms.
+    // The cutoff of 95 was chosen empirically for best performance.
     private static final int SMALL_PRIME_THRESHOLD = 95;
 
     // Certainty required to meet the spec of probablePrime
@@ -553,7 +591,6 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
         if (bitLength < 2)
             throw new ArithmeticException("bitLength < 2");
 
-        // The cutoff of 95 was chosen empirically for best performance
         return (bitLength < SMALL_PRIME_THRESHOLD ?
                 smallPrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd) :
                 largePrime(bitLength, DEFAULT_PRIME_CERTAINTY, rnd));
@@ -986,6 +1023,7 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
     private final static int MAX_CONSTANT = 16;
     private static BigInteger posConst[] = new BigInteger[MAX_CONSTANT+1];
     private static BigInteger negConst[] = new BigInteger[MAX_CONSTANT+1];
+
     static {
         for (int i = 1; i <= MAX_CONSTANT; i++) {
             int[] magnitude = new int[1];
@@ -1015,6 +1053,11 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
     private static final BigInteger TWO = valueOf(2);
 
     /**
+     * The BigInteger constant -1.  (Not exported.)
+     */
+    private static final BigInteger NEGATIVE_ONE = valueOf(-1);
+
+    /**
      * The BigInteger constant ten.
      *
      * @since   1.5
@@ -1290,17 +1333,29 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
     public BigInteger multiply(BigInteger val) {
         if (val.signum == 0 || signum == 0)
             return ZERO;
-        int resultSign = signum == val.signum ? 1 : -1;
-        if (val.mag.length == 1) {
-            return  multiplyByInt(mag,val.mag[0], resultSign);
-        }
-        if(mag.length == 1) {
-            return multiplyByInt(val.mag,mag[0], resultSign);
+
+        int xlen = mag.length;
+        int ylen = val.mag.length;
+
+        if ((xlen < KARATSUBA_THRESHOLD) || (ylen < KARATSUBA_THRESHOLD))
+        {
+            int resultSign = signum == val.signum ? 1 : -1;
+            if (val.mag.length == 1) {
+                return multiplyByInt(mag,val.mag[0], resultSign);
+            }
+            if(mag.length == 1) {
+                return multiplyByInt(val.mag,mag[0], resultSign);
+            }
+            int[] result = multiplyToLen(mag, xlen,
+                                         val.mag, ylen, null);
+            result = trustedStripLeadingZeroInts(result);
+            return new BigInteger(result, resultSign);
         }
-        int[] result = multiplyToLen(mag, mag.length,
-                                     val.mag, val.mag.length, null);
-        result = trustedStripLeadingZeroInts(result);
-        return new BigInteger(result, resultSign);
+        else
+            if ((xlen < TOOM_COOK_THRESHOLD) && (ylen < TOOM_COOK_THRESHOLD))
+                return multiplyKaratsuba(this, val);
+            else
+               return multiplyToomCook3(this, val);
     }
 
     private static BigInteger multiplyByInt(int[] x, int y, int sign) {
@@ -1402,6 +1457,272 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
     }
 
     /**
+     * Multiplies two BigIntegers using the Karatsuba multiplication
+     * algorithm.  This is a recursive divide-and-conquer algorithm which is
+     * more efficient for large numbers than what is commonly called the
+     * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
+     * multiplied have length n, the "grade-school" algorithm has an
+     * asymptotic complexity of O(n^2).  In contrast, the Karatsuba algorithm
+     * has complexity of O(n^(log2(3))), or O(n^1.585).  It achieves this
+     * increased performance by doing 3 multiplies instead of 4 when
+     * evaluating the product.  As it has some overhead, should be used when
+     * both numbers are larger than a certain threshold (found
+     * experimentally).
+     *
+     * See:  http://en.wikipedia.org/wiki/Karatsuba_algorithm
+     */
+    private static BigInteger multiplyKaratsuba(BigInteger x, BigInteger y)
+    {
+        int xlen = x.mag.length;
+        int ylen = y.mag.length;
+
+        // The number of ints in each half of the number.
+        int half = (Math.max(xlen, ylen)+1) / 2;
+
+        // xl and yl are the lower halves of x and y respectively,
+        // xh and yh are the upper halves.
+        BigInteger xl = x.getLower(half);
+        BigInteger xh = x.getUpper(half);
+        BigInteger yl = y.getLower(half);
+        BigInteger yh = y.getUpper(half);
+
+        BigInteger p1 = xh.multiply(yh);  // p1 = xh*yh
+        BigInteger p2 = xl.multiply(yl);  // p2 = xl*yl
+
+        // p3=(xh+xl)*(yh+yl)
+        BigInteger p3 = xh.add(xl).multiply(yh.add(yl));
+
+        // result = p1 * 2^(32*2*half) + (p3 - p1 - p2) * 2^(32*half) + p2
+        BigInteger result = p1.shiftLeft(32*half).add(p3.subtract(p1).subtract(p2)).shiftLeft(32*half).add(p2);
+
+        if (x.signum != y.signum)
+            return result.negate();
+        else
+            return result;
+    }
+
+    /**
+     * Multiplies two BigIntegers using a 3-way Toom-Cook multiplication
+     * algorithm.  This is a recursive divide-and-conquer algorithm which is
+     * more efficient for large numbers than what is commonly called the
+     * "grade-school" algorithm used in multiplyToLen.  If the numbers to be
+     * multiplied have length n, the "grade-school" algorithm has an
+     * asymptotic complexity of O(n^2).  In contrast, 3-way Toom-Cook has a
+     * complexity of about O(n^1.465).  It achieves this increased asymptotic
+     * performance by breaking each number into three parts and by doing 5
+     * multiplies instead of 9 when evaluating the product.  Due to overhead
+     * (additions, shifts, and one division) in the Toom-Cook algorithm, it
+     * should only be used when both numbers are larger than a certain
+     * threshold (found experimentally).  This threshold is generally larger
+     * than that for Karatsuba multiplication, so this algorithm is generally
+     * only used when numbers become significantly larger.
+     *
+     * The algorithm used is the "optimal" 3-way Toom-Cook algorithm outlined
+     * by Marco Bodrato.
+     *
+     *  See: http://bodrato.it/toom-cook/
+     *       http://bodrato.it/papers/#WAIFI2007
+     *
+     * "Towards Optimal Toom-Cook Multiplication for Univariate and
+     * Multivariate Polynomials in Characteristic 2 and 0." by Marco BODRATO;
+     * In C.Carlet and B.Sunar, Eds., "WAIFI'07 proceedings", p. 116-133,
+     * LNCS #4547. Springer, Madrid, Spain, June 21-22, 2007.
+     *
+     */
+    private static BigInteger multiplyToomCook3(BigInteger a, BigInteger b)
+    {
+        int alen = a.mag.length;
+        int blen = b.mag.length;
+
+        int largest = Math.max(alen, blen);
+
+        // k is the size (in ints) of the lower-order slices.
+        int k = (largest+2)/3;   // Equal to ceil(largest/3)
+
+        // r is the size (in ints) of the highest-order slice.
+        int r = largest - 2*k;
+
+        // Obtain slices of the numbers. a2 and b2 are the most significant
+        // bits of the numbers a and b, and a0 and b0 the least significant.
+        BigInteger a0, a1, a2, b0, b1, b2;
+        a2 = a.getToomSlice(k, r, 0, largest);
+        a1 = a.getToomSlice(k, r, 1, largest);
+        a0 = a.getToomSlice(k, r, 2, largest);
+        b2 = b.getToomSlice(k, r, 0, largest);
+        b1 = b.getToomSlice(k, r, 1, largest);
+        b0 = b.getToomSlice(k, r, 2, largest);
+
+        BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1, db1;
+
+        v0 = a0.multiply(b0);
+        da1 = a2.add(a0);
+        db1 = b2.add(b0);
+        vm1 = da1.subtract(a1).multiply(db1.subtract(b1));
+        da1 = da1.add(a1);
+        db1 = db1.add(b1);
+        v1 = da1.multiply(db1);
+        v2 = da1.add(a2).shiftLeft(1).subtract(a0).multiply(
+             db1.add(b2).shiftLeft(1).subtract(b0));
+        vinf = a2.multiply(b2);
+
+        /* The algorithm requires two divisions by 2 and one by 3.
+           All divisions are known to be exact, that is, they do not produce
+           remainders, and all results are positive.  The divisions by 2 are
+           implemented as right shifts which are relatively efficient, leaving
+           only an exact division by 3, which is done by a specialized
+           linear-time algorithm. */
+        t2 = v2.subtract(vm1).exactDivideBy3();
+        tm1 = v1.subtract(vm1).shiftRight(1);
+        t1 = v1.subtract(v0);
+        t2 = t2.subtract(t1).shiftRight(1);
+        t1 = t1.subtract(tm1).subtract(vinf);
+        t2 = t2.subtract(vinf.shiftLeft(1));
+        tm1 = tm1.subtract(t2);
+
+        // Number of bits to shift left.
+        int ss = k*32;
+
+        BigInteger result = vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
+
+        if (a.signum != b.signum)
+            return result.negate();
+        else
+            return result;
+    }
+
+
+    /**
+     * Returns a slice of a BigInteger for use in Toom-Cook multiplication.
+     *
+     * @param lowerSize The size of the lower-order bit slices.
+     * @param upperSize The size of the higher-order bit slices.
+     * @param slice The index of which slice is requested, which must be a
+     * number from 0 to size-1. Slice 0 is the highest-order bits, and slice
+     * size-1 are the lowest-order bits. Slice 0 may be of different size than
+     * the other slices.
+     * @param fullsize The size of the larger integer array, used to align
+     * slices to the appropriate position when multiplying different-sized
+     * numbers.
+     */
+    private BigInteger getToomSlice(int lowerSize, int upperSize, int slice,
+                                    int fullsize)
+    {
+        int start, end, sliceSize, len, offset;
+
+        len = mag.length;
+        offset = fullsize - len;
+
+        if (slice == 0)
+        {
+            start = 0 - offset;
+            end = upperSize - 1 - offset;
+        }
+        else
+        {
+            start = upperSize + (slice-1)*lowerSize - offset;
+            end = start + lowerSize - 1;
+        }
+
+        if (start < 0)
+            start = 0;
+        if (end < 0)
+           return ZERO;
+
+        sliceSize = (end-start) + 1;
+
+        if (sliceSize <= 0)
+            return ZERO;
+
+        // While performing Toom-Cook, all slices are positive and
+        // the sign is adjusted when the final number is composed.
+        if (start==0 && sliceSize >= len)
+            return this.abs();
+
+        int intSlice[] = new int[sliceSize];
+        System.arraycopy(mag, start, intSlice, 0, sliceSize);
+
+        return new BigInteger(trustedStripLeadingZeroInts(intSlice), 1);
+    }
+
+    /**
+     * Does an exact division (that is, the remainder is known to be zero)
+     * of the specified number by 3.  This is used in Toom-Cook
+     * multiplication.  This is an efficient algorithm that runs in linear
+     * time.  If the argument is not exactly divisible by 3, results are
+     * undefined.  Note that this is expected to be called with positive
+     * arguments only.
+     */
+    private BigInteger exactDivideBy3()
+    {
+        int len = mag.length;
+        int[] result = new int[len];
+        long x, w, q, borrow;
+        borrow = 0L;
+        for (int i=len-1; i>=0; i--)
+        {
+            x = (mag[i] & LONG_MASK);
+            w = x - borrow;
+            if (borrow > x)       // Did we make the number go negative?
+                borrow = 1L;
+            else
+                borrow = 0L;
+
+            // 0xAAAAAAAB is the modular inverse of 3 (mod 2^32).  Thus,
+            // the effect of this is to divide by 3 (mod 2^32).
+            // This is much faster than division on most architectures.
+            q = (w * 0xAAAAAAABL) & LONG_MASK;
+            result[i] = (int) q;
+
+            // Now check the borrow. The second check can of course be
+            // eliminated if the first fails.
+            if (q >= 0x55555556L)
+            {
+                borrow++;
+                if (q >= 0xAAAAAAABL)
+                    borrow++;
+            }
+        }
+        result = trustedStripLeadingZeroInts(result);
+        return new BigInteger(result, signum);
+    }
+
+    /**
+     * Returns a new BigInteger representing n lower ints of the number.
+     * This is used by Karatsuba multiplication and Karatsuba squaring.
+     */
+    private BigInteger getLower(int n) {
+        int len = mag.length;
+
+        if (len <= n)
+            return this;
+
+        int lowerInts[] = new int[n];
+        System.arraycopy(mag, len-n, lowerInts, 0, n);
+
+        return new BigInteger(trustedStripLeadingZeroInts(lowerInts), 1);
+    }
+
+    /**
+     * Returns a new BigInteger representing mag.length-n upper
+     * ints of the number.  This is used by Karatsuba multiplication and
+     * Karatsuba squaring.
+     */
+    private BigInteger getUpper(int n) {
+        int len = mag.length;
+
+        if (len <= n)
+            return ZERO;
+
+        int upperLen = len - n;
+        int upperInts[] = new int[upperLen];
+        System.arraycopy(mag, 0, upperInts, 0, upperLen);
+
+        return new BigInteger(trustedStripLeadingZeroInts(upperInts), 1);
+    }
+
+    // Squaring
+
+    /**
      * Returns a BigInteger whose value is {@code (this<sup>2</sup>)}.
      *
      * @return {@code this<sup>2</sup>}
@@ -1409,8 +1730,18 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
     private BigInteger square() {
         if (signum == 0)
             return ZERO;
-        int[] z = squareToLen(mag, mag.length, null);
-        return new BigInteger(trustedStripLeadingZeroInts(z), 1);
+        int len = mag.length;
+
+        if (len < KARATSUBA_SQUARE_THRESHOLD)
+        {
+            int[] z = squareToLen(mag, len, null);
+            return new BigInteger(trustedStripLeadingZeroInts(z), 1);
+        }
+        else
+            if (len < TOOM_COOK_SQUARE_THRESHOLD)
+                return squareKaratsuba();
+            else
+               return squareToomCook3();
     }
 
     /**
@@ -1481,6 +1812,83 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
     }
 
     /**
+     * Squares a BigInteger using the Karatsuba squaring algorithm.  It should
+     * be used when both numbers are larger than a certain threshold (found
+     * experimentally).  It is a recursive divide-and-conquer algorithm that
+     * has better asymptotic performance than the algorithm used in
+     * squareToLen.
+     */
+    private BigInteger squareKaratsuba()
+    {
+        int half = (mag.length+1) / 2;
+
+        BigInteger xl = getLower(half);
+        BigInteger xh = getUpper(half);
+
+        BigInteger xhs = xh.square();  // xhs = xh^2
+        BigInteger xls = xl.square();  // xls = xl^2
+
+        // xh^2 << 64  +  (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2
+        return xhs.shiftLeft(half*32).add(xl.add(xh).square().subtract(xhs.add(xls))).shiftLeft(half*32).add(xls);
+    }
+
+    /**
+     * Squares a BigInteger using the 3-way Toom-Cook squaring algorithm.  It
+     * should be used when both numbers are larger than a certain threshold
+     * (found experimentally).  It is a recursive divide-and-conquer algorithm
+     * that has better asymptotic performance than the algorithm used in
+     * squareToLen or squareKaratsuba.
+     */
+    private BigInteger squareToomCook3()
+    {
+        int len = mag.length;
+
+        // k is the size (in ints) of the lower-order slices.
+        int k = (len+2)/3;   // Equal to ceil(largest/3)
+
+        // r is the size (in ints) of the highest-order slice.
+        int r = len - 2*k;
+
+        // Obtain slices of the numbers. a2 is the most significant
+        // bits of the number, and a0 the least significant.
+        BigInteger a0, a1, a2;
+        a2 = getToomSlice(k, r, 0, len);
+        a1 = getToomSlice(k, r, 1, len);
+        a0 = getToomSlice(k, r, 2, len);
+        BigInteger v0, v1, v2, vm1, vinf, t1, t2, tm1, da1;
+
+        v0 = a0.square();
+        da1 = a2.add(a0);
+        vm1 = da1.subtract(a1).square();
+        da1 = da1.add(a1);
+        v1 = da1.square();
+        vinf = a2.square();
+        v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();
+
+        /* The algorithm requires two divisions by 2 and one by 3.
+           All divisions are known to be exact, that is, they do not produce
+           remainders, and all results are positive.  The divisions by 2 are
+           implemented as right shifts which are relatively efficient, leaving
+           only a division by 3.
+           The division by 3 is done by an optimized algorithm for this case.
+        */
+        t2 = v2.subtract(vm1).exactDivideBy3();
+        tm1 = v1.subtract(vm1).shiftRight(1);
+        t1 = v1.subtract(v0);
+        t2 = t2.subtract(t1).shiftRight(1);
+        t1 = t1.subtract(tm1).subtract(vinf);
+        t2 = t2.subtract(vinf.shiftLeft(1));
+        tm1 = tm1.subtract(t2);
+
+        // Number of bits to shift left.
+        int ss = k*32;
+
+        return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
+    }
+
+    // Division
+
+    /**
      * Returns a BigInteger whose value is {@code (this / val)}.
      *
      * @param  val value by which this BigInteger is to be divided.
@@ -1549,23 +1957,100 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
         if (signum==0)
             return (exponent==0 ? ONE : this);
 
-        // Perform exponentiation using repeated squaring trick
-        int newSign = (signum<0 && (exponent&1)==1 ? -1 : 1);
-        int[] baseToPow2 = this.mag;
-        int[] result = {1};
+        BigInteger partToSquare = this.abs();
+
+        // Factor out powers of two from the base, as the exponentiation of
+        // these can be done by left shifts only.
+        // The remaining part can then be exponentiated faster.  The
+        // powers of two will be multiplied back at the end.
+        int powersOfTwo = partToSquare.getLowestSetBit();
+
+        int remainingBits;
 
-        while (exponent != 0) {
-            if ((exponent & 1)==1) {
-                result = multiplyToLen(result, result.length,
-                                       baseToPow2, baseToPow2.length, null);
-                result = trustedStripLeadingZeroInts(result);
+        // Factor the powers of two out quickly by shifting right, if needed.
+        if (powersOfTwo > 0)
+        {
+            partToSquare = partToSquare.shiftRight(powersOfTwo);
+            remainingBits = partToSquare.bitLength();
+            if (remainingBits == 1)  // Nothing left but +/- 1?
+                if (signum<0 && (exponent&1)==1)
+                    return NEGATIVE_ONE.shiftLeft(powersOfTwo*exponent);
+                else
+                    return ONE.shiftLeft(powersOfTwo*exponent);
+        }
+        else
+        {
+            remainingBits = partToSquare.bitLength();
+            if (remainingBits == 1)  // Nothing left but +/- 1?
+                if (signum<0 && (exponent&1)==1)
+                    return NEGATIVE_ONE;
+                else
+                    return ONE;
+        }
+
+        // This is a quick way to approximate the size of the result,
+        // similar to doing log2[n] * exponent.  This will give an upper bound
+        // of how big the result can be, and which algorithm to use.
+        int scaleFactor = remainingBits * exponent;
+
+        // Use slightly different algorithms for small and large operands.
+        // See if the result will safely fit into a long. (Largest 2^63-1)
+        if (partToSquare.mag.length==1 && scaleFactor <= 62)
+        {
+            // Small number algorithm.  Everything fits into a long.
+            int newSign = (signum<0 && (exponent&1)==1 ? -1 : 1);
+            long result = 1;
+            long baseToPow2 = partToSquare.mag[0] & LONG_MASK;
+
+            int workingExponent = exponent;
+
+            // Perform exponentiation using repeated squaring trick
+            while (workingExponent != 0) {
+                if ((workingExponent & 1)==1)
+                    result = result * baseToPow2;
+
+                if ((workingExponent >>>= 1) != 0)
+                    baseToPow2 = baseToPow2 * baseToPow2;
             }
-            if ((exponent >>>= 1) != 0) {
-                baseToPow2 = squareToLen(baseToPow2, baseToPow2.length, null);
-                baseToPow2 = trustedStripLeadingZeroInts(baseToPow2);
+
+            // Multiply back the powers of two (quickly, by shifting left)
+            if (powersOfTwo > 0)
+            {
+                int bitsToShift = powersOfTwo*exponent;
+                if (bitsToShift + scaleFactor <= 62) // Fits in long?
+                    return valueOf((result << bitsToShift) * newSign);
+                else
+                    return valueOf(result*newSign).shiftLeft(bitsToShift);
             }
+            else
+                return valueOf(result*newSign);
+        }
+        else
+        {
+            // Large number algorithm.  This is basically identical to
+            // the algorithm above, but calls multiply() and square()
+            // which may use more efficient algorithms for large numbers.
+            BigInteger answer = ONE;
+
+            int workingExponent = exponent;
+            // Perform exponentiation using repeated squaring trick
+            while (workingExponent != 0) {
+                if ((workingExponent & 1)==1)
+                    answer = answer.multiply(partToSquare);
+
+                if ((workingExponent >>>= 1) != 0)
+                    partToSquare = partToSquare.square();
+            }
+            // Multiply back the (exponentiated) powers of two (quickly,
+            // by shifting left)
+            if (powersOfTwo > 0)
+                answer = answer.shiftLeft(powersOfTwo*exponent);
+
+            if (signum<0 && (exponent&1)==1)
+                return answer.negate();
+            else
+                return answer;
         }
-        return new BigInteger(result, newSign);
     }
 
     /**
@@ -2117,7 +2602,7 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
          * Perform exponentiation using repeated squaring trick, chopping off
          * high order bits as indicated by modulus.
          */
-        BigInteger result = valueOf(1);
+        BigInteger result = ONE;
         BigInteger baseToPow2 = this.mod2(p);
         int expOffset = 0;
 
@@ -2850,6 +3335,7 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
         return buf.toString();
     }
 
+
     /* zero[i] is a string of i consecutive zeros. */
     private static String zeros[] = new String[64];
     static {
@@ -3218,21 +3704,21 @@ public class BigInteger extends Number implements Comparable<BigInteger> {
      * little-endian binary representation of the magnitude (int 0 is the
      * least significant). If the magnitude is zero, return value is undefined.
      */
-     private int firstNonzeroIntNum() {
-         int fn = firstNonzeroIntNum - 2;
-         if (fn == -2) { // firstNonzeroIntNum not initialized yet
-             fn = 0;
-
-             // Search for the first nonzero int
-             int i;
-             int mlen = mag.length;
-             for (i = mlen - 1; i >= 0 && mag[i] == 0; i--)
-                 ;
-             fn = mlen - i - 1;
-             firstNonzeroIntNum = fn + 2; // offset by two to initialize
-         }
-         return fn;
-     }
+    private int firstNonzeroIntNum() {
+        int fn = firstNonzeroIntNum - 2;
+        if (fn == -2) { // firstNonzeroIntNum not initialized yet
+            fn = 0;
+
+            // Search for the first nonzero int
+            int i;
+            int mlen = mag.length;
+            for (i = mlen - 1; i >= 0 && mag[i] == 0; i--)
+                ;
+            fn = mlen - i - 1;
+            firstNonzeroIntNum = fn + 2; // offset by two to initialize
+        }
+        return fn;
+    }
 
     /** use serialVersionUID from JDK 1.1. for interoperability */
     private static final long serialVersionUID = -8287574255936472291L;