圆周率计算器

Ickiverar

完全由我自己编的一个圆周率计算器。控制台程序，运行在x64架构windows系统下，需要 vc redist 2012 支持（否则会提醒什么 MSVCR110.dll 找不到，找不到就下个 vc redist 2012 x86 & x64）。
至多可以计算圆周率到2.68亿（256M）位。其实更多也可以，但我的机子内存不够，没试过。
比不上现在比较有名的圆周率计算软件，但也算我的努力成果了。

Ickiverar

@KeyTo9_Fans
据说复杂度是O（N (lnN) ^3）,但事实是多少我没有太在意过。总之肯定优于O（N^a）, a>1 这样的。

liangbch

如果大数乘法的复杂度提 M(P）,p是精度。那么使用AGM(superPI使用的算法)算法，计算圆周率的算法的复杂度是ln(p)*M(p)

只是呼吸

我的是32位的，运行不了。

只是呼吸

现在，我终于用自己写的程序计算圆周率到一百万位了！
我自己造轮子，慢慢的写、慢慢的改，一步一步的，能把圆周率计算到计算到一百万位，还是蛮高兴的。

我用到的“最高级”的乘法是toom3,其次是karatusba乘法，基础乘法是comba。
除法是用Newton除法。
感谢数学研发论坛的支持。我的汇编就是学习了

只是呼吸

x86上128位二进制乘法最快速算法征解
https://bbs.emath.ac.cn/forum.ph ... 16&fromuid=2542

只是呼吸

入门的。

感谢liangbch 提供的加法汇编进位技术。

感谢liangbch 提供的开平方算法。

我的计算是采用agl算法。计算一百万位pi，迭代19次。我有一台笔记本、一台台式机，笔记本年代较久，速度慢，台式机速度快。这两台电脑的计算结果如下：（一百万位pi）

采用1000000000进制（10的九次方）

笔记本： 113秒   台式机：  38秒

采用1000000000000000000进制（10的18次方）

笔记本： 46秒   台式机：  16秒

只是呼吸

下面贴出计算pi的agl代码。这里面是有很大的坑，我把它指出来，再填了它，如果要用agl算pi的，那就借鉴去，以免掉坑里。

1. 
   
2. void re_give(Number& src, int b)
   
3. {
   
4.         src.pbig[0] = b;
   
5.         src.used = 1;
   
6.         src.decimal = 0;
   
7. }
   
8. static void re_give(Number& src, int a, int b)
   
9. {
   
10.         src.pbig[0] = 250000000;
    
11.         src.used = 1;
    
12.         src.decimal = 1;
    
13. }
    
14. int pi_agl(Number& pi, int pilen)
    
15. {
    
16.         int len;
    
17.         len = 2 * (pilen / 9 + 1);
    
18.         Number a(len), b(len), t(len), a1(len), b1(len), t1(len);//pi(100000);
    
19.         int w = 0, p = 1;
    
20.         ////a = 1.0;a,b,t,p赋初值
    
21.         ////b = 1 / sqrt(2);
    
22.         ////t = 1.0 / 4;
    
23.         ////p = 1;
    
24.         //
    
25.         re_give(a, 1);
    
26.         re_give(b1, 2);
    
27.         mp_sqrt(b, b1, pilen);///////
    
28.         div_one(b, b, 2);////////////////
    
29.         re_give(t, 1, 4);// t= 1/4
    
30.         std::cout << "初始化计算完成" << std::endl;
    
31.         while (w < 19)
    
32.         {
    
33.                 //a1 = (a + b) / 2;
    
34.                 //b1 = sqrt(a * b);
    
35.                 //t1 = t - p * ((a - (a + b) / 2) * (a - (a + b) / 2));
    
36.                 //p1 = 2 * p;
    
37.                 //pi = ((a1 + b1) * (a1 + b1)) / (4 * t1);
    
38.                 re_add(a1, a, b);//////---------------------------       
    
39.                 div_one(a1, a1, 2);////  a1
    
40.                 mp_mul(t1, a, b);/////-------------------------------乘法
    
41.                 re_shr(t1, a.used);///-----------------;右移
    
42.                 mp_sqrt(b1, t1, pilen);//  b1               
    
43.                 re_sub_abs(a, a, a1);///===
    
44.                 mp_sqr(b, a);/////====-------------------------------/// 平方
    
45.                 re_shr(b, a.used);///-------------------------------/////右移
    
46.                 re_one_mul(b, b, p);
    
47.                 re_sub_abs(t1, t, b);  ////  t1
    
48.                 p <<= 1;
    
49.                 re_mov(a, a1);
    
50.                 re_mov(b, b1);
    
51.                 re_mov(t, t1);
    
52.                 w++;
    
53.                 std::cout << "第" << w << "次计算完成" << std::endl;
    
54.         }////////
    
55.         std::cout << "最后计算π" << std::endl;
    
56.         //pi = ((a1 + b1) * (a1 + b1)) / (4 * t1);
    
57.         add_help(a1, b1, 0);
    
58.         div_one(b1, a1, 2);
    
59.         mp_sqr(a1, b1);/////-------------------------------
    
60.         //re_shr(a1, b1.used);///-------------------------------
    
61.         re_shr(t, t.used / 4);
    
62.         re_Newton(pi, a1, t, pilen / 9 + 1);////////
    
63.         std::cout << "π计算完成" << std::endl;
    
64.         return 0;
    
65. }
    

复制代码

只是呼吸

1）关于初始值：agl的初始值是
a=1;
b=1/sqrt(2);
t=1/4;
P=1;

这里面的b 的初始值（根号2的倒数）就决定了我们最终求的pi的精度。直白点说，如果b的精确值是一万位，那么我们最后得到的pi就是一万位，无论迭代多少次都不会增加一点精度，也不管后面有几万位参加迭代，也只能让pi精确到一万位。


2）关于中间运算的数据膨胀：在agl运算中，有一次高精度的乘法运算，一次高精度的平方运算。那么乘积的最终结果（长度）就会加倍。如果不处理，多次迭代数据迅速膨胀，占用内存。由1）可知，这些数据是没有用的，要在程序中删掉。这就是在agl步骤中有re_shr()的原因，函数 re_shr()是右移，在agl算法中，差不多要裁掉一半的数据，而且不影响精度，能大大加快运算速度。


3）关于平方：在agl算法中，有平方运算。尽量采用平方运算，比如toom3也要实现平方算法，它能加快运算速度，不要用乘法代替。


4）关于开平方算法：我用的开平方算法是在本论坛上，由liangbch
提供的。性能较好。不要用Newton的那个运算，就是用泰勒级数做除法逼近的那个，因为要用到除法，会很慢。另外，开平方运算也要注意中间数据膨胀问题，它会影响我们算pi的时间。


5）关于进制问题：计算机本身就是用二进制运算的，如果采用二进制来计算pi是不是就是最好的呢？从理论上来说，用二进制来计算pi是首选，但是，计算机算出来的pi就是我们看不懂的十六进制字符，需要转换为十进制。这种转换是比较耗时间的，如果把握不好，那就不如用十进制。
我自己就对十六进制转换为十进制没有信心，所以用十进制。


6）关于“高等级乘法”：如果用渐进指数小的乘法，那肯定是首选。我用的toom3乘法，渐进指数只比karatusba好一点点，用来求一百万位pi,有点勉强。
那你可能会问，我怎么不用更好的乘法呢?实际上，是我根本不会，我的能力不够。也许今后会慢慢的涉及到。

只是呼吸

这里我可能搞错了，是不是应该叫agm算法？我写成了agl了。