一个计算圆周率任意精度的spigot算法研究
最初接触到这算法是论坛的这个帖子：http://bbs.emath.ac.cn/viewthread.php?tid=1935
网上的那段C代码在这里贴下long a=10000,b,c=2800,d,e,f,g;
main(){for(;b-c;)f=a/5;
for(;d=0,g=c*2;c -=14,printf("%.4d",e+d/a),e=d%a)
for(b=c; d+=f*a,f=d%--g,d/=g--,--b; d*=b);}
用到的级数是pi/2= sum{k!/(2*k-1)!!,k=0,1,2,...}两年前我把它用批处理代码改写后，让它可以计算任意精度。@echo off&title 计算圆周率任意精度
:: windows系统直接保存为pi.bat，双击运行，默认计算100位精度
:: 命令行下执行，第一个参数（整数）为圆周率小数点的位数。
setlocal enabledelayedexpansion
if not %1.==. (set c=%1) else set c=100
set/a c=(c*100/3)+70,cc=c/10
for /l %%a in (1 1 %cc%)do set/a f_%%a=2000
for /l %%a in (%c% -132 100)do (set/a n=%%a/10,m=2*n-1
set/a "d=f_!n!*10000,f_!b!=d%%m,d=d/m,n-=1"
for /l %%b in (!n! -1 1)do (set/a n=%%b,m=2*n-1
set/a "d=d*n+f_!n!*10000,f_!n!=d%%m,d=d/m,n-=1"
)
set/a "an=e+d/10000,e=d%%10000"
if !an! lss 1000 set an=000!an!
set/p=!an:~-4!<nul
)
pause因为网上流传的那段的C代码没有任何可读性，消化那段C代码花了我好几个晚上的闲暇时间，但当时仅是为了在一个bat论坛娱乐，一点注释也没写。现在自己也是一点看不懂了。
对这样所谓的spigot算法，我仅仅知道这个代码，这对时间又对它感兴趣了，我想问，有没有计算整数开方的此类spigot算法，因为它无需建立大数据类型，实现起来很简洁。
我想对这类算法做深入的研究，不知论坛是否已有讨论过，或者有这方面的学习资料也好。
提前谢谢各位回帖。 建基在BBP公式上的，而BBP公式还没有理论上的突破。
http://zh.wikipedia.org/zh/%E8%B4%9D%E5%88%A9%EF%BC%8D%E6%B3%A2%E5%B0%94%E6%B8%A9%EF%BC%8D%E6%99%AE%E5%8A%B3%E5%A4%AB%E5%85%AC%E5%BC%8F
我也幻想所有基本超越函数能BBP化，这样在超越函数运算能转化到多项式计算。 1楼的代码太强了，强的我几乎完全看不懂，能否提供一些关于Dos批处理语法方面的参考资料，也好让我辈学习一下。 nice job 1楼的代码太强了，强的我几乎完全看不懂，能否提供一些关于Dos批处理语法方面的参考资料，也好让我辈学习一下。
liangbch 发表于 2011-4-22 20:12 http://bbs.emath.ac.cn/images/common/back.gif
版主这样的 回复 有点让我不好意思了，你的很多帖子我都跟不上思路，惭愧~~~
批处理方面的语法，你在cmd下输入help即可，然后就不断的help命令名（或者命令名/?），就可查看命令语法
如果嫌这样慢，想一次看完，cmd下粘贴如下代码，可以把他们全输出：pushd %temp%
(for /f %a in ('"help|findstr/b "')do help %a)>help.txt
start help.txt
popd因为代码页是默认中文的，而中文的翻译又有歧义，也许看英文语法更专业些：chcp 437
pushd %temp%
(for /f "skip=1" %a in ('"help|findstr/b "')do help %a)>help-english.txt
start help-english.txt
popd国内的批处理论坛，环境不好不说，版主的回复比较让人无语，新人去有时被“你怎么连这个都不会，”“你都是一年的会员了，怎么还问这个问题。。。”
这个外国批处理论坛很不错： www.dostips.com
至于这方面的资料，国内很少，非主流的东西，没有多少人写；有也只是些论坛的精华帖，但是确实不精华啊。还是少看了。 1# plp626
一楼代码有一处笔误f_!n!=d%%m写成了f_!b!=d%%m，导致后面精度丢失;两年了现在才发现。
现在改过来了：@echo off&title 计算圆周率任意精度
setlocal enabledelayedexpansion
if not %1.==. (set c=%1) else set c=100
set/ac=c*100/3+70,cc=c/10,n=cc
for /l %%a in (1 1 !cc!)do set/af_%%a=2000
for /l %%a in (!c! -132 100)do (
set/an=%%a/10,m=2*n-1,d=f_!n!*10000,d=d/m,n-=1
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
set/aan="e+d/10000,e=d%%10000
if !an! lss 1000 set an=000!an!
set/p=!an:~-4!<nul
)
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!<nul
)
-----------------------
报告一个论坛bug，我这里从editplus上粘贴代码到回复框，还没点击编辑帖子直接就回复了，并且是俩回帖 批处理无疑可以处理很复杂的任务。毕竟是解释执行，不适合作数值计算。 看不懂批处理，能改个BASIC的吗 本帖最后由 plp626 于 2011-4-24 11:26 编辑
因为所谓的spigot算法都是建立在 收敛级数 之上的，我想如何得到一个整数开方，开立方，或者开n次方的 spigot算法，就是找一个 通项相对简洁 的收敛级数。
找到一个资料，得知网上那个输出sqrt(2) 小数点2400位的4行C代码（还是一点看不懂，但还是和一楼一样的spigot算法）main(){int a=1000,b=0,c=1413,d,f,n=800,k;
for(;b<c;f=14);
for(;n--;d+=*f*a,printf("%.3d",d/a),*f=d%a)
for(d=0,k=c;--k;d/=b,d*=2*k-1)f=(d+=f*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<b/a<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+......))))
问题现在来了， 解(*)处的不定方程，并求出最佳解，难住了我，大家帮我看看，提前谢谢了。 若不限定b=1，更一般的情况就是，解关于正整数a,p,q,非0整数b 的不定方程：
a*q^m-(a-b)*N*p^m=0
其中N,m为正整数常数，p,q互质，a,b互质，|b|<a