java大數運算詳解【其十】大數除法之Burnikel-Ziegler除法演算法
目錄
java大數運算詳解【其三】大數乘法之平方演算法之按位二次展開式演算法
java大數運算詳解【其四】大數乘法之平方演算法之Karatsuba平方演算法
java大數運算詳解【其五】大數乘法之平方演算法之ToomCook3平方演算法
java大數運算詳解【其七】大數乘法之Karatsuba乘法和ToomCook3乘法
java大數運算詳解【其九】大數除法之試商法（Knuth除法）核心演算法
java大數運算詳解【其十】大數除法之Burnikel-Ziegler除法演算法
2、Burnikel-Ziegler除法演算法（分塊迴圈帶餘數試商法演算法，簡稱分塊試商法）
/**
*使用Burnikel-Ziegler演算法計算{@code this / val}.
* @param val 除數
* @return {@code this / val}
*/
private BigInteger divideBurnikelZiegler(BigInteger val) {
return divideAndRemainderBurnikelZiegler(val)[0];
}
/**
* 使用Burnikel-Ziegler演算法計算{@code this / val}和{@code this % val}.
* @param val 除數
* @return 一個包含商和餘數的陣列
*/
private BigInteger[] divideAndRemainderBurnikelZiegler(BigInteger val) {
MutableBigInteger q = new MutableBigInteger();
MutableBigInteger r = new MutableBigInteger(this).divideAndRemainderBurnikelZiegler(new MutableBigInteger(val), q);
BigInteger qBigInt = q.isZero() ? ZERO : q.toBigInteger(signum*val.signum);
BigInteger rBigInt = r.isZero() ? ZERO : r.toBigInteger(signum);
return new BigInteger[] {qBigInt, rBigInt};
}
核心演算法：
/**
* 使用Burnikel-Ziegler演算法計算{@code this/b}和{@code this%b}.
* <a href="http://cr.yp.to/bib/1998/burnikel.ps"> Burnikel-Ziegler algorithm</a>.
* 該方法的實現來自Burnikel-Ziegler論文第9頁中的演算法3。
* 引數β(beta)被選為2 ^ 32，所以幾乎所有移位都是32位的倍數。
* {@code this}和{@code b}必須是非負的。
* @param b 除數
* @param quotient 輸出{@code this/b}的引數
* @return 餘數
*/
MutableBigInteger divideAndRemainderBurnikelZiegler(MutableBigInteger b, MutableBigInteger quotient) {
int r = intLen;
int s = b.intLen;
// 清除商
quotient.offset = quotient.intLen = 0;
if (r < s) {
return this;
} else {
/**
* 分塊的大小是由s,即b.intLen來決定的。
* 在該演算法的執行過程中，會繼續將這些塊進行分割，
* 最終分割成長度（以整型長度為單位）不大於BURNIKEL_ZIEGLER_THRESHOLD的資料來進行計算。
* 因此，希望每一個塊都能分割為m段，m為2的冪次，每個段最多BURNIKEL_ZIEGLER_THRESHOLD個整型資料。
* 如此一來，需先確定m，才能確定分塊的大小。
*/
// 與Knuth除法不同，我們這裡不檢查2的公共冪，因為BZ已經執行得更快了，
// 如果兩個數都包含2的冪，然而取消它們沒有額外的好處。
// 步驟 1: 使 m = min{2^k | (2^k)*BURNIKEL_ZIEGLER_THRESHOLD > s}
int m = 1 << (32-Integer.numberOfLeadingZeros(s/BigInteger.BURNIKEL_ZIEGLER_THRESHOLD));//BURNIKEL_ZIEGLER_THRESHOLD=80
int j = (s+m-1) / m; // 步驟 2a: j = ceil(s/m)
int n = j * m; // 步驟 2b: 以無符號整型長度為單位，塊的長度
long n32 = 32L * n; // 以bit為單位，塊的長度
int sigma = (int) Math.max(0, n32 - b.bitLength()); // 步驟 3: sigma = max{T | (2^T)*B < beta^n}，n32必定不小於b.bitLength()，故有多餘操作
MutableBigInteger bShifted = new MutableBigInteger(b);
bShifted.safeLeftShift(sigma); // 步驟 4a: 左移b使它的長度是n的倍數，即有n32個bit
MutableBigInteger aShifted = new MutableBigInteger (this);
aShifted.safeLeftShift(sigma); // 步驟 4b: 同步左移a
// 步驟 5: t是容納加一個額外的位所需的塊數
int t = (int) ((aShifted.bitLength()+n32) / n32);
if (t < 2) {
t = 2;
}
// 步驟 6: 概念上分割成t個塊：a[t-1], ..., a[0]
MutableBigInteger a1 = aShifted.getBlock(t-1, t, n); // a中最重要的塊，在最高位上的塊
// 步驟 7: z[t-2] = [a[t-1], a[t-2]]
MutableBigInteger z = aShifted.getBlock(t-2, t, n); // 次重要的塊
z.addDisjoint(a1, n); // z[t-2]
// 做教科書式的除法，2個塊的資料除以1個塊的資料
MutableBigInteger qi = new MutableBigInteger();
MutableBigInteger ri;
for (int i=t-2; i > 0; i--) {
// 步驟 8a: 計算 (qi,ri) 以滿足 z=b*qi+ri
ri = z.divide2n1n(bShifted, qi);
// 步驟 8b: z = [ri, a[i-1]]
z = aShifted.getBlock(i-1, t, n); // a[i-1]
z.addDisjoint(ri, n);
quotient.addShifted(qi, i*n); // 更新商 (見論文的步驟9)
}
// 步驟8的最後一次迴圈：對i=0再做一次迴圈，但保持z不變
ri = z.divide2n1n(bShifted, qi);
quotient.add(qi);
ri.rightShift(sigma); // 步驟 9: a和b同步左移，因此，ri要右移回來以得到餘數
return ri;
}
}
該演算法通過以一個不小於除數長度（儘量接近）的分塊大小，將被除數分為多塊，利用試商法的原理，從高位的塊開始依次試商（使用divide2n1n計算），從而得到結果。
/** @see BigInteger#bitLength() */
long bitLength() {
if (intLen == 0)
return 0;
return intLen*32L - Integer.numberOfLeadingZeros(value[offset]);
}
/**
* 如{@link #leftShift(int)}但是{@code n}可以是零。
*/
void safeLeftShift(int n) {
if (n > 0) {
leftShift(n);
}
}
/**
* 返回{@code MutableBigInteger}，
* 來自於{@code this}中的開始於{@code index*blockLength}的包含{@code blockLength}長度的整型資料,
* 被用於Burnikel-Ziegler除法。
* @param index 塊的索引，從0到numBlocks-1，對應{@code this}中從低位到高位的資料。
* @param numBlocks 在{@code this}中的塊的總數
* @param blockLength 以無符號整型的32位長度為單位，一個塊的長度。
* @return
*/
private MutableBigInteger getBlock(int index, int numBlocks, int blockLength) {
int blockStart = index * blockLength;
if (blockStart >= intLen) {
return new MutableBigInteger();
}
int blockEnd;
if (index == numBlocks-1) {
blockEnd = intLen;
} else {
blockEnd = (index+1) * blockLength;
}
if (blockEnd > intLen) {
return new MutableBigInteger();
}
int[] newVal = Arrays.copyOfRange(value, offset+intLen-blockEnd, offset+intLen-blockStart);
return new MutableBigInteger(newVal);
}
/**
* 如{@link #addShifted(MutableBigInteger, int)}，
* 但是{@code this.intLen}必須不大於{@code n}，否則會從addend複製被截斷的資料到result中。
* 換句話說,用來連線{@code this}和{@code addend}.
*/
void addDisjoint(MutableBigInteger addend, int n) {
if (addend.isZero())
return;
int x = intLen;
int y = addend.intLen + n;
int resultLen = (intLen > y ? intLen : y);//max(x,y)
int[] result;
if (value.length < resultLen)
result = new int[resultLen];
else {
result = value;
Arrays.fill(value, offset+intLen, value.length, 0);//？多餘的操作
}
int rstart = result.length-1;
/**
* 把this.value中有效資料複製到result陣列的末尾，複製起始位置為result.length-this.intLen.
*/
// 如果需要則複製this中的資料
System.arraycopy(value, offset, result, rstart+1-x, x);
y -= x;
rstart -= x;
int len = Math.min(y, addend.value.length-addend.offset);
/**
* 把addend.value中從addend.offset開始的資料複製到result陣列的尾部，
* 複製起始位置為result.length-(addend.intLen+n)，且不覆蓋上一次複製的資料。
*/
System.arraycopy(addend.value, addend.offset, result, rstart+1-y, len);
// 兩次複製的資料的間隙置零
for (int i=rstart+1-y+len; i < rstart+1; i++)
result[i] = 0;
value = result;
intLen = resultLen;
offset = result.length - resultLen;
}
/**
* 該方法實現了來自Burnikel-Ziegler論文的第4頁的演算法1。
* 它將一個長度為2n的資料除以一個長度為n的資料。
* 引數β(beta)被選為2 ^ 32，所以幾乎所有移位都是32位的倍數。
* {@code this}必須是一個非負數，那樣還需滿足{@code this.bitLength() <= 2*b.bitLength()}.
* @param b 一個正數以使{@code b.bitLength()}足夠。
* @param quotient 用於輸出{@code this/b}的引數。
* @return {@code this%b}
*/
private MutableBigInteger divide2n1n(MutableBigInteger b, MutableBigInteger quotient) {
int n = b.intLen;
// 步驟 1: 基本情況
if (n%2 != 0 || n < BigInteger.BURNIKEL_ZIEGLER_THRESHOLD) {
return divideKnuth(b, quotient);
}
// 步驟 2: 把this視為[a1,a2,a3,a4]，其中，每個ai都是n/2個整型的資料
MutableBigInteger aUpper = new MutableBigInteger(this);
aUpper.safeRightShift(32*(n/2)); // aUpper = [a1,a2,a3]
keepLower(n/2); // this = a4
// 步驟 3: q1=aUpper/b, r1=aUpper%b
MutableBigInteger q1 = new MutableBigInteger();
MutableBigInteger r1 = aUpper.divide3n2n(b, q1);
// 步驟 4: quotient=[r1,this]/b, r2=[r1,this]%b
addDisjoint(r1, n/2); // this = [r1,this]
MutableBigInteger r2 = divide3n2n(b, quotient);
// 步驟 5: 使quotient=[q1,quotient]，並且返回r2
quotient.addDisjoint(q1, n/2);
return r2;
}
該方法利用了除法堅式的原理，將其分解為兩個divide3n2n的除法。
/**
* 如{@link #rightShift(int)}但是{@code n}可以大於value陣列的長度。
*/
void safeRightShift(int n) {
if (n/32 >= intLen) {
reset();
} else {
rightShift(n);
}
}
/**
* 將一個MutableBigInteger置為零，移除它的偏移量。
*/
void reset() {
offset = intLen = 0;
}
/**
* Discards all ints whose index is greater than {@code n}.
*/
private void keepLower(int n) {
if (intLen >= n) {
offset += intLen - n;
intLen = n;
}
}
/**
* 該方法實現了來自Burnikel-Ziegler論文的第5頁的演算法2。
* 它將一個長度為3n的資料除以一個長度為2n的資料。
* 引數β(beta)被選為2 ^ 32，所以幾乎所有移位都是32位的倍數。
* {@code this}必須是一個非負數，那樣還需滿足{@code 2*this.bitLength() <= 3*b.bitLength()}.
* @param quotient 用於輸出{@code this/b}的引數。
* @return {@code this%b}
*/
private MutableBigInteger divide3n2n(MutableBigInteger b, MutableBigInteger quotient) {
int n = b.intLen / 2; // b的整型資料的個數的一半
// 步驟 1: 把this視為[a1,a2,a3]，其中，每個ai都是n個（或更少）整型的資料；並使a12=[a1,a2]
MutableBigInteger a12 = new MutableBigInteger(this);
a12.safeRightShift(32*n);
// 步驟 2: 把b視為[b1,b2]，其中，每個bi都是n個（或更少）整型的資料
MutableBigInteger b1 = new MutableBigInteger(b);
b1.safeRightShift(n * 32);
BigInteger b2 = b.getLower(n);
MutableBigInteger r;
MutableBigInteger d;
if (compareShifted(b, n) < 0) {
// 步驟 3a: 若a1<b1，使quotient=a12/b1且r=a12%b1
r = a12.divide2n1n(b1, quotient);
// 步驟 4: d=quotient*b2
d = new MutableBigInteger(quotient.toBigInteger().multiply(b2));//注意此處，Burnikel-Ziegler除法演算法的時間複雜度主要取決於這一步所消耗的時間
} else {
// 步驟 3b: 若a1>=b1，使quotient=beta^n-1且r=a12-b1*2^n+b1（？註釋有誤，r=a12-b1*beta^n+b1）
quotient.ones(n);
a12.add(b1);
b1.leftShift(32*n);
a12.subtract(b1);
r = a12;
// 步驟 4: d=quotient*b2=(b2 << 32*n) - b2
d = new MutableBigInteger(b2);
d.leftShift(32 * n);
d.subtract(new MutableBigInteger(b2));
}
// 步驟 5: r = r*beta^n + a3 - d (論文說是a4)（？意義未明）
// 然而，不減去d直到迴圈結束，所以r不會變成負數
r.leftShift(32 * n);
r.addLower(this, n);
// 步驟 6: 加b，直到r>=d為止
while (r.compare(d) < 0) {
r.add(b);
quotient.subtract(MutableBigInteger.ONE);
}
r.subtract(d);
return r;
}
該方法利用了試商法的原理，使用divide2n1n得到議商，再修正議商。
/**
* 返回一個{@code BigInteger}，它等於this的低位上的{@code n}個整型資料構成的值。
*/
private BigInteger getLower(int n) {
if (isZero()) {
return BigInteger.ZERO;
} else if (intLen < n) {
return toBigInteger(1);
} else {
// 跳過前面的零
int len = n;
while (len > 0 && value[offset+intLen-len] == 0)
len--;
int sign = len > 0 ? 1 : 0;
return new BigInteger(Arrays.copyOfRange(value, offset+intLen-len, offset+intLen), sign);
}
}
/**
* 返回true，如果this MutableBigInteger值為零。
*/
boolean isZero() {
return (intLen == 0);
}
/**
* 返回一個等價於{@code b.leftShift(32*ints); return compare(b);}的值，
* 但是並不改變{@code b}的值。
*/
private int compareShifted(MutableBigInteger b, int ints) {
int blen = b.intLen;
int alen = intLen - ints;
if (alen < blen)
return -1;
if (alen > blen)
return 1;
// 加上Integer.MIN_VALUE使之作為無符號比較
// 從高位開始，迴圈比較
int[] bval = b.value;
for (int i = offset, j = b.offset; i < alen + offset; i++, j++) {
int b1 = value[i] + 0x80000000;
int b2 = bval[j] + 0x80000000;
if (b1 < b2)
return -1;
if (b1 > b2)
return 1;
}
return 0;
}
/**
* 使this成為{@code n}個整型資料組成的MutableBigInteger，它的所有位都是1。
* 被用於Burnikel-Ziegler除法。
* @param n {@code value}陣列中整型資料的個數。
* @return 一個值為{@code ((1<<(32*n)))-1}的MutableBigInteger
*/
private void ones(int n) {
if (n > value.length)
value = new int[n];
Arrays.fill(value, -1);
offset = 0;
intLen = n;
}
/**
* 將兩個MutableBigInteger物件的內容相加。
* 結果放在this MutableBigInteger內。加數的內容沒有更改。
*/
void add(MutableBigInteger addend) {
int x = intLen;
int y = addend.intLen;
int resultLen = (intLen > addend.intLen ? intLen : addend.intLen);
int[] result = (value.length < resultLen ? new int[resultLen] : value);
int rstart = result.length-1;
long sum;
long carry = 0;
// 公共部分相加
while(x > 0 && y > 0) {
x--; y--;
sum = (value[x+offset] & LONG_MASK) +
(addend.value[y+addend.offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
// 加到較長的部分上
while(x > 0) {
x--;
if (carry == 0 && result == value && rstart == (x + offset))
return;
sum = (value[x+offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
while(y > 0) {
y--;
sum = (addend.value[y+addend.offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
if (carry > 0) { // 長度上結果必須增長
resultLen++;
if (result.length < resultLen) {
int temp[] = new int[resultLen];
// 結果溢位，增長了一個整型的長度，把低位上的資料複製到新的結果中
System.arraycopy(result, 0, temp, 1, result.length);
temp[0] = 1;
result = temp;
} else {
result[rstart--] = 1;
}
}
value = result;
intLen = resultLen;
offset = result.length - resultLen;
}
/**
* 從較大的數中減去較小的數，並將結果放在this MutableBigInteger內。
*/
int subtract(MutableBigInteger b) {
MutableBigInteger a = this;
int[] result = value;
int sign = a.compare(b);
if (sign == 0) {
reset();
return 0;
}
if (sign < 0) {
MutableBigInteger tmp = a;
a = b;
b = tmp;
}
int resultLen = a.intLen;
if (result.length < resultLen)
result = new int[resultLen];
long diff = 0;
int x = a.intLen;
int y = b.intLen;
int rstart = result.length - 1;
// 公共部分相減
while (y > 0) {
x--; y--;
diff = (a.value[x+a.offset] & LONG_MASK) -
(b.value[y+b.offset] & LONG_MASK) - ((int)-(diff>>32));
result[rstart--] = (int)diff;
}
// 減到較長的部分上
while (x > 0) {
x--;
diff = (a.value[x+a.offset] & LONG_MASK) - ((int)-(diff>>32));
result[rstart--] = (int)diff;
}
value = result;
intLen = resultLen;
offset = value.length - resultLen;
normalize();
return sign;
}
/**
* 加上{@code addend}末尾的{@code n}個整型資料。
*/
void addLower(MutableBigInteger addend, int n) {
MutableBigInteger a = new MutableBigInteger(addend);
if (a.offset + a.intLen >= n) {
a.offset = a.offset + a.intLen - n;
a.intLen = n;
}
a.normalize();
add(a);
}
/**
* 比較兩個MutableBigIntegers的大小，返回 -1, 0 or 1，
* 作為this MutableBigInteger數值上小於、等於、大於{@code b}的值。
*/
final int compare(MutableBigInteger b) {
int blen = b.intLen;
if (intLen < blen)
return -1;
if (intLen > blen)
return 1;
// 加上Integer.MIN_VALUE使之作為無符號比較
// 從高位開始，迴圈比較
int[] bval = b.value;
for (int i = offset, j = b.offset; i < intLen + offset; i++, j++) {
int b1 = value[i] + 0x80000000;
int b2 = bval[j] + 0x80000000;
if (b1 < b2)
return -1;
if (b1 > b2)
return 1;
}
return 0;
}
/**
* 將{@code addend}的值左移{@code n}個整型位的加到this上。
* 具有與{@code addend.leftShift(32*ints); add(addend);}相同的效果；但是不改變{@code addend}的值。
*/
void addShifted(MutableBigInteger addend, int n) {
if (addend.isZero()) {
return;
}
int x = intLen;
int y = addend.intLen + n;
int resultLen = (intLen > y ? intLen : y);
int[] result = (value.length < resultLen ? new int[resultLen] : value);
int rstart = result.length-1;
long sum;
long carry = 0;
// 公共部分相加
while (x > 0 && y > 0) {
x--; y--;
int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
sum = (value[x+offset] & LONG_MASK) +
(bval & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
// 加到較長的部分上
while (x > 0) {
x--;
if (carry == 0 && result == value && rstart == (x + offset)) {
return;
}
sum = (value[x+offset] & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
while (y > 0) {
y--;
int bval = y+addend.offset < addend.value.length ? addend.value[y+addend.offset] : 0;
sum = (bval & LONG_MASK) + carry;
result[rstart--] = (int)sum;
carry = sum >>> 32;
}
if (carry > 0) { // 長度上結果必須增長
resultLen++;
if (result.length < resultLen) {
int temp[] = new int[resultLen];
// 結果溢位，增長了一個整型的長度，把低位上的資料複製到新的結果中
System.arraycopy(result, 0, temp, 1, result.length);
temp[0] = 1;
result = temp;
} else {
result[rstart--] = 1;
}
}
value = result;
intLen = resultLen;
offset = result.length - resultLen;
}
Burnikel-Ziegler除法演算法，又稱為分塊迴圈帶餘數試商法演算法，簡稱分塊試商法。它將被除數按除數的長度進行分塊，再採用試商法計算，使用 2個塊的資料除以1個塊的資料 得到議商，由於除數長度為一個塊，故不需修正議商。
對於 2n個整型的資料除以n個整型的資料 ，將之化為兩個 3n個整型的資料除以2n個整型的資料 ，再採用試商法計算，從而化為 2n個整型的資料除以n個整型的資料 得到議商，再修正。
Burnikel-Ziegler除法演算法的時間複雜度主要取決於其所使用的大數乘法的時間複雜度，因為分塊試商（不包括其採用的乘法）的時間複雜度是線性的。
相關文章
- java大數運算詳解【其九】大數除法之試商法（Knuth除法）核心演算法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之數值運算
- 巨量算數：2020十大消費新機遇
- 【演算法分析與設計】輾轉相除法演算法
- spark RDD textFile運算元 分割槽數量詳解Spark
- 輾轉相除法原理
- 乘法◆除法◆謎題(7)
- 乘法◆除法◆謎題(11)
- 乘法◆除法◆謎題(9)
- 乘法◆除法◆謎題(8)
- 乘法◆除法◆謎題(6)
- 乘法◆除法◆謎題(5)
- 乘法◆除法◆謎題(4)
- 乘法◆除法◆謎題(3)