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

目录

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

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

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

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

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

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

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

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

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

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

    核心算法：
    /**
     * 该MutableBigInteger除以除数div。
     * 商将被放置到提供的quotient对象中并将余数对象返回（当needRemainder值true）。
     */
    private MutableBigInteger divideMagnitude(MutableBigInteger div,
                                              MutableBigInteger quotient,
                                              boolean needRemainder ) {
        // 断言 div.intLen > 1
        // D1 规范化除数
        int shift = Integer.numberOfLeadingZeros(div.value[div.offset]);
        // 复制除数值以保护除数
        final int dlen = div.intLen;
        int[] divisor;
        MutableBigInteger rem; // 余数始于带着前导零的被除数
        if (shift > 0) {
            divisor = new int[dlen];
            copyAndShift(div.value,div.offset,dlen,divisor,0,shift);
            if (Integer.numberOfLeadingZeros(value[offset]) >= shift) {
                int[] remarr = new int[intLen + 1];
                rem = new MutableBigInteger(remarr);
                rem.intLen = intLen;
                rem.offset = 1;
                copyAndShift(value,offset,intLen,remarr,1,shift);
            } else {
                int[] remarr = new int[intLen + 2];
                rem = new MutableBigInteger(remarr);
                rem.intLen = intLen+1;
                rem.offset = 1;
                int rFrom = offset;
                // 相当于内联copyAndShift
                int c=0;
                int n2 = 32 - shift;
                for (int i=1; i < intLen+1; i++,rFrom++) {
                    int b = c;
                    c = value[rFrom];
                    remarr[i] = (b << shift) | (c >>> n2);
                }
                remarr[intLen+1] = c << shift;
            }
        } else {
            divisor = Arrays.copyOfRange(div.value, div.offset, div.offset + div.intLen);
            rem = new MutableBigInteger(new int[intLen + 1]);
            System.arraycopy(value, offset, rem.value, 1, intLen);
            rem.intLen = intLen;
            rem.offset = 1;
        }
        int nlen = rem.intLen;
        // 设置商的长度
        final int limit = nlen - dlen + 1;
        if (quotient.value.length < limit) {
            quotient.value = new int[limit];
            quotient.offset = 0;
        }
        quotient.intLen = limit;
        int[] q = quotient.value;
        // 如果长度不变，必须在rem中插入前导零
        if (rem.intLen == nlen) {//此处必定执行，故不必判断
            rem.offset = 0;
            rem.value[0] = 0;
            rem.intLen++;
        }
        int dh = divisor[0];
        long dhLong = dh & LONG_MASK;
        int dl = divisor[1];
        // D2 初始化 j（？意义未明）
        for (int j=0; j < limit-1; j++) {
            // D3 计算 qhat（？意义未明）
            // 评估 qhat（？意义未明）
            int qhat = 0;
            int qrem = 0;
            boolean skipCorrection = false;
            int nh = rem.value[j+rem.offset];
            int nh2 = nh + 0x80000000;
            int nm = rem.value[j+1+rem.offset];
            if (nh == dh) {
                qhat = ~0;//值为0xffffffff
                qrem = nh + nm;
                skipCorrection = qrem + 0x80000000 < nh2;//相当于(qrem & LONG_MASK) < (nh & LONG_MASK)
            } else {
                long nChunk = (((long)nh) << 32) | (nm & LONG_MASK);
                if (nChunk >= 0) {
                    qhat = (int) (nChunk / dhLong);
                    qrem = (int) (nChunk - (qhat * dhLong));
                } else {
                    long tmp = divWord(nChunk, dh);
                    qhat = (int) (tmp & LONG_MASK);
                    qrem = (int) (tmp >>> 32);
                }
            }
            if (qhat == 0)
                continue;
            if (!skipCorrection) { // 修正 qhat（？意义未明）
                long nl = rem.value[j+2+rem.offset] & LONG_MASK;
                long rs = ((qrem & LONG_MASK) << 32) | nl;
                long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
                if (unsignedLongCompare(estProduct, rs)) {
                    qhat--;
                    //((qrem&LONG_MASK)+dhLong)>>32==0则还需修正
                    qrem = (int)((qrem & LONG_MASK) + dhLong);
                    if ((qrem & LONG_MASK) >=  dhLong) {
                        estProduct -= (dl & LONG_MASK);
                        rs = ((qrem & LONG_MASK) << 32) | nl;
                        if (unsignedLongCompare(estProduct, rs))
                            qhat--;
                    }
                }
            }
            // D4 相乘并减去
            rem.value[j+rem.offset] = 0;
            int borrow = mulsub(rem.value, divisor, qhat, dlen, j+rem.offset);
            // D5 测试余数
            if (borrow + 0x80000000 > nh2) {
                // D6 相加回复
                divadd(divisor, rem.value, j+1+rem.offset);
                qhat--;
            }
            // 存储商的数据
            q[j] = qhat;
        } // D7 循环于 j
        // D3 计算 qhat
        // 评估 qhat
        int qhat = 0;
        int qrem = 0;
        boolean skipCorrection = false;
        int nh = rem.value[limit - 1 + rem.offset];
        int nh2 = nh + 0x80000000;
        int nm = rem.value[limit + rem.offset];
        if (nh == dh) {
            qhat = ~0;
            qrem = nh + nm;
            skipCorrection = qrem + 0x80000000 < nh2;
        } else {
            long nChunk = (((long) nh) << 32) | (nm & LONG_MASK);
            if (nChunk >= 0) {
                qhat = (int) (nChunk / dhLong);
                qrem = (int) (nChunk - (qhat * dhLong));
            } else {
                long tmp = divWord(nChunk, dh);
                qhat = (int) (tmp & LONG_MASK);
                qrem = (int) (tmp >>> 32);
            }
        }
        if (qhat != 0) {
            if (!skipCorrection) { // 修正 qhat
                long nl = rem.value[limit + 1 + rem.offset] & LONG_MASK;
                long rs = ((qrem & LONG_MASK) << 32) | nl;
                long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
                if (unsignedLongCompare(estProduct, rs)) {
                    qhat--;
                    //((qrem&LONG_MASK)+dhLong)>>32==0则还需修正
                    qrem = (int) ((qrem & LONG_MASK) + dhLong);
                    if ((qrem & LONG_MASK) >= dhLong) {
                        estProduct -= (dl & LONG_MASK);
                        rs = ((qrem & LONG_MASK) << 32) | nl;
                        if (unsignedLongCompare(estProduct, rs))
                            qhat--;
                    }
                }
            }
            // D4 相乘并减去
            int borrow;
            rem.value[limit - 1 + rem.offset] = 0;
            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);
            // D5 测试余数
            if (borrow + 0x80000000 > nh2) {
                // D6 相加回复
                if(needRemainder)
                    divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
                qhat--;
            }
            // 存储商的数据
            q[(limit - 1)] = qhat;
        }
        if (needRemainder) {
            // D8 反规范化
            if (shift > 0)
                rem.rightShift(shift);
            rem.normalize();
        }
        quotient.normalize();
        return needRemainder ? rem : null;
    }
    // 把src中从srcFrom开始的srcLen长度的数据（作为value数组）左移shift位再复制到dst的dstFrom处。
    private static void copyAndShift(int[] src, int srcFrom, int srcLen, int[] dst, int dstFrom, int shift) {
        int n2 = 32 - shift;
        int c=src[srcFrom];
        for (int i=0; i < srcLen-1; i++) {
            int b = c;
            c = src[++srcFrom];
            dst[dstFrom+i] = (b << shift) | (c >>> n2);
        }
        dst[dstFrom+srcLen-1] = c << shift;
    }
    /**
     * 作为无符号比较两个长整型数。
     * 若one大于two则返回true.
     */
    private boolean unsignedLongCompare(long one, long two) {
        return (one+Long.MIN_VALUE) > (two+Long.MIN_VALUE);
    }
    /**
     * 这种方法用于除法。
     * 它将数组a乘以整型数x，然后从q中offset处减去其结果。
     */
    private int mulsub(int[] q, int[] a, int x, int len, int offset) {
        long xLong = x & LONG_MASK;
        long carry = 0;
        offset += len;
        for (int j=len-1; j >= 0; j--) {
            long product = (a[j] & LONG_MASK) * xLong + carry;
            long difference = q[offset] - product;//qof=q[offset]
            q[offset--] = (int)difference;
            carry = (product >>> 32)
                     + (((difference & LONG_MASK) >
                         (((~(int)product) & LONG_MASK))) ? 1:0);
            /**
             * 讨论(((difference & LONG_MASK) >
                         (((~(int)product) & LONG_MASK))) ? 1:0)的值，如下：
             * 设product=(ph<<32)+pl，qof=q[offset]=difference+product，根据要求，我们希望其值为qof<pl.
             * 当qof>=pl时，difference=qof-product=-(ph<<32)+qof-pl,difference&LONG_MASK=qof-pl.
                    (~(int)product)&LONG_MASK=(~pl)&LONG_MASK=(-pl-1)&LONG_MASK=(-pl+0xffffffff)&LONG_MASK.
                    故(((difference & LONG_MASK) > (((~(int)product) & LONG_MASK))) ? 1:0)=(qof>0xffffffff?1:0)=0.
             * 当qof<pl时，difference=qof-product=-(ph<<32)+qof-pl,difference&LONG_MASK=qof-pl+0x100000000.
                    (~(int)product)&LONG_MASK=(~pl)&LONG_MASK=(-pl-1)&LONG_MASK=(-pl+0xffffffff)&LONG_MASK.
                    故(((difference & LONG_MASK) > (((~(int)product) & LONG_MASK))) ? 1:0)=(qof+0x100000000>0xffffffff?1:0)=1.
             * 因而计算无误。
             */
        }
        return (int)carry;
    }
    /**
     * 简单用于除法。
     * 该方法在指定的偏移量上将除数的一倍添加回被除数结果中。当qhat估计太大时使用，必须进行调整。
     * 即把数组a添加到数组result的offset处。
     */
    private int divadd(int[] a, int[] result, int offset) {
        long carry = 0;
        for (int j=a.length-1; j >= 0; j--) {
            long sum = (a[j] & LONG_MASK) +
                       (result[j+offset] & LONG_MASK) + carry;
            result[j+offset] = (int)sum;
            carry = sum >>> 32;
        }
        return (int)carry;
    }
    
    对于算法的流程，我们需要明确。该算法全称为循环带余数试商除法算法，简称试商法，其中，“试商”是算法的关键。
    所谓除法，“除”的意思即“减去”，从余数中减去“试商”与除数的乘积，这就构成了一层循环，因此试商法是二重偱环算法。
    要说明该算法流程，先建立几个概念，“余数”“议商”“积商”“余估值”“除估值”“定商”。
    余数，长度与除数相同，但小于除数，算法开始，余数是被除数的前几位。
    议商，被除数下一位到余数，再将余估值除以除估值的结果作为议商，修正议商。
    积商，议商经过修正后，与除数相乘的结果。
    余估值，余数的前两位，修正议商时，余数下一位到余估值，因此余估值变为余数的前三位。
    除估值，除数的首位，修正议商时，除数下一位到除估值，因此除估值娈为除数的前二位。
    定商，余数减去积商，值非负，则其结果作为下一次试商的余数，议商就是定商，否则，再加上除数作为下一次试商的余数，议商自减作为定商，这便是第二次修正议商。
    如此，算法流程十分清晰。
    ㈠为保护除数和被除数，将除数左移复制，使之没有前导零，被除数同步左移复制。
    ㈡确定余数，被除数的前几位，长度同除数。
    ㈢议商，被除数下一位到余数，再将余估值除以除估值的结果作为议商。
    ㈣修正议商，余数下一位到余估值，除数下一位到除估值，根据余估值和除估值修正议商。
    ㈤减去积商，议商与除数相乘得到积商，再从余数中减去。
    ㈥第二次修正议商，得到定商和下一次试商的余数。
    ㈦将定商存储到商数组中。
    ㈧跳转到步骤㈢，直到被除数末位也已经下到余数中。
    ㈨将余数右移，以得到真正的余数，至此算法结束。
    
    为证明divideMagnitude算法的正确性，只需证明其中两次修正议商无误，证明如下：
    先证明步骤㈢得到的议商一定不小于定商：
        步骤㈢得到的议商等价于将除数首位之外的位都置为零之后得到的定商，因而得证。
    第一次修正议商，代码如下：
        ⒈if (!skipCorrection) { // 修正 qhat（？意义未明）
        ⒉   long nl = rem.value[j+2+rem.offset] & LONG_MASK;
        ⒊   long rs = ((qrem & LONG_MASK) << 32) | nl;
        ⒋   long estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);
        ⒌   if (unsignedLongCompare(estProduct, rs)) {
        ⒍       qhat--;
                 //((qrem&LONG_MASK)+dhLong)>>32==0则还需修正
        ⒎       qrem = (int)((qrem & LONG_MASK) + dhLong);
        ⒏       if ((qrem & LONG_MASK) >=  dhLong) {
        ⒐           estProduct -= (dl & LONG_MASK);
        ⒑           rs = ((qrem & LONG_MASK) << 32) | nl;
        ⒒           if (unsignedLongCompare(estProduct, rs))
        ⒓               qhat--;
        ⒔       }
        ⒕   }
        ⒖}
        先明确修正议商的目的，是为了使议商等于余数的前三位除以除数的前两位。只需要证明该算法能达到此目的即可。
        除数的前两位分别为dh,dl，余数的前三位分别为nh,nm,nl，执行到该代码时，qhat=((nh<<32)+nm)/dh,qrem=((nh<<32)+nm)%dh=((nh<<32)+nm)-dh*qhat.
        设qu=((((nh<<32)+nm)<<32)+nl)/((dh<<32)+dl),qr=((((nh<<32)+nm)<<32)+nl)%((dh<<32)+dl)=((((nh<<32)+nm)<<32)+nl)-((dh<<32)+dl)*qu.
        ⒊   rs = ((qrem & LONG_MASK) << 32) | nl;=((((nh<<32)+nm)-dh*qhat)<<32)+nl.
        ⒋   estProduct = (dl & LONG_MASK) * (qhat & LONG_MASK);=dl*qhat.
        ((dh<<32)+dl)*qhat=((dh*qhat)<<32)+dl*qhat.
        ⒌   rs-estProduct=((((nh<<32)+nm)-dh*qhat)<<32)+nl-dl*qhat=(((nh<<32)+nm)<<32)+nl-(((dh*qhat)<<32)+dl*qhat)=((((nh<<32)+nm)<<32)+nl)-((dh<<32)+dl)*qhat.
            若rs-estProduct>=0，则rs-estProduct=((((nh<<32)+nm)<<32)+nl)-((dh<<32)+dl)*qhat<=((((nh<<32)+nm)<<32)+nl)-((dh<<32)+dl)*qu=qr.
                0<=rs-estProduct<=qr<(dh<<32)+dl，故qhat=qr，不需修正。
            若rs-estProduct<0，则rs-estProduct<0<=qr<(dh<<32)+dl，故qhat>qu，需要修正。
        ⒍   qhat--;=((nh<<32)+nm)/dh-1.
        ⒎   qrem = (int)((qrem & LONG_MASK) + dhLong);=(int)((((nh<<32)+nm)-dh*(qhat+1)+dh)&LONG_MASK)=(int)((((nh<<32)+nm)-dh*qhat)&LONG_MASK).
        ⒏   若还需修正议商，即qhat>qu，则((((nh<<32)+nm)<<32)+nl)-((dh<<32)+dl)*qhat<0，
            即((((nh<<32)+nm)-dh*qhat)<<32)+nl<dl*qhat<(1<<64)，故((nh<<32)+nm)-dh*qhat<=LONG_MASK，
            此时，qrem=((nh<<32)+nm)-dh*qhat>dh，故(qrem&LONG_MASK)>=dhLong.
        ⒐   estProduct -= (dl & LONG_MASK);=dl*(qhat+1)-dl=dl*qhat.
        ⒑   rs = ((qrem & LONG_MASK) << 32) | nl;=((((nh<<32)+nm)-dh*qhat)<<32)+nl.
        ⒒   此时情况与第五行相同。
        可以看到qhat最多自减两次，因此还需证明qhat最多自减两次就能达到修正议商的目的。
        执行到该代码时，qhat=((nh<<32)+nm)/dh,qrem=((nh<<32)+nm)%dh=((nh<<32)+nm)-dh*qhat,qu=((((nh<<32)+nm)<<32)+nl)/((dh<<32)+dl).
        qhat-qu=((nh<<32)+nm)/dh-((((nh<<32)+nm)<<32)+nl)/((dh<<32)+dl).
        令f=(decimal)((nh<<32)+nm)/dh-(decimal)((((nh<<32)+nm)<<32)+nl)/((dh<<32)+dl).
           <=(decimal)((nh<<32)+nm)/0x80000000-(decimal)((((nh<<32)+nm)<<32)+nl)/((0x80000000<<32)+dl)
           <=(decimal)((nh<<32)+nm)/0x80000000-(decimal)((((nh<<32)+nm)<<32)+nl)/((0x80000000<<32)+0xffffffff)
           <(decimal)((nh<<32)+nm)/0x80000000-(decimal)((((nh<<32)+nm)<<32)+nl)/(0x80000001<<32)
           <(decimal)((0x80000000<<32)+0xffffffff)/0x80000000-(decimal)(((0x80000000<<32)+0xffffffff)<<32)+nl)/(0x80000001<<32)
           <(decimal)((0x80000000<<32)+0xffffffff)/0x80000000-(decimal)((0x80000000<<32)+0xffffffff)/0x80000001
           <2
        故qhat-qu<=2.
        因此，第一次修正议商无误。
    第二次修正议商，代码如下：
        // D5 测试余数
        if (borrow + 0x80000000 > nh2) {
            // D6 相加回复
            if(needRemainder)
                divadd(divisor, rem.value, limit - 1 + 1 + rem.offset);
            qhat--;
        }
        可以看到qhat最多自减一次，因此还需证明qhat最多自减一次就能达到修正议商的目的。
        令f=(decimal)((((nh<<32)+nm)<<32)+nl)/((dh<<32)+dl)-(decimal)((((((nh<<32)+nm)<<32)+nl)<<32)+nr)/((((dh<<32)+dl)<<32)+dr).
           <1.
        因此，第二次修正议商无误。
    由此，该算法无误。