目录
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除法算法的时间复杂度主要取决于其所使用的大数乘法的时间复杂度，因为分块试商（不包括其采用的乘法）的时间复杂度是线性的。
511
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
