大一那年的C语言课程设计我选了个听起来很传统的题目高精度计算π的值。当时我以为写个循环套公式就行结果用double算到小数点后第17位后面的数字就开始“发疯”。后来才明白这根本不是数学公式的问题而是C语言自身的浮点精度扛不住——double用52位二进制尾数换算成十进制也就是15到17位有效数字想输出1000位π相当于要求一个人用肉眼数清一碗米的粒数。这个题目其实特别适合当C语言的综合练手项目它逼着你把数组、指针、函数、内存布局、运算顺序、文件读写这些零散知识点串起来还要自己做算法选型和边界测试。无论你是刚学完C语言基础还是在准备计算机二级、课程设计只要能完整跑通这个项目你对“C语言到底能干什么”的理解会提升一个档次。这篇文章就按我实际调试的路径从原理到代码再到坑位完整拆一遍。1. 高精度计算π到底在跟什么较劲1.1 double的精度墙16位之后只剩幻觉先做个直观实验。下面这段代码很多初学者都写过#include stdio.h int main() { double pi 0.0; int sign 1; for (int i 0; i 1000000; i) { pi sign * 4.0 / (2 * i 1); sign -sign; } printf(%.20f\n, pi); return 0; }这是用莱布尼茨级数算π跑100万项输出的结果前面几位是3.141591...后面全是错的。问题有两个层面第一个是这个级数收敛极慢百万项只能准到小数点后5位第二个更本质就算换一个收敛快的公式double的精度天花板也摆在那里。IEEE 754标准里double有53位二进制有效位大约只能精确表达15到17位十进制数字。你别看printf能打出20位小数的样子那是格式化输出在“补妆”真实计算时的舍入误差早就把后面污染了。所以“高精度计算π”这个需求天然就超出C语言内置数值类型的射程。想算到100位、1000位、1万位唯一的办法是自己造一套能表示任意多位数字的存储结构。1.2 高精度计算的核心用数组模拟任意精度既然一个double装不下1000位小数那就把1000位拆开用数组一个个装着。这是所有高精度计算的底层思想。最简单的方案是数组里每个元素存一位十进制数字0到9。比如3.14159可以存成一组int{3,1,4,1,5,9}。但这种方案有两个问题一是空间浪费一个int用4字节只存0到9内存利用率只有不到1%二是运算效率低模拟竖式加法时一位位处理1000位就是1000次循环乘法则需要双重循环达到百万次级别。更实用的做法是“压位”让一个数组元素尽量多存几位。我后面会详细讲10^9进制方案也就是一个int存9位十进制数字。这样算1000位π数组长度只需要1000/9约112个元素运算次数直接砍到原来的约1/81。你可以把数组想成一个超大整数把π的每一位数字当成这个整数的一部分所有计算都在这个自建的大数系统里完成。1.3 为什么说“C语言没有高精度库”反而是优点有人会问Python的decimal模块、Java的BigDecimal、甚至C的boost::multiprecision都能处理高精度为什么非要用C语言自己写搜索热词里也有“为什么高精度计算不自己写”这种问题。我的看法是如果你要的是工程结果当然可以直接调库但如果你要的是理解计算机如何处理数值、如何管理内存、如何做算法优化手写一遍就是最好的训练。C语言的数组是一个连续内存块你对每一位的操作本质上是在做内存地址偏移运算。这种“裸”的感觉在其他高级语言里很难体会到。另外高精度计算π是一个非常干净的算法题目数学公式固定、验证结果容易网上随便查π的前1000位来对比、性能瓶颈清晰乘法最多、也没有复杂的业务逻辑干扰。它就像算法界的“hello world”看着简单做透却要过好几关。2. 动手前先选好数学公式2.1 无穷级数、反三角函数、迭代法三条路线高精度计算π的数学公式非常多但归纳起来主要三条路线。第一条是无穷级数比如莱布尼茨级数π/4 1 - 1/3 1/5 - ...。优点是公式简单缺点也致命收敛太慢算1000位需要几百万项。第二条是反三角函数变换核心是Machin公式这类基于arctan加法定理的公式。它把π表示成几个反正切函数的线性组合每个arctan(1/x)里的x越大对应项衰减越快收敛速度远超莱布尼茨级数。第三条是迭代法代表是Gauss-Legendre算法和Borwein兄弟算法它们通过反复迭代逼近π收敛速度指数级但每一步迭代都涉及大数乘法、除法甚至开方代码复杂度高。如果目标只是课程设计或者学习演示Machin公式系列是性价比最高的数学上不晦涩实现上只需要高精度的“乘小整数”和“除小整数”两个操作就能稳定输出上千位。2.2 Machin公式收敛速度与实现难度最平衡Machin公式长这样π/4 4·arctan(1/5) - arctan(1/239)为什么选它不选别的关键在于两个底数5和239。先看arctan的泰勒展开arctan(x) x - x³/3 x⁵/5 - x⁷/7 ...把x1/5代入每一项大约以1/25的速度衰减也就是每算一项精度大约增加1.4位十进制数字。算到700多项就能得到1000位精度。再看x1/2391/239²≈1/57121每一项衰减极快每项贡献约4.8位十进制精度算200多项就够了。两者加起来大概900多次迭代每次迭代只是几次大数乘除整体计算量在毫秒级。这就是Machin公式的妙处把一个大问题拆成两个收敛速度都很快的小问题而且公式里只出现整数除法不涉及开方这类复杂运算。相比之下其他一些“看起来很美”的公式比如反正切函数组合公式π/4 arctan(1/2) arctan(1/3)收敛速度就慢很多不太适合做高位数计算。2.3 算到多少位需要多少项先做估算写代码之前先估算循环次数可以避免盲写。对Machin公式里的arctan(1/5)那一串第n项的绝对值大约是:1 / (5^(2n1) · (2n1))我要让最后一项小于10^(-1000)取个对数估算(2n1)·log10(5) log10(2n1) 10002n1大约等于1000/0.699 ≈ 1431考虑到第二项log10(2n1)约3.2n大约在714左右。对arctan(1/239)类似可得n大约在210左右。所以循环上限设个5000轮完全足够还可以提前判断当前项归零就退出。这个估算过程很有用它能让你提前知道程序要跑多久也方便后面测试不同位数时心里有数。3. C语言手写大数运算从一行int到一个数组3.1 存储方案设计为什么用10^9进制现在进入真正的C语言实现环节。我用的存储结构是#define DIGITS 1000 #define BASE 1000000000 #define LEN ((DIGITS 8) / 9 2) typedef int Big[LEN];一个int数组每个元素存0到999999999之间的数也就是10^9进制。数组下标0是最低位下标LEN-1是最高位。这个数和普通十进制数的关系是value a[0] a[1]·10^9 a[2]·10^18 ...为什么选用10^9进制而不是10进制三个原因。第一int能存的最大值是2147483647如果选10^9进制两个数相乘的中间结果可能接近10^18这已经超过int范围了但刚好在long long9.2×10^18之内。如果选10^10进制int就爆了选10^4进制又太浪费。第二压位能大幅减少循环次数。1000位小数在10进制下一个数组要1000个元素乘法要双重循环实现100万次操作换成10^9进制只需要112个元素乘法虽然理论上复杂度还是O(n²)但常数小得多。第三printf格式化这里占便宜。我后面输出时用%09d能自动补前导零正好把每个数组元素展开成9位十进制数。这里要留一个余量LEN比实际需要的多2个元素。原因稍后讲先记住这个细节后面调试“末尾位数不对”时你就知道多重要的。3.2 加、减、乘小整数、除小整数的实现高精度计算的四则运算我只需要四个基础函数就够了其中核心是乘小整数和除小整数。大数加小整数这里其实是两个大数相加void add(Big a, const Big b) { long long carry 0; for (int i 0; i LEN; i) { long long cur (long long)a[i] b[i] carry; a[i] (int)(cur % BASE); carry cur / BASE; } }大数减小整数void sub(Big a, const Big b) { long long borrow 0; for (int i 0; i LEN; i) { long long cur (long long)a[i] - b[i] - borrow; if (cur 0) { cur BASE; borrow 1; } else { borrow 0; } a[i] (int)cur; } }这两个函数和普通十进制竖式一模一样只不过满1000000000进一位。乘小整数比如大数乘以16void mul_small(Big a, int m) { long long carry 0; for (int i 0; i LEN; i) { long long cur (long long)a[i] * m carry; a[i] (int)(cur % BASE); carry cur / BASE; } }除小整数同时返回余数int div_small(Big a, int d) { long long rem 0; for (int i LEN - 1; i 0; i--) { long long cur (long long)rem * BASE a[i]; a[i] (int)(cur / d); rem cur % d; } return (int)rem; }注意除法的循环方向是从高位到低位和加法乘法相反。这是模拟手算长除法每次把上一位的余数乘上BASE再加到当前位再除以除数。如果你把方向写反结果会完全错误这是新手最容易犯的错。3.3 最容易翻车的三个边界情况第一个是乘法加法里的类型问题。a[i]本身是int乘以m后很可能超过int范围所以必须先把a[i]提升成long long再运算。我之前见过很多同学直接写cur a[i] * m carry结果m稍微大一点就溢出成负数程序输出乱七八糟。第二个是sub函数里borrow的处理。当cur 0时要加上BASE并设置借位但cur最小可能是-(BASE-1)-1-BASE加上BASE后变成0这没问题。不过如果你把a和b的关系搞反也就是a b时结果会变成负数而且数组里没有负数表示所以调用sub前要确认左边大于右边。第三个是数组越界。我故意设计了“乘小整数”的进位可能一直传到数组的最高位所以LEN要留余量。如果LEN刚好等于需要的位数乘法进位时可能把最高位挤出去数据就丢了。留2个元素的余量等于给你的进位和误差一个缓冲。4. 实操输出1000位π的完整代码4.1 完整C语言程序下面这个程序可以直接编译运行用Machin公式计算π的前1000位小数并输出到一个文本文件pi.txt。代码逻辑分成三部分大数运算基础函数、arctan(1/x)函数、主函数。#include stdio.h #include string.h #define DIGITS 1000 #define BASE 1000000000 #define LEN ((DIGITS 8) / 9 2) typedef int Big[LEN]; void clear(Big a) { memset(a, 0, sizeof(int) * LEN); } void set_one(Big a) { clear(a); a[LEN - 1] 1; } int is_zero(Big a) { for (int i 0; i LEN; i) if (a[i] ! 0) return 0; return 1; } void add(Big a, const Big b) { long long carry 0; for (int i 0; i LEN; i) { long long cur (long long)a[i] b[i] carry; a[i] (int)(cur % BASE); carry cur / BASE; } } void sub(Big a, const Big b) { long long borrow 0; for (int i 0; i LEN; i) { long long cur (long long)a[i] - b[i] - borrow; if (cur 0) { cur BASE; borrow 1; } else { borrow 0; } a[i] (int)cur; } } void mul_small(Big a, int m) { long long carry 0; for (int i 0; i LEN; i) { long long cur (long long)a[i] * m carry; a[i] (int)(cur % BASE); carry cur / BASE; } } int div_small(Big a, int d) { long long rem 0; for (int i LEN - 1; i 0; i--) { long long cur (long long)rem * BASE a[i]; a[i] (int)(cur / d); rem cur % d; } return (int)rem; } // 计算 dest One * arctan(1/x) void arctan(int x, Big dest) { Big term; set_one(term); div_small(term, x); clear(dest); add(dest, term); int odd 1; // 当前分母中的奇数项 int x2 x * x; int idx 0; // 第几项从0开始 while (!is_zero(term) idx 5000) { mul_small(term, odd); div_small(term, x2); odd 2; div_small(term, odd); idx; if (idx % 2 0) add(dest, term); else sub(dest, term); } } void print_pi(FILE *fp, const Big a) { fprintf(fp, %d., a[LEN - 1]); int printed 0; for (int i LEN - 2; i 0 printed DIGITS; i--) { if (printed 9 DIGITS) { fprintf(fp, %09d, a[i]); printed 9; } else { int need DIGITS - printed; int p 1; for (int j 0; j need; j) p * 10; fprintf(fp, %0*d, need, a[i] / (BASE / p)); printed need; } } fprintf(fp, \n); } int main() { Big pi, t1, t2; arctan(5, t1); arctan(239, t2); mul_small(t1, 16); mul_small(t2, 4); sub(t1, t2); FILE *fp fopen(pi.txt, w); if (fp NULL) { printf(无法创建文件\n); return 1; } print_pi(fp, t1); fclose(fp); printf(前1000位π已写入 pi.txt\n); return 0; }这段代码在GCC和VS Code配置好的C环境下都能直接运行Linux下用gcc pi.c -o pi编译Windows下在Visual Studio或MinGW都行。4.2 核心过程逐段拆解首先整体思路想把π算成高精度不能直接拿π的泰勒展开算而是先构造一个很大的“单位数”One它等于10的9×(LEN-1)次方。然后所有计算都基于这个One做整数定点运算。arctan函数里dest实际上存储的是One × arctan(1/x)从数学上看这是把小数放大成整数所有运算都是整数运算避免浮点误差。arctan(1/x)的每一项系数被一步一步算出来后放到dest里累加。看arctan里的循环逻辑。第一项是term One/x直接存进dest。接着进入循环mul_small(term, odd)这里的odd初始是1。这一步其实是把当前项的分母“去掉”一个奇数因子准备让它除以x²。div_small(term, x2)这一步是关键它把当前项除以x²让分母的指数从5的1次方变到3次方、从3次方到5次方依次推进。odd加2再div_small(term, odd)把新的奇数分母乘进去。交替加减完成符号变换第0项是正第1项是负第2项是正。你会发现这个循环里没有计算“当前项是否为0”的任何外部判断而是用is_zero(term)来判断。当term所有位都是0时说明当前项已经小于10的负几千次方它对目标位数没有影响了。最后主函数里面t1 arctan(1/5)t2 arctan(1/239)t1乘以16t2乘以4然后相减得到π的高精度定点表示。print_pi函数负责把t1按10^9进制展开成十进制数字输出。4.3 把结果写入文件验证精度在实际运行时1000位π直接打印到控制台滚动太快不方便核对所以我把它写到了pi.txt。你可以用任意文本编辑器打开和网上查到的π前1000位对比。验证方法我推荐这种你先找一段公认的前200位π存成标准字符串然后写一个小C程序逐字符比较或者直接把两个文件拖进Beyond Compare。如果只靠肉眼在终端里扫很容易漏掉某个数字的差别。我实测手写这段代码在本地跑1000位运行时间在几十毫秒级别几乎感觉不到卡顿。你可以在pi.txt里看到3.141592653589793238462643383279502884197169399375105820974944592...这串数字完全对得上。能对上意味着你的大数运算和迭代公式都写对了后面想加位数只要改DIGITS再重新编译就行。5. 更高精度的路线Gauss-Legendre与更快的算法5.1 Gauss-Legendre迭代原理与代码难点Machin公式算1000位很轻松但想算到10万位、100万位就不太够看了迭代次数随位数线性增长始终是个瓶颈。Gauss-Legendre算法则完全不同它的核心是四个递推式a₀ 1b₀ 1/√2t₀ 1/4p₀ 1然后重复迭代aₙ₊₁ (aₙ bₙ) / 2bₙ₊₁ √(aₙ · bₙ)tₙ₊₁ tₙ - pₙ · (aₙ - aₙ₊₁)²pₙ₊₁ 2 · pₙ最后π的近似值是π ≈ (aₙ₊₁ bₙ₊₁)² / (4 · tₙ₊₁)这个算法的恐怖之处在于收敛速度每迭代一次精度大约翻一倍。算1000位也许只需要10次迭代算100万位也就20多次。但难点全在“√”上面你必须自己实现一个大数开平方函数。开平方通常用牛顿迭代而牛顿迭代每轮又需要大数除法除法复杂度比乘法高所以整体实现难度直线上升。5.2 Chudnovsky公式为什么能破百万位如果目标是“百万位以上的π”目前的主流选择是Chudnovsky公式。它是基于模椭圆函数推导出来的级数每迭代一项大约能增加14位十进制精度。计算100万位π只需要约73000项配合高性能的大数乘法比如FFT乘法在个人电脑上也能几分钟内算完。Chudnovsky公式的表达式1/π 12 · Σₖ₌₀^∞ (-1)ᵏ (6k)! (13591409 545140134k) / ((3k)! (k!)³ 640320^(3k3/2))这个公式长得吓人而且涉及阶乘、大数乘大数、甚至开平方因为分母里有640320的3k3/2次方代码量比Machin高一个量级。所以我建议初学者先把Machin吃透理解了大数数组和高精度加减乘除以后再考虑挑战Chudnovsky。5.3 如果只是想快速出结果用现成库行不行如果你做这个项目不是为了学习而是为了“得到一个高精度π值”那完全可以不用自己造轮子。GNU MPFR库支持任意精度浮点libgmp提供任意精度整数在C语言里调它们可以一行算出几万位π。甚至有时候我调试自己的高精度代码也会先用Python的decimal模块算出标准答案再和我C程序输出的文件对比。工具没有高低贵贱知道什么时候该手写、什么时候该用库本身就是工程能力的一部分。不过话说回来这一篇的核心价值恰恰在“手写”这两个字上。用库当然方便但你自己模拟一遍大数运算之后再看计算机组成原理里的浮点数、内存对齐、溢出处理都会有“原来如此”的顿悟感。6. 常见问题与调试记录6.1 为什么末尾几位不稳定我第一次跑这段代码时输出的前990位都对但最后十位和标准值有差异。排查下来有两个原因一是LEN余量不足导致乘法进位时把最高位挤掉了二是定点除法的余数被丢弃末位会有1到2个单位的误差。解决办法很简单LEN保留2个元素的余量DIGITS只用来控制输出不参与存储分配。这样即使末尾有误差也被截断在输出之外了。6.2 程序跑得很慢瓶颈在哪里如果你把DIGITS改到10000会发现运行时间明显变长。瓶颈几乎都在arctan循环里的mul_small和div_small上它们每次都要完整扫一遍整个数组。想提速方向有两个一是把10^9进制改成更大的基比如用long long存更多位但要小心溢出需要拆成更细的乘法二是对大数乘法做优化简单的可以用分块乘法进阶的用Karatsuba或FFT。对我这个1000位的小规模项目来说这些优化暂时用不上但作为知识储备是很好的。6.3 这个题目还能怎么练算完π之后你可以顺手做几个变形把C语言基础打得再扎实一点。其一改成通过命令行参数接收位数而不是改宏重新编译这能练到main(int argc, char *argv[])和字符串转数字。其二用结构体封装大数把clear、add、sub、mul_small、div_small都变成操作结构体的函数。这等于从“面向过程”迈向“封装”为以后学C做铺垫。其三写一个自动校验函数把结果和标准π字符串逐字节比较如果发现不一致就定位到第一个错误位这个功能能让你改代码时心里踏实。我自己当年做完这个题目之后最大的收获不是“我会算π了”而是终于搞懂了内存和循环之间那种密切的关系一次乘法进位的传播方向、一次除法余数的传递顺序这些在课本上念一百遍不如自己跑一遍代码来得直观。如果你也正在被C语言的指针和数组折磨不妨试试这个项目它比背100个代码片段有用得多。