java大數運算詳解【其九】大數除法之試商法(Knuth除法)核心演算法

一念之卓發表於2020-12-17

目錄

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.
        因此,第二次修正議商無誤。
    由此,該演算法無誤。

 

相關文章