美文网首页编程语言爱好者
怎样计算圆周率到10000位以后

怎样计算圆周率到10000位以后

作者: aubell | 来源:发表于2020-04-02 09:08 被阅读0次

业余爱好之一,就是计算圆周率。
我把步骤记录下来,以后好重复,或者改进。

步骤一:

下载并安装GNU的GMP库,这个库提供了高精度的浮点数,精度之高,超出你的想象。

步骤二:

下载并安装GMP库的文档,用info命令查看,了解几个基本函数,初始化,赋值,开方,加,减,乘,除。

步骤三:

找一个计算pi的迭代算法,如AGM算法。


AGM.png

步骤四:

参照GMP文档和AGM算法编程,代码如下

#include <stdio.h>
#include <gmp.h>

#define P 65536

int main()
{
  int it_times = 20;
  int i;
  mpf_t a,b,c,x,y,pi,temp1,temp2,temp3;

  mpf_set_default_prec(P);
  
  mpf_inits(a,b,c,x,y,temp1,temp2,temp3,pi,NULL);

  mpf_set_ui(a,1);   //a:=1
  mpf_set_ui(x,1);   //x:=1

  mpf_sqrt_ui(temp1,2);
  
  mpf_ui_div(b,1,temp1); //set b:=1/sqrt(2)

  mpf_set_ui(temp2,4);
  mpf_ui_div(c,1,temp2); // set c:= 1/4
  
  //LOOP NOW
  for(i=0; i<it_times; i++){
    
    mpf_set(y,a);   // y = a
    
    mpf_add(temp1,a,b);
    mpf_div_ui(a,temp1,2); // a = (a+b)/2   
    
    mpf_mul(temp1,b,y);
    mpf_sqrt(b,temp1);  // b = sqrt(b*y);

    mpf_sub(temp1,a,y);
    mpf_mul(temp2,temp1,temp1); 
    mpf_mul(temp3,x,temp2);     

    mpf_sub(temp1,c,temp3);   // c = c-x(a-y)^2
    mpf_set(c,temp1);

    mpf_mul_ui(temp1,x,2);
    mpf_set(x,temp1);       //x = x*2
  }  
  
  mpf_add(temp1,a,b);
  mpf_mul(temp2,temp1,temp1);  // (a+b)^2
  mpf_mul_ui(temp3,c,4);       // 4c  
    
  mpf_div(pi,temp2,temp3);    // div to get pi
  printf("π = \n");
  mpf_out_str(NULL,10,0,pi);
  printf("\n");

  mpf_clears(a,b,c,x,y,temp1,temp2,temp3,pi,NULL);
  return 1;  
}

//gcc -lgmp mypi.c to compile
//search GMP methord

步骤五 编译运行

在终端输入:
gcc -lgmp mypi.c -o mypi
以编译,

输入
./mypi以运行

得到结果pi,结果太长了,就不贴图了。

相关文章

网友评论

    本文标题:怎样计算圆周率到10000位以后

    本文链接:https://www.haomeiwen.com/subject/odetphtx.html