解密那一段神奇的计算圆周率的程序(二)

这段程序可以更美观些，尽管 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的表示

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

三角形的面积

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

扇形

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

半径

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

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

函数

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

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

微分

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

多项式形式

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

多项式方程

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

系数通项

因此，

反正切展开

下面，进行关键的换元，

换元

得到

反正切公式

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

圆周率公式