【发布时间】：2015-03-04 06:45:33
【问题描述】：
我正在制作一个 BigInt 类作为编程练习。它使用以 65536 为基数的 2 的补码有符号整数向量（这样 32 位乘法就不会溢出。一旦它完全工作，我将增加基数）。
所有基本的数学运算都已编码，但有一个问题：使用我能够创建的基本算法，除法非常痛苦。 （它有点像商的每个数字的二进制除法......除非有人想看到它，否则我不会发布它......）
我想使用 Newton-Raphson 来找到（移位的）倒数，然后乘以（和移位），而不是我的慢速算法。我想我对基础知识有所了解：你给公式 (x1 = x0(2 - x0 * divisor)) 一个很好的初始猜测，然后经过一些迭代，x 收敛到互惠的。这部分似乎很容易......但是在尝试将此公式应用于大整数时遇到了一些问题：
问题一：
因为我正在使用整数...嗯...我不能使用分数。这似乎导致 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;
}
【问题讨论】：
-
your solution returns
1instead of2for2/1¶ 你认为你找到了解决方案，你可以post it as your own answer （答案应该作为答案发布，而不是问题更新）。 -
这是working (in my tests)
unsigned_div_newton()implementation in Python (text in Russian)。对于我尝试过的情况，基于长除法 (unsigned_div_long()) 的实现要快得多。
标签： c++ algorithm integer division bigint