企业官网建设流程全解析

学完C语言的函数和数组之后，我一直想找一道能把这些知识串起来的综合题。“用C语言高精度计算π的值”就是在这种心态下动手写的。一开始我以为这只是个数学题，查个公式循环累加就行，真动起手来才发现难点根本不在π，而在于你怎么在内存里表示一个超过double精度的小数、怎么对这样的小数做乘法和除法。这篇文章把我完整实现过程、选择Machin公式的理由以及调试时踩过的坑都写在这里，适合已经把C语言基础知识过了一遍、想通过一个综合练习加深理解的读者。

1. 为什么printf("%.20f", pi)打不出真正的π

1.1 浮点数精度：double到底能存多少位

先说一个很多人初学C语言时的误区：以为double能存很多位小数，至少小数点后20位没问题。我当时也这么想，直到我写了这样一段代码：

double pi = 3.14159265358979323846; printf("%.20f\n", pi);

输出结果是：

3.14159265358979311600

前15位对得上，从第16位开始就开始“胡说八道”了。原因在于double在计算机内部是用64位存储的，其中1位符号位、11位指数位、52位尾数位。52位二进制尾数换算成十进制，大约只有15到16位有效数字。

注意这里说的是“有效数字”，不是“小数点后的位数”。有效数字是整个数字从第一个非零数字开始算的位数。对π这种数字来说，double大概只能保证前15到16个数字是正确的，后面全是浮点表示误差。

换句话说，不是printf不舍得打印，是double本身就没存下那么多精确信息。你让一个只认识15位数字的容器去打印100位，它只能把内存里的二进制浮点近似值转换成十进制，再往后就是无法避免的噪声。

1.2 高精度计算的核心思想：把每一位都单独存起来

既然基本数据类型存不下，那就绕开它们。高精度计算的本质很简单：用数组模拟十进制数字的每一位，自己实现加法、减法、乘法和除法。

比如我想保存数值3.1415926535，就定义这样一个数组：

int a[11]; // a[0] = 3 // a[1] = 1 // a[2] = 4 // a[3] = 1 // ...

整数部分占一个元素，之后每一位小数占一个数组元素，每个元素的取值只能是0到9。这样一来，只要数组够长，想存多少位小数都行。

这种思路不只适用于C语言练习。很多大数运算库内部也是同样的想法，只不过它们不会奢侈到用0到9表示一位，而是用一个int表示0到9999甚至更大的“块”，以提高存储效率和运算速度。但对于学习阶段，用十进制数组更容易理解，调试时逐位打印也直观。

真正写代码时你会发现，乘法、除法、加法和减法在数组上的实现细节完全不同，尤其是进位的方向和借位的处理方式。这也是这个练习最有价值的地方：让你彻底搞懂“竖式运算”在程序里是怎么发生的。

2. 选公式比写代码更早：Machin公式凭什么能用

2.1 常见π公式的收敛速度对比

高精度算π，数学公式的选择决定了你后面要写的代码复杂度和运行时间。我能想到的公式大概有这么几类：

莱布尼茨级数虽然代码最简单，但收敛速度让人绝望：想算到100位，至少要迭代几百亿次。Chudnovsky公式是现代超级计算机算π的主流公式之一，但它的每一项都涉及巨大的阶乘和幂运算，对刚接触高精度计算的同学来说，光是把各项的大整数乘除处理好就够头疼了。

Machin公式是性价比很高的选择：数学上不复杂，收敛速度也足够快，用我们手写的数组乘除法就能轻松算出几百位甚至几千位。

2.2 arctan级数与Machin公式的推导

Machin公式长这样：

π/4 = 4·arctan(1/5) - arctan(1/239)

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

arctan(x) = x - x³/3 + x⁵/5 - x⁷/7 + ...

把x分别替换成1/5和1/239，就能得到两个收敛的级数，再按公式组合，就能得到π。

为什么选1/5和1/239而不是直接用arctan(1)？因为泰勒级数收敛的速度取决于|x|的大小，|x|越小，级数收敛越快。

• 如果直接算arctan(1)，x = 1，需要非常多项才能收敛到高精度，几乎不可用。
• 如果算arctan(1/5)，x = 0.2，收敛速度明显变快。
• 如果算arctan(1/239)，x约等于0.00418，收敛飞快。

Machin公式相当于把一个大任务拆成了两个小任务：一个收敛适中，一个收敛极快，组合起来效率远高于直接算。

2.3 迭代次数估算：算500位到底要循环多少轮

写代码之前，最好估算一下循环次数，不然你会不知道terms该设多大。

看arctan(1/5)这一项，第k项大约是：

(1/5)^(2k+1) / (2k+1)

项的大小随着k增大指数级衰减。要让它小于10^(-500)，粗略估算：

5^(2k+1) > 10^500

取对数后得到k大约需要160左右。也就是说，对于500位精度，arctan(1/5)循环160轮就足够了。arctan(1/239)收敛更快，同样500位只需更少的轮数。

实操中我不会卡着理论值算，直接把循环次数设成和数组长度一样，也就是terms = LEN。反正多算几轮不会带来多少额外开销，但能确保精度足够，省去反复调整的麻烦。

用数组长度作为循环次数还有个好处：当LEN是510时，循环510轮，时间开销依然很小。整个程序跑下来也就是毫秒级的事。

3. 一张数组搞定高精度四则运算

3.1 数组布局：下标越小的数位越重要

我采用这样的存储方案：a[0]存整数部分，a[1]存小数点后第1位，a[2]存小数点后第2位，依此类推。也就是数组下标越靠前，位权越高。

这种布局的打印太方便了，直接遍历输出，在a[0]和a[1]之间加一个小数点就行。但它的代价是乘法和加法的进位方向要特别小心，因为进位是从低位往高位走，在数组里就表示从下标大的一端往下标小的一端走。

为了避免每次分配数组空间的麻烦，我会用malloc动态申请，长度定义为LEN。在这个项目里，LEN取510，其中500位是最终要打印的，多出来的10位作为“冗余位”，作用后面会讲。

3.2 高精度除以小整数：竖式除法的还原

高精度除以小整数的核心是模拟我们在小学学过的竖式除法。假设要把一个数组表示的数除以d，做法是从最高位开始，一位一位往下处理：

void div_small(int a[], int len, int d) { int r = 0; for (int i = 0; i < len; i++) { int tmp = r * 10 + a[i]; a[i] = tmp / d; r = tmp % d; } }

每一步中，r是上一位除完剩下的余数，余数乘以10，再加上当前位的数字，构成当前“被除数”。商写入数组当前位置，新的余数继续留给下一位。

举个例子，计算1除以5：

• 数组初始为[1, 0, 0, 0, ...]。
• i=0：tmp = 0×10 + 1 = 1，商0，余1。
• i=1：tmp = 1×10 + 0 = 10，商2，余0。
• i=2及以后：全部为0。

最终数组变成[0, 2, 0, 0, ...]，正好是0.2。

特别注意，除法从高位开始处理，这和乘法正好相反。如果你用从低位往高位的顺序做除法，结果会一塌糊涂。

3.3 高精度乘以小整数：进位方向不能错

高精度乘以一个小整数m的算法也来自竖式乘法，但方向跟除法相反，要从最低位开始：

void mul_small(int a[], int len, int m) { int carry = 0; for (int i = len - 1; i >= 0; i--) { int tmp = a[i] * m + carry; a[i] = tmp % 10; carry = tmp / 10; } }

每位的计算结果要拆成两部分：tmp % 10留下当前位，tmp / 10作为进位传给下一位。因为数组下标越靠前位权越高，所以循环要从len - 1往0走。

新手很容易把顺序写反。如果从a[0]开始，那么进位会覆盖还没处理的低位，结果完全错乱。这点我在第一次写的时候踩过坑，调试时怎么都想不通为什么结果会多出一大截。

在这道题里，m的值都很小，最多也就是迭代次数相关的奇数，比如1000以内；tmp最多是9×1000+999的量级，int完全够用。但如果你把位数放大很多倍，就要重新评估carry会不会溢出。

3.4 加法和减法：最终组装用的两个工具

有了乘法和除法，最后还需要把各项累加起来。加法和减法同样要处理好进位方向：

void add_to(int dest[], int src[], int len) { int carry = 0; for (int i = len - 1; i >= 0; i--) { int tmp = dest[i] + src[i] + carry; dest[i] = tmp % 10; carry = tmp / 10; } } void sub_to(int dest[], int src[], int len) { int borrow = 0; for (int i = len - 1; i >= 0; i--) { int tmp = dest[i] - src[i] - borrow; if (tmp < 0) { tmp += 10; borrow = 1; } else { borrow = 0; } dest[i] = tmp; } }

加法从低位往高位逐位相加，满10就向高位进位；减法从低位往高位逐位相减，不够减就向高位借位。只要你确保dest整体大于src，减法的借位值在循环结束后一定是0。

这四个函数组合起来，就已经构成了一个简单的高精度小数运算工具箱。

4. 把Machin公式翻译成C语言

4.1 每一项都是上一项的变形：迭代更新

有了四则运算工具，接下来要解决的是怎么生成arctan级数的每一项。

arctan(x)的级数是：

x - x³/3 + x⁵/5 - x⁷/7 + ...

如果每一项都从零开始重新算x的幂次，会做很多重复工作。更好的办法是利用相邻两项的关系。

设第k项是：

term_k = x^(2k+1) / (2k+1)

那么第k+1项是：

term_(k+1) = term_k · x² · (2k+1) / (2k+3)

对x = 1/5来说，x² = 1/25，所以更新公式可以变成：

term = term × (2k+1) / (25 × (2k+3))

这样每轮迭代只需要一次乘法和一次除法，大大减少计算量。

这里有个顺序问题：到底是先乘后除，还是先除后乘？我的建议是先乘后除。如果先除，第一次除法就会产生截断误差，丢掉的信息后面乘回来也补不回来。先乘后除能尽量保留精度，代价只是中间结果稍微大一点，但远不会溢出int。

4.2 完整代码

把上面的思路组合起来，就能写出一份可运行的程序。我贴一个能直接编译运行的版本，输出小数点后500位：

#include <stdio.h> #include <stdlib.h> #define LEN 510 void set_int(int a[], int len, int v) { for (int i = 0; i < len; i++) { a[i] = 0; } a[0] = v; } void mul_small(int a[], int len, int m) { int carry = 0; for (int i = len - 1; i >= 0; i--) { int tmp = a[i] * m + carry; a[i] = tmp % 10; carry = tmp / 10; } } void div_small(int a[], int len, int d) { int r = 0; for (int i = 0; i < len; i++) { int tmp = r * 10 + a[i]; a[i] = tmp / d; r = tmp % d; } } void add_to(int dest[], int src[], int len) { int carry = 0; for (int i = len - 1; i >= 0; i--) { int tmp = dest[i] + src[i] + carry; dest[i] = tmp % 10; carry = tmp / 10; } } void sub_to(int dest[], int src[], int len) { int borrow = 0; for (int i = len - 1; i >= 0; i--) { int tmp = dest[i] - src[i] - borrow; if (tmp < 0) { tmp += 10; borrow = 1; } else { borrow = 0; } dest[i] = tmp; } } void arctan_series(int res[], int len, int denom, int terms) { int *term = (int *)malloc(len * sizeof(int)); set_int(term, len, 1); div_small(term, len, denom); int sign = 1; for (int k = 0; k < terms; k++) { if (sign == 1) { add_to(res, term, len); } else { sub_to(res, term, len); } sign = -sign; mul_small(term, len, 2 * k + 1); div_small(term, len, denom * denom * (2 * k + 3)); } free(term); } int main(void) { int *pi = (int *)malloc(LEN * sizeof(int)); int *b = (int *)malloc(LEN * sizeof(int)); set_int(pi, LEN, 0); set_int(b, LEN, 0); arctan_series(pi, LEN, 5, LEN); mul_small(pi, LEN, 4); arctan_series(b, LEN, 239, LEN); sub_to(pi, b, LEN); for (int i = 0; i < LEN; i++) { if (i == 1) { putchar('.'); } putchar('0' + pi[i]); } putchar('\n'); free(pi); free(b); return 0; }

运行后，前几位输出应该是3.141592653589793238462643383279...，和公认的π值完全一致。

4.3 几个容易写错的细节

第一，arctan_series函数里denom * denom * (2 * k + 3)在k较大时会变大。以239为例，当k = 500时，这个除数约等于239×239×1003，大约5700万，仍在int范围内。但如果你把LEN加大到十万级别，这里就要注意溢出了。

第二，mul_small(pi, LEN, 4)这段代码，很多人会问：为什么要先算arctan(1/5)再乘4，而不是把4乘到每一项里？其实两种做法数学上等价。我选择先算完再乘4，是为了让代码结构更清晰：arctan_series只负责算反三角级数，乘4是Machin公式层面的操作，分开放不容易出错。

第三，数组下标从0开始，所以打印时在i == 1处输出小数点。如果你数组长度不够，打印的时候小数点的位置会跟着错位。我建议选LEN时多给自己留一点冗余，后面验证的时候你会感谢这个决定。

5. 验证π：别让小数点后面全是幻觉

5.1 肉眼对照法与小规模调试

代码写完之后，最激动人心的一步就是运行。但是输出一大堆数字，怎么知道它们是对的？我推荐先算一个小规模的结果，比如LEN取20，把输出和已知的3.14159265358979323846逐位对照。这样人工检查不费力，出错也容易定位。

小规模验证通过后，再放大到500位。放大之后可以用网上公开的π小数位数据来核对前100位，至少能确认大方向没问题。

我还用过一个小技巧：把程序输出重定向到文件，然后用文本编辑器打开，和参考文件逐行对比。如果有差异，直接把文件拖到对比工具里看，能快速锁定是哪一位开始不对。

5.2 末位误差的根源：截断误差与冗余位

很多人会问：为什么我算出来的π，最后几位和参考数据对不上？

这是正常的，原因有两个。一个是循环次数不够，最后几项还没小到可以忽略的程度，它们的省略导致末位产生偏差。另一个是每一次除法都会在小数末位引入一点截断误差，这些误差在迭代中不断累积。

解决办法就是在计算时多留几位“冗余位”。比如你想输出500位，就把数组长度设为510。多出来的10位不会打印出来，但它们能吸收中间计算产生的误差，保证打印出来的前500位是稳定的。

LEN = 510就是我在这个目的下选的。如果你想输出1000位，数组长度至少设1010，安全起见1015或1020更好。

5.3 常见Bug清单

我在写和调试过程中，遇到过下面几个代表性的问题：

还有一个隐蔽的问题：在sub_to函数里，如果tmp < 0补10之后，dest[i]一定在0到9之间，不需要再做% 10。有些版本的代码会用tmp % 10，在这个场景下也能工作，但逻辑上不严谨。保留% 10反而容易在负数的取模行为上踩坑。

调试时建议多用小数组、小长度，把中间每一步的数组内容打印出来看。比如计算arctan_series(pi, 10, 5, 10)后，打印一下pi数组，手算验证一下前几位的值，就能迅速定位问题出在乘法还是除法。

6. 从500位到一万位：性能瓶颈与后续进阶

6.1 为什么越往后算越慢

这套实现跑500位非常轻松，但如果你把LEN改成10000，就会明显感到程序变慢了。

原因很简单：这个算法的时间复杂度是O(N²)级别。每一项的乘法和除法都要遍历整个长度为N的数组，而项数本身又和N成正比。于是总开销是N乘N，也就是N²。当N从500变到10000，计算量大约增加了400倍。

在我自己测试的机器上，算500位基本瞬间完成，算5000位大约需要几秒，算10000位就能感受到明显的等待。如果你只想算个一两千位，这个实现完全够用。

6.2 提速方向一：把进制从10改成10000

一个很容易理解的优化是：不要把每个数组元素存0到9，而是存0到9999。

这样做的好处是数组长度变成原来的四分之一，遍历次数大幅减少。与此同时，每一位上的取值不再是单个十进制数字，而是“四位一体”的块。打印的时候需要特殊处理，把一个元素格式化成4位数字输出，不够4位的前面补0。

这个优化能把速度提升好几倍，而且实现思路和十进制版本没有本质区别，非常适合作为下一个练习目标。

6.3 提速方向二：换算法与换策略

如果目标是算十万位甚至百万位，手写的十进制乘除法就不够看了。业界常用的方案包括：

• 用快速傅里叶变换（FFT）加速大整数乘法，把乘法复杂度从O(N²)降到O(N log N)。
• 换用Chudnovsky公式，配合二进制分裂算法，大幅减少级数迭代次数。
• 使用流式算法（比如Spigot算法），一边算一边输出π的每一位，不需要一次性维护整个数组。

这些方向每一个都值得单独写一篇长文。对我们这个练习来说，知道它们的存在、知道当前实现的上限在哪里，就已经很有价值了。

6.4 作为C语言综合练习题，它到底练了什么

这个项目让我对C语言的理解加深了不少，主要在三方面：

第一，数组作为函数参数传递时，长度信息会丢失。你必须在外面把长度传进去，不然函数里根本不知道数组有多大。这和我们平时写for (int i = 0; i < len; i++)的习惯不一样，一旦漏传长度，程序就会越界。

第二，内存分配一定要记得释放。代码里用了malloc，用完free，这是很多人初学时容易忽略的。虽然程序结束系统会回收内存，但养成释放的习惯对以后写长时间运行的程序很重要。

第三，算法选择和优化策略往往比代码细节更重要。同样算π，莱布尼茨级数和Machin公式的运行时间天差地别。多掌握几种实现思路，遇到实际项目时才能选出最合适的方案。

最后分享一个小技巧：算π的时候别急着打印全部位数，先把计算和打印分两步，把结果存到一个数组里，确认前50位完全正确后再把LEN放大。我当年就是没这么做，一上来直接算500位，结果是中间某一位开始全是错乱的数字，排查了很长时间才发现是乘法进位方向的问题。先验证小规模，再挑战大规模，能省下不少调试时间。把这个项目认真做完，再看那些动辄“输出π小数点后十万位”的代码，你会觉得它们不过是用更精巧的工具，解决同一个朴素的问题。