1. 先别急着敲代码这道题到底在考什么格里高利公式求 π 这道题几乎每个学 C 语言的人都会在在线判题系统里碰上一次也几乎每个人都会在循环边界上翻一次车。它的公式本体并不复杂π/4 1 − 1/3 1/5 − 1/7 1/9 − …一项正、一项负分母是连续的奇数把整个式子乘 4 就得到 π 的近似值。题目要求输入一个精度 eps程序不断累加直到某一项的绝对值小于 eps 为止最后输出累加结果乘 4 得到的 π 值。看起来就是求和 循环 判断可真正提交上去报错信息往往集中在几个很固定的位置。这道题的本质是把数学上的交替级数翻译成一段 C 语言的循环。它同时考了三件事浮点数的类型与精度、循环终止条件的语义、以及输出格式的严格匹配。三者任何一环出问题判题系统都会直接给个红叉而且反馈往往只有答案错误三个字不告诉你错在哪。所以我写这篇文章不是给你一份能 AC 的代码而是把这道题从数学推导、变量选型、循环结构、踩坑排查到进阶优化完整拆一遍让你下次遇到任何按精度截断的级数求和题目都能自己推出来。适合的读者范围很宽刚学完 for 循环和 if 判断的大一同学能拿到一份可以照着敲、照着调的实现已经写过几道判题题但总在浮点输出上失手的人能在第 4 章的问题速查表里找到对应症状对数值计算感兴趣的人第 5 章关于收敛速度和求和误差的分析应该会有点意思。全程用 C 语言标准库只用 stdio.h 和 math.h不依赖任何平台特性本地编译器一样能跑。2. 数学底子格里高利公式到底在算什么2.1 从反正切的泰勒展开推出来格里高利公式不是天上掉下来的它是反正切函数在 0 点泰勒展开的一个特例。arctan(x) x − x³/3 x⁵/5 − x⁷/7 …这个展开在 |x| ≤ 1 时成立x 1 时左边就是 arctan(1) π/4右边变成 1 − 1/3 1/5 − 1/7 …于是得到公式本身。明白这个来源有个直接好处你会知道这个级数在 x 1 这一点上收敛得极其勉强。泰勒展开的收敛半径是 1x 1 恰好落在收敛区间的端点上属于条件收敛——级数确实收敛到 π/4但它不是绝对收敛去掉绝对值之后 Σ1/(2n1) 是发散的。这一点直接决定了后面所有的性能问题收敛速度慢慢到你会怀疑电脑卡住了。生活化的类比想象你要从 1 走到 0.785每次迈的步子正负交替步子越来越小。前几步跨度巨大1.0、0.667、0.867、0.723……你能迅速逼近目标附近但一旦靠近剩下的就是无数个小碎步在目标两侧来回晃越晃越近却永远只是更近一点。这就是为什么 eps 从 1e-4 降到 1e-6项数会从五千项暴涨到五十万项。2.2 交替级数的误差界为什么最后一项小于 eps就够了教材上一般会讲莱布尼茨判别法对于各项绝对值单调递减趋于零的交替级数其前 n 项和与真实值之间的误差绝对值不超过第 n1 项的绝对值。这条定理是这道题能够用最后一项的大小作为停止条件的理论依据也是整个解法成立的根基。把它翻成大白话你算到第 n 项停手你欠下的那份误差最多也就和下一项差不多大。所以只要保证下一项小于 eps误差就一定小于 eps。这个界非常干净不需要做任何复杂的误差估计这也是出题人选择格里高利公式的原因——它让给定精度这个要求变得可判定、可编程。注意这个误差界是截断误差的界不含浮点运算的舍入误差。也就是说算到项小于 eps 就停只能保证数学意义上的误差受控实际程序里还有计算机本身的表示误差两者要分开看。第 5 章会专门聊这个。还有一个细节值得说定理说的是不超过下一项但没说误差一定等于下一项。实际上真值总是夹在相邻两个部分和之间所以误差大约是下一项的一半。这给了我们一个免费的精度提升技巧第 5 章会讲怎么用。2.3 收敛速度的残酷现实先算一笔账写代码之前先估算一下工作量这是我养成的习惯能避免很多程序跑不完的冤枉路。设累加到第 k 项k 从 0 开始计数这一项是 1/(2k1)。题目要求这一项小于 eps 就停也就是 1/(2k1) eps反解得到 2k1 1/eps即 k (1/eps − 1)/2。忽略常数项所需项数大约是 1/(2·eps)。按这个估算列一张表心里就有底了。下表的 π 近似值是按误差约 1/(4N) 估算的参考量级实际运行结果会在这个数附近浮动输入精度 eps需要累加的项数π 的截断误差量级π 近似值参考1e-1约 5 项约 1e-13.1416 附近抖动大1e-2约 50 项约 1e-2约 3.12181e-3约 500 项约 1e-3约 3.13961e-4约 5000 项约 4e-4约 3.14141e-6约 50 万项约 4e-6约 3.14161e-8约 5000 万项约 4e-8约 3.1416从这张表能读出两个关键信息。第一eps 1e-4 时只有几千项任何现代机器都是瞬间完成不用担心超时eps 1e-8 时是五千万项单次跑下来通常要零点几秒到一两秒具体取决于编译优化等级某些判题系统会对运行时间比较敏感。第二当题目要求输出四位小数时eps 取到 1e-4 这一档截断误差约 4e-4已经可能影响到小数点后第三、第四位——这就是很多人本地看着对、提交就错一位的根源不是代码写错是公式本身的精度不够。3. 方案设计循环怎么搭、变量怎么定3.1 通项怎么生成直接算还是递推生成每一项有两种写法各有取舍我两种都写过建议你都试试。第一种是直接计算维护一个分母变量 d初值 1每轮 d 2当前项就是 sign / (double)dsign 每轮取反。写法直白可读性最好缺点是每轮要做一次浮点除法而且 sign 需要额外维护一个变量或者用 i%2 判断。第二种是递推已知 term_k (−1)^k/(2k1)那么 term_{k1} −term_k × (2k1)/(2k3)。只做乘法和除法不需要单独维护符号但引入了当前项 × 比例因子的误差传递项数多了以后累积误差略大一点点实际上在本量级可以忽略。我平时首选直接计算因为它把符号和分母两个概念分得很清楚出错时一眼能看出来是哪半边的问题。递推版适合项数上千万的场景省一次除法但对这道题来说意义不大。3.2 循环终止时机先判断还是先累加这是本题最容易出语义错误的地方必须掰开说。题目的原话是直到最后一项的值小于给定精度为止。这句话有两种理解一是这一项要加进去二是这一项不要加。两种写法都能让误差受控但具体数值会差一项在某些 eps 取值下足以改变第四位小数。我一般采用先判断、后累加算出当前项如果它的绝对值已经小于 eps直接跳出循环不加。这样保证所有被加进去的项都 ≥ eps而没被加的第一项 eps误差严格小于 eps与莱布尼茨界完全吻合。写出来是这样sum 0 term 1.0 while |term| eps: sum term 生成下一项 term对应的另一种是先累加、后判断sum 0 term 1.0 loop: sum term 生成下一项 term if |term| eps: break两种写法的部分和差异就是最后那一项大小不超过 eps。如果 eps 1e-4这个差异乘以 4 之后是 4e-4 量级在四位小数输出时有可能造成末位不同。所以严格来说你得按题目的原话来定。我的经验是绝大多数这类题目的标准答案采用不加最后那一项的写法因为误差界更干净。你在本地两种都跑一下比对判题结果一次就能确定下来。提示如果题目描述写成直到某项小于 eps 为止该项也计入那就用第二种写成直到某项小于 eps不再累加就用第一种。中文题面里这两个表述经常混着用别凭感觉猜。3.3 变量类型double 够不够什么时候要 long double这道题的标准答案用 double 完全够。double 在主流平台上有效十进制位约 15~16 位累加几千项到几百万项累加误差远小于截断误差不影响结果。我做过对比测试同样的 eps 1e-6double 和 long double 算出来的四位小数结果完全一致。long double 的位置有点尴尬。它在 x86 上是 80 位扩展精度有效位约 18~19 位但在某些平台上和 double 一模一样跨平台行为不一致。而且用 long double 就要把 scanf 的格式符从 %lf 改成 %Lfprintf 改成 %Lf属于纸上多写几笔收益极小的操作。除非题目明确要求 eps 小到 1e-12 这种量级否则我不用。真正需要关注的是计数器类型。如果你用 int 做分母eps 1e-9 时 2k1 会跑到 2e9 以上逼近 int 上限 2147483647一不小心就溢出变成负数然后分母为 0程序直接崩。用 long long 保平安代价只是几个字节没有任何理由省。4. 核心实现从最小可用版到工程加强版4.1 最小可用版本与逐行解释先给一份能直接提交的版本代码不长但每一行都有说法#include stdio.h #include math.h int main(void) { double eps; if (scanf(%lf, eps) ! 1) { return 0; } double sum 0.0; long long k 0; /* 第几项从 0 开始 */ double term 1.0; /* 当前项第一项是 1/1 */ while (fabs(term) eps) { sum term; k; /* 第 k 项的分母是 2k1符号正负交替 */ double sign (k % 2 1) ? -1.0 : 1.0; term sign / (double)(2 * k 1); } printf(pi %.4f\n, 4.0 * sum); return 0; }逐行看几个要点。scanf(%lf, eps)里的 l 不能省。double 必须用 %lf 读入float 才用 %f。虽然 printf 里 double 用 %f 就行因为可变参数默认提升但 scanf 必须严格区分。这是新手常犯的第一类错误。fabs(term) eps里的 fabs 来自 math.h作用是取绝对值。不能写成 term eps因为负项的绝对值才是我们关心的量。也不能写成 fabs(term) eps浮点数永远不要用等号比较。sign / (double)(2 * k 1)这个写法看起来啰嗦其实是为了防止整数除法。如果写成sign / (2 * k 1)由于 sign 已经是 double表达式会自动转成浮点除法结果其实是对的。但我还是习惯显式转换因为一旦哪天把 sign 改成 int 类型的 ±1这一行立刻变成整数除法结果全是 0而且编译器不报错——这种错误最难查。printf(pi %.4f\n, 4.0 * sum)里两个细节。一是 4.0 写成浮点字面量虽然 4 * sum 因为 sum 是 double 也会自动提升但显式写 4.0 更保险意图也更清楚。二是换行符\n判题系统对行尾是否换行通常宽容但加上更安全某些严格的比较器会逐字符比对。4.2 参数计算与实测记录拿 eps 1e-4 走一遍k 0 时 term 1.01.0 ≥ 1e-4累加k 变成 1term 更新为 −1/3 ≈ −0.3333继续直到 k 5000 时 term sign/(2×50001) −1/10001 ≈ −9.999e-5此时 fabs(term) 9.999e-5 1e-4循环退出。总共累加了 k 0 到 4999 共 5000 项与第 2.3 节的估算完全吻合。我在本地实测了不同 eps 下的项数和耗时用 gcc -O2 编译结果如下eps实际累加项数累加和 sum4×sum输出四位小数1e-150.83493.33973.33971e-2500.78043.12163.12161e-35000.78493.13963.13961e-450000.78533.14143.14141e-5500000.785393.141573.14161e-65000000.7853973.141593.1416这张表里有个值得注意的现象eps 1e-4 输出的 3.1414 和真值 3.1416 差了两位。这不是 bug是 5000 项截断必然的结果。如果题目给定的测试用例恰好是 eps 1e-4 并且要求输出 3.1416那只能说明题目期望的是累加到下一项小于 eps 才停或者eps 判据用更严格的方式此时需要调整循环边界。所以我建议大家本地测的时候把 eps 从 1e-1 一直扫到 1e-6把六组输出都记下来再去比对判题反馈能快速定位语义分歧。4.3 加强版输入保护、Kahan 求和与倒序累加实际跑 eps 很小的场景时普通版本会露出两个尾巴一是 eps ≤ 0 会死循环二是项数过亿后累加误差开始抬头。下面这版加了保护#include stdio.h #include math.h int main(void) { double eps; if (scanf(%lf, eps) ! 1) { return 0; } if (eps 0.0 || eps 1.0) { /* 非法精度或过大精度直接给个合理兜底 */ printf(pi %.4f\n, 4.0 * (eps 1.0 ? 1.0 : 0.0)); return 0; } double sum 0.0; double comp 0.0; /* Kahan 补偿项 */ long long k 0; double term 1.0; while (fabs(term) eps) { /* Kahan 补偿求和 */ double y term - comp; double t sum y; comp (t - sum) - y; sum t; k; double sign (k % 2 1) ? -1.0 : 1.0; term sign / (double)(2 * k 1); } printf(pi %.4f\n, 4.0 * sum); return 0; }Kahan 补偿求和的核心思路是每次加法都会丢掉一点低位信息那就把丢掉的那部分单独存起来下一次加的时候补回去。代码里 comp 就是这个被丢掉的部分。它在项数少的时候毫无作用在项数上千万的时候能明显压住误差增长。关于 eps 的边界题目一般保证 0 eps 1但真实判题数据偶尔会有 1e-9 甚至更小的值。eps 太小会导致项数爆炸五分钟都跑不完这是公式本身的极限不是你能优化的。遇到这种输入比较现实的选择是确认判题系统的时限或者改用第 6 章的加速方案。5. 判题系统上的高频翻车点与排查表5.1 整数除法这道题的头号杀手写了这么多年代码我见过最多的错误类型就是这个。当你写下term 1 / (2 * k 1)时编译器会老老实实地做整数除法1/3 得 01/5 得 0所有项全是 0sum 永远是 0输出 0.0000。更气人的是如果第一项你写成 1.0后面全 0输出会变成 4.0000比 0 还误导人。修法只有一条确保分子或分母有一个是浮点。写成1.0 / (2 * k 1)或者(double)1 / (2 * k 1)或者把分母强转成 double。我个人的习惯是把常量写成浮点字面量1.0比(double)1短也更容易一眼看出意图。顺带说一个变体错误term sign / (2 * k 1)里如果 sign 是 int 类型的 ±1同样会退化成整数除法。这类错误编译器不会报警需要你自己在写完循环后把第一两项手算一遍和程序打印的中间值对一对。我在调试期通常会在循环里加一句if (k 5) printf(debug: k%lld term%.10f sum%.10f\n, k, term, sum);跑完马上删掉这个习惯救过我无数次。5.2 输出格式差一个空格就是答案错误判题系统的输出比对通常是很严格的pi 3.1416和pi3.1416和pi 3.1416末尾多空格在机器眼里是三个不同的字符串。题目要求的格式必须逐字符照抄包括等号两边的空格、小写 pi、保留几位小数。我整理了一张格式对照表方便你按题目要求直接套题目要求printf 写法pi 3.1416printf(pi %.4f\n, v);pi3.1416printf(pi%.4f\n, v);PI 3.141593六位printf(PI %.6f\n, v);3.1416只输出数值printf(%.4f\n, v);关于保留位数还有一点%.4f采用四舍五入更准确地说是就近舍入、平局取偶而有些人用(int)(v * 10000 0.5) / 10000.0手动处理这种做法在负数时会翻车也没必要。直接用格式符最稳。5.3 死循环与边界输入三个常见的死循环来源。第一eps 输入为 0 或负数fabs(term) eps恒为真循环永不退出。第二符号忘记翻转term 变成常数或者绝对值不递减循环条件永远成立。第三有时判断写成了term eps而不是fabs(term) eps当 term 变成负数后条件立刻为假循环提前退出输出一个明显偏小的值——这个错误不会死循环但结果全错更隐蔽。排查死循环的土办法在循环体里加一个计数器超过一百万次就强制 break 并打印当前状态一眼就能看出是条件写错了还是 eps 本身有问题。long long guard 0; while (fabs(term) eps) { /* ... */ if (guard 100000000LL) { fprintf(stderr, loop guard triggered, k%lld term%g\n, k, term); break; } }5.4 常见问题速查表把上面散落的症状集中成一张表出问题时按症状反查比一行行读代码快得多症状最可能的根因修法输出 pi 0.0000整数除法所有项为 0分子写成 1.0输出 pi 4.0000只有第一项参与累加检查 term 的更新是否被短路输出约为真值 1/4忘了最后乘 4输出时写4.0 * sum输出为负的 π第一项符号写反确认首项是 1程序卡住不返回eps ≤ 0或符号不翻转加 eps 校验检查 sign 更新前两位对、末位差 1截断误差或停止条件差一项换用另一套边界写法试一次大 eps 时结果离谱eps 1程序按规则提前停属于题目数据问题不用改代码提交是格式错误而非答案错误空格、大小写、换行不符逐字符对照题面经验遇到格式错误时最快的定位办法是把题面要求的输出字符串原样复制到 printf 的格式串里只把数值部分换成 %.4f不要自己凭记忆敲。6. 想算得更快更准几条实用的优化路线6.1 成对相减把交替级数变成正项级数既然交替求和容易产生相消那就把相邻两项先配对(1 − 1/3) (1/5 − 1/7) …每一对化简成 2/((4j1)(4j3))全是正数累加过程中不再有正负抵消。改写后的代码double sum 0.0; long long j 0; double pair; do { double d1 (double)(4 * j 1); double d2 (double)(4 * j 3); pair 2.0 / (d1 * d2); sum pair; j; } while (pair eps); printf(pi %.4f\n, 4.0 * sum);这个写法有两个好处。一是累加全为正项不存在大数吃小数的相消损失数值上更稳。二是每轮处理两项循环次数减半代码里少一个符号变量逻辑负担更轻。缺点是停止条件变得不那么直观——pair 对应的是两小项之和量级约等于原级数两项的合并误差与题面某一项小于 eps的字面要求略有出入。所以这个版本适合你自己练习和对答案正式提交时还是要按题目字面语义来。6.2 平均值修正免费的额外一位精度前面提过交替级数的真值总是夹在相邻两个部分和之间。也就是说算到 S_n 时真值在 S_n 和 S_{n1} 之间取二者的平均 (S_n S_{n1})/2 会显著减小误差通常能把误差从 O(1/n) 压到 O(1/n²) 的量级。实现上很简单在循环退出时你手里已经有 sum对应 S_n和 term就是下一项 a_{n1}但还没加。于是可以写double better sum term / 2.0;这个 better 的精度一般比 sum 高一个数量级。注意这个技巧在给定精度求值的题目里不能随便用因为判题系统比的是标准程序的结果你用了修正反而可能对不上。它更适合你自己做数值实验、验证公式正确性的时候用。6.3 换个公式拿马青公式做交叉验证格里高利公式最大的问题是收敛太慢慢到 eps 1e-8 就要五千万项。如果你只是想验证自己的循环逻辑对不对可以临时换一个收敛快的公式来对照比如马青类公式π/4 4·arctan(1/5) − arctan(1/239)其中 arctan(x) 用同样的泰勒级数展开但 x 很小每加一项数值掉得飞快几十项就能到十几位精度。我常拿它当参考答案格里高利版算出来的值和马青版算出来的值对比如果差在 eps 的量级说明循环逻辑没问题如果差得离谱那就是代码错了不是精度问题。这种交叉验证的思路在做任何数值计算题时都值得用。/* 仅供参考的交叉验证用 arctan 小参数展开算 pi */ static double arctan_small(double x) { double sum 0.0, term x, x2 x * x; long long n 0; while (fabs(term) 1e-18) { sum term; n; term * -x2; term / (double)(2 * n 1) / (double)(2 * n - 1); } return sum; }顺带说一句这种递推写法里那个除法因子很容易写错我第一版就写反了结果算出来是 0.785 的两倍。所以写完一定要拿已知值验证arctan(1) 应该接近 0.785398。6.4 Kahan 与倒序累加什么时候才真的需要前面给的加强版里用了 Kahan但我要老实说一句对这道题Kahan 基本是杀鸡用牛刀。算一下就知道eps 1e-4 时五千项双精度机器 epsilon 约 2.2e-16部分和约 0.785单次加法的舍入误差上限约 1.7e-16五千项线性累加最坏也就 8.6e-13和 4e-4 的截断误差差了九个数量级。就算 eps 降到 1e-8五千万项累加的最坏误差约 8.6e-9和 4e-8 的截断误差仍有一个数量级的余量。真正需要 Kahan 或者倒序累加的场合是 eps 小到 1e-12 以下、项数上十亿级别的场景或者你要求结果精确到最后一位有效数字。倒序累加先数项数再从最小项往回加能彻底消除大数吃小数但需要两遍扫描代码量翻倍对这道题来说是纯粹的过度设计。我的建议是先把最朴素的版本写对、提交通过再谈优化。很多人一上来就琢磨高精度求和结果连整数除法都没改过来本末倒置。7. 我在这类题目上踩过的几个坑第一次做这道题的时候我把停止条件写成了term eps因为当时脑子里想的是项小于精度就停完全没意识到 term 变成负数后条件立刻为真。程序输出 4.0000我盯着屏幕看了十几分钟直到手动打印了前五项才反应过来——第一项 1.0 加上去第二项 −0.3333 一判断−0.3333 0.0001 成立直接退出。这个错我现在都还记得所以每次看到有人问为什么输出 4我都先让他检查有没有 fabs。第二个坑更隐蔽我用int存分母eps 取 1e-9 时2k1 一路加到 2e9 以上溢出成负数然后除出来一个巨大的数循环条件瞬间失效程序吐出个荒唐的结果。那次之后我把所有和项数、分母相关的计数器全换成 long long再没出过问题。这种溢出错误编译器不会提醒运行时也不报错只给你一个静默的错误答案属于最难查的一类。第三个坑是输出格式。题目要求pi 3.1416我写成pi %.4f但忘了等号左边的空格连续三次格式错误。后来学乖了直接把题面的输出示例整段复制进格式串只改数字部分。这个习惯我保持到现在凡是涉及固定格式输出的题一律不手打。最后分享一个我自己整理的小流程遇到这类按精度求值的题按顺序走一遍基本不会翻车先手推三项确认公式和符号再写最小循环只打印中间结果不判题然后用 eps 1e-1 到 1e-3 三组数据验证趋势是否合理值应该逐步逼近 π接着对照题面逐字符核对输出格式最后才提交。这套流程看起来多花五分钟实际上省下的调试时间远不止五分钟。