32 import java.io.IOException;

    33 import java.io.ObjectInputStream;

    34 import java.io.ObjectOutputStream;

    35 import java.io.ObjectStreamField;

    36 import java.util.Arrays;

   100  * @author  Alan Eliasen

   180 

   181     /**

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

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

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

   185      * experimentally to work well.

   186      */

   187     private static final int KARATSUBA_THRESHOLD = 50;

   188 

   189     /**

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

   191      * If the number of ints in each mag array is greater than the

   192      * Karatsuba threshold, and the number of ints in at least one of

   193      * the mag arrays is greater than this threshold, then Toom-Cook

   194      * multiplication will be used.

   195      */

   196     private static final int TOOM_COOK_THRESHOLD = 75;

   197 

   198     /**

   199      * The threshold value for using Karatsuba squaring.  If the number

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

   201      * Karatsuba squaring will be used.   This value is found

   202      * experimentally to work well.

   203      */

   204     private static final int KARATSUBA_SQUARE_THRESHOLD = 90;

   205 

   206     /**

   207      * The threshold value for using Toom-Cook squaring.  If the number

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

   209      * Toom-Cook squaring will be used.   This value is found

   210      * experimentally to work well.

   211      */

   212     private static final int TOOM_COOK_SQUARE_THRESHOLD = 140;

  1026 

  1054 

  1055     /**

  1056      * The BigInteger constant -1.  (Not exported.)

  1057      */

  1058     private static final BigInteger NEGATIVE_ONE = valueOf(-1);

  1460      * Multiplies two BigIntegers using the Karatsuba multiplication

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

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

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

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

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

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

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

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

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

  1470      * experimentally).

  1471      *

  1472      * See:  http://en.wikipedia.org/wiki/Karatsuba_algorithm

  1473      */

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

  1475     {

  1476         int xlen = x.mag.length;

  1477         int ylen = y.mag.length;

  1478 

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

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

  1481 

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

  1483         // xh and yh are the upper halves.

  1484         BigInteger xl = x.getLower(half);

  1485         BigInteger xh = x.getUpper(half);

  1486         BigInteger yl = y.getLower(half);

  1487         BigInteger yh = y.getUpper(half);

  1488 

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

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

  1491 

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

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

  1494 

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

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

  1497 

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

  1499             return result.negate();

  1500         else

  1501             return result;

  1502     }

  1503 

  1504     /**

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

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

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

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

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

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

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

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

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

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

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

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

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

  1518      * only used when numbers become significantly larger.

  1519      *

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

  1521      * by Marco Bodrato.

  1522      *

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

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

  1525      *

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

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

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

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

  1530      *

  1531      */

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

  1533     {

  1534         int alen = a.mag.length;

  1535         int blen = b.mag.length;

  1536 

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

  1538 

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

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

  1541 

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

  1543         int r = largest - 2*k;

  1544 

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

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

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

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

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

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

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

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

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

  1554 

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

  1556 

  1557         v0 = a0.multiply(b0);

  1558         da1 = a2.add(a0);

  1559         db1 = b2.add(b0);

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

  1561         da1 = da1.add(a1);

  1562         db1 = db1.add(b1);

  1563         v1 = da1.multiply(db1);

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

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

  1566         vinf = a2.multiply(b2);

  1567 

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

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

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

  1571            implemented as right shifts which are relatively efficient, leaving

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

  1573            linear-time algorithm. */

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

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

  1576         t1 = v1.subtract(v0);

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

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

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

  1580         tm1 = tm1.subtract(t2);

  1581 

  1582         // Number of bits to shift left.

  1583         int ss = k*32;

  1584 

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

  1586 

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

  1588             return result.negate();

  1589         else

  1590             return result;

  1591     }

  1592 

  1593 

  1594     /**

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

  1596      *

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

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

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

  1600      * number from 0 to size-1. Slice 0 is the highest-order bits, and slice

  1601      * size-1 are the lowest-order bits. Slice 0 may be of different size than

  1602      * the other slices.

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

  1604      * slices to the appropriate position when multiplying different-sized

  1605      * numbers.

  1606      */

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

  1608                                     int fullsize)

  1609     {

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

  1611 

  1612         len = mag.length;

  1613         offset = fullsize - len;

  1614 

  1615         if (slice == 0)

  1616         {

  1617             start = 0 - offset;

  1618             end = upperSize - 1 - offset;

  1619         }

  1620         else

  1621         {

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

  1623             end = start + lowerSize - 1;

  1624         }

  1625 

  1626         if (start < 0)

  1627             start = 0;

  1628         if (end < 0)

  1629            return ZERO;

  1630 

  1631         sliceSize = (end-start) + 1;

  1632 

  1633         if (sliceSize <= 0)

  1634             return ZERO;

  1635 

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

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

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

  1639             return this.abs();

  1640 

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

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

  1643 

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

  1645     }

  1646 

  1647     /**

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

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

  1650      * multiplication.  This is an efficient algorithm that runs in linear

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

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

  1653      * arguments only.

  1654      */

  1655     private BigInteger exactDivideBy3()

  1656     {

  1657         int len = mag.length;

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

  1659         long x, w, q, borrow;

  1660         borrow = 0L;

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

  1662         {

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

  1664             w = x - borrow;

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

  1666                 borrow = 1L;

  1667             else

  1668                 borrow = 0L;

  1669 

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

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

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

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

  1674             result[i] = (int) q;

  1675 

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

  1677             // eliminated if the first fails.

  1678             if (q >= 0x55555556L)

  1679             {

  1680                 borrow++;

  1681                 if (q >= 0xAAAAAAABL)

  1682                     borrow++;

  1683             }

  1684         }

  1685         result = trustedStripLeadingZeroInts(result);

  1686         return new BigInteger(result, signum);

  1687     }

  1688 

  1689     /**

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

  1691      * This is used by Karatsuba multiplication and Karatsuba squaring.

  1692      */

  1693     private BigInteger getLower(int n) {

  1694         int len = mag.length;

  1695 

  1696         if (len <= n)

  1697             return this;

  1698 

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

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

  1701 

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

  1703     }

  1704 

  1705     /**

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

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

  1708      * Karatsuba squaring.

  1709      */

  1710     private BigInteger getUpper(int n) {

  1711         int len = mag.length;

  1712 

  1713         if (len <= n)

  1714             return ZERO;

  1715 

  1716         int upperLen = len - n;

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

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

  1719 

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

  1721     }

  1722 

  1723     // Squaring

  1724 

  1725     /**

  1815      * Squares a BigInteger using the Karatsuba squaring algorithm.  It should

  1816      * be used when both numbers are larger than a certain threshold (found

  1817      * experimentally).  It is a recursive divide-and-conquer algorithm that

  1818      * has better asymptotic performance than the algorithm used in

  1819      * squareToLen.

  1820      */

  1821     private BigInteger squareKaratsuba()

  1822     {

  1823         int half = (mag.length+1) / 2;

  1824 

  1825         BigInteger xl = getLower(half);

  1826         BigInteger xh = getUpper(half);

  1827 

  1828         BigInteger xhs = xh.square();  // xhs = xh^2

  1829         BigInteger xls = xl.square();  // xls = xl^2

  1830 

  1831         // xh^2 << 64  +  (((xl+xh)^2 - (xh^2 + xl^2)) << 32) + xl^2

  1832         return xhs.shiftLeft(half*32).add(xl.add(xh).square().subtract(xhs.add(xls))).shiftLeft(half*32).add(xls);

  1833     }

  1834 

  1835     /**

  1836      * Squares a BigInteger using the 3-way Toom-Cook squaring algorithm.  It

  1837      * should be used when both numbers are larger than a certain threshold

  1838      * (found experimentally).  It is a recursive divide-and-conquer algorithm

  1839      * that has better asymptotic performance than the algorithm used in

  1840      * squareToLen or squareKaratsuba.

  1841      */

  1842     private BigInteger squareToomCook3()

  1843     {

  1844         int len = mag.length;

  1845 

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

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

  1848 

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

  1850         int r = len - 2*k;

  1851 

  1852         // Obtain slices of the numbers. a2 is the most significant

  1853         // bits of the number, and a0 the least significant.

  1854         BigInteger a0, a1, a2;

  1855         a2 = getToomSlice(k, r, 0, len);

  1856         a1 = getToomSlice(k, r, 1, len);

  1857         a0 = getToomSlice(k, r, 2, len);

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

  1859 

  1860         v0 = a0.square();

  1861         da1 = a2.add(a0);

  1862         vm1 = da1.subtract(a1).square();

  1863         da1 = da1.add(a1);

  1864         v1 = da1.square();

  1865         vinf = a2.square();

  1866         v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();

  1867 

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

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

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

  1871            implemented as right shifts which are relatively efficient, leaving

  1872            only a division by 3.

  1873            The division by 3 is done by an optimized algorithm for this case.

  1874         */

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

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

  1877         t1 = v1.subtract(v0);

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

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

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

  1881         tm1 = tm1.subtract(t2);

  1882 

  1883         // Number of bits to shift left.

  1884         int ss = k*32;

  1885 

  1886         return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);

  1887     }

  1888 

  1889     // Division

  1890 

  1891     /**

  3338