java大数运算详解【其五】大数乘法之平方算法之ToomCook3平方算法

目录

java大数运算详解【其一】大数加减法

java大数运算详解【其二】大数乘法

java大数运算详解【其三】大数乘法之平方算法之按位二次展开式算法

java大数运算详解【其四】大数乘法之平方算法之Karatsuba平方算法

java大数运算详解【其五】大数乘法之平方算法之ToomCook3平方算法

java大数运算详解【其六】大数乘法之单位乘法和经典乘法

java大数运算详解【其七】大数乘法之Karatsuba乘法和ToomCook3乘法

java大数运算详解【其八】大数除法

java大数运算详解【其九】大数除法之试商法（Knuth除法）核心算法

java大数运算详解【其十】大数除法之Burnikel-Ziegler除法算法

    1.3、ToomCook3平方算法（三分展开式算法）
    /**
     * 使用3路Toom-Cook平方算法平方一个大整数。
     * 当两个数字都大于某一阈值时(实验发现)，应该使用它。
     * 它是一种递归分治算法，其渐近性优于squareToLen和squareKaratsuba算法。
     */
    private BigInteger squareToomCook3() {
        int len = mag.length;
        // k是低阶切片的大小(以int表示)。
        int k = (len+2)/3;   // 等于 ceil(largest/3)
        // r是最高阶片的大小(以int为单位)。
        int r = len - 2*k;
        // 获取这些数字的切片。a2是数字中最重要的位（高位），a0是最不重要的位（低位）。
        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();//v0=a0.s
        da1 = a2.add(a0);//da1=a0+a2
        vm1 = da1.subtract(a1).square();//vm1=(a0+a2-a1).s
        da1 = da1.add(a1);//da1=a0+a1+a2
        v1 = da1.square();//v1=(a0+a1+a2).s
        vinf = a2.square();//vinf=a2.s
        v2 = da1.add(a2).shiftLeft(1).subtract(a0).square();//v2=(2*(a0+a1+2*a2)-a0).s
        // 这个算法需要两个2的除法和一个3的除法。
        // 所有的除法都是精确的，也就是说，它们不产生余数，所有的结果都是肯定的。
        // 以2为单位的除法是按右移执行的，这样比较有效，只剩下以3为单位的除法。
        // 3的除法是由这种情况下的优化算法完成的。
        t2 = v2.subtract(vm1).exactDivideBy3();//t2=((2*(a0+a1+2*a2)-a0).s-(a0+a2-a1).s)/3
        tm1 = v1.subtract(vm1).shiftRight(1);//tm1=((a0+a1+a2).s-(a0+a2-a1).s)/2
        t1 = v1.subtract(v0);//t1=(a0+a1+a2).s-a0.s
        t2 = t2.subtract(t1).shiftRight(1);//t2=(((2*(a0+a1+2*a2)-a0).s-(a0+a2-a1).s)/3-((a0+a1+a2).s-a0.s))/2
        t1 = t1.subtract(tm1).subtract(vinf);//t1=(a0+a1+a2).s-a0.s-((a0+a1+a2).s-(a0+a2-a1).s)/2-a2.s
        t2 = t2.subtract(vinf.shiftLeft(1));//t2=(((2*(a0+a1+2*a2)-a0).s-(a0+a2-a1).s)/3-((a0+a1+a2).s-a0.s))/2-2*a2.s
        tm1 = tm1.subtract(t2);//tm1=((a0+a1+a2).s-(a0+a2-a1).s)/2-((((2*(a0+a1+2*a2)-a0).s-(a0+a2-a1).s)/3-((a0+a1+a2).s-a0.s))/2-2*a2.s)
        // 向左移位的位数。
        int ss = k*32;
        //return=(((((((a2.s<<ss)+t2)<<ss)+t1)<<ss)+tm1)<<ss)+a0.s
        return vinf.shiftLeft(ss).add(t2).shiftLeft(ss).add(t1).shiftLeft(ss).add(tm1).shiftLeft(ss).add(v0);
    }
    现证明算法的正确性：
    this=(((a2<<ss)+a1)<<ss)+a0,记：1<<(4*ss)=G,1<<(3*ss)=H,1<<(2*ss)=M,1<<ss=L，则：
    this*this=this<2>=(a2*M+a1*L+a0)<2>=a2<2>*G+2*a2*a1*H+(a1<2>+2*a2*a0)*M+2*a1*a0*L+a0<2>.
    t2=(((2*(a0+a1+2*a2)-a0)<2>-(a0+a2-a1)<2>)/3-((a0+a1+a2)<2>-a0<2>))/2-2*a2<2>
      =(((a0+2*a1+4*a2)<2>-(a0+a2-a1)<2>)/3-(a1<2>+a2<2>+2*a0*a1+2*a1*a2+2*a0*a2))/2-2*a2<2>
      =((3*a1<2>+15*a2<2>+6*a0*a1+18*a1*a2+6*a2*a0)/3-(a1<2>+a2<2>+2*a0*a1+2*a1*a2+2*a0*a2))/2-2*a2<2>
      =(4*a2<2>+4*a1*a2)/2-2*a2<2>
      =2*a1*a2.
    t1=(a0+a1+a2)<2>-a0<2>-((a0+a1+a2)<2>-(a0+a2-a1)<2>)/2-a2<2>
      =(a1<2>+a2<2>+2*a0*a1+2*a1*a2+2*a0*a2)-4*a1(a0+a2)/2-a2<2>
      =a1<2>+2*a0*a2.
    tm1=((a0+a1+a2)<2>-(a0+a2-a1)<2>)/2-((((2*(a0+a1+2*a2)-a0)<2>-(a0+a2-a1)<2>)/3-((a0+a1+a2)<2>-a0<2>))/2-2*a2<2>)
       =2*a1*(a0+a2)-t2
       =2*a0*a1.
    展开式的每一项都正确，故该算法正确无误。
    这里把this分成了三部分，用到了五次平方，之所以不用更为简单的算式：a1<2>+2*a0*a2=(a0+a1+a2)<2>-a0<2>-a2<2>-2*a1*(a0+a2)，是因为平方运算要比乘法快。
    /**
     * 返回一个用于Toom-Cook乘法的大整数切片。
     *
     * @param lowerSize 低阶位片的长度。
     * @param upperSize 高阶位片的长度。
     * @param slice 请求哪个切片的索引，该索引必须是从0到size-1的数字。
     * 切片0是最高阶位，切片size-1是最低阶位。切片0的大小可能与其他切片不同。
     * @param fullsize 较大整数数组的大小，用于在乘不同大小的数时将片对齐到适当的位置。
     */
    private BigInteger getToomSlice(int lowerSize, int upperSize, int slice,
                                    int fullsize) {
        int start, end, sliceSize, len, offset;
        len = mag.length;
        offset = fullsize - len;//offset>=0
        /**
         * 注：将数组切分时会依照较长的数组的长度确定lowerSize和upperSize的值，
         * 若需切分较短的数组，此时offset值将为正，得到的高阶位片的长度将小于upperSize.
         */
        if (slice == 0) {
            start = 0 - offset;
            end = upperSize - 1 - offset;
        } else {
            start = upperSize + (slice-1)*lowerSize - offset;
            end = start + lowerSize - 1;
        }
        if (start < 0) {//执行到此处，说明不正常切分，切片长度值不为upperSize或lowerSize.
            start = 0;
        }
        if (end < 0) {
           return ZERO;
        }
        sliceSize = (end-start) + 1;
        if (sliceSize <= 0) {
            /**
             * 这里是多余的，因为执行到此处，start>=0，有：
             * （1）start>0,切片长度值为upperSize或lowerSize.
             * （2）start==0,分以下几种情况：
             *     1、offset==0，则slice==0,end =upperSize-1,sliceSize=upperSize>0,正常切分。
             *     2、offset>0，则只需考虑不正常切分的情况：
             *         当end>=0时，sliceSize>=1；
             *         当end<0时，在此前将会返回，不会执行到此处。
             */
            return ZERO;
        }
        // 在执行Toom-Cook时，所有的切片都是正数，当最终数字组成时，将调整符号。
        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);
    }
    /**
     * 对指定的数字进行3的精确除法(即，余数已知为零)。
     * 这是在图姆-库克乘法中使用的。这是一个在线性时间内有效运行的算法。
     * 如果参数不能被3整除，则结果是未定义的。请注意，这只会以符号为正的参数调用。
     */
    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;
            /**
             * 注意，这里进行了减法，w可能为负（此时，产生借位1）。
             * 这种情况下，borrow > x.
             */
            if (borrow > x) {      // 我们把数字变成负数了吗?
                borrow = 1L;
            } else {
                borrow = 0L;
            }
            // 0xAAAAAAAB是3 (mod 2^32)的模块化倒数。
            // 因此,这是除以3的效果(mod 2^32)。
            // 这比大多数架构上的除法要快得多。
            q = (w * 0xAAAAAAABL) & LONG_MASK;
            result[i] = (int) q;
            // 现在检查借位。如果第一个检查失败，第二个检查当然可以取消。
            if (q >= 0x55555556L) {
                borrow++;//借位1
                if (q >= 0xAAAAAAABL)//借位2
                    borrow++;
            }
        }
        result = trustedStripLeadingZeroInts(result);
        return new BigInteger(result, signum);
    }
    该算法的原理是除以一个数等价于乘以一个数的倒数，而0x100000000/3=0x55555555，0xAAAAAAAB=0x55555555*2+1.
    现证明该算法的正确性：
    考虑计算整数w/3的情况（w<2^32），我们假设可以向更高的位借位，那么我们只需计算得到末32位的整型值和借位值即可，分别设为q,borrow，有：
    q=(borrow*2^32+w)/3，根据对3整除的性质，可以知道borrow=0或1或2即可。
    当borrow=0时，w可以整除3，那么：
        w*0xAAAAAAAB=w*(0x55555555*2+1)=(w/3)*(0x55555555*2*3+3)=(w/3)*(0xFFFFFFFF*2+3)=(w/3)*((0xFFFFFFFF+1)*2+1)=(w/3)*(0x200000000+1).
        (int)(w*0xAAAAAAAB)=w/3.
        w/3<0x100000000/3=0x55555555+1/3，即：w/3<=0x55555555.
    当borrow=1时，w mod 3=2，那么：
        w*0xAAAAAAAB=w*(0x55555555*2+1)=(((w+1)/3)*3-1)*(0x55555555*2+1)=((w+1)/3)*(0x200000000+1)-(0x55555555*2+1).
        (int)(w*0xAAAAAAAB)=(w+1)/3-(0x55555555*2+1)+0x100000000=(w+1)/3+0x55555555.
        而：(w+0x100000000)/3=(w+1+0xFFFFFFFF)/3=(w+1)/3+0x55555555>=0x55555556.
        且：(w+1)/3+0x55555555<0xAAAAAAAA+2/3，即：(w+1)/3+0x55555555<=0xAAAAAAAA.
    当borrow=2时，w mod 3=1，那么：
        w*0xAAAAAAAB=w*(0x55555555*2+1)=(((w-1)/3)*3+1)*(0x55555555*2+1)=((w-1)/3)*(0x200000000+1)+(0x55555555*2+1).
        (int)(w*0xAAAAAAAB)=(w-1)/3+(0x55555555*2+1)=(w-1)/3+0xAAAAAAAB.
        而：(w+0x200000000)/3=(w-1)/3+0xAAAAAAAB>=0xAAAAAAAB.
        且：(w-1)/3+0xAAAAAAAB<0x55555555+0xAAAAAAAB=0x100000000，即：(w-1)/3+0xAAAAAAAB<=0xFFFFFFFF.
    故该程序计算无误。
    /**
     * 返回值为{@code (this >> n)}的大型整数。
     * 右移距离{@code n}，可能为负, 那种情况下，要做正确的左移运算。
     */
    public BigInteger shiftRight(int n) {
        if (signum == 0)
            return ZERO;
        if (n > 0) {
            return shiftRightImpl(n);
        } else if (n == 0) {
            return this;
        } else {
            // (-n)中可能的int溢出不是问题，
            // 因为shiftRightImpl认为它的参数是无符号的。
            return new BigInteger(shiftLeft(mag, -n), signum);
        }
    }
    ToomCook3平方算法，将mag数组分为三份，再进行乘法运算。因而，又称为三分展开式算法。