目录
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.
因此，第二次修正议商无误。
由此，该算法无误。
扫一扫
2252
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
