企业官网建设流程全解析

算圆周率这种事，听起来特别像大学期末编程作业。但真把需求摆到台面上,目标是“小数点后10000位”，性质就完全变了——它不再是你用Math.PI糊弄一下、截图交差的小练习，而是实打实的高精度计算问题。你在网上搜一圈就会发现，做这件事的人大致分三类：一类是刚学完Python循环，想看看自己能算到多少位；一类是被面试题问到了，比如“不用库，怎么算pi到任意精度”；还有一类是纯粹想验证某台机器或者某个语言的数值计算能力，拿10000位当benchmark用。

不管你是哪种，这篇文章都打算把这件事彻底讲透：从算法选型到浮点数陷阱，从代码实现到末位校验，全部走一遍。我会用Python为主，顺便提一下如果是用C/C++或者JavaScript，思路要怎么变。

1. 方案选型：先想清楚，再写循环

1.1 为什么不能直接double乘到底

先说一个最常见的坑。很多人的第一反应是：圆周率不是有现成公式吗，pi = 4 * arctan(1)，然后对着arctan的泰勒展开一项一项算不就行了？理论上是这样，但只要你真的用Python的float去跑，算到小数点后十几位就开始不对劲，到四五十位的时候已经完全不可信了。

原因很简单：IEEE 754的双精度浮点数只有53位有效二进制位，换算成十进制大约是15到17位有效数字。你连第50位都算不出来，更别提10000位。所以这个问题的第一课是：如果你需要的小数位数超过15位，一切浮点类型都要扔到一边，必须用“高精度数”来做运算。Python里的decimal.Decimal、fractions.Fraction，或者干脆用整数去模拟小数点，都是路子。

1.2 自己算，还是用现成的高精度库

很多人会问，那直接用mpmath这种库不就完了，一行mpmath.mp.dps = 20000; print(mpmath.pi)搞定。这确实是最快的办法。但如果你正在准备面试、正在做算法作业，或者只是想把计算过程彻底搞明白，那用库反而会丢掉核心的乐趣。我自己倾向于把这当成两种梯度：

• 如果用库，那么重点在于理解库背后用了什么公式、为什么快。
• 如果自己写，那么重点在于精度控制、存储设计和迭代收敛。

我在这篇文章里会把两种路线都演示一遍：先用Python标准库decimal配合一个收敛很快的公式拿到结果，然后手写一个用整数运算实现的高精度计算版本，最后对比校验。

1.3 公式选哪个，决定了你要写多少循环

计算pi的公式非常多，但并不是所有公式都适合算到10000位。

• Leibniz公式：收敛慢到令人绝望，算到几千项才几个有效位，基本只能用来讲概念。
• Machin公式：pi/4 = 4*arctan(1/5) - arctan(1/239)，收敛快，手工推导友好，适合教学。每算一项大约能多出1.4位十进制精度。
• Chudnovsky算法：收敛极快，每项能多出约14位精度，是目前各种破纪录计算的标准选项。

打个比方，如果Leibniz公式是走路，Machin公式是骑自行车，那Chudnovsky就是高铁。所以这篇文章里，用decimal实现时我会直接用Chudnovsky算法，自己手写整数版本时则用Machin公式——原因后面会说，手写整数高精度除法和开方是有代价的，Chudnovsky涉及大整数平方根，初学者容易栽进去。

2. Chudnovsky算法：用Python标准库算出10000位

2.1 核心原理

Chudnovsky算法的公式长这样：

1/pi = 12 * sum_{k=0..∞} (-1)^k * (6k)! * (13591409 + 545140134k) / ((3k)! * (k!)^3 * 640320^(3k+1.5))

别被这串东西吓到。它本质上是一个收敛极快的级数。你每多算一项，正确数字的位数大概能多出14位。这意味着，算到10000位，理论上只需要迭代到10000 / 14次左右，也就是大约714项。实际计算时考虑到末尾截断误差，你多迭代几十项就够了。

这个公式能快速收敛的原因，和模形式的深层数学性质有关，不是这篇博客的重点。我们只需要把它当作一个“黑箱但透明”的迭代公式来用：知道它快，知道怎么实现，知道每一项的计算方式即可。

2.2 关键参数与预处理

当我们要用decimal.Decimal来算时，需要提前做几件事：

1. 设置全局精度。如果要输出10000位，不能只把精度设为10000，因为计算过程中中间值有舍入误差，一般要留出额外的保护位。我通常会设置到10000 + 50，多的50位留给中间过程的舍入误差。

2. 把公式里的常数预先算出来，不要放在循环里反复用高精度除法。比如640320**1.5、545140134、13591409这些，循环外先算好。

3. 尽量用整数运算保存迭代过程中的递推关系，而不是每次都从头计算阶乘。比如我们可以维护当前k对应的(6k)!、(3k)!、(k!)^3，通过递推公式从k推进到k+1，这样能省掉大量乘法。

2.3 代码实现

下面这段代码是我实际测试过可以直接跑的版本：

from decimal import Decimal, getcontext # 计算位数 DIGITS = 10000 # 多算一些保护位，减少中间舍入误差 getcontext().prec = DIGITS + 50 def chudnovsky_pi(n_terms): """ 利用 Chudnovsky 级数计算圆周率 n_terms: 迭代项数，每项约贡献 14 位精度 """ C = 426880 * Decimal(10005).sqrt() M = Decimal(1) L = Decimal(13591409) X = Decimal(1) K = Decimal(6) S = Decimal(13591409) for k in range(1, n_terms): M = M * (K**3 - 16*K) / (Decimal(k)**3) L += Decimal(545140134) X *= Decimal(-262537412640768000) S += Decimal(M * L) / X K += Decimal(12) pi = C / S return pi if __name__ == "__main__": pi = chudnovsky_pi(800) print(str(pi)[:DIGITS + 2]) # 包含整数部分的 3 和小数点

这里n_terms我直接传了800。刚才说了大约714项够用，多个几十项是为了让末尾更稳定。输出的时候，因为getcontext().prec设置成了DIGITS + 50，得到的Decimal会多出一些尾部数字，直接截断到我们需要的位数即可。

2.4 为什么精度要设置成DIGITS + 50

很多第一次接触高精度计算的人会问，10000位就10000位，为什么要多算50位？

因为.这个公式虽然理论上每项贡献14位，但每一步都涉及乘除法，Decimal在运算中会对中间结果做舍入。如果精度卡在10000位整，最后一次迭代的误差可能会污染最后几位。多出50位保护位，等于给整个计算过程留了缓冲带。最终输出前再从第一位开始截取，就不用担心末位抖动。

有一个经验规律是：保护位至少是迭代次数的两位数级别。我一般习惯取DIGITS / 100 + 20，对于10000位来说就是120，比较保守。但这里为了简洁，用50也够。要注意的是，保护位越多计算耗时越长，因为高精度乘法耗时大约是位数的一次方到二次方之间，没必要无限放大。

2.5 输出后怎么验证

输出一串数字很容易，难的是确认这串数字是正确的。

第一步，先确认前几位是3.14159265。这能排除大部分错误。 第二步，截取一段已知的正确圆周率片段做比对。比如我知道圆周率小数点后第1到20位是897932384626433832，可以直接用字符串查找验证。 第三步，更具说服力的办法是：用两种完全不同的算法分别计算10000位，再逐位比较。比如这一节的Chudnovsky算法算一遍，下一节手写的基于整数和Machin公式的版本再算一遍，两者比对。如果完全一致，基本可以认定结果正确。

3. 不用浮点数：手写高精度整数运算算pi

3.1 核心思路

如果说上一部分你只是写了一个调用Decimal的脚本，那这一部分才算真正理解高精度计算的底层原理。你要丢掉Decimal，丢掉float，只用Python原生的整数类型，自己定义数组来模拟“任意精度实数”。

一个很自然的想法是：把pi表示成“整数部分 + 一个小数点后N位的整数”。比如我要算100位，那我实际计算的是int(pi * 10^100)，一个巨大的整数，最后需要输出时，我再在合适的位置插入小数点。这样一来，所有高精度计算都转化成整数运算，完全绕开浮点数精度限制。

3.2 用整数实现高精度算术

什么叫“用整数实现高精度算术”？比如你要算两个“拥有10000位小数”的数之和，最简单的方式是：

• 每个数用Python列表存每一位数字（0-9）。
• 从最低位开始逐位相加，记下进位。

同理，高精度乘法和除法都可以用类似手算的方法实现。Python的整数类型本身可以极大，所以还有个取巧的做法：不逐位存0-9，而是把每9位数字拼成一个int作为“一坨”，然后像数组元素一样处理，从而大幅减少进位处理逻辑。

对于10000位圆周率，位数并不算夸张，所以我直接用“整个数放大成一个大整数”的策略。比如我计算int(pi * 10^10000)，这个数大约有10001位十进制数字，Python的int处理起来毫无压力。

3.3 基于Machin公式的整数实现

Machin公式是：

pi/4 = 4*arctan(1/5) - arctan(1/239)

而arctan(x)的泰勒展开是：

arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + ...

当x = 1/5或x = 1/239时，每一项都是有理数。我们把整个计算乘以一个足够大的缩放因子SCALE = 10^(DIGITS + 10)，就能把小数运算变成整数运算。

举个例子，arctan(1/5)的第一项1/5乘以SCALE就是SCALE // 5。第二项1/(3*5^3)乘以SCALE就是SCALE // (3 * 5**3)。逐项累加，直到新增项小到对DIGITS位不再有影响。

实现代码如下：

DIGITS = 10000 SCALE = 10 ** (DIGITS + 10) # 多出10位保护 def arctan_inverse(x_inv, scale): """ 计算 scale * arctan(1 / x_inv) 利用级数: arctan(1/x) = 1/x - 1/(3x^3) + 1/(5x^5) - ... """ x = x_inv * x_inv # x^2 term = scale // x_inv # 第一项: scale / x result = term n = 1 while term != 0: term //= x n += 2 if n % 4 == 1: # n = 1, 5, 9... 这一项系数为正？需注意符号 pass if n % 4 == 3: result -= term // n else: result += term // n return result

等等，上面的符号判断有点绕。按标准写法，更清晰的做法是：

def arctan_inverse(x_inv, scale): x = x_inv * x_inv term = scale // x_inv result = term n = 1 while term: term //= x n += 2 # 符号交替：n=3 为减，n=5 为加，n=7 为减... if n % 4 == 3: result -= term // n else: result += term // n return result

然后：

pi_scaled = 4 * (4 * arctan_inverse(5, SCALE) - arctan_inverse(239, SCALE))

这里pi_scaled就是pi * 10^(DIGITS + 10)的整数表示。最后要输出时，只需要先除以10^10得到放大10^DIGITS倍的结果，再转成字符串插入小数点。

3.4 为什么这个版本更“手工”

这个版本里没有Decimal.sqrt，没有任意精度库，全程只用了整数加、减、乘、除（整除）。每一步你都可以在纸上验算，逻辑完全透明。对于想搞懂原理的人来说，这是很好的学习材料。唯一的“性能问题”是：SCALE = 10^10010，每一项的整数都很大，Python的大整数乘法虽然做过优化，但反复除法还是有点耗时。对于10000位这个量级，实测在几秒到几十秒之间，完全能接受。

这也解释了一个常见误区：“高精度计算一定要用高精度库”。其实库本质上是帮你封装了这些整数运算。当需求只是算10000位圆周率时，直接用大整数模拟就够了。

4. 常见问题与排查技巧实录

4.1 为什么算出来的pi前面是3.14159，后面全错了

这是最常见的翻车现场。原因通常是两类：

• SCALE不够大。如果你只放大10^(DIGITS)，但中间运算需要再多几位来保存进位误差，末尾就会错。解决办法就是留出保护位。
• 迭代次数不够。比如Machin级数中，term //= x会越变越小，当term变成0时，说明当前项对最终结果已经没有贡献。如果你在while term != 0的循环条件里退出了，按理说是够的。但如果某个实现里用固定range(N)，而N太小，就会漏项。排查方法是临时打印一下循环结束时term的大小，确认它已经完全为0或足够小。

4.2 字符串插入小数点的位置错了

10000位pi算出来是一个大整数，表示的是pi * 10^(DIGITS+extra)。打印前一定要把这个“额外放大倍率”还原。

一个我自己经常用的做法是：最终整数记为pi_big，它等于floor(pi * 10^DIGITS * 10^extra)。输出时先pi_big //= 10^extra，得到floor(pi * 10^DIGITS)，再转字符串，第一位是整数部分3，后面DIGITS位是小数部分。如果直接拿着最原始的整数去插小数点，你会得到一堆莫名其妙的数字。

4.3 算得太慢怎么办

10000位其实不算大，慢多半是因为实现太粗糙。两个最简单的优化手段，一个是在算法层面：用Chudnovsky而非Machin。Machin公式每项约贡献1.4位，算10000位大约需要7000多次大整数除法，虽然能跑完但明显更吃力。Chudnovsky每项约14位，迭代次数少了一个数量级。

另一个是在数据结构层面：如果你真的在列表中逐位存数字，那么每次加法都要遍历整个列表，10000位还不算大，但如果你哪天想算100万位，这种写法会慢到你怀疑人生。改进办法就是“压位”：每个列表元素存9位或18位整数，而不是存0到9的一位数字，这样列表长度缩短9到18倍，循环次数相应减少。

4.4 末位数字正确吗：误差检验方法

高精度计算最怕的就是“看起来对”但末尾错。这里给出三个实用的校验手段：

1. 查表法：从网上找已知的pi前几万位文本文件，把自己算出的结果的前100位、中段100位、末尾100位拿出来比对。
2. 双算法法：用Chudnovsky和Machin两个不同级数分别计算，比对结果是否完全一致。因为两者逻辑完全不同，同一位数同时出错的概率几乎不存在。
3. 增量法：把保护位从10改成50，重新跑一遍，看前10000位是否保持不变。如果变了，说明原结果的保护位不足，需增加保护位。

有一个很实用的自查点：pi的小数点后前100位是公开的，可以快速验证：

3.14159265358979323846264338327950288419716939937510 58209749445923078164062862089986280348253421170679

拿自己算出来的结果从第1位到第100位对一遍，如果全中，说明大方向对了。

5. 实操总结：从10000位到更大规模的思考

5.1 10000位之后还能玩什么

算完10000位之后，很多人会觉得“就这？”其实这一步的意义在于，你已经完整走了一遍高精度数值计算的流程：公式选型、精度控制、整数化运算、结果校验。

下一步可以尝试几件事：

• 把算法改成并行或分块，挑战100万位。
• 用C语言写一个大整数乘法，把FFT（快速傅里叶变换）引入高精度乘法，体验一下从O(n^2)到O(n log n)的飞跃。
• 把输出结果做可视化、做艺术化展示。

在这些方向上，Chudnovsky算法的地位会越来越高，因为极限计算场景下收敛速度决定一切。

5.2 我个人的实操心得

我最早自己写这个程序时，犯过一个特别低级的错误：忘了给中间结果保留保护位，导致算出来的数最后20位全是乱的。当时我拿标准值一对比，整个人都懵了，排查了大半天才意识到，原来不是公式问题，而是精度不够。后来我养成了一个习惯，就是所有高精度计算在开头先问自己三句话：

1. 我的中间计算需要多少额外保护位？
2. 我的最终结果会放大多少倍？
3. 我的结果校验用什么标准？

把这三个问题想清楚，基本不会翻车。

最后再分享一个小技巧：如果只是临时验证一下某段高精度结果，可以在计算时故意把DIGITS设成100和10000各跑一次，看前100位是否完全一致。如果一致，说明算法稳定；如果不一致，说明存在精度或边界问题。这个方法几乎能定位90%以上的隐藏bug。