如何快速将一个大整数分成一个词？
How to quickly divide a big uinteger into a word?
我目前正在开发一个 class 来处理大的无符号整数。但是，我需要不完整的功能，即：
bi_uint+=bi_uint - 已经实施。没有投诉。bi_uint*=std::uint_fast64_t - 已经实施。没有投诉。bi_uint/=std::uint_fast64_t - 已实现但工作起来非常 缓慢 ，还需要一个两倍于 uint_fast64_t 的类型。在测试用例中，除法比乘法慢 35 倍
接下来给出我的除法实现，基于简单的长除法算法：
#include <climits>
#include <cstdint>
#include <limits>
#include <vector>
class bi_uint
{
public:
using u64_t = std::uint_fast64_t;
constexpr static std::size_t u64_bits = CHAR_BIT * sizeof(u64_t);
using u128_t = unsigned __int128;
static_assert(sizeof(u128_t) >= sizeof(u64_t) * 2);
//little-endian
std::vector<u64_t> data;
//User should guarantee data.size()>0 and val>0
void self_div(const u64_t val)
{
auto it = data.rbegin();
if(data.size() == 1) {
*it /= val;
return;
}
u128_t rem = 0;
if(*it < val) {
rem = *it++;
data.pop_back();
}
u128_t r = rem % val;
while(it != data.rend()) {
rem = (r << u64_bits) + *it;
const u128_t q = rem / val;
r = rem % val;
*it++ = static_cast<u64_t>(q);
}
}
};
您可以看到使用了 unsigned __int128 类型，因此，此选项不可移植并且绑定到单个编译器 - GCC 并且还需要 x64 平台。
看了page关于除法算法的文章后，我觉得合适的算法应该是“Newton-Raphson division”。然而，“Newton-Raphson division”算法对我来说似乎很复杂。我想有一个更简单的算法来划分类型“big_uint/uint”，它的性能几乎相同。
问:如何快速将bi_uint分成u64_t?
我有大约 10^6 次迭代，每次迭代都使用列出的所有操作
如果这很容易实现，那么我希望具有可移植性而不使用unsigned __int128。否则，我宁愿放弃可移植性而选择更简单的方法。
编辑 1：
这是一个学术项目，我无法使用第三方库。
在大多数现代 CPU 中，除法确实比乘法慢得多。
参考
https://agner.org/optimize/instruction_tables.pdf
在 Intel Skylake 上 MUL/IMUL 有 3-4 个周期的延迟；而 DIV/IDIV 可能需要 26-90 个周期；比 MUL 慢 7 - 23 倍；所以您的初始基准测试结果并不令人意外。
如果您恰好在 x86 CPU 上，如下面的答案所示，如果这确实是瓶颈，您可以尝试使用 AVX/SSE 指令。基本上你需要依赖特殊说明而不是像 DIV/IDIV.
这样的一般说明
第 1 部分（请参阅下面的第 2 部分）
我设法在我的旧笔记本电脑上将你的分区代码加速 5x 次（甚至在 GodBolt 服务器上加速 7.5x 次）使用 Barrett Reduction，这是一种允许用多次乘法和加法代替一次除法的技术。今天刚刚从 sctracth 实现了整个代码。
如果您愿意，可以直接跳转到我的 post 末尾的代码位置，而无需阅读冗长的描述，因为代码完全 运行nable 无需任何知识或依赖。
下面的代码仅适用于 Intel x64，因为我只使用了 Intel instructions 并且只使用了它们的 64 位变体。当然对于 x32 和其他处理器也可以是 re-written，因为 Barrett 算法是通用的。
为了简短地解释整个 Barrett Reduction pseudo-code 我会用 Python 来写它，因为它是最简单的语言，适合于理解 pseudo-code:
# https://www.nayuki.io/page/barrett-reduction-algorithm
def BarrettR(n, s):
return (1 << s) // n
def BarrettDivMod(x, n, r, s):
q = (x * r) >> s
t = x - q * n
return (q, t) if t < n else (q + 1, t - n)
基本上在上面的伪代码中 BarrettR() 只对同一个除数完成一次（你对整个大整数除法使用相同的 single-word 除数）。 BarrettDivMod() 每次要进行除法或模运算时都使用，基本上给定输入 x 和除数 n 它 returns 元组 (x / n, x % n)，没有别的，但它比常规除法指令更快。
在下面的 C++ 代码中，我实现了 Barrett 的两个相同函数，但做了一些 C++ 特定的优化以使其更快。优化是可能的，因为除数 n 总是 64 位，x 是 128 位，但较高的一半总是小于 n（最后一个假设发生是因为较高的一半在你的大整数除法始终是余数模 n).
Barrett 算法使用非 2 的幂的除数 n，因此不允许使用 0, 1, 2, 4, 8, 16, ... 这样的除数。你可以通过对大整数做正确的 bit-shift 来涵盖除数的这种简单情况，因为除以 2 的幂只是 bit-shift。允许任何其他除数，甚至包括不是 2 的幂的除数。
同样重要的是要注意我的 BarrettDivMod() 只接受严格小于 divisor * 2^64 的股息 x，换句话说，128 位股息的高一半 [=17] =] 应该小于除数。对于您的大整数除法函数，这始终是正确的，因为较高的一半始终是余数模数除数。 x 的这条规则应该由你检查，它在我的 BarrettDivMod() 中仅作为 DEBUG 断言检查，在发布时被删除。
你可以注意到 BarrettDivMod() 有两个大分支，这是同一算法的两个变体，第一个使用 CLang/GCC 仅类型 unsigned __int128，第二个仅使用 64 位指令和因此适用于 MSVC。
我试图针对三个编译器 CLang/GCC/MSVC，但有些 MSVC 版本在使用 Barrett 时速度只提高了 2 倍，而 CLang/GCC 速度提高了 5 倍。也许我在 MSVC 代码中做了一些错误。
您可以看到我使用您的 class bi_uint 对两个版本的代码进行时间测量 - 使用常规除法和使用 Barrett。 重要 请注意，我对您的代码进行了相当大的更改，首先不使用 u128（以便 MSVC 版本编译时没有 u128），其次不修改 data 向量，所以它只读除法并且不存储最终结果（我需要这个 read-only 来 运行 非常快地加速测试，而无需在每次测试迭代时复制 data 向量） .所以你的代码在我的代码片段中被破坏了，它可以't-be复制粘贴直接使用，我只用你的代码来测量速度。
Barrett 归约工作得更快，不仅因为除法比乘法慢，而且因为乘法和加法在现代 CPUs 上都很好地流水线化，现代 CPU 可以执行多个 mul 或 add一个周期内的指令，但前提是这几个 mul/add 不依赖于彼此的结果，换句话说 CPU 可以 运行 在一个周期内并行执行多个指令。据我所知除法不能 运行 并行，因为 CPU 中只有一个模块可以进行除法，但它仍然有点流水线，因为在第一次除法的 50% 之后是可以在 CPU 流水线的开头并行启动完成的第二部分。
在某些计算机上，您可能会注意到常规除法版本有时会慢得多，这是因为 CLang/GCC 回退到 library-based 除法的软实现，即使对于 128 位除数也是如此。在这种情况下，我的 Barrett 甚至可以显示 7-10x 倍的加速，因为它不使用库函数。
为了克服上述问题，关于软除法，最好直接使用 DIV 指令添加 Assembly code，或者在你的编译器中找到一些实现这个的内部函数（我认为 CLang/GCC 有这样的内在）。如果需要，我也可以编写这个 Assembly 实现，请在评论中告诉我。
更新。正如承诺的那样，为 CLang/GCC、函数 UDiv128Asm() 实现了 128 位除法的汇编变体。之后是更改它用作 CLang/GCC 128 分区的主要实现，而不是常规 u128(a) / b。您可以通过在 UDiv128() 函数体内将 #if 0 替换为 #if 1 来恢复常规的 u128 实现。
#include <cstdint>
#include <bit>
#include <stdexcept>
#include <string>
#include <immintrin.h>
#if defined(_MSC_VER) && !defined(__clang__)
#define IS_MSVC 1
#else
#define IS_MSVC 0
#endif
#if IS_MSVC
#include <intrin.h>
#endif
#define ASSERT_MSG(cond, msg) { if (!(cond)) throw std::runtime_error("Assertion (" #cond ") failed at line " + std::to_string(__LINE__) + "! Msg: '" + std::string(msg) + "'."); }
#define ASSERT(cond) ASSERT_MSG(cond, "")
#ifdef _DEBUG
#define DASSERT_MSG(cond, msg) ASSERT_MSG(cond, msg)
#else
#define DASSERT_MSG(cond, msg)
#endif
#define DASSERT(cond) DASSERT_MSG(cond, "")
using u16 = uint16_t;
using u32 = uint32_t;
using i64 = int64_t;
using u64 = uint64_t;
using UllPtr = unsigned long long *;
inline int GetExp(double x) {
return int((std::bit_cast<uint64_t>(x) >> 52) & 0x7FF) - 1023;
}
inline size_t BitSizeWrong(uint64_t x) {
return x == 0 ? 0 : (GetExp(x) + 1);
}
inline size_t BitSize(u64 x) {
size_t r = 0;
if (x >= (u64(1) << 32)) {
x >>= 32;
r += 32;
}
while (x >= 0x100) {
x >>= 8;
r += 8;
}
while (x) {
x >>= 1;
++r;
}
return r;
}
#if !IS_MSVC
inline u64 UDiv128Asm(u64 h, u64 l, u64 d, u64 * r) {
u64 q;
asm (R"(
.intel_syntax
mov rdx, %V[h]
mov rax, %V[l]
div %V[d]
mov %V[r], rdx
mov %V[q], rax
)"
: [q] "=r" (q), [r] "=r" (*r)
: [h] "r" (h), [l] "r" (l), [d] "r" (d)
: "rax", "rdx"
);
return q;
}
#endif
inline std::pair<u64, u64> UDiv128(u64 hi, u64 lo, u64 d) {
#if IS_MSVC
u64 r, q = _udiv128(hi, lo, d, &r);
return std::make_pair(q, r);
#else
#if 0
using u128 = unsigned __int128;
auto const dnd = (u128(hi) << 64) | lo;
return std::make_pair(u64(dnd / d), u64(dnd % d));
#else
u64 r, q = UDiv128Asm(hi, lo, d, &r);
return std::make_pair(q, r);
#endif
#endif
}
inline std::pair<u64, u64> UMul128(u64 a, u64 b) {
#if IS_MSVC
u64 hi, lo = _umul128(a, b, &hi);
return std::make_pair(hi, lo);
#else
using u128 = unsigned __int128;
auto const x = u128(a) * b;
return std::make_pair(u64(x >> 64), u64(x));
#endif
}
inline std::pair<u64, u64> USub128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
u64 r_hi, r_lo;
_subborrow_u64(_subborrow_u64(0, a_lo, b_lo, (UllPtr)&r_lo), a_hi, b_hi, (UllPtr)&r_hi);
return std::make_pair(r_hi, r_lo);
}
inline std::pair<u64, u64> UAdd128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
u64 r_hi, r_lo;
_addcarry_u64(_addcarry_u64(0, a_lo, b_lo, (UllPtr)&r_lo), a_hi, b_hi, (UllPtr)&r_hi);
return std::make_pair(r_hi, r_lo);
}
inline int UCmp128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
if (a_hi != b_hi)
return a_hi < b_hi ? -1 : 1;
return a_lo == b_lo ? 0 : a_lo < b_lo ? -1 : 1;
}
std::pair<u64, size_t> BarrettRS64(u64 n) {
// https://www.nayuki.io/page/barrett-reduction-algorithm
ASSERT_MSG(n >= 3 && (n & (n - 1)) != 0, "n " + std::to_string(n))
size_t const nbits = BitSize(n);
// 2^s = q * n + r; 2^s = (2^64 + q0) * n + r; 2^s - n * 2^64 = q0 * n + r
u64 const dnd_hi = (nbits >= 64 ? 0ULL : (u64(1) << nbits)) - n;
auto const q0 = UDiv128(dnd_hi, 0, n).first;
return std::make_pair(q0, nbits);
}
template <bool Use128 = true, bool Adjust = true>
std::pair<u64, u64> BarrettDivMod64(u64 x_hi, u64 x_lo, u64 n, u64 r, size_t s) {
// ((x_hi * 2^64 + x_lo) * (2^64 + r)) >> (64 + s)
DASSERT(x_hi < n);
#if !IS_MSVC
if constexpr(Use128) {
using u128 = unsigned __int128;
u128 const xf = (u128(x_hi) << 64) | x_lo;
u64 q = u64((u128(x_hi) * r + xf + u64((u128(x_lo) * r) >> 64)) >> s);
if (s < 64) {
u64 t = x_lo - q * n;
if constexpr(Adjust) {
u64 const mask = ~u64(i64(t - n) >> 63);
q += mask & 1;
t -= mask & n;
}
return std::make_pair(q, t);
} else {
u128 t = xf - u128(q) * n;
return t < n ? std::make_pair(q, u64(t)) : std::make_pair(q + 1, u64(t) - n);
}
} else
#endif
{
auto const w1a = UMul128(x_lo, r).first;
auto const [w2b, w1b] = UMul128(x_hi, r);
auto const w2c = x_hi, w1c = x_lo;
u64 w1, w2 = _addcarry_u64(0, w1a, w1b, (UllPtr)&w1);
w2 += _addcarry_u64(0, w1, w1c, (UllPtr)&w1);
w2 += w2b + w2c;
if (s < 64) {
u64 q = (w2 << (64 - s)) | (w1 >> s);
u64 t = x_lo - q * n;
if constexpr(Adjust) {
u64 const mask = ~u64(i64(t - n) >> 63);
q += mask & 1;
t -= mask & n;
}
return std::make_pair(q, t);
} else {
u64 const q = w2;
auto const [b_hi, b_lo] = UMul128(q, n);
auto const [t_hi, t_lo] = USub128(x_hi, x_lo, b_hi, b_lo);
return t_hi != 0 || t_lo >= n ? std::make_pair(q + 1, t_lo - n) : std::make_pair(q, t_lo);
}
}
}
#include <random>
#include <iomanip>
#include <iostream>
#include <chrono>
void TestBarrett() {
std::mt19937_64 rng{123}; //{std::random_device{}()};
for (size_t i = 0; i < (1 << 11); ++i) {
size_t const nbits = rng() % 63 + 2;
u64 n = 0;
do {
n = (u64(1) << (nbits - 1)) + rng() % (u64(1) << (nbits - 1));
} while (!(n >= 3 && (n & (n - 1)) != 0));
auto const [br, bs] = BarrettRS64(n);
for (size_t j = 0; j < (1 << 6); ++j) {
u64 const hi = rng() % n, lo = rng();
auto const [ref_q, ref_r] = UDiv128(hi, lo, n);
u64 bar_q = 0, bar_r = 0;
for (size_t k = 0; k < 2; ++k) {
bar_q = 0; bar_r = 0;
if (k == 0)
std::tie(bar_q, bar_r) = BarrettDivMod64<true>(hi, lo, n, br, bs);
else
std::tie(bar_q, bar_r) = BarrettDivMod64<false>(hi, lo, n, br, bs);
ASSERT_MSG(bar_q == ref_q && bar_r == ref_r, "i " + std::to_string(i) + ", j " + std::to_string(j) + ", k " + std::to_string(k) +
", nbits " + std::to_string(nbits) + ", n " + std::to_string(n) + ", bar_q " + std::to_string(bar_q) +
", ref_q " + std::to_string(ref_q) + ", bar_r " + std::to_string(bar_r) + ", ref_r " + std::to_string(ref_r));
}
}
}
}
class bi_uint
{
public:
using u64_t = std::uint64_t;
constexpr static std::size_t u64_bits = 8 * sizeof(u64_t);
//little-endian
std::vector<u64_t> data;
static auto constexpr DefPrep = [](auto n){
return std::make_pair(false, false);
};
static auto constexpr DefDivMod = [](auto dnd_hi, auto dnd_lo, auto dsr, auto const & prep){
return UDiv128(dnd_hi, dnd_lo, dsr);
};
//User should guarantee data.size()>0 and val>0
template <typename PrepT = decltype(DefPrep), typename DivModT = decltype(DefDivMod)>
void self_div(const u64_t val, PrepT const & Prep = DefPrep, DivModT const & DivMod = DefDivMod)
{
auto it = data.rbegin();
if(data.size() == 1) {
*it /= val;
return;
}
u64_t rem_hi = 0, rem_lo = 0;
if(*it < val) {
rem_lo = *it++;
//data.pop_back();
}
auto const prep = Prep(val);
u64_t r = rem_lo % val;
u64_t q = 0;
while(it != data.rend()) {
rem_hi = r;
rem_lo = *it;
std::tie(q, r) = DivMod(rem_hi, rem_lo, val, prep);
//*it++ = static_cast<u64_t>(q);
it++;
auto volatile out = static_cast<u64_t>(q);
}
}
};
void TestSpeed() {
auto Time = []{
static auto const gtb = std::chrono::high_resolution_clock::now();
return std::chrono::duration_cast<std::chrono::duration<double>>(
std::chrono::high_resolution_clock::now() - gtb).count();
};
std::mt19937_64 rng{123};
std::vector<u64> limbs, divisors;
for (size_t i = 0; i < (1 << 17); ++i)
limbs.push_back(rng());
for (size_t i = 0; i < (1 << 8); ++i) {
size_t const nbits = rng() % 63 + 2;
u64 n = 0;
do {
n = (u64(1) << (nbits - 1)) + rng() % (u64(1) << (nbits - 1));
} while (!(n >= 3 && (n & (n - 1)) != 0));
divisors.push_back(n);
}
std::cout << std::fixed << std::setprecision(3);
double div_time = 0;
{
bi_uint x;
x.data = limbs;
auto const tim = Time();
for (auto dsr: divisors)
x.self_div(dsr);
div_time = Time() - tim;
std::cout << "Divide time " << div_time << " sec" << std::endl;
}
{
bi_uint x;
x.data = limbs;
for (size_t i = 0; i < 2; ++i) {
if (IS_MSVC && i == 0)
continue;
auto const tim = Time();
for (auto dsr: divisors)
x.self_div(dsr, [](auto n){ return BarrettRS64(n); },
[i](auto dnd_hi, auto dnd_lo, auto dsr, auto const & prep){
return i == 0 ? BarrettDivMod64<true>(dnd_hi, dnd_lo, dsr, prep.first, prep.second) :
BarrettDivMod64<false>(dnd_hi, dnd_lo, dsr, prep.first, prep.second);
});
double const bar_time = Time() - tim;
std::cout << "Barrett" << (i == 0 ? "128" : "64 ") << " time " << bar_time << " sec, boost " << div_time / bar_time << std::endl;
}
}
}
int main() {
try {
TestBarrett();
TestSpeed();
return 0;
} catch (std::exception const & ex) {
std::cout << "Exception: " << ex.what() << std::endl;
return -1;
}
}
输出：
Divide time 3.171 sec
Barrett128 time 0.675 sec, boost 4.695
Barrett64 time 0.642 sec, boost 4.937
第 2 部分
因为你有一个非常有趣的问题，在我第一次发布这个 post 几天后，我决定从头开始实现所有大整数数学。
下面的代码对自然大数（正整数）实现数学运算 +, -, *, /, <<, >>，对浮动大数实现 +, -, *, /。两种类型的数字都是任意大小的（甚至数百万位）。除了您要求的那些之外，我还完全实现了 Newton-Raphson (both square and cubic variants) and Goldschmidt 快速除法算法。
这里是 Newton-Raphson/Golschmidt 函数的代码片段，剩余的代码非常大，链接在下面的外部服务器上：
BigFloat & DivNewtonRaphsonSquare(BigFloat b) {
// https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division
auto a = *this;
a.exp_ += b.SetScale(0);
if (b.sign_) {
a.sign_ = !a.sign_;
b.sign_ = false;
}
thread_local BigFloat two, c_48_17, c_32_17;
thread_local size_t static_prec = 0;
if (static_prec != BigFloat::prec_) {
two = 2;
c_48_17 = BigFloat(48) / BigFloat(17);
c_32_17 = BigFloat(32) / BigFloat(17);
static_prec = BigFloat::prec_;
}
BigFloat x = c_48_17 - c_32_17 * b;
for (size_t i = 0, num_iters = std::ceil(std::log2(double(static_prec + 1)
/ std::log2(17.0))) + 0.1; i < num_iters; ++i)
x = x * (two - b * x);
*this = a * x;
return BitNorm();
}
BigFloat & DivNewtonRaphsonCubic(BigFloat b) {
// https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division
auto a = *this;
a.exp_ += b.SetScale(0);
if (b.sign_) {
a.sign_ = !a.sign_;
b.sign_ = false;
}
thread_local BigFloat one, c_140_33, c_m64_11, c_256_99;
thread_local size_t static_prec = 0;
if (static_prec != BigFloat::prec_) {
one = 1;
c_140_33 = BigFloat(140) / BigFloat(33);
c_m64_11 = BigFloat(-64) / BigFloat(11);
c_256_99 = BigFloat(256) / BigFloat(99);
static_prec = BigFloat::prec_;
}
BigFloat e, y, x = c_140_33 + b * (c_m64_11 + b * c_256_99);
for (size_t i = 0, num_iters = std::ceil(std::log2(double(static_prec + 1)
/ std::log2(99.0)) / std::log2(3.0)) + 0.1; i < num_iters; ++i) {
e = one - b * x;
y = x * e;
x = x + y + y * e;
}
*this = a * x;
return BitNorm();
}
BigFloat & DivGoldschmidt(BigFloat b) {
// https://en.wikipedia.org/wiki/Division_algorithm#Goldschmidt_division
auto a = *this;
a.exp_ += b.SetScale(0);
if (b.sign_) {
a.sign_ = !a.sign_;
b.sign_ = false;
}
BigFloat one = 1, two = 2, f;
for (size_t i = 0;; ++i) {
f = two - b;
a *= f;
b *= f;
if (i % 3 == 0 && (one - b).GetScale() < -i64(prec_) + i64(bit_sizeof(Word)))
break;
}
*this = a;
return BitNorm();
}
请参阅下面的 Output:，它将显示 Newton-Raphson 和 Goldschmidt 方法实际上比常规 School-grade 慢 10x 倍（称为 Reference 在输出中）算法。这 3 种高级算法彼此之间的速度大致相同。如果使用 Fast Fourier Transform 进行乘法运算，可能 Raphson/Goldschmidt 会更快，因为两个大数的乘法运算占用了这些算法 95% 的时间。在下面的代码中，Raphson/Goldschmidt 算法的所有结果不仅是 time-measured，而且还检查了与 School-grade（参考）算法相比结果的正确性（请参阅控制台输出中的 diff 2^...，这显示结果与学校成绩相比有多大差异）。
完整源代码在这里。完整的代码非常庞大，由于每个 post 的 30 000 个字符的 SO 限制，它不适合这个 Whosebug，尽管我从 scracth 专门为此 post 编写了这段代码。这就是为什么提供 external download link（PasteBin 服务器），同时单击下面的 Try it online! 链接，它是 运行 在 GodBolt 的 linux 服务器上运行的相同代码副本：
输出：
========== 1 K bits ==========
Reference 0.000029 sec
Raphson2 0.000066 sec, boost 0.440x, diff 2^-8192
Raphson3 0.000092 sec, boost 0.317x, diff 2^-8192
Goldschmidt 0.000080 sec, boost 0.365x, diff 2^-1022
========== 2 K bits ==========
Reference 0.000071 sec
Raphson2 0.000177 sec, boost 0.400x, diff 2^-16384
Raphson3 0.000283 sec, boost 0.250x, diff 2^-16384
Goldschmidt 0.000388 sec, boost 0.182x, diff 2^-2046
========== 4 K bits ==========
Reference 0.000319 sec
Raphson2 0.000875 sec, boost 0.365x, diff 2^-4094
Raphson3 0.001122 sec, boost 0.285x, diff 2^-32768
Goldschmidt 0.000881 sec, boost 0.362x, diff 2^-32768
========== 8 K bits ==========
Reference 0.000484 sec
Raphson2 0.002281 sec, boost 0.212x, diff 2^-65536
Raphson3 0.002341 sec, boost 0.207x, diff 2^-65536
Goldschmidt 0.002432 sec, boost 0.199x, diff 2^-8189
========== 16 K bits ==========
Reference 0.001199 sec
Raphson2 0.009042 sec, boost 0.133x, diff 2^-16382
Raphson3 0.009519 sec, boost 0.126x, diff 2^-131072
Goldschmidt 0.009047 sec, boost 0.133x, diff 2^-16380
========== 32 K bits ==========
Reference 0.004311 sec
Raphson2 0.039151 sec, boost 0.110x, diff 2^-32766
Raphson3 0.041058 sec, boost 0.105x, diff 2^-262144
Goldschmidt 0.045517 sec, boost 0.095x, diff 2^-32764
========== 64 K bits ==========
Reference 0.016273 sec
Raphson2 0.165656 sec, boost 0.098x, diff 2^-524288
Raphson3 0.210301 sec, boost 0.077x, diff 2^-65535
Goldschmidt 0.208081 sec, boost 0.078x, diff 2^-65534
========== 128 K bits ==========
Reference 0.059469 sec
Raphson2 0.725865 sec, boost 0.082x, diff 2^-1048576
Raphson3 0.735530 sec, boost 0.081x, diff 2^-1048576
Goldschmidt 0.703991 sec, boost 0.084x, diff 2^-131069
========== 256 K bits ==========
Reference 0.326368 sec
Raphson2 3.007454 sec, boost 0.109x, diff 2^-2097152
Raphson3 2.977631 sec, boost 0.110x, diff 2^-2097152
Goldschmidt 3.363632 sec, boost 0.097x, diff 2^-262141
========== 512 K bits ==========
Reference 1.138663 sec
Raphson2 12.827783 sec, boost 0.089x, diff 2^-524287
Raphson3 13.799401 sec, boost 0.083x, diff 2^-524287
Goldschmidt 15.836072 sec, boost 0.072x, diff 2^-524286
我目前正在开发一个 class 来处理大的无符号整数。但是，我需要不完整的功能，即：
bi_uint+=bi_uint- 已经实施。没有投诉。bi_uint*=std::uint_fast64_t- 已经实施。没有投诉。bi_uint/=std::uint_fast64_t- 已实现但工作起来非常 缓慢 ，还需要一个两倍于uint_fast64_t的类型。在测试用例中，除法比乘法慢 35 倍
接下来给出我的除法实现，基于简单的长除法算法：
#include <climits>
#include <cstdint>
#include <limits>
#include <vector>
class bi_uint
{
public:
using u64_t = std::uint_fast64_t;
constexpr static std::size_t u64_bits = CHAR_BIT * sizeof(u64_t);
using u128_t = unsigned __int128;
static_assert(sizeof(u128_t) >= sizeof(u64_t) * 2);
//little-endian
std::vector<u64_t> data;
//User should guarantee data.size()>0 and val>0
void self_div(const u64_t val)
{
auto it = data.rbegin();
if(data.size() == 1) {
*it /= val;
return;
}
u128_t rem = 0;
if(*it < val) {
rem = *it++;
data.pop_back();
}
u128_t r = rem % val;
while(it != data.rend()) {
rem = (r << u64_bits) + *it;
const u128_t q = rem / val;
r = rem % val;
*it++ = static_cast<u64_t>(q);
}
}
};
您可以看到使用了 unsigned __int128 类型，因此，此选项不可移植并且绑定到单个编译器 - GCC 并且还需要 x64 平台。
看了page关于除法算法的文章后，我觉得合适的算法应该是“Newton-Raphson division”。然而，“Newton-Raphson division”算法对我来说似乎很复杂。我想有一个更简单的算法来划分类型“big_uint/uint”，它的性能几乎相同。
问:如何快速将bi_uint分成u64_t?
我有大约 10^6 次迭代，每次迭代都使用列出的所有操作
如果这很容易实现，那么我希望具有可移植性而不使用unsigned __int128。否则，我宁愿放弃可移植性而选择更简单的方法。
编辑 1： 这是一个学术项目，我无法使用第三方库。
在大多数现代 CPU 中，除法确实比乘法慢得多。
参考
https://agner.org/optimize/instruction_tables.pdf
在 Intel Skylake 上 MUL/IMUL 有 3-4 个周期的延迟；而 DIV/IDIV 可能需要 26-90 个周期；比 MUL 慢 7 - 23 倍；所以您的初始基准测试结果并不令人意外。
如果您恰好在 x86 CPU 上，如下面的答案所示，如果这确实是瓶颈，您可以尝试使用 AVX/SSE 指令。基本上你需要依赖特殊说明而不是像 DIV/IDIV.
这样的一般说明第 1 部分（请参阅下面的第 2 部分）
我设法在我的旧笔记本电脑上将你的分区代码加速 5x 次（甚至在 GodBolt 服务器上加速 7.5x 次）使用 Barrett Reduction，这是一种允许用多次乘法和加法代替一次除法的技术。今天刚刚从 sctracth 实现了整个代码。
如果您愿意，可以直接跳转到我的 post 末尾的代码位置，而无需阅读冗长的描述，因为代码完全 运行nable 无需任何知识或依赖。
下面的代码仅适用于 Intel x64，因为我只使用了 Intel instructions 并且只使用了它们的 64 位变体。当然对于 x32 和其他处理器也可以是 re-written，因为 Barrett 算法是通用的。
为了简短地解释整个 Barrett Reduction pseudo-code 我会用 Python 来写它，因为它是最简单的语言，适合于理解 pseudo-code:
# https://www.nayuki.io/page/barrett-reduction-algorithm
def BarrettR(n, s):
return (1 << s) // n
def BarrettDivMod(x, n, r, s):
q = (x * r) >> s
t = x - q * n
return (q, t) if t < n else (q + 1, t - n)
基本上在上面的伪代码中 BarrettR() 只对同一个除数完成一次（你对整个大整数除法使用相同的 single-word 除数）。 BarrettDivMod() 每次要进行除法或模运算时都使用，基本上给定输入 x 和除数 n 它 returns 元组 (x / n, x % n)，没有别的，但它比常规除法指令更快。
在下面的 C++ 代码中，我实现了 Barrett 的两个相同函数，但做了一些 C++ 特定的优化以使其更快。优化是可能的，因为除数 n 总是 64 位，x 是 128 位，但较高的一半总是小于 n（最后一个假设发生是因为较高的一半在你的大整数除法始终是余数模 n).
Barrett 算法使用非 2 的幂的除数 n，因此不允许使用 0, 1, 2, 4, 8, 16, ... 这样的除数。你可以通过对大整数做正确的 bit-shift 来涵盖除数的这种简单情况，因为除以 2 的幂只是 bit-shift。允许任何其他除数，甚至包括不是 2 的幂的除数。
同样重要的是要注意我的 BarrettDivMod() 只接受严格小于 divisor * 2^64 的股息 x，换句话说，128 位股息的高一半 [=17] =] 应该小于除数。对于您的大整数除法函数，这始终是正确的，因为较高的一半始终是余数模数除数。 x 的这条规则应该由你检查，它在我的 BarrettDivMod() 中仅作为 DEBUG 断言检查，在发布时被删除。
你可以注意到 BarrettDivMod() 有两个大分支，这是同一算法的两个变体，第一个使用 CLang/GCC 仅类型 unsigned __int128，第二个仅使用 64 位指令和因此适用于 MSVC。
我试图针对三个编译器 CLang/GCC/MSVC，但有些 MSVC 版本在使用 Barrett 时速度只提高了 2 倍，而 CLang/GCC 速度提高了 5 倍。也许我在 MSVC 代码中做了一些错误。
您可以看到我使用您的 class bi_uint 对两个版本的代码进行时间测量 - 使用常规除法和使用 Barrett。 重要 请注意，我对您的代码进行了相当大的更改，首先不使用 u128（以便 MSVC 版本编译时没有 u128），其次不修改 data 向量，所以它只读除法并且不存储最终结果（我需要这个 read-only 来 运行 非常快地加速测试，而无需在每次测试迭代时复制 data 向量） .所以你的代码在我的代码片段中被破坏了，它可以't-be复制粘贴直接使用，我只用你的代码来测量速度。
Barrett 归约工作得更快，不仅因为除法比乘法慢，而且因为乘法和加法在现代 CPUs 上都很好地流水线化，现代 CPU 可以执行多个 mul 或 add一个周期内的指令，但前提是这几个 mul/add 不依赖于彼此的结果，换句话说 CPU 可以 运行 在一个周期内并行执行多个指令。据我所知除法不能 运行 并行，因为 CPU 中只有一个模块可以进行除法，但它仍然有点流水线，因为在第一次除法的 50% 之后是可以在 CPU 流水线的开头并行启动完成的第二部分。
在某些计算机上，您可能会注意到常规除法版本有时会慢得多，这是因为 CLang/GCC 回退到 library-based 除法的软实现，即使对于 128 位除数也是如此。在这种情况下，我的 Barrett 甚至可以显示 7-10x 倍的加速，因为它不使用库函数。
为了克服上述问题，关于软除法，最好直接使用 DIV 指令添加 Assembly code，或者在你的编译器中找到一些实现这个的内部函数（我认为 CLang/GCC 有这样的内在）。如果需要，我也可以编写这个 Assembly 实现，请在评论中告诉我。
更新。正如承诺的那样，为 CLang/GCC、函数 UDiv128Asm() 实现了 128 位除法的汇编变体。之后是更改它用作 CLang/GCC 128 分区的主要实现，而不是常规 u128(a) / b。您可以通过在 UDiv128() 函数体内将 #if 0 替换为 #if 1 来恢复常规的 u128 实现。
#include <cstdint>
#include <bit>
#include <stdexcept>
#include <string>
#include <immintrin.h>
#if defined(_MSC_VER) && !defined(__clang__)
#define IS_MSVC 1
#else
#define IS_MSVC 0
#endif
#if IS_MSVC
#include <intrin.h>
#endif
#define ASSERT_MSG(cond, msg) { if (!(cond)) throw std::runtime_error("Assertion (" #cond ") failed at line " + std::to_string(__LINE__) + "! Msg: '" + std::string(msg) + "'."); }
#define ASSERT(cond) ASSERT_MSG(cond, "")
#ifdef _DEBUG
#define DASSERT_MSG(cond, msg) ASSERT_MSG(cond, msg)
#else
#define DASSERT_MSG(cond, msg)
#endif
#define DASSERT(cond) DASSERT_MSG(cond, "")
using u16 = uint16_t;
using u32 = uint32_t;
using i64 = int64_t;
using u64 = uint64_t;
using UllPtr = unsigned long long *;
inline int GetExp(double x) {
return int((std::bit_cast<uint64_t>(x) >> 52) & 0x7FF) - 1023;
}
inline size_t BitSizeWrong(uint64_t x) {
return x == 0 ? 0 : (GetExp(x) + 1);
}
inline size_t BitSize(u64 x) {
size_t r = 0;
if (x >= (u64(1) << 32)) {
x >>= 32;
r += 32;
}
while (x >= 0x100) {
x >>= 8;
r += 8;
}
while (x) {
x >>= 1;
++r;
}
return r;
}
#if !IS_MSVC
inline u64 UDiv128Asm(u64 h, u64 l, u64 d, u64 * r) {
u64 q;
asm (R"(
.intel_syntax
mov rdx, %V[h]
mov rax, %V[l]
div %V[d]
mov %V[r], rdx
mov %V[q], rax
)"
: [q] "=r" (q), [r] "=r" (*r)
: [h] "r" (h), [l] "r" (l), [d] "r" (d)
: "rax", "rdx"
);
return q;
}
#endif
inline std::pair<u64, u64> UDiv128(u64 hi, u64 lo, u64 d) {
#if IS_MSVC
u64 r, q = _udiv128(hi, lo, d, &r);
return std::make_pair(q, r);
#else
#if 0
using u128 = unsigned __int128;
auto const dnd = (u128(hi) << 64) | lo;
return std::make_pair(u64(dnd / d), u64(dnd % d));
#else
u64 r, q = UDiv128Asm(hi, lo, d, &r);
return std::make_pair(q, r);
#endif
#endif
}
inline std::pair<u64, u64> UMul128(u64 a, u64 b) {
#if IS_MSVC
u64 hi, lo = _umul128(a, b, &hi);
return std::make_pair(hi, lo);
#else
using u128 = unsigned __int128;
auto const x = u128(a) * b;
return std::make_pair(u64(x >> 64), u64(x));
#endif
}
inline std::pair<u64, u64> USub128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
u64 r_hi, r_lo;
_subborrow_u64(_subborrow_u64(0, a_lo, b_lo, (UllPtr)&r_lo), a_hi, b_hi, (UllPtr)&r_hi);
return std::make_pair(r_hi, r_lo);
}
inline std::pair<u64, u64> UAdd128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
u64 r_hi, r_lo;
_addcarry_u64(_addcarry_u64(0, a_lo, b_lo, (UllPtr)&r_lo), a_hi, b_hi, (UllPtr)&r_hi);
return std::make_pair(r_hi, r_lo);
}
inline int UCmp128(u64 a_hi, u64 a_lo, u64 b_hi, u64 b_lo) {
if (a_hi != b_hi)
return a_hi < b_hi ? -1 : 1;
return a_lo == b_lo ? 0 : a_lo < b_lo ? -1 : 1;
}
std::pair<u64, size_t> BarrettRS64(u64 n) {
// https://www.nayuki.io/page/barrett-reduction-algorithm
ASSERT_MSG(n >= 3 && (n & (n - 1)) != 0, "n " + std::to_string(n))
size_t const nbits = BitSize(n);
// 2^s = q * n + r; 2^s = (2^64 + q0) * n + r; 2^s - n * 2^64 = q0 * n + r
u64 const dnd_hi = (nbits >= 64 ? 0ULL : (u64(1) << nbits)) - n;
auto const q0 = UDiv128(dnd_hi, 0, n).first;
return std::make_pair(q0, nbits);
}
template <bool Use128 = true, bool Adjust = true>
std::pair<u64, u64> BarrettDivMod64(u64 x_hi, u64 x_lo, u64 n, u64 r, size_t s) {
// ((x_hi * 2^64 + x_lo) * (2^64 + r)) >> (64 + s)
DASSERT(x_hi < n);
#if !IS_MSVC
if constexpr(Use128) {
using u128 = unsigned __int128;
u128 const xf = (u128(x_hi) << 64) | x_lo;
u64 q = u64((u128(x_hi) * r + xf + u64((u128(x_lo) * r) >> 64)) >> s);
if (s < 64) {
u64 t = x_lo - q * n;
if constexpr(Adjust) {
u64 const mask = ~u64(i64(t - n) >> 63);
q += mask & 1;
t -= mask & n;
}
return std::make_pair(q, t);
} else {
u128 t = xf - u128(q) * n;
return t < n ? std::make_pair(q, u64(t)) : std::make_pair(q + 1, u64(t) - n);
}
} else
#endif
{
auto const w1a = UMul128(x_lo, r).first;
auto const [w2b, w1b] = UMul128(x_hi, r);
auto const w2c = x_hi, w1c = x_lo;
u64 w1, w2 = _addcarry_u64(0, w1a, w1b, (UllPtr)&w1);
w2 += _addcarry_u64(0, w1, w1c, (UllPtr)&w1);
w2 += w2b + w2c;
if (s < 64) {
u64 q = (w2 << (64 - s)) | (w1 >> s);
u64 t = x_lo - q * n;
if constexpr(Adjust) {
u64 const mask = ~u64(i64(t - n) >> 63);
q += mask & 1;
t -= mask & n;
}
return std::make_pair(q, t);
} else {
u64 const q = w2;
auto const [b_hi, b_lo] = UMul128(q, n);
auto const [t_hi, t_lo] = USub128(x_hi, x_lo, b_hi, b_lo);
return t_hi != 0 || t_lo >= n ? std::make_pair(q + 1, t_lo - n) : std::make_pair(q, t_lo);
}
}
}
#include <random>
#include <iomanip>
#include <iostream>
#include <chrono>
void TestBarrett() {
std::mt19937_64 rng{123}; //{std::random_device{}()};
for (size_t i = 0; i < (1 << 11); ++i) {
size_t const nbits = rng() % 63 + 2;
u64 n = 0;
do {
n = (u64(1) << (nbits - 1)) + rng() % (u64(1) << (nbits - 1));
} while (!(n >= 3 && (n & (n - 1)) != 0));
auto const [br, bs] = BarrettRS64(n);
for (size_t j = 0; j < (1 << 6); ++j) {
u64 const hi = rng() % n, lo = rng();
auto const [ref_q, ref_r] = UDiv128(hi, lo, n);
u64 bar_q = 0, bar_r = 0;
for (size_t k = 0; k < 2; ++k) {
bar_q = 0; bar_r = 0;
if (k == 0)
std::tie(bar_q, bar_r) = BarrettDivMod64<true>(hi, lo, n, br, bs);
else
std::tie(bar_q, bar_r) = BarrettDivMod64<false>(hi, lo, n, br, bs);
ASSERT_MSG(bar_q == ref_q && bar_r == ref_r, "i " + std::to_string(i) + ", j " + std::to_string(j) + ", k " + std::to_string(k) +
", nbits " + std::to_string(nbits) + ", n " + std::to_string(n) + ", bar_q " + std::to_string(bar_q) +
", ref_q " + std::to_string(ref_q) + ", bar_r " + std::to_string(bar_r) + ", ref_r " + std::to_string(ref_r));
}
}
}
}
class bi_uint
{
public:
using u64_t = std::uint64_t;
constexpr static std::size_t u64_bits = 8 * sizeof(u64_t);
//little-endian
std::vector<u64_t> data;
static auto constexpr DefPrep = [](auto n){
return std::make_pair(false, false);
};
static auto constexpr DefDivMod = [](auto dnd_hi, auto dnd_lo, auto dsr, auto const & prep){
return UDiv128(dnd_hi, dnd_lo, dsr);
};
//User should guarantee data.size()>0 and val>0
template <typename PrepT = decltype(DefPrep), typename DivModT = decltype(DefDivMod)>
void self_div(const u64_t val, PrepT const & Prep = DefPrep, DivModT const & DivMod = DefDivMod)
{
auto it = data.rbegin();
if(data.size() == 1) {
*it /= val;
return;
}
u64_t rem_hi = 0, rem_lo = 0;
if(*it < val) {
rem_lo = *it++;
//data.pop_back();
}
auto const prep = Prep(val);
u64_t r = rem_lo % val;
u64_t q = 0;
while(it != data.rend()) {
rem_hi = r;
rem_lo = *it;
std::tie(q, r) = DivMod(rem_hi, rem_lo, val, prep);
//*it++ = static_cast<u64_t>(q);
it++;
auto volatile out = static_cast<u64_t>(q);
}
}
};
void TestSpeed() {
auto Time = []{
static auto const gtb = std::chrono::high_resolution_clock::now();
return std::chrono::duration_cast<std::chrono::duration<double>>(
std::chrono::high_resolution_clock::now() - gtb).count();
};
std::mt19937_64 rng{123};
std::vector<u64> limbs, divisors;
for (size_t i = 0; i < (1 << 17); ++i)
limbs.push_back(rng());
for (size_t i = 0; i < (1 << 8); ++i) {
size_t const nbits = rng() % 63 + 2;
u64 n = 0;
do {
n = (u64(1) << (nbits - 1)) + rng() % (u64(1) << (nbits - 1));
} while (!(n >= 3 && (n & (n - 1)) != 0));
divisors.push_back(n);
}
std::cout << std::fixed << std::setprecision(3);
double div_time = 0;
{
bi_uint x;
x.data = limbs;
auto const tim = Time();
for (auto dsr: divisors)
x.self_div(dsr);
div_time = Time() - tim;
std::cout << "Divide time " << div_time << " sec" << std::endl;
}
{
bi_uint x;
x.data = limbs;
for (size_t i = 0; i < 2; ++i) {
if (IS_MSVC && i == 0)
continue;
auto const tim = Time();
for (auto dsr: divisors)
x.self_div(dsr, [](auto n){ return BarrettRS64(n); },
[i](auto dnd_hi, auto dnd_lo, auto dsr, auto const & prep){
return i == 0 ? BarrettDivMod64<true>(dnd_hi, dnd_lo, dsr, prep.first, prep.second) :
BarrettDivMod64<false>(dnd_hi, dnd_lo, dsr, prep.first, prep.second);
});
double const bar_time = Time() - tim;
std::cout << "Barrett" << (i == 0 ? "128" : "64 ") << " time " << bar_time << " sec, boost " << div_time / bar_time << std::endl;
}
}
}
int main() {
try {
TestBarrett();
TestSpeed();
return 0;
} catch (std::exception const & ex) {
std::cout << "Exception: " << ex.what() << std::endl;
return -1;
}
}
输出：
Divide time 3.171 sec
Barrett128 time 0.675 sec, boost 4.695
Barrett64 time 0.642 sec, boost 4.937
第 2 部分
因为你有一个非常有趣的问题，在我第一次发布这个 post 几天后，我决定从头开始实现所有大整数数学。
下面的代码对自然大数（正整数）实现数学运算 +, -, *, /, <<, >>，对浮动大数实现 +, -, *, /。两种类型的数字都是任意大小的（甚至数百万位）。除了您要求的那些之外，我还完全实现了 Newton-Raphson (both square and cubic variants) and Goldschmidt 快速除法算法。
这里是 Newton-Raphson/Golschmidt 函数的代码片段，剩余的代码非常大，链接在下面的外部服务器上：
BigFloat & DivNewtonRaphsonSquare(BigFloat b) {
// https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division
auto a = *this;
a.exp_ += b.SetScale(0);
if (b.sign_) {
a.sign_ = !a.sign_;
b.sign_ = false;
}
thread_local BigFloat two, c_48_17, c_32_17;
thread_local size_t static_prec = 0;
if (static_prec != BigFloat::prec_) {
two = 2;
c_48_17 = BigFloat(48) / BigFloat(17);
c_32_17 = BigFloat(32) / BigFloat(17);
static_prec = BigFloat::prec_;
}
BigFloat x = c_48_17 - c_32_17 * b;
for (size_t i = 0, num_iters = std::ceil(std::log2(double(static_prec + 1)
/ std::log2(17.0))) + 0.1; i < num_iters; ++i)
x = x * (two - b * x);
*this = a * x;
return BitNorm();
}
BigFloat & DivNewtonRaphsonCubic(BigFloat b) {
// https://en.wikipedia.org/wiki/Division_algorithm#Newton%E2%80%93Raphson_division
auto a = *this;
a.exp_ += b.SetScale(0);
if (b.sign_) {
a.sign_ = !a.sign_;
b.sign_ = false;
}
thread_local BigFloat one, c_140_33, c_m64_11, c_256_99;
thread_local size_t static_prec = 0;
if (static_prec != BigFloat::prec_) {
one = 1;
c_140_33 = BigFloat(140) / BigFloat(33);
c_m64_11 = BigFloat(-64) / BigFloat(11);
c_256_99 = BigFloat(256) / BigFloat(99);
static_prec = BigFloat::prec_;
}
BigFloat e, y, x = c_140_33 + b * (c_m64_11 + b * c_256_99);
for (size_t i = 0, num_iters = std::ceil(std::log2(double(static_prec + 1)
/ std::log2(99.0)) / std::log2(3.0)) + 0.1; i < num_iters; ++i) {
e = one - b * x;
y = x * e;
x = x + y + y * e;
}
*this = a * x;
return BitNorm();
}
BigFloat & DivGoldschmidt(BigFloat b) {
// https://en.wikipedia.org/wiki/Division_algorithm#Goldschmidt_division
auto a = *this;
a.exp_ += b.SetScale(0);
if (b.sign_) {
a.sign_ = !a.sign_;
b.sign_ = false;
}
BigFloat one = 1, two = 2, f;
for (size_t i = 0;; ++i) {
f = two - b;
a *= f;
b *= f;
if (i % 3 == 0 && (one - b).GetScale() < -i64(prec_) + i64(bit_sizeof(Word)))
break;
}
*this = a;
return BitNorm();
}
请参阅下面的 Output:，它将显示 Newton-Raphson 和 Goldschmidt 方法实际上比常规 School-grade 慢 10x 倍（称为 Reference 在输出中）算法。这 3 种高级算法彼此之间的速度大致相同。如果使用 Fast Fourier Transform 进行乘法运算，可能 Raphson/Goldschmidt 会更快，因为两个大数的乘法运算占用了这些算法 95% 的时间。在下面的代码中，Raphson/Goldschmidt 算法的所有结果不仅是 time-measured，而且还检查了与 School-grade（参考）算法相比结果的正确性（请参阅控制台输出中的 diff 2^...，这显示结果与学校成绩相比有多大差异）。
完整源代码在这里。完整的代码非常庞大，由于每个 post 的 30 000 个字符的 SO 限制，它不适合这个 Whosebug，尽管我从 scracth 专门为此 post 编写了这段代码。这就是为什么提供 external download link（PasteBin 服务器），同时单击下面的 Try it online! 链接，它是 运行 在 GodBolt 的 linux 服务器上运行的相同代码副本：
输出：
========== 1 K bits ==========
Reference 0.000029 sec
Raphson2 0.000066 sec, boost 0.440x, diff 2^-8192
Raphson3 0.000092 sec, boost 0.317x, diff 2^-8192
Goldschmidt 0.000080 sec, boost 0.365x, diff 2^-1022
========== 2 K bits ==========
Reference 0.000071 sec
Raphson2 0.000177 sec, boost 0.400x, diff 2^-16384
Raphson3 0.000283 sec, boost 0.250x, diff 2^-16384
Goldschmidt 0.000388 sec, boost 0.182x, diff 2^-2046
========== 4 K bits ==========
Reference 0.000319 sec
Raphson2 0.000875 sec, boost 0.365x, diff 2^-4094
Raphson3 0.001122 sec, boost 0.285x, diff 2^-32768
Goldschmidt 0.000881 sec, boost 0.362x, diff 2^-32768
========== 8 K bits ==========
Reference 0.000484 sec
Raphson2 0.002281 sec, boost 0.212x, diff 2^-65536
Raphson3 0.002341 sec, boost 0.207x, diff 2^-65536
Goldschmidt 0.002432 sec, boost 0.199x, diff 2^-8189
========== 16 K bits ==========
Reference 0.001199 sec
Raphson2 0.009042 sec, boost 0.133x, diff 2^-16382
Raphson3 0.009519 sec, boost 0.126x, diff 2^-131072
Goldschmidt 0.009047 sec, boost 0.133x, diff 2^-16380
========== 32 K bits ==========
Reference 0.004311 sec
Raphson2 0.039151 sec, boost 0.110x, diff 2^-32766
Raphson3 0.041058 sec, boost 0.105x, diff 2^-262144
Goldschmidt 0.045517 sec, boost 0.095x, diff 2^-32764
========== 64 K bits ==========
Reference 0.016273 sec
Raphson2 0.165656 sec, boost 0.098x, diff 2^-524288
Raphson3 0.210301 sec, boost 0.077x, diff 2^-65535
Goldschmidt 0.208081 sec, boost 0.078x, diff 2^-65534
========== 128 K bits ==========
Reference 0.059469 sec
Raphson2 0.725865 sec, boost 0.082x, diff 2^-1048576
Raphson3 0.735530 sec, boost 0.081x, diff 2^-1048576
Goldschmidt 0.703991 sec, boost 0.084x, diff 2^-131069
========== 256 K bits ==========
Reference 0.326368 sec
Raphson2 3.007454 sec, boost 0.109x, diff 2^-2097152
Raphson3 2.977631 sec, boost 0.110x, diff 2^-2097152
Goldschmidt 3.363632 sec, boost 0.097x, diff 2^-262141
========== 512 K bits ==========
Reference 1.138663 sec
Raphson2 12.827783 sec, boost 0.089x, diff 2^-524287
Raphson3 13.799401 sec, boost 0.083x, diff 2^-524287
Goldschmidt 15.836072 sec, boost 0.072x, diff 2^-524286