翻日历的时候发现2026年3月14日正好是周六这个日期放在其他年份只是三月的普通一天但在玩数学的人眼里0314几个数字一出现大脑就会自动接上3.1415926535。于是第11期《数学周刊》从3月9日到3月15日这一周想绕开π都绕不开。这一期我干脆做成一个π特辑不打算从定义重讲π而是把社区里这周最热闹的玩法、算π的历史和算法、一道能想一整天的证明题以及在业余条件下怎么亲手复现一遍π按自己的经验整理一遍。看完这篇文章你可以直接抄代码去跑也可以只挑感兴趣的章节读不会影响理解。1. 本周社区在折腾什么π Day把围观者分成了几拨1.1 跑分党把计算π当硬件压力测试圆周率日前后最热闹的一拨人往往不是数学系出身而是跑分党。原因很简单π的位数计算是高度重复的浮点运算对CPU指令集、内存带宽、存储速度都特别敏感而且结果不容易造假——两个程序跑出来几十亿位以后只要失配一位就能立刻发现问题。所以你会看到一些人互相较劲用y-cruncher这类工具把π算到几亿位、几十亿位顺便观察机器的温度曲线。这个玩法我不打算劝退毕竟我自己也干过。但要提醒一句为了压测而算π和为了体验算法而算π是两码事。用普通笔记本直接挑战十亿位以上内存和散热大概率先顶不住。想在小机子上体验跑分乐趣把目标定在一亿位以内记录一下耗时感受算法和硬件在I/O上的拉扯就够了。真正冲纪录的人用的是服务器级内存阵列和专门调优的存储那已经属于发烧友赛道不适合本文展开。1.2 教学派把3.14变成动手日另一拨是老师。不少中学数学老师把3月14日设计成“测量日”让学生用绳子、圆规和A4纸分别量出圆周率。这类活动比讲一节课有效得多因为学生很快会意识到手工量出来的数据总在3.1到3.2之间晃想要更精确就必须换思路。真正好玩的地方不是记住3.14而是回答一个问题为什么所有圆的周长和直径之比都恰好相同这个问题听起来简单实际上已经把几何学逼到了极限。它是从“特殊圆”到“所有圆”的跳跃也是后来微积分登场的一个推动力。我见过不少学生从测量日之后开始对数学产生兴趣原因不是他们记住了π的值而是他们第一次发现自己动手得到的数字和传说中的常数对不上从而产生了“到底为什么”的好奇心。1.3 数字侦探在π里找生日和电话还有一拨比较文艺在π的前几亿位里找自己的生日、毕业年份或者某句号码。因为π看起来分布均匀甚至有人直接把它当随机数源用。这里是本周最容易产生误解的地方π不是随机产生的它是确定性计算出来的无限不循环小数只是目前没有证据证明它是一个正规数。也就是说虽然前几十亿位的数字看起来均匀但在数学上我们并不知道每一个有限数字序列最终会不会以理想频率出现。沿着这个思路往深里走其实是“伪随机性”问题在作祟。计算机生成的所有随机数也都是确定性的但这种确定性和统计意义上的随机性可以共存。π只是把这件事推到了极致一个由精确公式产生的数通过了大量统计随机性检验却仍然无法被证明满足最自然的均匀分布定义。这个反差很有嚼头适合在周日夜里慢慢想。2. 重新认识π从内接多边形到十万亿位2.1 几何时代的暴力与巧劲在无穷级数出现之前算π主要靠几何逼近。阿基米德用内外切正96边形夹出了一个范围这个办法直观但笨重每加一次边运算量就涨一截。同一思路在刘徽、祖冲之手里继续精进祖冲之给出的约率22/7和密率355/113尤其值得一说。22/7只精确到小数点后两位但355/113能精确到小数点后六位误差只有约两千万分之一。这个分数在欧洲直到16世纪才重新出现中间隔了近千年。这里有一个容易被忽略的点355/113并不是靠数万名工人手工量出来的而是从不等式反推出来的有理逼近。用现代语言说这是一个“已知目标精度反求分母不太大的最佳分数”的问题。这种思想比分数本身更值钱因为它讲的不是“怎么算得更准”而是“在不会算很多位数的情况下怎么用一个小巧的分数把事情办妥”。2.2 级数把π从几何题变成了加法题真正的转折点是无穷级数。Machin公式把π/4表示成两个反正切函数的线性组合只需要把无穷级数一项一项加下去就能稳定推进小数位。1706年前后人们用这类公式把π算到了一百位以上。这个时代的标志性变化是计算π不再依赖画图或几何不等式而是依赖“算得快的好公式”。这背后是一个值得记住的道理同一件事换个表示方式难度可以天差地别。用几何逼近要提心吊胆地处理误差界用级数求和却可以按需多取一项误差结构清楚得多。数学史上很多突破本质上都是“换表示形式”从算π到解方程、展开函数全都在干同一件事。2.3 计算机时代与一个漂亮到不合理的公式进入计算机时代后计算π的里程碑变成1949年ENIAC用七十小时算出两千多位。此后纪录一路飙升真正让我觉得“算法比机器更值钱”的是1995年的BBP公式def bbp_pi(terms20): total 0.0 power 1.0 for k in range(terms): total (4 / (8*k 1) - 2 / (8*k 4) - 1 / (8*k 5) - 1 / (8*k 6)) * power power / 16.0 return total这段代码只有几行跑完二十项就让你看到3.141592653589793——因为浮点数精度到这里就饱和了。BBP公式真正惊人的地方在于它允许跳过前面的位直接提取π在十六进制下的某一段数字不需要先算前面所有位。这违背了大多数人对“计算后一位必须依赖前一位”的直觉也是当年发表后立刻引起轰动的原因。现代计算纪录用的是另一个武器Chudnovsky级数。它的完整形式写成这样1/π 12/640320^(3/2) * Σ (-1)^k (6k)! (13591409545140134k) / ((3k)! (k!)^3 640320^(3k))每一个新项大约多贡献14位十进制精度。同样跑二十项Chudnovsky已经远超双精度能表示的范围所以实际的大规模计算必须配合高精度整数库和极其聪明的存储策略。看到世界纪录从几千位涨到十万亿位很多人以为是电脑变快了其实算法换代贡献了同样大的力量。我把几种算法放进一张表里方便对比算法时代收敛特点适合场景多边形逼近古典每多一位几乎要重复一遍几何计算教学演示莱布尼茨级数17世纪极慢误差约1/n理解级数收敛Machin公式18世纪中等速度手工可算几百位历史编程练习BBP公式1995年可跳到十六进制任意位附近理论趣味与分布式计算Chudnovsky级数现代每项约14位十进制世界纪录与高精度计算3. 我写的第一段算π代码蒙特卡洛方法及它的收敛性格3.1 原理与代码蒙特卡洛方法大概是程序员最先接触的算π方式。在边长为1的正方形里画一个半径0.5的内切圆随机丢点落在圆内的比例应该趋近于圆面积与正方形面积之比也就是π/4。把比例乘4就得到π的估计值。我常用的写法是直接用NumPy批量生成随机点import numpy as np def monte_carlo_pi(n): points np.random.rand(n, 2) inside (points[:, 0] ** 2 points[:, 1] ** 2) 1.0 return 4.0 * inside.sum() / n print(monte_carlo_pi(1_000_000))一次跑一百万点结果往往在3.139和3.146之间晃不会稳定落在3.14159附近。第一次跑出来的读者不用怀疑代码写错了蒙特卡洛就是这个性格。3.2 为什么误差只按根号下降这里的关键是误差的数学模型。每次随机点相当于一次独立试验落在圆内的概率是π/4估计量是4乘以样本比例。根据二项分布估计值的标准差大约等于1.64除以根号下的样本量N。想要把误差缩小一十进制位样本量要扩大约一百倍这就是蒙特卡洛方法的典型代价简单、可并行但收敛慢。我用表格记录过理论误差样本量 N标准差大约直观效果10^40.016只能保证两位小数10^60.0016有时能在第三位上碰对10^80.00016可以稳定到三位小数10^100.000016勉强摸到四位小数所以蒙特卡洛算π在数字精度上根本打不过级数法。它的价值不在这里而在于一个方法论示范当问题复杂到无法直接积分比如高维积分、期权定价、物理模拟这套“随机抽样估计期望”的思路依然成立。π只是那颗最好用的教学弹珠。3.3 实操时最容易踩的两个坑第一是随机数生成器质量。现代numpy的默认生成器已经足够但如果你在某个老环境里遇到古老的线性同余生成器点分布可能出现可察觉的网格结构估算结果会系统偏移。判断方法很简单把生成的二维点画成散点图肉眼看有没有规律排列。第二是一次性生成大数组的内存问题。N取十亿时np.random.rand(n, 2)会占用大量内存普通机器很可能直接卡死。更稳妥的做法是分块采样比如每次生成十万个点累加命中次数最后统一除以总样本数。这和训练神经网络时用mini-batch的道理一样都是为了在有限内存里跑更大规模的计算。4. 一道与π有关的证明题从二重积分到巴塞尔问题4.1 题目与第一直觉这周给读者留的作业是一道证明题证明∫₀¹ ∫₀¹ 1/(1-xy) dxdy π²/6我第一次看到这个式子时第一反应是“这跟π有什么关系”。等把积分区域画出来才发现难点在于分母里的xy把x和y耦合在一起不能直接拆成两个独立积分。真正的突破口来自一个中学级的知识当|xy|1时1/(1-xy)可以展开成等比级数。因为积分区域是[0,1]×[0,1]xy始终小于1展开完全合理。4.2 用级数交换积分与求和把分母展开后原积分变成对(xy)^n从n0到无穷求和再积分的结构。由于被积函数非负这里可以放心交换求和与积分的顺序∫₀¹ ∫₀¹ 1/(1-xy) dxdy Σ_{n0}^∞ ∫₀¹ ∫₀¹ (xy)^n dxdy Σ_{n0}^∞ (∫₀¹ x^n dx)(∫₀¹ y^n dy) Σ_{n0}^∞ 1/(n1)² Σ_{k1}^∞ 1/k²到这一步原积分被完整转化成了巴塞尔级数。剩下的核心结论是欧拉在1735年前后证明的著名结果所有正整数平方的倒数之和等于π²/6。于是题目得证。4.3 欧拉当年的跳跃正弦函数的无穷乘积如果读者问我“巴塞尔级数为什么等于π²/6”最直观的现代证明来自正弦函数的无穷乘积展开sin x / x ∏_{n1}^∞ (1 - x²/(n²π²))把左边按泰勒级数展开x²项的系数是-1/6把右边展开x²项的系数是所有-1/(n²π²)的总和。让这两个系数相等立刻得到Σ1/n²π²/6。欧拉当年用多项式类比推出这个结论时无穷乘积的严格性还没有完全建立但结果正确到令人惊叹。后来数学分析逐渐成熟这一步步被补全成了经典的严整证明。这里的启示很实际很多困难问题不是直接算出来的而是换了一个表示形式让隐藏结构自己跳出来。正弦函数和正方形上的积分表面上看毫无关系但它们在π²/6这个数字上汇合说明数学里的“表面无关联”常常只是视角错位。4.4 顺手用蒙特卡洛验证既然上一节写了蒙特卡洛干脆用同样的随机方法来验证这道积分。把被积函数1/(1-xy)当作用随机点计算均值的对象import numpy as np samples np.random.rand(100_000, 2) values 1.0 / (1.0 - samples[:, 0] * samples[:, 1]) print(values.mean()) # 大约 1.6449真实值π²/6约等于1.64493406685。十万次随机采样通常能给你1.64到1.65之间的结果误差在两三位小数上。用同一套随机工具既算了π又验证了一道证明题这大概就是“常数之间的暗号”最好玩的注脚。4.5 概率视角的趣味补充如果换成概率语言还能把这个积分说成一句话设X和Y是[0,1]上的独立均匀随机变量那么E[1/(1-XY)]π²/6。一个涉及两个均匀变量的期望居然落在一个几何常数上这种反直觉是数学常态。它也让这个题目在讲解时多了很多入口——你可以从重积分讲从级数讲从概率讲结果是同一个数字。5. 本周菜单按自己的目标接着往下钻5.1 想读点故事推荐《π的历史》如果对算π的历史感兴趣Petr Beckmann写的《A History of Pi》很值得翻。这本书不堆公式重点在讲清每个时代的计算者是怎么在有限条件里想办法的。看完它你会明白今天跑一次代码就能得到百万位是靠几百年来无数次“换表示形式”才做到的。它虽然成书较早但历史部分没有过时。5.2 想无脑算任意精度mpmath一行搞定在Python里想快速得到高精度的π不必自己实现Chudnovsky级数mpmath库已经封装好了from mpmath import mp mp.dps 50 print(mp.pi)把mp.dps改大一些就能输出更多小数位。mpmath内部用的正是现代快速算法适合做验证和教学不建议用它冲击世界纪录因为纯Python在高精度运算上的速度远不如C语言优化过的程序。5.3 想挑战极限y-cruncher的注意事项我自己在折腾y-cruncher时踩过不少坑最大的教训是别一开始就贪大。y-cruncher几乎是最流行的π计算工具常被用来压测硬件但它对内存和临时存储的要求随位数暴涨。普通笔记本先跑一千万到一亿位热热身再逐步加量比直接挑战万亿位明智得多。真正冲纪录的人会花几天时间准备内存盘和磁盘阵列这个门槛已经不只是数学问题更是硬件工程问题。5.4 一句经验总结把算法、硬件和题目都过一遍之后我最深的体会是算法之间的差距往往比机器之间的差距更悬殊。蒙特卡洛跑一亿次采样只能勉强逼近三五位小数Chudnovsky级数跑几十项已经超出浮点数表达能力前者适合培养直觉后者才适合书写纪录。下周第12期会恢复常规选题但π这笔账记下之后再看到某个积分里蹦出3.14159你就不会觉得奇怪了。最后分享一个小技巧想快速心算π的近似值时355/113比22/7更值得记它精确到小数点后六位日常生活和工程设计里基本够用。3月14日就算过完了这个分数可以一直带在身边。