数学研发网设为首页收藏本站
切换到窄版

 找回密码
 欢迎注册
数学研发论坛»论坛 › 【算法研发】 › 算法交流 › 一个计算圆周率任意精度的spigot算法研究 ...
12345678下一页
返回列表 发新帖
楼主: plp626

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

[复制链接]
liangbch
发表于 2011-4-26 11:44:26 | 显示全部楼层
以开平方为例，解释一下我程序的功能 楼主的本质是想找到这样一个表达式，对于给定的数N，求正整数p,q和a. 使得p,q,a 满足 sqrt(N)=(q/p)sqrt(1-1/a), 显然p,q是sqrt(N)的一个有理数逼近。下面的代码中的函数print_pq打印出p小于limit且满足以上公式的表达式
  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <math.h>
  5. typedef unsigned long DWORD;
  6. typedef __int64 INT64;
  8. void print_pq(INT64 N,INT64 limit)
  9. {
  10. INT64 low_p,low_q,high_p,high_q,mid_p,mid_q;
  11. double f=sqrt(N);
  13. //初始化low_q,low_q,high_p,high_q
  14. low_q=DWORD(floor(f));
  15. high_q=DWORD(ceil(f));
  16. low_p=high_p=1;
  18. mid_p=low_p+high_p;
  19. mid_q=low_q+high_q;
  21. while ( mid_p<limit )
  22. {
  23. INT64 q=mid_q;
  24. INT64 p=mid_p;
  26. if ( double(q)/double(p)>f)
  27. {
  28. INT64 t1=q*q-N*p*p;
  29. INT64 t2=q*q;
  30. if ( t2 % N ==0)
  31. {
  32. t1/=N;
  33. t2/=N;
  34. }
  35. printf("\$sqrt(%I64d)=(%I64d/%I64d)sqrt(1-%I64d/%I64d)\$\n",N,q,p,t1,t2);
  36. }
  39. if ( double(mid_q)/double(mid_p) > f)
  40. {
  41. high_p=mid_p; high_q=mid_q;
  42. }
  43. else
  44. {
  45. low_p=mid_p; low_q=mid_q;
  46. }
  47. mid_p=low_p+high_p;
  48. mid_q=low_q+high_q;
  49. }
  50. }
  52. int main(int argc, char* argv[])
  53. {
  54. print_pq(2,1000000I64);
  55. print_pq(3,1000000I64);
  56. return 0;
  57. }
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

plp626
 楼主| 发表于 2011-4-26 11:46:01 | 显示全部楼层
31# liangbch XX？ 没看明白。。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

plp626
 楼主| 发表于 2011-4-26 11:48:59 | 显示全部楼层
$sqrt(n)=(q/p)sqrt(1-1/a)$ liangbch 发表于 2011-4-26 11:44
就是啊，我怎么受启发选了那个级数呢，这个更简单些。。 我再研究下。。
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2011-4-26 11:50:03 | 显示全部楼层
贴出31楼的运行结果
  1. $sqrt(2)=(3/2)sqrt(1-1/9)$
  2. $sqrt(2)=(10/7)sqrt(1-1/50)$
  3. $sqrt(2)=(17/12)sqrt(1-1/289)$
  4. $sqrt(2)=(58/41)sqrt(1-1/1682)$
  5. $sqrt(2)=(99/70)sqrt(1-1/9801)$
  6. $sqrt(2)=(338/239)sqrt(1-1/57122)$
  7. $sqrt(2)=(577/408)sqrt(1-1/332929)$
  8. $sqrt(2)=(1970/1393)sqrt(1-1/1940450)$
  9. $sqrt(2)=(3363/2378)sqrt(1-1/11309769)$
  10. $sqrt(2)=(11482/8119)sqrt(1-1/65918162)$
  11. $sqrt(2)=(19601/13860)sqrt(1-1/384199201)$
  12. $sqrt(2)=(66922/47321)sqrt(1-1/2239277042)$
  13. $sqrt(2)=(114243/80782)sqrt(1-1/13051463049)$
  14. $sqrt(2)=(390050/275807)sqrt(1-1/76069501250)$
  15. $sqrt(2)=(665857/470832)sqrt(1-1/443365544449)$
  16. $sqrt(3)=(7/4)sqrt(1-1/49)$
  17. $sqrt(3)=(26/15)sqrt(1-1/676)$
  18. $sqrt(3)=(97/56)sqrt(1-1/9409)$
  19. $sqrt(3)=(362/209)sqrt(1-1/131044)$
  20. $sqrt(3)=(1351/780)sqrt(1-1/1825201)$
  21. $sqrt(3)=(5042/2911)sqrt(1-1/25421764)$
  22. $sqrt(3)=(18817/10864)sqrt(1-1/354079489)$
  23. $sqrt(3)=(70226/40545)sqrt(1-1/4931691076)$
  24. $sqrt(3)=(262087/151316)sqrt(1-1/68689595569)$
  25. $sqrt(3)=(978122/564719)sqrt(1-1/956722646884)$
复制代码

评分

参与人数 1威望 +1 金币 +1 贡献 +1 鲜花 +1 收起 理由
plp626 + 1 + 1 + 1 + 1 谢谢分享

查看全部评分

毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

plp626
 楼主| 发表于 2011-4-26 11:54:29 | 显示全部楼层
31# liangbch 不过，$(q/p)sqrt(1-b/a)$相比较于$(q/p)*1/sqrt(1-b/a)$那个级数通项稍复杂点点.
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2011-4-26 11:59:42 | 显示全部楼层
贴出 (1-x) n次方根的展开式： 公式4.gif
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

plp626
 楼主| 发表于 2011-4-26 12:47:16 | 显示全部楼层
这个是$1/root{m}{1-x}$的展示
$1/root{m}{1-x}=1+x/m+(1+m)/{m^2*2!}*x^2+{(1+m)(1+2m)}/{m^3*3!}*x^3+......$
$=(1+1/mx(1+{1+m}/{2m}x(1+{1+2m}/{3m}x(1+...))))$
|x|<1 ------------------------------------------------
$root{m}{1-x}=1-x/m-(m-1)/{m^2*2!}*x^2-{(m-1)(2m-1)}/{m^3*3!}*x^3+......$
$=(1-x/m(1+{m-1}/{2m}x(1-{2m-1}/{3m}x(1+...(1-... )))))$
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2011-4-26 12:53:10 | 显示全部楼层
OK, 换成 $p/q*1/sqrt(1-1/a)$的形式
  1. #include <stdlib.h>
  2. #include <stdio.h>
  3. #include <math.h>
  5. typedef unsigned long DWORD;
  6. typedef __int64 INT64;
  8. void print_pq(INT64 N,INT64 limit)
  9. {
  10. INT64 low_p,low_q,high_p,high_q,mid_p,mid_q;
  11. double f=sqrt(N);
  13. //初始化low_q,low_q,high_p,high_q
  14. low_q=DWORD(floor(f));
  15. high_q=DWORD(ceil(f));
  16. low_p=high_p=1;
  18. mid_p=low_p+high_p;
  19. mid_q=low_q+high_q;
  21. while ( mid_p<limit )
  22. {
  23. INT64 q=mid_q;
  24. INT64 p=mid_p;
  26. if ( double(q)/double(p)<f)
  27. {
  28. INT64 t1=p*p*N-q*q;
  29. INT64 t2=p*p*N;
  30. if ( t1 % N ==0)
  31. {
  32. t1/=N;
  33. t2/=N;
  34. }
  35. printf("sqrt(%I64d)=\$(%I64d/%I64d)1/sqrt(1-%I64d/%I64d)\$\n",N,q,p,t1,t2);
  36. }
  39. if ( double(mid_q)/double(mid_p) > f)
  40. {
  41. high_p=mid_p; high_q=mid_q;
  42. }
  43. else
  44. {
  45. low_p=mid_p; low_q=mid_q;
  46. }
  47. mid_p=low_p+high_p;
  48. mid_q=low_q+high_q;
  49. }
  50. }
  52. int main(int argc, char* argv[])
  53. {
  54. print_pq(2,1000000I64);
  55. print_pq(3,1000000I64);
  56. return 0;
  57. }
复制代码
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

liangbch
发表于 2011-4-26 13:05:16 | 显示全部楼层
贴出楼上代码的运行结果 sqrt(2)=$(4/3)1/sqrt(1-1/9)$ sqrt(2)=$(7/5)1/sqrt(1-1/50)$ sqrt(2)=$(24/17)1/sqrt(1-1/289)$ sqrt(2)=$(41/29)1/sqrt(1-1/1682)$ sqrt(2)=$(140/99)1/sqrt(1-1/9801)$ sqrt(2)=$(239/169)1/sqrt(1-1/57122)$ sqrt(2)=$(816/577)1/sqrt(1-1/332929)$ sqrt(2)=$(1393/985)1/sqrt(1-1/1940450)$ sqrt(2)=$(4756/3363)1/sqrt(1-1/11309769)$ sqrt(2)=$(8119/5741)1/sqrt(1-1/65918162)$ sqrt(2)=$(27720/19601)1/sqrt(1-1/384199201)$ sqrt(2)=$(47321/33461)1/sqrt(1-1/2239277042)$ sqrt(2)=$(161564/114243)1/sqrt(1-1/13051463049)$ sqrt(2)=$(275807/195025)1/sqrt(1-1/76069501250)$ sqrt(2)=$(941664/665857)1/sqrt(1-1/443365544449)$ sqrt(3)=$(3/2)1/sqrt(1-1/4)$ sqrt(3)=$(5/3)1/sqrt(1-2/27)$ sqrt(3)=$(12/7)1/sqrt(1-1/49)$ sqrt(3)=$(19/11)1/sqrt(1-2/363)$ sqrt(3)=$(45/26)1/sqrt(1-1/676)$ sqrt(3)=$(71/41)1/sqrt(1-2/5043)$ sqrt(3)=$(168/97)1/sqrt(1-1/9409)$ sqrt(3)=$(265/153)1/sqrt(1-2/70227)$ sqrt(3)=$(627/362)1/sqrt(1-1/131044)$ sqrt(3)=$(989/571)1/sqrt(1-2/978123)$ sqrt(3)=$(2340/1351)1/sqrt(1-1/1825201)$ sqrt(3)=$(3691/2131)1/sqrt(1-2/13623483)$ sqrt(3)=$(8733/5042)1/sqrt(1-1/25421764)$ sqrt(3)=$(13775/7953)1/sqrt(1-2/189750627)$ sqrt(3)=$(32592/18817)1/sqrt(1-1/354079489)$ sqrt(3)=$(51409/29681)1/sqrt(1-2/2642885283)$ sqrt(3)=$(121635/70226)1/sqrt(1-1/4931691076)$ sqrt(3)=$(191861/110771)1/sqrt(1-2/36810643323)$ sqrt(3)=$(453948/262087)1/sqrt(1-1/68689595569)$ sqrt(3)=$(716035/413403)1/sqrt(1-2/512706121227)$ sqrt(3)=$(1694157/978122)1/sqrt(1-1/956722646884)$
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

plp626
 楼主| 发表于 2011-4-26 13:13:33 | 显示全部楼层
本帖最后由 plp626 于 2011-4-26 14:07 编辑 我现在发现把b/a强制换成1/a会 使得p,q数值过大，p超过(2^32/10000)，做分母的时候，得到的商就不够4位了，这样不好统一输出。。。 当然b的绝对值越小越好，还是假定b的绝对值小于10把， 这样不会使得p太大，又适合于32位整数的直接除法
毋因群疑而阻独见　　毋任己意而废人言
毋私小惠而伤大体　　毋借公论以快私情
回复 支持 反对

使用道具 举报

12345678下一页
返回列表 发新帖
高级模式
B Color Image Link Quote Code Smilies
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

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

GMT+8, 2024-11-22 05:04 , Processed in 0.030150 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

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