计算π数位的涓流算法

作者：Stanley Rabinowitz 和 Stan Wagon

令人惊叹的是，表 1 所示的这个算法完全不使用浮点运算，却能算出的各位数字。算法开始时，在以其上方所示分数为列标题的各列中，都填入若干个 2。接着将每一项乘以 10；然后从最右边一列开始，把各项对分母取模（其中该列的标题为分子/分母），从而得到商与余数。余数留在原位置不动，而把进位到左边相邻的一列。这一"取模—进位"的操作一直持续到最左列为止。最左列结果的十位数字，就是的下一位数。此后过程不断重复：余数乘以 10，对各个分母取模，并累加进位。

表 1. 一个生成数位的算法的演算过程。虚线标出了关键步骤：从右边开始，将各格的数对该列列头的分母（依次是 25、23、21……）取模，再把所得的商乘上分子（依次是 12、11、10……），进位到左边一列。例如，在那一列中，20 除以 19 得到余数 1，并向左进位。经过最左列的进位之后，依次得到的十位数字是 3、1、4、1。若想获得更多的数位，就必须从更长的一串 2 开始。

译自《美国数学月刊》（The American Mathematical Monthly）第 102 卷（1995 年），第 195–203 页

这是一个"涓流"（spigot）算法：它每次泵出一位数字，而且数字一经算出便不再使用。此外，这些数字在生成过程中完全不需要对浮点实数做高精度（或低精度）运算；整个算法只对很小的整数做普通的整数运算。例如，要得到的前 5000 位数字，只需要对小于 600,000,000 的整数进行运算即可。尽管高精度浮点例程本身也是由整数运算搭建起来的，但本文的算法相当简单，并不去模拟浮点计算。

为了引出的算法，我们先讨论简单得多的的情形——Sale 为发明了这样一个涓流算法［Sale］。他的算法是本文第 1 节讨论的基础。

• e的数字在此数制中呈循环（周期）性质的计数系统。

一个实数的十进制表示可以理解为一种无穷嵌套的表达式；例如：

如果我们改变"进制序列"——在上面的例子中，这个序列是就能得到一些既有趣又有用的表示形式。例如，采用混合递增的进制，便可得到如下形式，称为混合进制表示（mixed-radix representation，参见［Knu, §4.1］）：

其中（即各位"数字"）是非负整数。若对都有，则称该表示是正规的（regular）。混合进制表示记作。对进制而言，每个正实数都有正规表示；并且只要我们排除以"最大数字"结尾的表示，表示就是唯一的（否则例如）；此后，对一切进制我们都排除此类表示。下面这条引理的证明见附录 1。

引理 1
（a）若，则；特别地，。
（b）使用混合进制基数 b 的表示是唯一的。
（c）的整数部分是，小数部分是。

在这个计数系统中，某些无理数会变成循环（周期）的。例如；这正好等价于把无穷级数重新表述为

在这个系统中，有理数对应于有限的（会终止的）数字序列（见附录 1 引理 2）。

求区间内实数 x 的十进制各位数字，可以这样做：取 x 的整数部分，把 x 的小数部分乘以 10，再取所得结果的整数部分，再把新的余下小数部分乘以 10，如此继续下去。在某些混合进制下，这个过程格外简单。若，那么。后者未必是一个正规表示：有些数字可能过大。但我们可以对各数字"减量"——把它们对取模，其中 i 是进制序列 b 中对应元素的分母。从最右端开始做这些取模运算，并把所得的商逐次向左进位，最终便得到 10x 的正规表示。因此，乘以 10 在算法上是直接明了的。而借助引理 1(c)，对 b-表示取整数部分与小数部分同样容易。

现在我们可以给出求的前位十进制数字的算法了。其正确性证明——即说明个混合进制数字足以得到位十进制数字的误差分析——见附录 1 的引理 3。

算法 e-spigot（e 涓流算法）

1. 初始化：首位数字取 2，并把长度为的数组 A 初始化为
2. 重复次：

乘以 10：将 A 的每个元素乘以 10。

取小数部分：从右端开始，将 A 的第个元素对取模，并把商向左进一位。

输出下一位数字：最后得到的那个商就是的下一位数字。

从含 10 个 1 的数组出发（这相当于 11 个混合进制数字，可保证输出的 9 位数字；表中只显示了前 5 位），该算法的前几步演算显示在表 2 中。

表 2. e 的涓流算法的演算过程（的数字以黑体标出）。列标题为的那一列，其归算是对 i 取模进行的。各行最左边的十进制实数，是把该行各列数字视作混合进制表示时所得的值。由于算法只从 11 个混合进制数字开始，第一个十进制数仅是的一个近似值。

[1]任何为一个假定为正规数（normal number）的数 x 产生数字的算法，都存在一个虽不太可能出现、却可能影响结果的缺陷。若 x 介于 1 与 10 之间，而算法说 x 的前 100 位数字是，比如说，4, 6, 5, 0, 7, …, 3, 9, 9, 9, 9, 9——那么就无法确定最后 6 位数字是否正确。它们实际上是 x 的某个近似值的数字，该近似值与真值之差在以内。你不能简单地往下多算几步、直到出现一个非 9 为止，因为内存分配必须事先确定。使用者必须意识到：一串以 9 结尾的数字是一个“危险信号”（red flag），意味着这些位数字可疑；而即使完全没有 9，最后一位数字也有可能是错的。实践中，人们可以（比如说）多要 6 位数字，从而把出这种问题的概率降到百万分之一。

• 数位的涓流算法。第 1节的思路可以导出一个针对的涓流算法，但与相比，还多出一些复杂性，也引出一些额外的有趣问题。我们的出发点是下面这个算是颇为知名的级数：

这个级数可以由的 Wallis（沃利斯）乘积导出；另一种做法是对级数施加一种称为"欧拉变换"（Euler's transform）的加速技巧。这些证明连同另外三个证明以及对早期文献的指引，都可以在［Li］中找到。我们用表示奇数的连乘积；于是该级数等价于，

可以展开为，

最后这个表达式引出了混合进制基数：其各列为，即第位的基数为，相对于它，干脆就是。在进制的正规表示中，第位的数字必须落在区间内。遗憾的是，基数不如基数那样方便好用。

引理 4（证明见附录 1）。进制中各位取最大值的那个数表示 2；因此形如的正规表示落在 0 与 2 之间。引理 4 意味着进制表示不唯一。

例如

由此，

与算法更相关的是：在进制下取整数部分和小数部分并不像进制那样直截了当。数的整数部分要么是、要么是，取决于是落在还是之内。这个麻烦可以这样克服：在下一轮迭代中把的个位数字留在原处，而把的十位数字称为一个"预位数字"（predigit）。预位数字必须暂时扣住不发，因为偶尔（大约每 20 次迭代出现一次）下一个预位数字会是 10；当进位值（介于 0 到 19 之间）大于 10，且同时留下的个位数字是 9（它在"乘以 10"那一步变成 90）时，这种情况就会发生。此时必须先把扣住的那个数加 1，再把它释放出来。算法的具体细节见下文；本文开头的讲述曾回避了这个"偶尔出现 10"的问题。关于个混合进制数字足以给出的位数字的证明，见附录 1（引理 5）。附录 2 给出了该算法的一个 Pascal 实现。

涓流算法

1. 初始化：为一个长度为的数组。
2. 重复次:

乘以 10：将 A 的每个元素乘以 10。

把 A 化为正规形式：从右端开始，将 A 的第个元素（对应基数中的第项）对取模，得到商和余数。余数留在原位，把向左进一位。最后一次被进出的整数（在的那个位置产生）可能大到 19。

取得下一个预位数字：将 A 的最左元素（最大不超过）对 10 取模。所得的商就是的新预位数字，余数留在原位。

调整预位数字：若既不是 9 也不是 10，就把扣住的各预位数字作为 π 的真数字释放输出，并扣住。若是 9，就把加入扣住的预位数字队列。若是 10，则：

把当前预位数字置为 0 并继续扣住；

把所有其他扣住的预位数字加 1（9 变成 0）；

除当前扣住的预位数字外，把其余全部作为的真数字释放输出。

这个算法只使用整数运算，很容易编程实现。论文开头的那张表（表 1）演示了它的运行过程，从
的 13 个混合进制数字起步（足以给出 4 位十进制数字）。为帮助理解算法的运作，请注意：表 1 的第一行（有限行）是 3.1414796… 的一个混合进制表示，第二行表示 31.414796…，第五行表示 1.414796…，第六行是 14.14796…，第九行是 4.14796…，依此类推。表 3 展示了用更大初始数组所做的计算结果；"扣住不发"（holding）的机制要到第 32 位数字时才变得相关。

表 3.的实际数字（最底行）与用 116 个 2 组成的初始数组（足以给出 35 位数字）迭代 35 次所得的最左列进制数字序列的对照。在第 32 次迭代时出现了一个 102，由此产生一个值为 10 的预位数字。

我们再次强调，该算法只使用整数运算。要得到的 5000 位数字，只需要对小于 600,000,000 的整数做算术运算。这个算法自然地引出了这样一个问题：能否把它改进成一个在本质上与 e-spigot 一样简单的算法。

问题. 是否存在一个由有理数构成的基数，使得的表示是周期的、或构成等差数列，并且对而言，恒等于其整数部分？

Gosper［Gos, 第 32 页］发现了一个
的级数，它使我们离"涓流算法的理想境界"仅一步之遥、令人急切的想要得到完美的结果：

他通过使用关于的欧拉变换的一个改进形式，得到了这个级数。

Gosper 的级数引出了基数，相对于它而言，。经过计算可以发现，比起进制情形下出现的 2，这已是实质性的改进。在通常关于数字分布的随机性假设下，进制中出一个"坏预位数字"的概率是 1/20，而在进制中它降到不足 1/110；这是因为进制的预位数字等于 10，只有当余数是 9（乘以 10 后变成 90）且进位是 10 时才会发生。前者的出现概率是 10%，而后者最多 11 次迭代才出现一次，因为进位等于一个介于 0 与 10.93 之间的实数的整数部分。所以说，基数 d 距离"涓流算法的完美境界"只差 1% 以内。由于 Gosper 的级数比我们使用的级数收敛得更快，它对内存的需求也更少：的位数字只需长度为的初始数组；不过，同样大小的数组在进制下调用的整数，会比进制下的更大。

改进 Gosper 级数方案的一条途径，是把中的分数约成最简分数。此时各位取最大值的正规数是，其值等于 1.0000476468…。不难看出，在这个新基数下的正规表示保持不变。然而，约分所花的功夫超过了"预位数字出现 10 的次数减少"所带来的收益。因此，上面那个问题的肯定答案，恐怕理论意义大于实用价值。

的涓流算法绝不能与最近发现的快速算法（博温兄弟、丘德诺夫斯基兄弟等人提出的那些）相提并论——后者已被用来计算的数亿位数字（见［BBB］）。但涓流算法的确有一个优点：完全避开所有浮点运算；因此它很容易在家用电脑上实现，几分钟内就能产出几千位数字。而且它直接以十进制给出结果（其他大多数算法输出的是二进制或某种内部格式，还得再跑一遍换算才能取得十进制数字）。

本文给出的算法可以通过一次输出多位数字来提速。例如，要一次输出五位十进制数字，只需用基数 100,000 来计算的数字即可：在主步骤中把"乘以 10"换成"乘以 100,000"就行。这时整数部分就是十万进制下的下一个"数字"。

如果在十万进制下工作，并且事先知道将要计算的那段数字中不含 "00000" 这个字串，那么算法中调整预位数字的那冗长部分就可以省去。这能带来一个极其短小的计算机程序。例如，Rabinowitz［Rab］就利用这个想法展示过一个仅 14 行的 Fortran 程序，能输出 π 的 1000 位十进制数字。最后我们还提一句：这个算法可以并行化，并行之后它在大约 10,000 位以内会快得惊人。关于其他函数的涓流算法的例子，可参见［Abd］。

附录1.五个引理

引理1

（a）如果；特别地，。

（b）使用混合进制基数的表示是唯一的。

（c）的整数部分是，小数部分是。

证明：

（a）只需证明，而这可由如下事实推出：该级数可以叠缩为

（b）设与表示同一个实数。于是对某个，有，其中且。但这样一来就有，与（a）矛盾。

（c）这由（a）即得。

引理 2.一个正数是有理数，当且仅当它在使用混合进制基数 b时的数字最终全为 0（即表示是有限的）。

证明： 反向（充分性）显然。对于正向（必要性），我们使用一个辅助引理（子引理）。

子引理. 对任意整数和
，只要
，就存在落在区间内的整数
，使得

。

证明： 对 t 作归纳。若
，把写成，其中。由归纳假设，存在一个序列(，它是关于各项等的表示；由此就是关于各项等的表示。

回到引理 2 的证明：设给定一个正有理数。用子引理把表成

再除以，就得到
的一个表示：它是若干阶乘倒数之和，且系数都足够小——这与混合进制基数下的一个有限（终止）表示是同一回事。

引理3.求数位的算法是正确的。

证明：必须证明：的个混合进制数字足以得到的位十进制数字。我们先证明：当时，个混合进制数字就够用于得到位十进制数字。使用个混合进制数字，意味着我们实际算出的其实是下面这个数的十进制数字：。因此我们只需证明（见论文开头的脚注）。对级数尾部作几何级数估计可得，再由 Stirling（斯特林）公式即可推出，

若，则直接计算各位数字便可验证
个混合进制数字已经够用。

引理4.c 进制中各位取最大值的那个数表示 2；因此形如的正规表示落在 0 与 2 之间。

证明：我们不做形式化的证明，而是展示一些 Mathematica 计算如何引导出这个结果（以及一个证明）。用级数来表述，该引理断言

一个粗略的计算显示这个和接近 2。接着对余项——部分和与 2 之差——做精确的有理数运算，得到序列，

这些余项的规律由"后一项除以前一项"揭示出来，得到

即比值通项为。用归纳法可证明这个规律普遍成立；由此余项具有闭式

它趋于 0，结论得证。

引理 5. 求数位的算法是正确的。

证明： 与的情形一样，我们考察，其中。这一误差正是我们主级数的尾部。由于其后每一项都小于前一项的一半，整个尾部小于其首项的两倍，于是问题化为考察。把分母拆成偶数因子与奇数因子，它便化为，而后者小于或者。容易看出（利用）这一事实），当时，最后这个量小于，结论得证。

附录 2.PASCAL代码

下面这段程序实现了涓流算法，我们要感谢麦克莱斯特学院的学生 Simeon Simeonov 提供了它。这段代码利用了这样一个事实：在"预位数字"（predigit）队列中，最左成员右边总跟着一串 9，因此只需记住最左边那个预位数字以及 9 的个数即可。程序计算 π 的前 1000 位数字，要求使用带有 longint 数据类型（32 位整数）的 Pascal 版本。

program Pi_Spigot; const n = 1000; len = 10 * n div 3; var i, j, k, q, x, nines, predigit: integer; a: array[0..(len - 1)] of longint; begin for j := 0 to (len - 1) do a[j] := 2; {Start with 2s} nines := 0; predigit := 0;{First predigit is a 0} for j := 0 to (n - 1) do begin q := 0; for i := (len - 1) downto 0 do {Work backwards} begin x := 10 * a[i] + q * (i + 1); a[i] := x mod (2 * i + 1); q := x div (2 * i + 1); end; a[0] := q mod 10; q := q div 10; if q = 9 then nines := nines + 1 else if q = 10 then begin write(predigit + 1); for k := 0 to (nines - 1) do write(0); {zeros} predigit := 0; nines := 0; end else begin write(predigit); predigit := q; if nines <> 0 then begin for k := 0 to (nines - 1) do write(9); nines := 0; end; end; end; writeln(predigit); end.

(翻译注解：增加一段同等功能的C#代码，以作参考）

static void Pi_Spigot() { const int n = 1000; var len = 10 * n / 3; long i, j, k, q, x, nines = 0, predigit = 0; long[] a = new long[len]; for (j = 0; j < len; j++) a[j] = 2; for (j = 0; j < n; j++) { q = 0; for (i = len - 1; i >= 0; i--) { x = 10 * a[i] + q * (i + 1); a[i] = x % (2 * (i + 1) - 1); q = x / (2 * (i + 1) - 1); } a[0] = q % 10; q /= 10; if (q == 9) { nines++; } else if (q == 10) { Console.Write(predigit + 1); for (k = 0; k < nines; k++) { Console.Write(0); } predigit = 0; nines = 0; } else { //正常输出位 Console.Write(predigit); predigit = q; //如果存在多个9，一起输出 if (nines != 0) { for (k = 0; k < nines; k++) Console.Write('9'); nines = 0; } } } Console.WriteLine(predigit); }

校样补记。Mathematica 的最新版本（2.3）能够求出本文中出现的许多级数的和。只需大约一秒钟，它就能算出第 2 节开头那个关键级数的和为
，算出引理 1 证明中级数的和为，还能算出引理 4 证明中级数的和为 2。

参考文献

[Abd] S. Kamal Abdali（萨迈尔·卡迈勒·阿卜达利），《算法 393——任意精度的特殊级数求和》（Algorithm 393 — Special series summation with arbitrary precision），《ACM 通讯》（Comm. ACM）第 13 卷（1970），第 570 页。

[BBB] J. M. Borwein（乔纳森·博温）、P. B. Borwein（彼得·博温）、D. H. Bailey（戴维·贝利），《拉马努金、模方程与 π 的近似，或：如何计算 π 的十亿位数》（Ramanujan, modular equations, and approximations to pi, or How to compute one billion digits of pi），《美国数学月刊》第 96 卷（1989），第 201–219 页。

[Gos] R. W. Gosper（拉尔夫·戈斯珀），《级数加速》（Acceleration of series），报告第 304 号，麻省理工学院人工智能实验室，剑桥（马萨诸塞州），1974 年。

[Knu] D. E. Knuth（高德纳），《计算机程序设计艺术》（The Art of Computer Programming）第 2 卷，里丁（马萨诸塞州），Addison-Wesley 出版社，1981 年。

[Li] J. C. R. Li问题 E854，《美国数学月刊》第 56 卷（1949），第 633–635 页。

[Rab] S. Rabinowitz（斯坦·拉宾诺维茨），摘要 863-11-482：《一个 π 的涓流算法》（A spigot algorithm for pi），《美国数学会摘要》（Abstracts Amer. Math. Society）第 12 卷（1991），第 30 页。

[Sale] A. H. J. Sale（阿利斯泰尔·塞尔），《e 的多位有效数字计算》（The calculation of e to many significant digits），《计算机杂志》（Comput. J.）第 11 卷（1968），第 229–230 页。

MathPro出版社

信箱 713

韦斯特福德（西福德），马萨诸塞州 01886

72717.3515@compuserve.com

数学系

麦克莱斯特学院（Macalester College）

圣保罗，明尼苏达州 55105

wagon@macalstr.edu（原版笔误，实为：wagon@macalester.edu）