一个计算圆周率任意精度的spigot算法研究

plp626

最初接触到这算法是论坛的这个帖子： http://bbs.emath.ac.cn/viewthread.php?tid=1935 网上的那段C代码在这里贴下

1. long a=10000,b,c=2800,d,e,f[2801],g;
2. main(){for(;b-c;)f[b++]=a/5;
3. for(;d=0,g=c*2;c -=14,printf("%.4d",e+d/a),e=d%a)
4. for(b=c; d+=f*a,f=d%--g,d/=g--,--b; d*=b);}
5. 用到的级数是pi/2= sum{k!/(2*k-1)!!,k=0,1,2,...}

复制代码

两年前我把它用批处理代码改写后，让它可以计算任意精度。

1. @echo off&title 计算圆周率任意精度
2. :: windows系统直接保存为pi.bat，双击运行，默认计算100位精度
3. :: 命令行下执行，第一个参数（整数）为圆周率小数点的位数。
4. setlocal enabledelayedexpansion
5. if not %1.==. (set c=%1) else set c=100
6. set/a c=(c*100/3)+70,cc=c/10
7. for /l %%a in (1 1 %cc%)do set/a f_%%a=2000
8. for /l %%a in (%c% -132 100)do (set/a n=%%a/10,m=2*n-1
9. set/a "d=f_!n!*10000,f_!b!=d%%m,d=d/m,n-=1"
10. for /l %%b in (!n! -1 1)do (set/a n=%%b,m=2*n-1
11. set/a "d=d*n+f_!n!*10000,f_!n!=d%%m,d=d/m,n-=1"
12. )
13. set/a "an=e+d/10000,e=d%%10000"
14. if !an! lss 1000 set an=000!an!
15. set/p=!an:~-4!<nul
16. )
17. pause

复制代码

因为网上流传的那段的C代码没有任何可读性，消化那段C代码花了我好几个晚上的闲暇时间，但当时仅是为了在一个bat论坛娱乐，一点注释也没写。现在自己也是一点看不懂了。 对这样所谓的spigot算法，我仅仅知道这个代码，这对时间又对它感兴趣了，我想问，有没有计算整数开方的此类spigot算法，因为它无需建立大数据类型，实现起来很简洁。 我想对这类算法做深入的研究，不知论坛是否已有讨论过，或者有这方面的学习资料也好。 提前谢谢各位回帖。

zeroieme

建基在BBP公式上的，而BBP公式还没有理论上的突破。 http://zh.wikipedia.org/zh/%E8%B ... B%E5%85%AC%E5%BC%8F 我也幻想所有基本超越函数能BBP化，这样在超越函数运算能转化到多项式计算。

liangbch

1楼的代码太强了，强的我几乎完全看不懂，能否提供一些关于Dos批处理语法方面的参考资料，也好让我辈学习一下。

〇〇

nice job

plp626

1楼的代码太强了，强的我几乎完全看不懂，能否提供一些关于Dos批处理语法方面的参考资料，也好让我辈学习一下。 liangbch 发表于 2011-4-22 20:12

版主这样的 回复 有点让我不好意思了，你的很多帖子我都跟不上思路，惭愧~~~ 批处理方面的语法，你在cmd下输入help即可，然后就不断的help命令名（或者命令名/?），就可查看命令语法 如果嫌这样慢，想一次看完，cmd下粘贴如下代码，可以把他们全输出：

1. pushd %temp%
2. (for /f %a in ('"help|findstr/b [A-Z]"')do help %a)>help.txt
3. start help.txt
4. popd

复制代码

因为代码页是默认中文的，而中文的翻译又有歧义，也许看英文语法更专业些：

1. chcp 437
2. pushd %temp%
3. (for /f "skip=1" %a in ('"help|findstr/b [A-Z]"')do help %a)>help-english.txt
4. start help-english.txt
5. popd

复制代码

国内的批处理论坛，环境不好不说，版主的回复比较让人无语，新人去有时被“你怎么连这个都不会，”“你都是一年的会员了，怎么还问这个问题。。。” 这个外国批处理论坛很不错： www.dostips.com 至于这方面的资料，国内很少，非主流的东西，没有多少人写；有也只是些论坛的精华帖，但是确实不精华啊。还是少看了。

plp626

1# plp626 一楼代码有一处笔误f_!n!=d%%m写成了f_!b!=d%%m，导致后面精度丢失;两年了现在才发现。 现在改过来了：

1. @echo off&title 计算圆周率任意精度
2. setlocal enabledelayedexpansion
3. if not %1.==. (set c=%1) else set c=100
4. set/ac=c*100/3+70,cc=c/10,n=cc
5. for /l %%a in (1 1 !cc!)do set/af_%%a=2000
6. for /l %%a in (!c! -132 100)do (
7. set/an=%%a/10,m=2*n-1,d=f_!n!*10000,d=d/m,n-=1
8. for /l %%b in (!n! -1 1)do set/an=%%b,m=2*n-1,d=d*n+f_!n!*10000,f_!n!=d%%m,d=d/m,n-=1
9. set/aan="e+d/10000,e=d%%10000
10. if !an! lss 1000 set an=000!an!
11. set/p=!an:~-4!<nul
12. )
13. pause

复制代码

其实还可以直接粘贴在cmd下，这样更方便些

:: 粘贴本代码在cmd命令行下回车即可 cmd /Q/V: ON set "c=200" :: 这里设置的是输出200位（含整数位3） set/ac=c*100/3+70,cc=c/10,n=cc>nul (for /l %a in (1 1 !cc!)do set/af%a=2000 >nul)&for /l %a in (!c! -132 100)do ( set/an=%a/10,m=2*n-1,d=f!n!*10000,f!n!=d%m,d=d/m,n-=1 >nul for /l %b in (!n! -1 1)do set/an=%b,m=2*n-1,d=d*n+f!n!*10000,f!n!=d%m,d=d/m,n-=1 >nul set/aan=e+d/10000,e=d%10000 >nul if !an! lss 1000 set an=000!an! set/p=!an:~-4!

----------------------- 报告一个论坛bug，我这里从editplus上粘贴代码到回复框，还没点击编辑帖子直接就回复了，并且是俩回帖

zeroieme

批处理无疑可以处理很复杂的任务。毕竟是解释执行，不适合作数值计算。

mjs1wh

看不懂批处理，能改个BASIC的吗

plp626

因为所谓的spigot算法都是建立在 收敛级数 之上的，我想如何得到一个整数开方，开立方，或者开n次方的 spigot算法，就是找一个 通项相对简洁 的收敛级数。 找到一个资料，得知网上那个输出$sqrt(2)$ 小数点2400位的4行C代码（还是一点看不懂，但还是和一楼一样的spigot算法）

1. main(){int a=1000,b=0,c=1413,d,f[1414],n=800,k;
2. for(;b<c;f[b++]=14);
3. for(;n--;d+=*f*a,printf("%.3d",d/a),*f=d%a)
4. for(d=0,k=c;--k;d/=b,d*=2*k-1)f[k]=(d+=f[k]*a)%(b=100*k);}

复制代码

基于如下公式：

$sqrt(2)=7/{5*sqrt(1-1/50)}$

用到的是这个麦克劳林级数：

$1/(1-x)^n=1+nx+{n(n+1)}/{2!}x^2+... $

受到这个启发，我们把它扩展一下； 给上面级数两边同乘以一个有理数$q/p$得到

$q/p*(1/(1-x)^n)=q/p*(1+nx+{n(n+1)}/{2!}x^2+...) $

我们令一个整数N的开m次方$root{m}{N}$为：

$root{m}{N}=q/p*(1/(1-b/a)^{1/m}), -1 我们给他两边再做m次方得到

$N={q^m*a}/{p^m*(a-b)}$

为了收敛速度快和以后的通项简洁些，我们可以令b=1，a为正整数； 所以问题可以转化为求下面关于正整数a,p,q的不定方程（*）：

$a*q^m-(a-1)*N*p^m=0$

这里N，m为正整数常数，p,q互质. 对于这个不定方程的解，我们取a较大者，p,q较小者； 这样便可以得到任意一个正整数N开m次方的spigot算法的级数:

$root{m}{N}=q/p*(1/(1-1/a)^{1/m})$

$=q/p*(1+1/{ma}+1/{m(m+1)*2!}*1/{a^2}+...+1/{m(m+1)...(m+n-1)*n!}*1/a^n+...)$

$=q/p*(1+1/{ma}(1+1/{(m+1)a}(1+...1/{(m+n-1)a}(1+......))))$

问题现在来了， 解(*)处的不定方程，并求出最佳解，难住了我，大家帮我看看，提前谢谢了。

plp626

若不限定b=1，更一般的情况就是，解关于正整数a,p,q,非0整数b 的不定方程：

$a*q^m-(a-b)*N*p^m=0$

其中N,m为正整数常数，p,q互质，a,b互质，|b|