割圆术求圆周率程序

忽然对数学起了兴趣。这个圆周率，神秘的圆周率到底怎么割出来的？
以我们现在的知识和条件，解它应该是小菜了。不过还是可以试试以探风景。
关键是找到递推公式，已知边长，求得下一次的边长。
这个好办。 图我不会画，就靠文字表达了。
根据勾股定理(平面几何都受勾股定理的约束,好在我用平面几何会证明勾股定理)，
根据当前弦长，先求得当前弦高，再求得弧高，再求得新弦长.
来一段c语言:
hjj@hjj-Inspiron:~/MyTest/Test$ cat hello1.cpp
#include <stdio.h>
#include <math.h>
double split(double a);
int main(int argc, char *argv[])
{
    double a=1.0;
    int n=6;
    double s1 = a * n;
    double b=0.0;
    double s2 = b * n * 2;
    int count = 0;
    while(count<20 && s1 != s2)
    {
        b=split(a);
        s1 = a * n;
        s2 = b * 2*n;
        n *=2;
        a = b;
        count ++;
        printf("<count:%d> s1:%18.16lf, s2:%18.16lf,pi:%18.16lf,a:%18.16lf\n",count,s1,s2,s2/2,a);
    }
    printf("pi:%18.16lf\n",s1/2);
    return 0;
}

double split(double a)
{
    // 如此计算误差甚大
    return sqrt(2-2*sqrt(1-(a/2)*(a/2)));
}
看看结果:
~/MyTest/Test$ cat txt1
<count:1> s1:6.0000000000000000, s2:6.2116570824604995,pi:3.1058285412302498,a:0.5176380902050416
<count:2> s1:6.2116570824604995, s2:6.2652572265624737,pi:3.1326286132812369,a:0.2610523844401031
<count:3> s1:6.2652572265624737, s2:6.2787004060937441,pi:3.1393502030468721,a:0.1308062584602863
<count:4> s1:6.2787004060937441, s2:6.2820639017810596,pi:3.1410319508905298,a:0.0654381656435527
<count:5> s1:6.2820639017810596, s2:6.2829049445706886,pi:3.1414524722853443,a:0.0327234632529723
<count:6> s1:6.2829049445706886, s2:6.2831152158232442,pi:3.1415576079116221,a:0.0163622792078730
<count:7> s1:6.2831152158232442, s2:6.2831677842978717,pi:3.1415838921489359,a:0.0081812080524712
<count:8> s1:6.2831677842978717, s2:6.2831809264735234,pi:3.1415904632367617,a:0.0040906125823395
<count:9> s1:6.2831809264735234, s2:6.2831842120860966,pi:3.1415921060430483,a:0.0020453073607051
<count:10> s1:6.2831842120860966, s2:6.2831850331763093,pi:3.1415925165881546,a:0.0010226538139935
<count:11> s1:6.2831850331763093, s2:6.2831852372815788,pi:3.1415926186407894,a:0.0005113269236069
<count:12> s1:6.2831852372815788, s2:6.2831852906424315,pi:3.1415926453212157,a:0.0002556634639747
<count:13> s1:6.2831852906424315, s2:6.2831852906424315,pi:3.1415926453212157,a:0.0001278317319874

靠，误差好大，才做了13次循环，就做到头了，而且结果也不准确，我希望精度是15-16位，可是结果才7位，
3.1415926 - 3.1415927 之间，后面就不对了。
计算误差也能害死人，换种方法, 打印其计算过程。

include <stdio.h>
#include <math.h>
void split(int count,double a1, int n,double &a2);
int main(int argc, char *argv[])
{
    // 初始值边长1，边数6
    double a1=1.0,a2=0.0;
    int n=6;
    int count = 0;
    printf("%10s,%22s,%22s,%22s,%22s,%22s\n","<次数>", "折半弦长", "弦高","弧高","新弦长","pi值");
    while(count<100)
    {
        split(count,a1,n,a2);
        count++;
        n*=2;
        a1=a2;
    }
    return 0;
}

void split(int count,double a1, int n,double &a2)
{
    double b,c,d,S2,pi;
    b = a1/2;            // 折半弦长
    c = sqrt( 1 - b * b);  // 弦高
    d = 1-c;                // 弧高
    a2 = sqrt(b*b + d*d);    // 新弦长
    S2 = a2 * 2 * n ;    // 周长
    pi = S2 / 2;        // pi 值
    printf("<count:%2d> b:%18.16lf, c:%18.16lf,d:%18.16lf,a2:%18.16lf,pi:%18.16lf\n",count,b,c,d,a2,pi);
}
结果：
~/MyTest/Test$ ./hello2
  <次数>,          折半弦长,                弦高,                弧高,             新弦长,                 pi值
<count: 0> b:0.5000000000000000, c:0.8660254037844386,d:0.1339745962155614,a2:0.5176380902050415,pi:3.1058285412302489
<count: 1> b:0.2588190451025207, c:0.9659258262890683,d:0.0340741737109317,a2:0.2610523844401031,pi:3.1326286132812378
<count: 2> b:0.1305261922200516, c:0.9914448613738104,d:0.0085551386261896,a2:0.1308062584602861,pi:3.1393502030468667
<count: 3> b:0.0654031292301431, c:0.9978589232386035,d:0.0021410767613965,a2:0.0654381656435523,pi:3.1410319508905098
<count: 4> b:0.0327190828217761, c:0.9994645874763657,d:0.0005354125236343,a2:0.0327234632529736,pi:3.1414524722854624
<count: 5> b:0.0163617316264868, c:0.9998661379095618,d:0.0001338620904382,a2:0.0163622792078743,pi:3.1415576079118579
<count: 6> b:0.0081811396039371, c:0.9999665339174011,d:0.0000334660825989,a2:0.0081812080524696,pi:3.1415838921483186
<count: 7> b:0.0040906040262348, c:0.9999916334443506,d:0.0000083665556494,a2:0.0040906125823282,pi:3.1415904632280505
<count: 8> b:0.0020453062911641, c:0.9999979083589001,d:0.0000020916410999,a2:0.0020453073606766,pi:3.1415921059992717
<count: 9> b:0.0010226536803383, c:0.9999994770895884,d:0.0000005229104116,a2:0.0010226538140274,pi:3.1415925166921577
<count:10> b:0.0005113269070137, c:0.9999998692723886,d:0.0000001307276114,a2:0.0005113269237248,pi:3.1415926193653840
<count:11> b:0.0002556634618624, c:0.9999999673180966,d:0.0000000326819034,a2:0.0002556634639513,pi:3.1415926450336911
<count:12> b:0.0001278317319757, c:0.9999999918295241,d:0.0000000081704759,a2:0.0001278317322368,pi:3.1415926514507682
<count:13> b:0.0000639158661184, c:0.9999999979573810,d:0.0000000020426190,a2:0.0000639158661510,pi:3.1415926530550373
<count:14> b:0.0000319579330755, c:0.9999999994893453,d:0.0000000005106547,a2:0.0000319579330796,pi:3.1415926534561045
<count:15> b:0.0000159789665398, c:0.9999999998723362,d:0.0000000001276638,a2:0.0000159789665403,pi:3.1415926535563719
<count:16> b:0.0000079894832702, c:0.9999999999680841,d:0.0000000000319159,a2:0.0000079894832702,pi:3.1415926535814380
<count:17> b:0.0000039947416351, c:0.9999999999920209,d:0.0000000000079791,a2:0.0000039947416351,pi:3.1415926535877046
<count:18> b:0.0000019973708176, c:0.9999999999980053,d:0.0000000000019947,a2:0.0000019973708176,pi:3.1415926535892713
<count:19> b:0.0000009986854088, c:0.9999999999995013,d:0.0000000000004987,a2:0.0000009986854088,pi:3.1415926535896630
<count:20> b:0.0000004993427044, c:0.9999999999998753,d:0.0000000000001247,a2:0.0000004993427044,pi:3.1415926535897611
<count:21> b:0.0000002496713522, c:0.9999999999999688,d:0.0000000000000312,a2:0.0000002496713522,pi:3.1415926535897856
<count:22> b:0.0000001248356761, c:0.9999999999999922,d:0.0000000000000078,a2:0.0000001248356761,pi:3.1415926535897913
<count:23> b:0.0000000624178380, c:0.9999999999999980,d:0.0000000000000020,a2:0.0000000624178380,pi:3.1415926535897931
<count:24> b:0.0000000312089190, c:0.9999999999999994,d:0.0000000000000006,a2:0.0000000312089190,pi:3.1415926535897936
<count:25> b:0.0000000156044595, c:0.9999999999999999,d:0.0000000000000001,a2:0.0000000156044595,pi:3.1415926535897936
<count:26> b:0.0000000078022298, c:0.9999999999999999,d:0.0000000000000001,a2:0.0000000078022298,pi:3.1415926535897936
<count:27> b:0.0000000039011149, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000039011149,pi:3.1415926535897936
<count:28> b:0.0000000019505574, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000019505574,pi:3.1415926535897936
<count:29> b:0.0000000009752787, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000009752787,pi:-1.0471975511965979
<count:30> b:0.0000000004876394, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000004876394,pi:-1.0471975511965979
<count:31> b:0.0000000002438197, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000002438197,pi:0.0000000000000000
<count:32> b:0.0000000001219098, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000001219098,pi:0.0000000000000000
<count:33> b:0.0000000000609549, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000000609549,pi:0.0000000000000000
<count:34> b:0.0000000000304775, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000000304775,pi:0.0000000000000000
<count:35> b:0.0000000000152387, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000000152387,pi:0.0000000000000000
<count:36> b:0.0000000000076194, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000000076194,pi:0.0000000000000000
<count:37> b:0.0000000000038097, c:1.0000000000000000,d:0.0000000000000000,a2:0.0000000000038097,pi:0.0000000000000000

最后精度为3.141592653589793是正确的，最后一位不准确。 这次是对的，其中弧高是0.0000000000000001, 不能比这更精确了。
而且计算也没有引入过大误差。

进一步想看看第一次计算误差到底是什么样的，下面是我的比较，看出第一次计算方法，误差是忽大忽小的，至于为什么引入如此误差，
就没有再继续研究了，我只是一个业余爱好者。

排序为上面的大，下面的小，以第二次正确值为ref 参考值
5176380902050416 第一次:
5176380902050415 ref:

2610523844401031 第一次:
2610523844401031 ref:

1308062584602863 第一次:
1308062584602861 ref:

0654381656435527 第一次:
0654381656435523 ref:

0327234632529736 ref:
0327234632529723 第一次:

0163622792078743 ref:
0163622792078730 第一次:

0081812080524712 第一次:
0081812080524696 ref:

0040906125823395 第一次:
0040906125823282 ref:

0020453073607051 第一次:
0020453073606766 ref:

0010226538140274 ref:
0010226538139935 第一次:

0005113269237248 ref:
0005113269236069 第一次:

0002556634639747 第一次:
0002556634639513 ref:

0001278317322368 ref:
0001278317319874 第一次:

总之，凡事需得研究，才能明白，纵使计算精度很高，不当的计算方法，也可能达不到目的。