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

gracias

http://user.qzone.qq.com/414353346/blog/1288275811
http://numbers.computation.free. ... grams/tinycodes.pdf

plp626

11# gracias


恩，后面那个资料我看过，所以才如此感兴趣，只是两年了，对于细节我又忘记了；但大概的思路我还有；

我是要问任意整数开m次方的那个收敛级数，这个问题我已经有了想法，就是如9楼所讨论的

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

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

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

给定一个不大的正整数N，当n=2，3，4，5，。。。时，上面那个不定方程的解（a,b,p,q）的通项怎么求出，并给出最佳解。

plp626

http://user.qzone.qq.com/414353346/blog/1288275811
http://numbers.computation.free. ... grams/tinycodes.pdf

gracias 发表于 2011-4-24 11:07

谢谢你提供的关于这个算法的详细资料，认真看过后有种相见恨晚的感觉；

http://www.comlab.ox.ac.uk/resea ... lgprog-20031107.pdf
http://www.cs.umb.edu/~offner/files/pi.pdf
http://www.mathpropress.com/stan/bibliography/spigot.pdf
http://numbers.computation.free. ... grams/tinycodes.pdf

plp626

言归正传，继续思考这个不定方程的解法。。

liangbch

5# plp626

谢谢楼主，有时间好好学下批处理。

▄︻┻═┳一‥

当m=2时，也就是开方的时候，用程序暴力搜索1000以内的所有整数解得到：

只是俺的这个程序算法太笨蛋了，运行了半天才得到这些

1. 
   
2. +------+------+------+------+
   
3. |    N |    a |    p |    q |
   
4. +------+------+------+------+
   
5. |    2 |    9 |    3 |    4 |
   
6. |    2 |   50 |    5 |    7 |
   
7. |    2 |  289 |   17 |   24 |
   
8. |    3 |    4 |    2 |    3 |
   
9. |    3 |   49 |    7 |   12 |
   
10. |    3 |  676 |   26 |   45 |
    
11. |    5 |   81 |    9 |   20 |
    
12. |    6 |   25 |    5 |   12 |
    
13. |    6 |  243 |    9 |   22 |
    
14. |    7 |   64 |    8 |   21 |
    
15. |    8 |    9 |    3 |    8 |
    
16. |    8 |   50 |    5 |   14 |
    
17. |    8 |  289 |   17 |   48 |
    
18. |   10 |  361 |   19 |   60 |
    
19. |   11 |  100 |   10 |   33 |
    
20. |   12 |   49 |    7 |   24 |
    
21. |   12 |  676 |   13 |   45 |
    
22. |   13 |  325 |    5 |   18 |
    
23. |   14 |    8 |    2 |    7 |
    
24. |   14 |  225 |   15 |   56 |
    
25. |   15 |   16 |    4 |   15 |
    
26. |   15 |  961 |   31 |  120 |
    
27. |   18 |   50 |    5 |   21 |
    
28. |   18 |  289 |   17 |   72 |
    
29. |   20 |   81 |    9 |   40 |
    
30. |   21 |   28 |    2 |    9 |
    
31. |   22 |   99 |    3 |   14 |
    
32. |   23 |  576 |   24 |  115 |
    
33. |   24 |   25 |    5 |   24 |
    
34. |   24 |  243 |    9 |   44 |
    
35. |   27 |    4 |    2 |    9 |
    
36. |   27 |   49 |    7 |   36 |
    
37. |   27 |  676 |   26 |  135 |
    
38. |   28 |   64 |    4 |   21 |
    
39. |   30 |  121 |   11 |   60 |
    
40. |   32 |    9 |    3 |   16 |
    
41. |   32 |   50 |    5 |   28 |
    
42. |   32 |  289 |   17 |   96 |
    
43. |   33 |   12 |    2 |   11 |
    
44. |   33 |  529 |   23 |  132 |
    
45. |   34 |   18 |    3 |   17 |
    
46. |   35 |   36 |    6 |   35 |
    
47. |   39 |  625 |   25 |  156 |
    
48. |   40 |  361 |   19 |  120 |
    
49. |   42 |  169 |   13 |   84 |
    
50. |   44 |  100 |    5 |   33 |
    
51. |   45 |   81 |    3 |   20 |
    
52. |   48 |   49 |    7 |   48 |
    
53. |   48 |  676 |   13 |   90 |
    
54. |   50 |    9 |    3 |   20 |
    
55. |   50 |  289 |   17 |  120 |
    
56. |   52 |  325 |    5 |   36 |
    
57. |   54 |   25 |    5 |   36 |
    
58. |   54 |  243 |    3 |   22 |
    
59. |   54 |  755 |  349 |  938 |
    
60. |   55 |   45 |    3 |   22 |
    
61. |   56 |  225 |   15 |  112 |
    
62. |   57 |   76 |    2 |   15 |
    
63. |   58 |  887 |  959 |  134 |
    
64. |   60 |   16 |    2 |   15 |

复制代码

▄︻┻═┳一‥

上面输出表面有些N在1000内还无解

猜想某些N也许在计算机的32位整数内无解。

看来这个级数对于N开m次方不具有通用性。

liangbch

// 对于一个给定的浮点数，求一个一定范围内的p/q,使得p/q近似等于这个浮点数，有1种2分迭代法，可迅速逼近，下面给出c语言代码

1. #include <stdlib.h>
   
2. #include <stdio.h>
   
3. #include <math.h>
   
4. 
   
5. typedef unsigned long DWORD;
   
6. 
   
7. void print_pq(double f,DWORD limit)
   
8. {
   
9.         DWORD low_p,low_q,high_p,high_q,mid_p,mid_q;
   
10.        
    
11.         //初始化low_q,low_q,high_p,high_q
    
12.         low_q=DWORD(floor(f));
    
13.         high_q=DWORD(ceil(f));
    
14.         low_p=high_p=1;
    
15.        
    
16.         mid_p=low_p+high_p;
    
17.         mid_q=low_q+high_q;
    
18.        
    
19.         while ( mid_p<limit )
    
20.         {
    
21.                 printf("%u/%u=%.12f\n",mid_q,mid_p,double(mid_q)/double(mid_p) );
    
22. 
    
23.                 if ( double(mid_q)/double(mid_p) > f)
    
24.                 {
    
25.                         high_p=mid_p; high_q=mid_q;
    
26.                 }
    
27.                 else
    
28.                 {
    
29.                         low_p=mid_p; low_q=mid_q;
    
30.                 }
    
31.                 mid_p=low_p+high_p;
    
32.                 mid_q=low_q+high_q;
    
33.         }
    
34. }
    
35. 
    
36. int main(int argc, char* argv[])
    
37. {
    
38.         print_pq(sqrt(2),1000000); //求分子在1000000以内的sqrt(2)有理数逼近
    
39.         return 0;
    
40. }

复制代码

liangbch

这里给出楼上代码的运行结果

1. 3/2=1.500000000000
   
2. 4/3=1.333333333333
   
3. 7/5=1.400000000000
   
4. 10/7=1.428571428571
   
5. 17/12=1.416666666667
   
6. 24/17=1.411764705882
   
7. 41/29=1.413793103448
   
8. 58/41=1.414634146341
   
9. 99/70=1.414285714286
   
10. 140/99=1.414141414141
    
11. 239/169=1.414201183432
    
12. 338/239=1.414225941423
    
13. 577/408=1.414215686275
    
14. 816/577=1.414211438475
    
15. 1393/985=1.414213197970
    
16. 1970/1393=1.414213926777
    
17. 3363/2378=1.414213624895
    
18. 4756/3363=1.414213499851
    
19. 8119/5741=1.414213551646
    
20. 11482/8119=1.414213573100
    
21. 19601/13860=1.414213564214
    
22. 27720/19601=1.414213560533
    
23. 47321/33461=1.414213562057
    
24. 66922/47321=1.414213562689
    
25. 114243/80782=1.414213562427
    
26. 161564/114243=1.414213562319
    
27. 275807/195025=1.414213562364
    
28. 390050/275807=1.414213562382
    
29. 665857/470832=1.414213562375
    
30. 941664/665857=1.414213562372

复制代码

▄︻┻═┳一‥

18# liangbch

你说的可是最佳有理逼近？

那个是连分数，只是当分母超过 21亿的时候，就不好算了，要自己构造大数据类型做除法了。。。。

还是不如spigot算法 省事。。。