java大數運算詳解【其九】大數除法之試商法（Knuth除法）核心演算法
目錄
java大數運算詳解【其三】大數乘法之平方演算法之按位二次展開式演算法
java大數運算詳解【其四】大數乘法之平方演算法之Karatsuba平方演算法
java大數運算詳解【其五】大數乘法之平方演算法之ToomCook3平方演算法
java大數運算詳解【其七】大數乘法之Karatsuba乘法和ToomCook3乘法
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.
因此，第二次修正議商無誤。
由此，該演算法無誤。
相關文章
- java大數運算詳解【其十】大數除法之Burnikel-Ziegler除法演算法Java演算法
- 大數運算—大數加法、減法、乘法、除法詳解
- 計算機計算小數除法的陷阱計算機
- PAT-B 1017 A除以B【模擬 大數除法】
- 利用位移法進行除法運算
- java大整數四則運算Java
- JavaScript / 除法運算子JavaScript
- C++除法運算 // 靜態斷言C++
- [除數函式+除法分塊]函式
- 01—Hbuilder—js—除法、取餘數UIJS
- 【leetcode】29. Divide Two Integers 不能使用乘除法的整數除法LeetCodeIDE
- 例項九— 除法器設計
- (C語言) int型之間除法運算，向零取整C語言
- 求最大公約數（輾轉相除法）
- 定點數的原碼補碼除法
- Java 基礎 之 算數運算子Java
- C語言高效程式設計的四大祕技之使用位操作，減少除法和取模的運算C語言程式設計
- NumPy 舍入小數、對數、求和和乘積運算詳解
- 除法與GCD演算法的相關分析GC演算法
- 輾轉相除法
- 除法謎題(2)
- 除法謎題(1)
- 多項式除法
- android kotlin 安全除法AndroidKotlin
- scss不能用除法?CSS
- SHELL之數值運算
- 【演算法分析與設計】輾轉相除法演算法
- oracle rac 核心引數詳解Oracle
- spark RDD textFile運算元 分割槽數量詳解Spark
- 輾轉相除法原理
- 乘法◆除法◆謎題(7)
- 乘法◆除法◆謎題(11)
- 乘法◆除法◆謎題(9)
- 乘法◆除法◆謎題(8)
- 乘法◆除法◆謎題(6)
- 乘法◆除法◆謎題(5)
- 乘法◆除法◆謎題(4)
- 乘法◆除法◆謎題(3)