infinite-pi，用C++算任意位数的π

为什么做

这是TUM C++课程的期末项目，题目是用Chudnovsky公式和binary splitting算任意位数的π，成绩全部来自这个项目。

课程要求的，和我多做的

下面按课程展示时的任务清单整理。

我做了什么

• 在 include/chudnovsky.h（329行）里写完整个算法：单项、合并、递归、并行、收尾公式、项数估计、√π。全部是模板，数字类型 T 是参数。
• 写了 main.cpp（命令行演示）和 bench/（计时程序加两个Python驱动脚本、两张图）。
• 写了19个单元测试和一个很小的测试框架（137行，没有用Catch2之类的库）。
• 写了CMake、.clang-tidy、.clang-format、Doxygen配置，和 README.md、FEATURES.md。

没做的：自己写BigInt/BigFloat、自己写除法和开方、FFT。

怎么做的

binary splitting

Chudnovsky级数的第k项是个分数。逐项相加要做很多次高精度除法。binary splitting的做法是不除，把一段区间 [lo, hi) 的和表示成三个整数 (P, Q, T)，这段的和等于 T / Q。

• 叶子（只有1项）：leafTerm(k) 用递推式直接算。k=0时是 (1, 1, 13591409)；k≥1时 p = (6k-5)(2k-1)(6k-1)，q = k³ × 10939058860032000，t = p × (13591409 + 545140134k)，k为奇数时 t 取负。
• 拆分：mid = (lo + hi) / 2，左右各递归一次，终止条件是 hi - lo == 1。
• 合并：两段相邻区间用同一个函数 mergePQT，串行和并行共用。

T P = left.P * right.P;
T Q = left.Q * right.Q;
T t = right.Q * left.T_ + left.P * right.T_;   // 不对称，左右不能换

• 收尾：整棵树里没有除法。只有最后一步 π = 426880 · √10005 · Q / T 做一次除法、一次开方。
• 项数：每项约14.18位，termsNeededFor(d) = ceil(d / 14.1816) + 1。30,000位要2,117项，100,000位要7,053项。多算1项是故意的，宁多勿少。
• √π：对算好的π再调一次 sqrt。

数字类型：Boost的cpp_dec_float，不是自己写的

算法只要求 T 有乘、加、取负，最后一步还要除和开方。我用两个C++20 concept写这些要求：BigNumber 管树，BigReal 管最后一步。BigReal 多加了一条 !numeric_limits<T>::is_integer，因为 cpp_int 也有 / 和 sqrt（整数版），只查语法的话会放它进来，然后悄悄算出错误结果。测试里用 static_assert 钉住这一点。

实际跑的类型是 cpp_dec_float<位数+10>：十进制浮点，每个元素存8位十进制（基数10⁸，32位整数）。我没有自己写limb、进位、Karatsuba、Newton迭代，除法和开方都是Boost的。

两个后果：

• cpp_dec_float的精度是编译期常数，运行时才知道要几位。所以 main.cpp 预先实例化1,000/3,000/10,000/30,000/100,000位五档，用fold expression从小到大试，选第一个够用的。
• 大于等于20,000位的档，数字存储从栈上的 std::array 换成堆分配（std::conditional_t）。不换的话，递归里每层要放好几个几十KB的数，栈会爆。

Boost在项目里的用途有两个：运行时的数字类型；测试里用 cpp_int 做精确整数对照。它不是拿来给我的实现做对拍的，对拍用的参考位数是用mpmath单独生成的。

并行：std::async，按子树分，只分最上面几层

binarySplittingParallel(lo, hi, depth)：

• 在最上面 depth 层，右半段交给 std::async(std::launch::async, ...)，当前线程自己算左半段，再 get() 取回右半段并合并。
• depth = splitDepthFor(threads)，即让 2^depth 不小于线程数的最小整数。8线程是depth 3，切成8棵子树，开7个async任务。30,000位时每棵子树约265项，100,000位时约882项。
• 低于这个深度，或区间不足2项，退回串行的 binarySplitting。
• 分的粒度由深度决定，没有按项数设阈值。线程数不是2的幂时向上取整，所以传5到8个线程，实际都是8棵子树。

并行版的树和串行版完全相同，用同一个 mid、同一个 mergePQT，只是兄弟子树谁先算不同。所以结果逐位相等。测试里测了两次：精确整数 (P, Q, T) 和最终的π值。

构建和工具

• C++23，CMake里可以改成20。需要CMake 3.25以上，find_package(Boost) 和 Threads。
• CMake文件里没有设自己的编译选项，也没有默认的构建类型。README第一段的构建命令没带 -DCMAKE_BUILD_TYPE=Release，只有跑基准的命令带了。不带的话演示程序是没开优化的。
• .clang-tidy 开了bugprone、performance、modernize、readability，关掉了一批（比如魔法数字检查）。.clang-format 是Google风格、4格缩进、每行100列。这两个只是配置文件，没有接进CMake。头文件里有一处 NOLINT，说明我跑过。
• Doxygen有配置，生成的 docs/ 在 .gitignore 里，没提交。

测试和基准怎么做的

• 19个用例分7组：leafTerm 4个、binarySplitting 4个、mergePQT 1个、并行5个、computePi 3个、termsNeededFor 1个、sqrtOfPi 1个。整数部分用 cpp_int 和手算的精确值逐位对。π部分和独立生成的参考位数比，容差最小到1e-100。CTest里只登记了一个测试 pi_tests，19个用例在它里面跑。
• 基准程序 bench_mine 用 steady_clock 只计时 computePi，每次一行输出。每个位数用刚好够的精度档，不多给。结果的前10位打到stderr，一是自查，二是防止编译器把计算删掉。
• 两个Python脚本调用它：位数小于20,000跑3次取最小，大于等于20,000只跑1次。GMP一侧是整个进程的墙钟时间。

结果

并行加速比

下表来自 bench/speedup.csv，加速比 = 1线程耗时 ÷ N线程耗时。

30,000位、8线程的实测耗时是6.71 s。100,000位算π加√π，8线程用了201 s，这个数字只在README里，原始数据不在 bench/，我也没有重跑。

对比gmp-chudnovsky

下表来自 bench/results.csv，比值是我用两列算的。

100位和300位不放：GMP一侧含约2 ms进程启动，比本项目的计算时间还长，数字没有意义。

本项目的曲线在log-log图上斜率约2.6到2.7（3,000到10,000位是2.57，10,000到30,000位是2.73，我用CSV算的），GMP的曲线几乎是平的。差距随位数越拉越大，不是常数倍。

README里还有一个对照实验：同样算30,000位，工作宽度开成100,010位要316 s，开成30,010位要43.7 s，慢7.2倍。这个实验没有脚本和CSV，只有README里的一句话。

没做好的地方

• 和GMP差得很远，原因我还没证实，见下一节。
• 没有自己写大数类。「从零实现」只限算法层，乘法、除法、开方都是Boost的。
• 6.40倍比的是自己的串行版本，而这个串行版本本身比GMP慢几千倍。
• 30,000位的每个数字都只跑了1次，没有取最小值，也没有重复。两次独立运行差1.7%，可以当作噪声的粗略估计。
• GMP在30,000位的7.67 ms是单次整个进程的墙钟时间，含约2 ms启动，噪声很大。比值5,697是低估，不是高估。
• 大位数的正确性：测试最大只到约130位。基准程序只打印结果的前10位。30,000位和100,000位的结果有没有和参考值逐位比过，仓库里看不出来。
• 「树里全是精确整数」这个说法，只对测试里的 cpp_int 成立。我按公式估算，30,000位时根节点的Q约有52,000位，比30,010位的工作宽度还宽，所以树的上半部分在做十进制浮点的舍入。估算没有实测，结果的有效位数仍然够，但README和展示稿里「No rounding」的说法需要改。
• 8线程不等于8个同样快的核。这台Mac有性能核和能效核，所以80%的并行效率只能当参考。
• README第一段构建命令没带Release。
• clang-tidy「0 warnings」是我展示稿里的说法，这次我没有重新跑。

为什么慢5,700倍：待profile验证的假设

这一节是读代码和Boost头文件得到的怀疑，都没有实测，不是结论。每条后面写了怎么证明它是错的。

H1：每次运算都付全宽的价。 我读了Boost 1.90的 cpp_dec_float.hpp：乘法按 min(prec_elem, v.prec_elem) 个元素做，而 prec_elem 默认等于编译期宽度。也就是说叶子上的小数和根节点上的大数，一次乘法的成本差不多。数一下代码里的 T × T：每个叶子6次，每次合并4次，30,000位约21,166次。43.7 s ÷ 21,166 ≈ 2.1 ms一次，量级上说得通。README里那个316 s对43.7 s的实验（宽度放大3.33倍，慢7.2倍，指数约1.64）也指向同一个方向。 怎么证伪：对30,000位的串行运行做Time Profiler，看时间是不是集中在乘法上；单独测一次宽度30,010的乘法是不是约2 ms；把叶子附近的子树换成 cpp_int 跑一遍，看总时间降多少。

H2：乘法算法，README的说法可能不对。 README和展示稿都写 cpp_dec_float 用schoolbook乘法，O(n²)。但Boost 1.90的头文件里，元素数达到129个（约1,032位）以上时会走Karatsuba。30,010位宽是3,755个元素，会走Karatsuba。如果这是真的，「乘法是O(n²)」这个解释要改。数字上也对得上：项数（一次方）乘Karatsuba（1.585次方）是2.585，实测斜率2.57到2.73；纯schoolbook应该接近3。吻合只是算术上吻合，不是证明。Karatsuba也比不过GMP的Toom-Cook和FFT，所以差距还是在，只是原因的分配要变。 怎么证伪：看profile的调用栈里有没有 eval_multiply_kara_n_by_n_to_2n；确认测量时用的Boost版本。

H3：十进制表示。 每个元素存8位十进制，约26.6个二进制位，GMP的limb是64位。元素数多约2.4倍；如果都是schoolbook，乘法次数多约5.8倍。每次乘完还要在十进制下进位整理。这是常数因子，解释不了斜率。 怎么证伪：只看H1和H2修正后剩下的差距有多大。

H4：内存分配。 大于等于20,000位的档用堆分配。我读头文件，走Karatsuba时每次乘法要建一个临时池，8 × 4,096个32位元素，约128 KB（用的是 dynamic_array，4,096这个数我没有细查）。如果每次乘法都分配一次，约21,000次乘法就是约21,000次128 KB的分配。README里「每个中间结果一次堆分配」说的是另一件事。 怎么证伪：用Instruments的Allocations或dtrace数 malloc 次数和大小；看malloc在profile里占多少。

H5：没有约分公因子。 README说gmp-chudnovsky会把P/Q/T里的公因子筛掉，让操作数更短，我的实现没做。我没有读过gmp-chudnovsky.c，这条只是转述README。影响多大不知道。 怎么证伪：读gmp-chudnovsky.c，在我的实现里加上，看操作数位数和耗时变不变。

H6：最后的除法和开方。 全宽的Boost除法和开方各做一次，算√π时再多一次开方。我猜占比小，但没有测过。 怎么证伪：在 finalizePi 前后各计一次时。

H7：并行为什么只有6.40倍，README的解释可能也要改。 README归因于Amdahl：树顶的合并最大，再加最后的除法和开方。但如果H1成立，每次乘法成本差不多，树顶那三层的合并（关键路径上约12次乘法）对总共21,166次来说几乎可以忽略，理想加速比应该接近8。更像样的候选是：8个线程里有一部分落在能效核上（4线程时是×3.90，到8线程才弯下来）；大于等于20,000位的档下多线程同时 malloc 的争用；内存带宽。 怎么证伪：分别记录每个线程的忙时间；换成小于20,000位的档（栈上存储）看弯曲是否还在；注意线程数传5到8都是8棵子树，扫线程数看不出能效核的影响，要换成固定大小的线程池。

怎么复现

仓库：papersii/infinite-pi，commit b6810e7。

cmake -S . -B build -DCMAKE_BUILD_TYPE=Release -DBUILD_BENCH=ON
cmake --build build
ctest --test-dir build --output-on-failure   # 19个用例
python3 bench/run_speedup.py                  # 加速比，30,000位串行约43 s

和GMP对比要多加 -DBUILD_REFERENCE=ON，需要 brew install gmp，CMake会自动下载gmp-chudnovsky.c，再跑 python3 bench/run_bench.py。

演示程序：./build/pi 1000、./build/pi 10000 4（4线程）、./build/pi -t 72（正好72项）。上限100,000位，超过30,000位会提示要等几分钟。

接下来

先搞清楚为什么慢5,700倍，再决定改哪里。

1. 在30,000位上做串行的profile，先验证H1、H2、H4，看时间花在哪。
2. 确认测量时用的Boost版本，再回头改README里「schoolbook」的说法。
3. 把30,000位的基准改成跑多次取最小，补上100,000位的原始数据。
4. 用 cpp_int 做叶子附近的子树，看H1对不对。