计算π数位的涓流算法作者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。引理 1a若则特别地。b使用混合进制基数 b 的表示是唯一的。c的整数部分是小数部分是。在这个计数系统中某些无理数会变成循环周期的。例如这正好等价于把无穷级数重新表述为在这个系统中有理数对应于有限的会终止的数字序列见附录 1 引理 2。求区间内实数 x 的十进制各位数字可以这样做取 x 的整数部分把 x 的小数部分乘以 10再取所得结果的整数部分再把新的余下小数部分乘以 10如此继续下去。在某些混合进制下这个过程格外简单。若那么。后者未必是一个正规表示有些数字可能过大。但我们可以对各数字减量——把它们对取模其中 i 是进制序列 b 中对应元素的分母。从最右端开始做这些取模运算并把所得的商逐次向左进位最终便得到 10x 的正规表示。因此乘以 10 在算法上是直接明了的。而借助引理 1(c)对 b-表示取整数部分与小数部分同样容易。现在我们可以给出求的前位十进制数字的算法了。其正确性证明——即说明个混合进制数字足以得到位十进制数字的误差分析——见附录 1 的引理 3。算法 e-spigote 涓流算法初始化首位数字取 2并把长度为的数组 A 初始化为重复次乘以 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沃利斯乘积导出另一种做法是对级数施加一种称为欧拉变换Eulers transform的加速技巧。这些证明连同另外三个证明以及对早期文献的指引都可以在Li中找到。我们用表示奇数的连乘积于是该级数等价于可以展开为最后这个表达式引出了混合进制基数其各列为即第位的基数为相对于它干脆就是。在进制的正规表示中第位的数字必须落在区间内。遗憾的是基数不如基数那样方便好用。引理 4证明见附录 1。进制中各位取最大值的那个数表示 2因此形如的正规表示落在 0 与 2 之间。引理 4 意味着进制表示不唯一。例如由此与算法更相关的是在进制下取整数部分和小数部分并不像进制那样直截了当。数的整数部分要么是、要么是取决于是落在还是之内。这个麻烦可以这样克服在下一轮迭代中把的个位数字留在原处而把的十位数字称为一个预位数字predigit。预位数字必须暂时扣住不发因为偶尔大约每 20 次迭代出现一次下一个预位数字会是 10当进位值介于 0 到 19 之间大于 10且同时留下的个位数字是 9它在乘以 10那一步变成 90时这种情况就会发生。此时必须先把扣住的那个数加 1再把它释放出来。算法的具体细节见下文本文开头的讲述曾回避了这个偶尔出现 10的问题。关于个混合进制数字足以给出的位数字的证明见附录 1引理 5。附录 2 给出了该算法的一个 Pascal 实现。涓流算法初始化为一个长度为的数组。重复次:乘以 10将 A 的每个元素乘以 10。把 A 化为正规形式从右端开始将 A 的第个元素对应基数中的第项对取模得到商和余数。余数留在原位把向左进一位。最后一次被进出的整数在的那个位置产生可能大到 19。取得下一个预位数字将 A 的最左元素最大不超过对 10 取模。所得的商就是的新预位数字余数留在原位。调整预位数字若既不是 9 也不是 10就把扣住的各预位数字作为 π 的真数字释放输出并扣住。若是 9就把加入扣住的预位数字队列。若是 10则把当前预位数字置为 0 并继续扣住把所有其他扣住的预位数字加 19 变成 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 一样简单的算法。问题. 是否存在一个由有理数构成的基数使得的表示是周期的、或构成等差数列并且对而言恒等于其整数部分GosperGos, 第 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 这个字串那么算法中调整预位数字的那冗长部分就可以省去。这能带来一个极其短小的计算机程序。例如RabinowitzRab就利用这个想法展示过一个仅 14 行的 Fortran 程序能输出 π 的 1000 位十进制数字。最后我们还提一句这个算法可以并行化并行之后它在大约 10,000 位以内会快得惊人。关于其他函数的涓流算法的例子可参见Abd。附录1.五个引理引理1a如果特别地。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韦斯特福德西福德马萨诸塞州 0188672717.3515compuserve.com数学系麦克莱斯特学院Macalester College圣保罗明尼苏达州 55105wagonmacalstr.edu原版笔误实为wagonmacalester.edu 
长河君 · 2026/10/31 20:19:13 · 阅读 1.2k 计算 π 数位的涓流算法 （翻译）
计算π数位的涓流算法作者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。引理 1a若则特别地。b使用混合进制基数 b 的表示是唯一的。c的整数部分是小数部分是。在这个计数系统中某些无理数会变成循环周期的。例如这正好等价于把无穷级数重新表述为在这个系统中有理数对应于有限的会终止的数字序列见附录 1 引理 2。求区间内实数 x 的十进制各位数字可以这样做取 x 的整数部分把 x 的小数部分乘以 10再取所得结果的整数部分再把新的余下小数部分乘以 10如此继续下去。在某些混合进制下这个过程格外简单。若那么。后者未必是一个正规表示有些数字可能过大。但我们可以对各数字减量——把它们对取模其中 i 是进制序列 b 中对应元素的分母。从最右端开始做这些取模运算并把所得的商逐次向左进位最终便得到 10x 的正规表示。因此乘以 10 在算法上是直接明了的。而借助引理 1(c)对 b-表示取整数部分与小数部分同样容易。现在我们可以给出求的前位十进制数字的算法了。其正确性证明——即说明个混合进制数字足以得到位十进制数字的误差分析——见附录 1 的引理 3。算法 e-spigote 涓流算法初始化首位数字取 2并把长度为的数组 A 初始化为重复次乘以 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沃利斯乘积导出另一种做法是对级数施加一种称为欧拉变换Eulers transform的加速技巧。这些证明连同另外三个证明以及对早期文献的指引都可以在Li中找到。我们用表示奇数的连乘积于是该级数等价于可以展开为最后这个表达式引出了混合进制基数其各列为即第位的基数为相对于它干脆就是。在进制的正规表示中第位的数字必须落在区间内。遗憾的是基数不如基数那样方便好用。引理 4证明见附录 1。进制中各位取最大值的那个数表示 2因此形如的正规表示落在 0 与 2 之间。引理 4 意味着进制表示不唯一。例如由此与算法更相关的是在进制下取整数部分和小数部分并不像进制那样直截了当。数的整数部分要么是、要么是取决于是落在还是之内。这个麻烦可以这样克服在下一轮迭代中把的个位数字留在原处而把的十位数字称为一个预位数字predigit。预位数字必须暂时扣住不发因为偶尔大约每 20 次迭代出现一次下一个预位数字会是 10当进位值介于 0 到 19 之间大于 10且同时留下的个位数字是 9它在乘以 10那一步变成 90时这种情况就会发生。此时必须先把扣住的那个数加 1再把它释放出来。算法的具体细节见下文本文开头的讲述曾回避了这个偶尔出现 10的问题。关于个混合进制数字足以给出的位数字的证明见附录 1引理 5。附录 2 给出了该算法的一个 Pascal 实现。涓流算法初始化为一个长度为的数组。重复次:乘以 10将 A 的每个元素乘以 10。把 A 化为正规形式从右端开始将 A 的第个元素对应基数中的第项对取模得到商和余数。余数留在原位把向左进一位。最后一次被进出的整数在的那个位置产生可能大到 19。取得下一个预位数字将 A 的最左元素最大不超过对 10 取模。所得的商就是的新预位数字余数留在原位。调整预位数字若既不是 9 也不是 10就把扣住的各预位数字作为 π 的真数字释放输出并扣住。若是 9就把加入扣住的预位数字队列。若是 10则把当前预位数字置为 0 并继续扣住把所有其他扣住的预位数字加 19 变成 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 一样简单的算法。问题. 是否存在一个由有理数构成的基数使得的表示是周期的、或构成等差数列并且对而言恒等于其整数部分GosperGos, 第 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 这个字串那么算法中调整预位数字的那冗长部分就可以省去。这能带来一个极其短小的计算机程序。例如RabinowitzRab就利用这个想法展示过一个仅 14 行的 Fortran 程序能输出 π 的 1000 位十进制数字。最后我们还提一句这个算法可以并行化并行之后它在大约 10,000 位以内会快得惊人。关于其他函数的涓流算法的例子可参见Abd。附录1.五个引理引理1a如果特别地。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韦斯特福德西福德马萨诸塞州 0188672717.3515compuserve.com数学系麦克莱斯特学院Macalester College圣保罗明尼苏达州 55105wagonmacalstr.edu原版笔误实为wagonmacalester.edu