我正在编写一个 BigInt 类作为编程练习。它在 base-65536 中使用 2 的补码有符号整数的 vector (这样 32 位乘法就不会溢出。一旦我完全正常工作,我将增加基数)。
所有基本数学运算都经过编码,但有一个问题:使用我能够创建的基本算法,除法痛苦地很慢。 (它有点像商的每个数字的二进制除法......我不会发布它,除非有人想看到它......)
我想使用 Newton-Raphson 来找到(移位的)倒数,然后相乘(和移位),而不是我的慢速算法。我想我已经掌握了基础知识:你给公式 (x1 = x0(2 - x0 * divisor)) 一个很好的初始猜测,然后经过一些迭代后,x 收敛到互惠的。这部分看起来很简单......但是在尝试将此公式应用于大整数时我遇到了一些问题:
问题 1:
因为我使用的是整数……好吧……我不能使用分数。这似乎导致 x 总是发散(x0 * 除数似乎必须 <2?)。我的直觉告诉我应该对方程式进行一些修改,使其可以计算整数(达到一定的准确性),但我真的很难找出它是什么。 (我缺乏数学技能在这里打败了我......)我想我需要找到一些等价的方程式,而不是 d 有 d*[base^somePower]强>?是否有像 (x1 = x0(2 - x0 * d)) 这样适用于整数的方程式?
问题 2:
当我使用牛顿公式求出一些数字的倒数时,结果最终只是比答案应该是的小一小部分……例如。当试图找到 4 的倒数(十进制)时:
x0 = 0.3
x1 = 0.24
x2 = 0.2496
x3 = 0.24999936
x4 = 0.2499999999983616
x5 = 0.24999999999999999999998926258176
如果我以 10 为基数表示数字,我会想要 25 的结果(并记住将乘积右移 2)。对于一些倒数,例如 1/3,您可以在知道自己有足够的准确性后简单地截断结果。但是如何从上面的结果中提取出正确的倒数呢?
抱歉,如果这一切都太模糊或者我要求太多。我浏览了 Wikipedia 和我可以在 Google 上找到的所有研究论文,但我觉得我正在用头撞墙。我感谢任何人能给我的帮助!
...
编辑:使算法正常工作,尽管它比我预期的要慢得多。与我的旧算法相比,我实际上失去了很多速度,即使是在有数千位数字的数字上......我仍然遗漏了一些东西。这不是乘法的问题,它非常快。 (我确实在使用 Karatsuba 的算法)。
对于任何感兴趣的人,这是我当前的 Newton-Raphson 算法迭代:
bigint operator/(const bigint& lhs, const bigint& rhs) {
if (rhs == 0) throw overflow_error("Divide by zero exception");
bigint dividend = lhs;
bigint divisor = rhs;
bool negative = 0;
if (dividend < 0) {
negative = !negative;
dividend.invert();
}
if (divisor < 0) {
negative = !negative;
divisor.invert();
}
int k = dividend.numBits() + divisor.numBits();
bigint pow2 = 1;
pow2 <<= k + 1;
bigint x = dividend - divisor;
bigint lastx = 0;
bigint lastlastx = 0;
while (1) {
x = (x * (pow2 - x * divisor)) >> k;
if (x == lastx || x == lastlastx) break;
lastlastx = lastx;
lastx = x;
}
bigint quotient = dividend * x >> k;
if (dividend - (quotient * divisor) >= divisor) quotient++;
if (negative)quotient.invert();
return quotient;
}
这是我的(非常丑陋的)旧算法,它更快:
bigint operator/(const bigint& lhs, const bigint & rhs) {
if (rhs == 0) throw overflow_error("Divide by zero exception");
bigint dividend = lhs;
bigint divisor = rhs;
bool negative = 0;
if (dividend < 0) {
negative = !negative;
dividend.invert();
}
if (divisor < 0) {
negative = !negative;
divisor.invert();
}
bigint remainder = 0;
bigint quotient = 0;
while (dividend.value.size() > 0) {
remainder.value.insert(remainder.value.begin(), dividend.value.at(dividend.value.size() - 1));
remainder.value.push_back(0);
remainder.unPad();
dividend.value.pop_back();
if (divisor > remainder) {
quotient.value.push_back(0);
} else {
int count = 0;
int i = MSB;
bigint value = 0;
while (i > 0) {
bigint increase = divisor * i;
bigint next = value + increase;
if (next <= remainder) {
value = next;
count += i;
}
i >>= 1;
}
quotient.value.push_back(count);
remainder -= value;
}
}
for (int i = 0; i < quotient.value.size() / 2; i++) {
int swap = quotient.value.at(i);
quotient.value.at(i) = quotient.value.at((quotient.value.size() - 1) - i);
quotient.value.at(quotient.value.size() - 1 - i) = swap;
}
if (negative)quotient.invert();
quotient.unPad();
return quotient;
}
最佳答案
首先,您可以在时间 O(n^2)
中实现除法并使用合理的常量,因此它不会比原始乘法慢(很多)。但是,如果您使用 Karatsuba -like 算法,甚至 FFT -based 乘法算法,那么你确实可以使用 Newton-Raphson 加速你的除法算法。
用于计算 x
倒数的 Newton-Raphson 迭代是 q[n+1]=q[n]*(2-q[n]*x)
.
假设我们要计算floor(2^k/B)
其中B
是一个正整数。 WLOG, B≤2^k
;否则,商为 0
。 x=B/2^k
的 Newton-Raphson 迭代产生 q[n+1]=q[n]*(2-q[n]*B/2^k)
。我们可以将其重新排列为
q[n+1]=q[n]*(2^(k+1)-q[n]*B) >> k
这种类型的每次迭代只需要整数乘法和位移。它是否收敛于 floor(2^k/B)
?不必要。然而,在最坏的情况下,它最终会在 floor(2^k/B)
和 ceiling(2^k/B)
之间交替(证明它!)。因此,您可以使用一些不太聪明的测试来查看您是否属于这种情况,然后提取 floor(2^k/B)
。 (这个“不太聪明的测试”应该比每次迭代中的乘法快很多;但是,优化这个东西会很好)。
事实上,计算 floor(2^k/B)
足以计算任何正整数 floor(A/B)
A,B
。取 k
使得 A*B≤2^k
,并验证 floor(A/B)=A*ceiling(2^k/B) >> k
.
最后,此方法的一个简单但重要的优化是在 Newton-Raphson 方法的早期迭代中截断乘法(即仅计算乘积的较高位)。之所以这样做,是因为早期迭代的结果与商数相去甚远,执行不准确也没有关系。 (改进这个论点并证明如果你适本地做这件事,你可以在 O(M(2n))
时间内划分两个 ≤n
位整数,假设你可以在时间M(k)
上乘以两个≤k
位的整数,M(x)
是一个递增的凸函数。
关于c++ - 具有大整数的牛顿-拉夫森除法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27801397/