具有大整数的牛顿-拉夫森除法
Newton-Raphson Division With Big Integers
我正在制作一个 BigInt class 作为编程练习。它在 base-65536 中使用 2 的补码有符号整数的向量（这样 32 位乘法就不会溢出。一旦我完全正常工作，我将增加基数）。
所有基本数学运算都经过编码，但有一个问题：使用我能够创建的基本算法，除法非常痛苦 很慢。 （它有点像商的每个数字的二进制除法......我不会post它除非有人想看到它......）
我想使用 Newton-Raphson 来找到（移位的）倒数，然后乘法（和移位），而不是我的慢速算法。我想我对基础知识有所了解：你给公式 (x1 = x0(2 - x0 * divisor)) 一个很好的初始猜测，然后经过一些迭代后，x收敛于倒数。这部分看起来很简单......但是我在尝试将这个公式应用于大整数时遇到了一些问题：
问题一：
因为我正在处理整数...嗯...我不会使用分数。这似乎导致 x 总是发散（x0 * 除数似乎必须 <2？）。我的直觉告诉我应该对方程式进行一些修改，使其可以计算整数（达到一定的准确性），但我真的很难找出它是什么。 （我缺乏数学技能在这里打败了我......）我想我需要找到一些等价的方程式而不是 d 有 d*[base ^somePower]？是否可以使用像 (x1 = x0(2 - x0 * d)) 这样的等式来处理整数？
问题二：
当我使用牛顿公式求一些数字的倒数时，结果只是比答案应该是下面的一小部分...例如。当试图找到 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的算法）。
对于任何感兴趣的人，这是我当前的牛顿-拉夫森算法迭代：
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;
}
Newton-Raphson 是一种近似算法 - 不适用于整数数学。您将得到舍入误差，这将导致您所看到的问题。你可以用浮点数来做这道题，然后看看你是否得到一个整数，精确到指定的位数（见下一段）
关于第二个问题，选择您想要的精度（小数位数）并四舍五入到该精度。如果您在问题中选择了 20 位精度，您将四舍五入为 0.25。您只需要迭代，直到所需的精度位数稳定为止。一般来说，在计算机上表示无理数通常会导致不精确。
首先，您可以在时间上实现除法 O(n^2) 并使用合理的常数，因此它不会比简单的乘法慢（很多）。但是，如果您使用基于 Karatsuba-like algorithm, or even FFT 的乘法算法，那么您确实可以使用 Newton-Raphson 加速除法算法。
用于计算 x 的倒数的 Newton-Raphson 迭代是 q[n+1]=q[n]*(2-q[n]*x).
假设我们要计算 floor(2^k/B) 其中 B 是一个正整数。博客，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)。 （这个 "not-so-clever test" 应该比每次迭代中的乘法快很多；但是，优化这个东西会很好）。
实际上，计算 floor(2^k/B) 足以计算任何正整数 A,B 的 floor(A/B)。取 k 使得 A*B≤2^k，并验证 floor(A/B)=A*ceiling(2^k/B) >> k.
最后，此方法的一个简单但重要的优化是在 Newton-Raphson 方法的早期迭代中截断乘法（即仅计算乘积的较高位）。之所以这样做，是因为早期迭代的结果与商数相去甚远，执行不准确也没有关系。 （细化这个论点并表明，如果你适当地做这件事，你可以在时间 O(M(2n)) 中除以两个 ≤n 位整数，假设你可以在时间中乘以两个 ≤k 位整数 M(k)，而M(x)是一个递增的凸函数）。
如果我没看错的话，一个重大改进就是为 x 选择一个好的起始值。知道除数有多少位你就知道逆的最高有效位在哪里，如
1/x = pow(2,log2(1/x))
1/x = pow(2,-log2(x))
1/x >= pow(2,-floor(log2(x)))
floor(log2(x)) 只是最高有效位集的索引。
正如 op 在评论中所建议的那样，使用 256 位查找 table 将进一步加快收敛速度，因为每一步大约使正确数字的数量翻倍。从 8 个正确的数字开始比从 1 开始好，甚至比从更少的数字开始要好得多。
constexpr fixpoint_integer_inverse(const T& d) {
uint8_t lut[256] = { 255u,254u,253u,252u,251u,250u,249u,248u,247u,246u,245u,244u,243u,242u,241u,
240u,240u,239u,238u,237u,236u,235u,234u,234u,233u,232u,231u,230u,229u,229u,228u,
227u,226u,225u,225u,224u,223u,222u,222u,221u,220u,219u,219u,218u,217u,217u,216u,
215u,214u,214u,213u,212u,212u,211u,210u,210u,209u,208u,208u,207u,206u,206u,205u,
204u,204u,203u,202u,202u,201u,201u,200u,199u,199u,198u,197u,197u,196u,196u,195u,
195u,194u,193u,193u,192u,192u,191u,191u,190u,189u,189u,188u,188u,187u,187u,186u,
186u,185u,185u,184u,184u,183u,183u,182u,182u,181u,181u,180u,180u,179u,179u,178u,
178u,177u,177u,176u,176u,175u,175u,174u,174u,173u,173u,172u,172u,172u,171u,171u,
170u,170u,169u,169u,168u,168u,168u,167u,167u,166u,166u,165u,165u,165u,164u,164u,
163u,163u,163u,162u,162u,161u,161u,161u,160u,160u,159u,159u,159u,158u,158u,157u,
157u,157u,156u,156u,156u,155u,155u,154u,154u,154u,153u,153u,153u,152u,152u,152u,
151u,151u,151u,150u,150u,149u,149u,149u,148u,148u,148u,147u,147u,147u,146u,146u,
146u,145u,145u,145u,144u,144u,144u,144u,143u,143u,143u,142u,142u,142u,141u,141u,
141u,140u,140u,140u,140u,139u,139u,139u,138u,138u,138u,137u,137u,137u,137u,136u,
136u,136u,135u,135u,135u,135u,134u,134u,134u,134u,133u,133u,133u,132u,132u,132u,
132u,131u,131u,131u,131u,130u,130u,130u,130u,129u,129u,129u,129u,128u,128u,128u,
127u
};
const auto l = log2(d);
T x;
if (l<8) {
x = T(1)<<(digits(d)-1-l);
} else {
if (digits(d)>(l+8)) x = T(lut[(d>>(l-8))-256])<<(digits(d)-l-8);
else x = T(lut[(d>>(l-8))-256])>>(l+8-digits(d));
}
if (x==0) x=1;
while(true) {
const auto lm = long_mul(x,T(1)-x*d);
const T i = get<0>(lm);
if (i) x+=i;
else return x;
}
return x;
}
// calculate a * b = r0r1
template<typename T>
typename std::enable_if<std::is_unsigned<T>::value,tuple<T,T>>::type
constexpr long_mul(const T& a, const T& b){
const T N = digits<T>()/2;
const T t0 = (a>>N)*(b>>N);
const T t1 = ((a<<N)>>N)*(b>>N);
const T t2 = (a>>N)*((b<<N)>>N);
const T t3 = ((a<<N)>>N)*((b<<N)>>N);
const T t4 = t3+(t1<<N);
const T r1 = t4+(t2<<N);
const T r0 = (r1<t4)+(t4<t3)+(t1>>N)+(t2>>N)+t0;
return {r0,r1};
}
我正在制作一个 BigInt class 作为编程练习。它在 base-65536 中使用 2 的补码有符号整数的向量（这样 32 位乘法就不会溢出。一旦我完全正常工作，我将增加基数）。
所有基本数学运算都经过编码，但有一个问题：使用我能够创建的基本算法，除法非常痛苦 很慢。 （它有点像商的每个数字的二进制除法......我不会post它除非有人想看到它......）
我想使用 Newton-Raphson 来找到（移位的）倒数，然后乘法（和移位），而不是我的慢速算法。我想我对基础知识有所了解：你给公式 (x1 = x0(2 - x0 * divisor)) 一个很好的初始猜测，然后经过一些迭代后，x收敛于倒数。这部分看起来很简单......但是我在尝试将这个公式应用于大整数时遇到了一些问题：
问题一：
因为我正在处理整数...嗯...我不会使用分数。这似乎导致 x 总是发散（x0 * 除数似乎必须 <2？）。我的直觉告诉我应该对方程式进行一些修改，使其可以计算整数（达到一定的准确性），但我真的很难找出它是什么。 （我缺乏数学技能在这里打败了我......）我想我需要找到一些等价的方程式而不是 d 有 d*[base ^somePower]？是否可以使用像 (x1 = x0(2 - x0 * d)) 这样的等式来处理整数？
问题二：
当我使用牛顿公式求一些数字的倒数时，结果只是比答案应该是下面的一小部分...例如。当试图找到 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的算法）。
对于任何感兴趣的人，这是我当前的牛顿-拉夫森算法迭代：
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;
}
Newton-Raphson 是一种近似算法 - 不适用于整数数学。您将得到舍入误差，这将导致您所看到的问题。你可以用浮点数来做这道题，然后看看你是否得到一个整数，精确到指定的位数（见下一段）
关于第二个问题，选择您想要的精度（小数位数）并四舍五入到该精度。如果您在问题中选择了 20 位精度，您将四舍五入为 0.25。您只需要迭代，直到所需的精度位数稳定为止。一般来说，在计算机上表示无理数通常会导致不精确。
首先，您可以在时间上实现除法 O(n^2) 并使用合理的常数，因此它不会比简单的乘法慢（很多）。但是，如果您使用基于 Karatsuba-like algorithm, or even FFT 的乘法算法，那么您确实可以使用 Newton-Raphson 加速除法算法。
用于计算 x 的倒数的 Newton-Raphson 迭代是 q[n+1]=q[n]*(2-q[n]*x).
假设我们要计算 floor(2^k/B) 其中 B 是一个正整数。博客，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)。 （这个 "not-so-clever test" 应该比每次迭代中的乘法快很多；但是，优化这个东西会很好）。
实际上，计算 floor(2^k/B) 足以计算任何正整数 A,B 的 floor(A/B)。取 k 使得 A*B≤2^k，并验证 floor(A/B)=A*ceiling(2^k/B) >> k.
最后，此方法的一个简单但重要的优化是在 Newton-Raphson 方法的早期迭代中截断乘法（即仅计算乘积的较高位）。之所以这样做，是因为早期迭代的结果与商数相去甚远，执行不准确也没有关系。 （细化这个论点并表明，如果你适当地做这件事，你可以在时间 O(M(2n)) 中除以两个 ≤n 位整数，假设你可以在时间中乘以两个 ≤k 位整数 M(k)，而M(x)是一个递增的凸函数）。
如果我没看错的话，一个重大改进就是为 x 选择一个好的起始值。知道除数有多少位你就知道逆的最高有效位在哪里，如
1/x = pow(2,log2(1/x))
1/x = pow(2,-log2(x))
1/x >= pow(2,-floor(log2(x)))
floor(log2(x)) 只是最高有效位集的索引。
正如 op 在评论中所建议的那样，使用 256 位查找 table 将进一步加快收敛速度，因为每一步大约使正确数字的数量翻倍。从 8 个正确的数字开始比从 1 开始好，甚至比从更少的数字开始要好得多。
constexpr fixpoint_integer_inverse(const T& d) {
uint8_t lut[256] = { 255u,254u,253u,252u,251u,250u,249u,248u,247u,246u,245u,244u,243u,242u,241u,
240u,240u,239u,238u,237u,236u,235u,234u,234u,233u,232u,231u,230u,229u,229u,228u,
227u,226u,225u,225u,224u,223u,222u,222u,221u,220u,219u,219u,218u,217u,217u,216u,
215u,214u,214u,213u,212u,212u,211u,210u,210u,209u,208u,208u,207u,206u,206u,205u,
204u,204u,203u,202u,202u,201u,201u,200u,199u,199u,198u,197u,197u,196u,196u,195u,
195u,194u,193u,193u,192u,192u,191u,191u,190u,189u,189u,188u,188u,187u,187u,186u,
186u,185u,185u,184u,184u,183u,183u,182u,182u,181u,181u,180u,180u,179u,179u,178u,
178u,177u,177u,176u,176u,175u,175u,174u,174u,173u,173u,172u,172u,172u,171u,171u,
170u,170u,169u,169u,168u,168u,168u,167u,167u,166u,166u,165u,165u,165u,164u,164u,
163u,163u,163u,162u,162u,161u,161u,161u,160u,160u,159u,159u,159u,158u,158u,157u,
157u,157u,156u,156u,156u,155u,155u,154u,154u,154u,153u,153u,153u,152u,152u,152u,
151u,151u,151u,150u,150u,149u,149u,149u,148u,148u,148u,147u,147u,147u,146u,146u,
146u,145u,145u,145u,144u,144u,144u,144u,143u,143u,143u,142u,142u,142u,141u,141u,
141u,140u,140u,140u,140u,139u,139u,139u,138u,138u,138u,137u,137u,137u,137u,136u,
136u,136u,135u,135u,135u,135u,134u,134u,134u,134u,133u,133u,133u,132u,132u,132u,
132u,131u,131u,131u,131u,130u,130u,130u,130u,129u,129u,129u,129u,128u,128u,128u,
127u
};
const auto l = log2(d);
T x;
if (l<8) {
x = T(1)<<(digits(d)-1-l);
} else {
if (digits(d)>(l+8)) x = T(lut[(d>>(l-8))-256])<<(digits(d)-l-8);
else x = T(lut[(d>>(l-8))-256])>>(l+8-digits(d));
}
if (x==0) x=1;
while(true) {
const auto lm = long_mul(x,T(1)-x*d);
const T i = get<0>(lm);
if (i) x+=i;
else return x;
}
return x;
}
// calculate a * b = r0r1
template<typename T>
typename std::enable_if<std::is_unsigned<T>::value,tuple<T,T>>::type
constexpr long_mul(const T& a, const T& b){
const T N = digits<T>()/2;
const T t0 = (a>>N)*(b>>N);
const T t1 = ((a<<N)>>N)*(b>>N);
const T t2 = (a>>N)*((b<<N)>>N);
const T t3 = ((a<<N)>>N)*((b<<N)>>N);
const T t4 = t3+(t1<<N);
const T r1 = t4+(t2<<N);
const T r0 = (r1<t4)+(t4<t3)+(t1>>N)+(t2>>N)+t0;
return {r0,r1};
}