学完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 常见π公式的收敛速度对比高精度算π数学公式的选择决定了你后面要写的代码复杂度和运行时间。我能想到的公式大概有这么几类公式收敛速度实现难度适合高精度吗π/4 1 - 1/3 1/5 - 1/7 ...极慢每项只增加约0.4位低不适合算100位要循环上亿次π/4 4·arctan(1/5) - arctan(1/239)较快每轮约增加1.4位中很适合Chudnovsky公式极快每项约14位极高需要大整数阶乘适合超高位但不适合入门莱布尼茨级数虽然代码最简单但收敛速度让人绝望想算到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)^(2k1) / (2k1)项的大小随着k增大指数级衰减。要让它小于10^(-500)粗略估算5^(2k1) 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, ...]。i0tmp 0×10 1 1商0余1。i1tmp 1×10 0 10商2余0。i2及以后全部为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×1000999的量级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^(2k1) / (2k1)那么第k1项是term_(k1) term_k · x² · (2k1) / (2k3)对x 1/5来说x² 1/25所以更新公式可以变成term term × (2k1) / (25 × (2k3))这样每轮迭代只需要一次乘法和一次除法大大减少计算量。这里有个顺序问题到底是先乘后除还是先除后乘我的建议是先乘后除。如果先除第一次除法就会产生截断误差丢掉的信息后面乘回来也补不回来。先乘后除能尽量保留精度代价只是中间结果稍微大一点但远不会溢出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清单我在写和调试过程中遇到过下面几个代表性的问题问题现象原因乘法结果多出很多位数字比预期大几倍乘法进位方向写反从高位往低位计算除法结果全是0算完的小数全是0除法从低位开始处理高位余数没有正确传递结果打印到一半停下来程序崩溃数组越界比如循环用了len而不是len最后几位经常跳动两次运行结果不稳定循环次数不足或者冗余位太少减完结果是负数输出变成“2.99999...”dest和src顺序反了应该用大的减小的还有一个隐蔽的问题在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位结果是中间某一位开始全是错乱的数字排查了很长时间才发现是乘法进位方向的问题。先验证小规模再挑战大规模能省下不少调试时间。把这个项目认真做完再看那些动辄“输出π小数点后十万位”的代码你会觉得它们不过是用更精巧的工具解决同一个朴素的问题。