数学研发论坛

 找回密码
 欢迎注册
查看: 15880|回复: 72

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

[复制链接]
发表于 2011-4-22 11:10:32 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
最初接触到这算法是论坛的这个帖子:
/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算法,因为它无需建立大数据类型,实现起来很简洁。

我想对这类算法做深入的研究,不知论坛是否已有讨论过,或者有这方面的学习资料也好。

提前谢谢各位回帖。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-4-22 12:13:28 | 显示全部楼层
建基在BBP公式上的,而BBP公式还没有理论上的突破。

http://zh.wikipedia.org/zh/%E8%B ... B%E5%85%AC%E5%BC%8F

我也幻想所有基本超越函数能BBP化,这样在超越函数运算能转化到多项式计算。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-4-22 20:12:26 | 显示全部楼层
1楼的代码太强了,强的我几乎完全看不懂,能否提供一些关于Dos批处理语法方面的参考资料,也好让我辈学习一下。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-4-23 10:41:10 | 显示全部楼层
nice job
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2011-4-23 15:50:08 | 显示全部楼层
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

至于这方面的资料,国内很少,非主流的东西,没有多少人写;有也只是些论坛的精华帖,但是确实不精华啊。还是少看了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2011-4-23 15:51:03 | 显示全部楼层
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!<nul
)

-----------------------
报告一个论坛bug,我这里从editplus上粘贴代码到回复框,还没点击编辑帖子直接就回复了,并且是俩回帖
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-4-23 16:51:54 | 显示全部楼层
批处理无疑可以处理很复杂的任务。毕竟是解释执行,不适合作数值计算。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2011-4-23 17:02:18 | 显示全部楼层
看不懂批处理,能改个BASIC的吗
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2011-4-24 08:47:19 | 显示全部楼层
本帖最后由 plp626 于 2011-4-24 11:26 编辑

因为所谓的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<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+......))))$



问题现在来了, 解(*)处的不定方程,并求出最佳解,难住了我,大家帮我看看,提前谢谢了。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2011-4-24 10:57:36 | 显示全部楼层
若不限定b=1,更一般的情况就是,解关于正整数a,p,q,非0整数b 的不定方程:
$a*q^m-(a-b)*N*p^m=0$

其中N,m为正整数常数,p,q互质,a,b互质,|b|<a
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2019-1-20 08:23 , Processed in 0.055013 second(s), 17 queries .

Powered by Discuz! X3.4

© 2001-2017 Comsenz Inc.

快速回复 返回顶部 返回列表