tan (pi/4) = 1

pi/4 = arctan 1

pi = 4arctan 1

积分其导数并固定在一点上的值给出反三角函数作为定积分的表达式:

当 x 等于 1 时，在有极限的域上的积分是瑕积分，但仍是良好定义的。

一个看起来很牛的计算Pi的程序: int a=10000,b,c=2800,d,e,f[2801],g;   main() { for(;b-c;)     f[b++]=a/5; for(;d=0,g=c*2;c -=14,printf("%.4d",e+d/a),e=d%a)       for(b=c; d+=f[b]*a,f[b]=d%--g,d/=g--,--b; d*=b); } 

这还不是最牛的，源码所在的这个网站里牛程序多的是：

http://numbers.computation.free.fr/Constants/TinyPrograms/tinycodes.html

下面是用C写的Pi计算程序：

The following tiny C code computes digits of p. Boris Gourevitch kindly sent me the information that this program is from Dik T. Winter (cwi institute, Holland).

The program has 158 characters long. Finding the algorithm used is not so easy. Sebastian Wedeniws worked hard to find a shorter program (142 characters) :

Jean-Charles Meyrignac pointed me out a tiny code from the page http://www.isr.umd.edu/~jasonp/pigjerr. The program is by Gjerrit Meinsma, of 141 characters long, and computes 1000 digits of pi.

二、矩阵拟合法计算Pi - 可用OpenMP并行

Sequetial Version:  #include   #include   void main ()  {  time_t start, finish;  static long num_steps = 1000000000;  double step;  int i;  double x, pi, sum = 0.0;  step = 1.0/(double) num_steps;  start = clock();  for (i=0;i < num_steps; i++)  {  x = (i+0.5)*step;  sum = sum + 4.0/(1.0+x*x);  }  pi = step * sum;  finish = clock();  printf( "Pi = %16.15f (%d steps), %ld ms\n ", pi, num_steps, finish-start );  return;  }

--------------------------------------------------------------------------------------------

Parallel Version:  #include   #include   #include   void main ()  {  time_t start, finish;  static long num_steps = 1000000000;  double step;  int i;  double x, pi, sum = 0.0;  step = 1.0/(double) num_steps;  start = clock();  #pragma omp parallel for reduction(+:sum) private(x) /*只加了这一句，其他不变*/  for (i=0;i < num_steps; i++)  {  x = (i+0.5)*step;  sum = sum + 4.0/(1.0+x*x);  }  pi = step * sum;  finish = clock();  printf( "Pi = %16.15f (%d steps), %ld ms\n ", pi, num_steps, finish-start );  return;  }  result:  Sequential version:  Pi = 3.141592653589792 (1000000000 steps), 13900000 ms  Parallel version for 8 threads:  Pi = 3.141592653589794 (1000000000 steps), 1820000 ms

从结果可以看到8线程的speedup=7.64，接近线性。因为这个程序本身具有良好的并发性，循环间几乎没有数据依赖，除了sum，但是用reduction(+:sum)把对于sum的相关也消除了。而且实验平台本身就有8个处理器核。

三、概率方法

Name: 蒙特卡罗法求解 pi    Description:         边长为r的正方形，内接1/4圆周，半径也为r.         正方形面积s1 = r*r         1/4圆周面积s2 = pi*r*r / 4          假设在正方形中随机选择许多点，那么这些点中有一部分也落在1/4圆周内。         如果所选择的点是均匀分布的，那么我们可以期望落于1/4圆周内的点的概率是：         f = s2 / s1 = pi / 4         因此，通过测定f , 就可以计算出 pi = 4*f . 

程序如下：