猫史档案馆

网络上广为流传的一段程序，计算圆周率的程序，看起来是这样的：

#include<stdio.h>

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);
}

把它重新写过：

#include "stdio.h"

int a=10000, bit, cur=2800, sum, rem, fill[2801], q;

void cal_to_next(){
  sum += fill[bit]*a;
  --q;
  fill[bit] = sum % q;
  sum /=q;
  --q;
  --bit;
}

main()   
{   
  while(bit-cur){    
    fill[bit++]=a/5;  
  }

  do{
    sum=0;
    q=cur*2;
    bit = cur;
    cal_to_next();    
     
    while(bit){  
      sum *= bit;
      cal_to_next();
    };
    
    printf("%.4d",rem + sum/a);
    
    rem  = sum % a;
    cur  = cur-14;    
  } while(cur);
return 0;   
}  

这样应该好懂了。如果还是不明白，那么，lisp版本就是这样的：

(defun frac-pi(n)
  (do ((s 0 (+ s v))
       (a 1 (+ a 1))
       (v 2 (* v (/ a b)))
       (b 3 (+ b 2)))
      ((= a n) s)))

如果还是不明白，那么，excel版本应该是这样的：

   

其中，唯一的公式是 D3=D2*C3，向下复制。D列所有数字的和就是要求的圆周率的近似值。

学数学的人是这样写的：

 

（虽然数学符号很混乱的说，但果然看起来紧凑，比lisp看起来还要紧凑。）
这么多解释了，你一定能看懂这个公式了。

 

这段程序可以更美观些，尽管 c 语言写的程序在美观程度上永远比不上 lisp 的。

#include "stdio.h"
int i,j;
long result,ret=0,tmp,tab[8401];
void main() {
  for(i=0;i<8401;i++) tab[i] = 2000;
  for(j=8400/14;j>0;j--){
    result = 0;
    for (i=j*14;i>0;i--){
      tmp = ( result*i + tab[i]*10000);
      result = tmp/(2*i-1);
      tab[i] = tmp%(2*i-1);
    }
    printf("%.4d",ret+result/10000);
    ret = result%10000;
  }
}

使用的公式其实就是：

 

那么这个公式是怎样推导/证明得来的呢？
这个公式不是普通的公式，推导起来相当复杂。因为这个公式是大数学家欧拉无聊的时候推导出来的。（姓欧的人都很厉害，用归纳法可以证明：欧阳修，欧阳询，欧阳锋，欧阳克，欧阳中石，欧几里德，欧拉......）

欧拉推导的时候，可以用纯粹的公式，也可用纯粹的图形，也可以什么都不用，直接写结果。
为了方便学习和记忆，我结合几何图形来推导一遍。

 

设有这样一个直角三角形ABC，三边的长度如图所标。则角B为直角。
那么，角A可以有多种表示方法，可以用反正弦，也可以用反正切。

 

这个三角形的面积可以写成：

 

假设要在角A上做一个扇形，使得扇形的面积正好等于这个三角形。很显然这个扇形的半径小于1,大于边AB。因此，可以在AC线段上取一个分点Y，以AY为半径做圆。

 

那么这个扇形的半径应该是多少呢？
假设为R,那么

 

A就是上面的角A，弧度表示。

为了方便运算，设R的平方的倒数为y，则

 

现在就考察这个函数。
在牛顿之前，人们就喜欢把函数展开成多项式的形式。因为多项式方便运算，在允许的误差内，还可以随时截断。那么，这个函数该如何展开呢？

展开一个函数，通常都会用到导数。那么，把这个函数求导，可以得到：

 

假设已经得到y的多项式展开，那么，形式上应该是：

 

把这个形式，连同它的导数，代入上面的微分方程，经过整理以后，就是：

 

通过比较方程两边的系数，可以得到：

 

因此，

 

下面，进行关键的换元，

 

得到

 

令t=1，就可以得到想要的结果了：

 

注： 本文并非原创，原文见链接： https://www.jianshu.com/p/3ef2cd036780 https://www.jianshu.com/p/8b70d73cde80