先把结论说在前面如果你只是想用C语言拿到3.14159265358979那直接调用math.h就够了double精度能给你大约15位有效数字。可一旦你想把π算到几百位、上千位比如验证“山巅一寺一壶酒尔乐苦煞吾”这句谐音诗是否真的对应圆周率小数点后22位double立刻就不够用了。这事只能靠高精度计算来解决而用纯C语言手写高精度计算恰恰是理解数组、进位、余数的绝佳练习题。这篇文章我打算完整拆解一个只靠标准C库实现的高精度π计算程序不依赖GMP、不依赖任何第三方大数库从存储方案、公式选型到完整代码最后再聊几个实际调试中踩过的坑。不管你是正在刷翁恺老师练习题的初学者还是想把“大数运算”这件事彻底搞明白的C语言玩家这份笔记都值得你跟着敲一遍。1. 别急着写代码先搞懂为什么要高精度1.1 double的天花板在哪里不少同学第一次看到这道题时第一反应是直接用double循环累加不就行了比如Leibniz级数π/41-1/31/5-1/7...写个for循环几千万次最后printf(%.15f)好像也能出个差不多的值。但问题在于double在内存里是IEEE 754双精度浮点数只有52位有效尾数换算成十进制大约是15到17位有效数字。这16位里整数部分“3”已经占了一位后面最多还能给出15位左右的小数。你想看第100位是几double根本表示不了。所以这时候就必须换思路不依赖浮点数用整数数组一位一位地存小数。这就是高精度计算的核心思想。1.2 高精度本质上是把竖式搬进数组你还记得小学列竖式做加减法吗个位满十向十位进一十位满十向百位进一。高精度计算其实完全是同一件事只是把“满十进一”改成了“满一万进一”“满一亿进一”把进制的基数设成10000或者更大。数组的下标就是数位每一位存一个0到BASE-1之间的整数进位、借位规则和竖式一模一样。这么说可能有点抽象举个最简单的例子。用十进制数组a[4]存1234就是a[0]1,a[1]2,a[2]3,a[3]4。如果我想给这个数加上9000先从个位算a[3]04十位303百位202千位1910超过9所以要向“万位”进位。如果数组没有留万位那就得扩容或者提前预留空间。高精度计算的所有麻烦本质上都围绕“进位往哪传、借位往哪借、数组留多少位”这三个问题展开。1.3 这道题到底在练什么“高精度计算π的值”是C语言里非常经典的一道综合练习你可能会在翁恺老师的练习题列表里见到它也可能在“C语言必背100代码”之类的合集里看到类似题目。它表面上考的是圆周率实际上考的是三件事第一数组边界的控制能力遍历时多一个等号可能就数组越界第二进位和借位的逻辑方向写反直接全盘崩溃第三把一个数学公式拆成程序的模块化能力。这三样东西是你后面写链表、写大整数运算、写嵌入式协议解析都逃不掉的基本功。2. 公式选型为什么偏偏选Machin公式2.1 三个常见公式的对比计算π的方法从古到今有几十种但适合我们这种手写高精度库的并不多。我拿三个最典型的公式对比过直接看表格公式收敛速度实现难度适合规模Leibniz级数很慢每项大约只增加0.3位小数极其简单适合演示不适合算高精度Machin公式主项每项增加约1.4位副项每项增加约4.7位中等只需要高精度除法几百位到几千位非常舒服Chudnovsky级数每项增加约14位高需要实现高精度开方和乘除万位以上性能才有优势你可能听说过Leibniz公式毕竟代码只有三行但它的收敛速度实在太感人想算到1000位需要循环大约数千亿次这在单机上是灾难。Chudnovsky公式是目前纪录级别的算法每项能多出14位精度但它的每一项里都带着复杂的平方根这意味着你还要先写一个高精度开平方函数题目难度瞬间飙升。Machin公式就卡在中间位速足够手写实现时只需要“高精度除以小整数”不用做高精度乘除法正好适合我们这篇博文的定位。2.2 Machin公式到底是个什么来头Machin公式长这样π/4 4·arctan(1/5) − arctan(1/239)也就是说π 16·arctan(1/5) − 4·arctan(1/239)这个公式看着像变魔术实际可以用正切的和角公式验证。先设tan α 1/5用倍角公式连续翻两次能得到tan 4α 120/119。再设tan β 1/239用差角公式算tan(4α − β)算出来的结果恰好等于tan(π/4)1。由于4α−β落在合适范围内所以4α−βπ/4。整个过程只需要高中三角函数知识但是数学家John Machin在1706年发现它时可是直接用到了手算π位数的顶峰——算到100位考虑到当时没有计算机这个公式的厉害之处你应该能感觉到了。为什么这个公式收敛快关键在于两项的分母。arctan(1/5)的级数展开中第k项分母里有5的奇数次幂每往后一项分母大约变成原来的25倍换算成小数位就是每位约增加1.4位而arctan(1/239)的分母是239的奇数次幂每项增加约4.7位。这就意味着要算1000位小数只需要把arctan(1/5)展开到大约720项arctan(1/239)展开到大约220项完全在可接受范围内。2.3 泰勒展开怎么和数组结合起来arctan(x)的泰勒展开是arctan(x) x − x³/3 x⁵/5 − x⁷/7 ...把x1/5代入得到arctan(1/5) 1/5 − 1/(3·5³) 1/(5·5⁵) − 1/(7·5⁷) ...注意到没有每一项都是“某个整数分之一”不存在浮点数。如果我们用数组表示一个高精度小数那么这个级数的每一项都能通过“高精度整数除一个普通整数”来获得。这正是Machin公式适合手写高精度库的核心原因整个程序里最复杂的运算只是高精度的除以单精度整数而不需要实现完整的高精度除法。3. 高精度运算库核心就四个函数3.1 数字怎么存为什么选万进制我见过有人用十进制数组一个元素存0到9的一位。这个方案写起来直观但效率低得离谱同样表示一个1000位的小数十进制数组需要1000个元素而万进制数组只需要250个元素。我这里选的基数是10000也就是万进制。数组的每个元素取值0到9999正好对应十进制里的4位。这相当于把十进制数字每4个压缩成一组比如0.14159265在万进制数组里存成a[0]0, a[1]1415, a[2]9265。这样做的第一个好处是省内存第二个好处是输出时也不用一位位拼printf(%04d, a[i])就能直接打出一组4位数字非常方便。数组下标约定也很关键a[0]存整数部分a[1]存小数点后第1到第4位a[2]存第5到第8位依此类推。比如最终π的整数部分是3那么a[0]3你直接输出3然后从a[1]开始按四位一组输出就行了。3.2 高精度除以单精度整数整个程序的核心这个函数是整个项目里最值得反复琢磨的一段。我们的目标是让一个高精度数组a除以一个普通int整数x也就是a a / x。怎么做到回忆一下小学竖式除法从最高位开始每一位除以除数得到商把余数乘上进制传到下一位再加上下一位原有的数字继续除。翻译成代码就是void div_int(int a[], int x) { int i, rem 0; for (i 0; i PREC; i) { int cur rem * BASE a[i]; a[i] cur / x; rem cur % x; } }这里的rem就是上一轮除完剩下的余数。因为万进制下满10000进一所以余数传给下一位时要乘上BASE也就是乘10000。我最开始写这个函数时总是习惯从低位往高位循环因为加法和乘法都是从低位开始的。但除法恰恰相反余数永远是从高位往低位流动的所以必须从高位下标0往低位去。你要是把方向写反了算出来的结果会完全不对而且很难查。这个函数的复杂度是O(n)n是数组长度。后面整个π的计算本质上就是在反复调用它所以它的效率基本决定程序的整体表现。3.3 加法、减法、乘法的实现要点加减乘虽然简单但进位方向同样是容易翻车的地方。加法从数组末尾最低位往开头最高位遍历每一位相加再加上上一位的进位。如果某个位置满了BASE就向高一位进1。void add(int a[], int b[]) { int i, carry 0; for (i PREC - 1; i 0; i--) { int cur a[i] b[i] carry; a[i] cur % BASE; carry cur / BASE; } }减法同样从低往高遍历每一位相减再减去借位。如果不够减就向高一位借1注意借来的1在本位相当于BASE。void sub(int a[], int b[]) { int i, borrow 0; for (i PREC - 1; i 0; i--) { int cur a[i] - b[i] - borrow; if (cur 0) { cur BASE; borrow 1; } else { borrow 0; } a[i] cur; } }乘一个小整数比如最后把结果乘4的方向和加法一样也是从低位往高位乘出来的carry一路向上传void mul_small(int a[], int x) { int i, carry 0; for (i PREC - 1; i 0; i--) { int cur a[i] * x carry; a[i] cur % BASE; carry cur / BASE; } }这四个函数加起来不超过50行但你把它们组合好就拥有了一台“高精度小数计算器”。3.4 固定精度数组的两个硬性提醒使用固定长度数组时有两点必须提前想好。第一数组要比你最终要输出的位数多留几个元素。因为级数计算会产生末尾误差如果只准备刚好够用的位数最后几位很可能是不正确的。我一般会在需要的位数基础上多加4个万进制元素相当于多算16位然后最后做一次四舍五入再截断。第二乘法进位可能会让数组最前面多出一个超出预期的元素。比如乘法进位一路传到a[0]而a[0]已经接近9999那么a[0]可能溢出。在π这道题里整数部分最大也就是3或4不会出现这种情况但如果你后续拿这个库算别的数一定要考虑进位到数组越界的问题。4. 完整实现从arctan到输出1000位π4.1 arctan函数怎么利用递推前面说过arctan(1/x)的级数每一项都是分数但直接每一轮都重新从1开始除以x的幂会做大量重复除法。更聪明的做法是用递推。设当前项为term_k 1 / ((2k1) · x^(2k1))我们可以维护一个只关于x的p 1/x^(2k1)。每一轮开始p其实就等于term_k乘以分母中的奇数部分。所以流程是先把p复制到term里让term除以当前奇数odd得到真正的当前项然后把这一项加到最终结果里最后让p除以x²得到下一轮需要的1/x^(2k3)同时odd加2。写成伪代码就是p 1先除以x得到p 1/x当前奇数odd 1term p然后term除以odd如果sign为正result term否则result - termp除以x²odd加2sign取反重复第3到第5步直到p的所有位都是0这个递推比从头逐项计算快很多而且逻辑清楚。整个while循环跑多少轮呢前面估算过算1000位时arctan(1/5)大约要跑720轮arctan(1/239)大约跑220轮每轮内部是做两次高精度除法总运算量非常小。4.2 完整代码可以直接抄走下面这段代码我用固定数组实现目标输出1000位小数#include stdio.h #include string.h #define DIGITS 1000 // 要输出的十进制小数位数 #define BASE 10000 // 万进制 #define PREC (DIGITS / 4 4) // 数组总长度多留4组防止末尾误差 void div_int(int a[], int x) { int i, rem 0; for (i 0; i PREC; i) { int cur rem * BASE a[i]; a[i] cur / x; rem cur % x; } } void add(int a[], int b[]) { int i, carry 0; for (i PREC - 1; i 0; i--) { int cur a[i] b[i] carry; a[i] cur % BASE; carry cur / BASE; } } void sub(int a[], int b[]) { int i, borrow 0; for (i PREC - 1; i 0; i--) { int cur a[i] - b[i] - borrow; if (cur 0) { cur BASE; borrow 1; } else { borrow 0; } a[i] cur; } } void mul_small(int a[], int x) { int i, carry 0; for (i PREC - 1; i 0; i--) { int cur a[i] * x carry; a[i] cur % BASE; carry cur / BASE; } } int is_zero(int a[]) { int i; for (i 0; i PREC; i) { if (a[i] ! 0) { return 0; } } return 1; } void copy_arr(int dst[], int src[]) { int i; for (i 0; i PREC; i) { dst[i] src[i]; } } // 保留到第 keep 个小数位组并根据后一组做四舍五入 void round_array(int a[], int keep) { int i; if (a[keep 1] BASE / 2) { for (i keep; i 0; i--) { a[i]; if (a[i] BASE) { break; } a[i] 0; } } for (i keep 1; i PREC; i) { a[i] 0; } } // 计算 arctan(1/x)结果累加到 result 数组 void arctan(int x, int result[]) { int temp[PREC], term[PREC]; int odd 1, sign 1; memset(temp, 0, sizeof(temp)); temp[0] 1; div_int(temp, x); while (!is_zero(temp)) { copy_arr(term, temp); div_int(term, odd); if (sign 0) { add(result, term); } else { sub(result, term); } div_int(temp, x * x); odd 2; sign -sign; } } int main(void) { int atan5[PREC] {0}; int atan239[PREC] {0}; int i; arctan(5, atan5); arctan(239, atan239); // pi 16*arctan(1/5) - 4*arctan(1/239) mul_small(atan5, 4); sub(atan5, atan239); mul_small(atan5, 4); round_array(atan5, DIGITS / 4); printf(pi %d., atan5[0]); for (i 1; i DIGITS / 4; i) { printf(%04d, atan5[i]); } printf(\n); return 0; }这段代码在VS Code配好的C环境中可以直接编译运行也兼容大多数主流IDE。输出会是一大串数字你在屏幕上可能看不到全部建议把结果重定向到文件里再对比。4.3 怎么确认算出来的值没问题跑完代码以后先别急着相信结果。最靠谱的验证方式是拿已知的π前100位对比3.1415926535 8979323846 2643383279 5028841971 6939937510 5820974944 5923078164 0628620899 8628034825 3421170679如果你的程序输出的前100位和这个完全一致说明至少前100位没问题。如果你想验证更多位可以和在线圆周率网站比对或者用Python的decimal库算一份出来对照。我自己跑上面这段代码1000位输出里前999位都是稳定的最后一位因为截断方式不同可能出现很小的偏差这属于正常现象。4.4 性能实测印象以我笔记本上的实测结果来说1000位几乎是瞬间完成肉眼看不到等待把DIGITS改成5000也能在一两秒内跑完到20000位耗时就会明显上来大约需要十几秒甚至更久毕竟这个实现的核心循环是O(n²)级别的。如果只是做练习题1000位完全够用了。5. 调试实录这几个坑我全都踩过5.1 算到一半全是0问题出在除法方向第一次独立写div_int时我潜意识里抄了加法的遍历方向从数组末尾往开头除。结果就是最高位除完后的余数被丢到了“下一位”但实际上“下一位”在数组中已经遍历过了余数根本没有传递下去。最后输出的数字从某一位开始全是0看起来就像级数突然断掉了。如果你遇到类似现象先检查div_int里是不是从高位到低位循环也就是for(i 0; i PREC; i)。5.2 减法算出负数借位逻辑写错了sub函数里的borrow处理有个容易忽视的细节如果cur小于0要先加上BASE然后把borrow置1。很多初学者会在else分支里忘记把borrow清零导致高一位莫名其妙减了2。这个问题在单独测试sub函数时可能看不出来因为结果是小负数时printf很难察觉但两个大数相减时就会错得离谱。建议写完以后专门造两组数比如a 10000b 1手动心算一下再对拍程序输出。5.3 最后一位总是差1别慌别忘了级数是无限项而我们只算到数组能表示的精度就截断了。这就意味着最后一位的舍入误差几乎不可能避免。解决办法很简单程序里PREC DIGITS / 4 4也就是多算16位然后用round_array对前DIGITS位做四舍五入。这个处理不是可有可无如果你直接截断1000位的结果里第1000位经常比标准值小1这是我在调试时反复确认过的。5.4 数组越界是越改越乱固定数组的边界问题非常阴险。特别是arctan里用到了temp[PREC]、term[PREC]result数组在main里是[PREC]{0}这几个数组长度都是一样的。但是如果把PREC定义成DIGITS/4而PREC本身要包括a[0]整数位那么小数位组数就少了一组输出最后几位会乱。我的建议是所有用PREC做循环边界的地方都要写成i PREC而不是i PREC。多一个等号数组最后一位就溢出了。5.5 想提速的话这几个方向值得试如果你以后要算更多位可以从三个方向优化。第一把万进制改成更大进制比如用1e9进制每一位存9位十进制数组长度进一步缩短但这时候remBASE可能会逼近int上限需要改用long long来存cur。第二在arctan里每轮都要执行div_int(temp, xx)而x*x是固定不变的可以考虑把这个除数作为参数传进去避免每次重复计算乘法。第三输出阶段用fwrite代替printf的逐组格式化大批量输出时能省不少时间。真正到万位以上还有更多基于FFT的快速算法但那已经超出这道练习题的讨论范围了。最后再说点实在的我个人写这个程序最大的体会是高精度计算并不神秘它就是把小学竖式重新翻译成了C语言的数组操作。如果你在学数组和指针时总觉得“会背概念但不会用”这道题值得你亲手做一遍尤其是先拿纸笔手动算一遍竖式除法再回来写div_int你会发现原本容易写错的方向问题一下子变得特别自然。另外也分享一个我自己的习惯写这种多函数程序时每写完一个函数立刻单独测试不要等所有函数都写完再整体调试。比如先写add和sub随便填两组数看结果再写div_int算一个1/3看输出是不是0.3333。每一步都验证最后组合起来时出错概率会低很多。把这份代码跑通后你手里相当于攒了一套最基础的高精度小工具以后再做任何大整数运算、高精度开方、算法竞赛题都能直接拿来打底。