(转)一个让人发狂的PI求解C程序

记得是大四的时候,我第一次看到了这个程序。那段时间，我对各种稀奇的C程序很感兴趣，也写了很多程序，但是当时看到这个程序时，着实是吓了一大跳。
今天无意着又想起这个程序，所以上网找了找。觉得找到了一个不错文章解析这个程序的。现在转载过来与大家再次分析。

在网上也找到个另类变态的类似程序，但是我没有执行成功。不过还是把地址贴出来与大家分享下，还请高手指教如何运行。：）

原文地址：

http://blogger.org.cn/blog/more.asp?name=njucs&id=10151

http://blog.est.im/archives/877

详见三楼的回复。
1988年IOCCC的Best layout奖，作者Merlyn LeRoy (Brian Westley).

http://ioccc.org/years.html#1988_westley

注：以上两个地址中，说的都是同一个程序。

本文转载自以下地址:
http://hi.baidu.com/jumbo/blog/item/155cb01b13cfba198618bf89.html

#include <stdio.h>
long 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);
}

简短的4行代码，就可以精确计算机出800位的PI（圆周率）值。
实在太震撼人心了。这样的程序也能运行，竟然还能完成这样让人难以置信的任务，真是太神了。

这是某一年The International Obfuscated C Code Contest（国际模糊C代码大赛）上的获奖作品（努力了，但是没有找到一个确切的时间）。这是属于C大师的盛会，因为这是一件极具挑战的活儿。

还有一个更酷的：

#define _ F-->00 || F-OO--;
long F=00,OO=00;
main(){F_OO();printf("%1.3f\n", 4.*-F/OO/OO);}F_OO()
{
             _-_-_
       _-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
_-_-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
 _-_-_-_-_-_-_-_-_-_-_-_-_-_-_
  _-_-_-_-_-_-_-_-_-_-_-_-_-_
    _-_-_-_-_-_-_-_-_-_-_-_
       _-_-_-_-_-_-_-_-_
            _-_-_-_
}

一、源程序
本文分析下面这个很流行的计算 PI 的小程序。下面这个程序初看起来似乎摸不到头脑，不过不用担心，当你读完本文的时候就能够基本读懂它了。
程序一：很牛的计算 PI 的程序

#include <stdio.h>
long 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);
}

二、数学公式
数学家们研究了数不清的方法来计算 PI ,这个程序所用的公式如下：

1
2
3

          1          2          3                    k
pi = 2 + --- * (2 + --- * (2 + --- * (2 + ...  (2 + --- * (2 + ... ))...)))
          3          5          7                   2k+1

至于这个公式为什么能够计算出 PI，已经超出了本文的能力范围。
下面要做的事情就是要分析清楚程序是如何实现这个公式的。
我们先来验证一下这个公式：
程序二：PI 公式验证程序

#include "stdio.h"
void main()
{
   float pi=2;
   int  i;
   for(i=100;i>=1;i--) pi=pi*(float)i/(2*i+1)+2;
   printf("%f\n",pi);
   getchar();
}

上面这个程序的结果是3.141593。

三、程序展开
在正式分析程序之前，我们需要对程序一进行一下展开。我们可以看出程序一都是使用 for 循环来完成计算的，这样做虽然可以使得程序短小，但是却很难读懂。根据for循环的运行顺序，我们可以把它展开为如下while循环的程序：
程序三：for转换为while之后的程序

#include <stdio.h>
int a=10000,b,c=2800,d,e,f[2801],g;
main() {
  int i;
  for(i=0;i<c;i++) f[i]=a/5;
  while(c!=0)
  {
    d=0;
    g=c*2;
    b=c;
    while(1)
    {
      d=d+f[b]*a;
      g--;
      f[b]=d%g;
      d=d/g;
      g--;
      b--;
      if(b==0) break;
      d=d*b;
    }
    c=c-14;
    printf("%.4d",e+d/a);
    e=d%a;
  }
}

注：for([1];[2];[3]) {[4];}的运行顺序是[1],[2],[4],[3]。如果有逗号操作符，例如：d=0,g=c*2，则先运行 d=0, 然后运行g=c*2,并且最终的结果是最后一个表达式的值，也就是这里的 c*2。

下面我们就针对展开后的程序来分析。

四、程序分析
要想计算出无限精度的PI，我们需要上述的迭代公式运行无数次，并且其中每个分数也是完全精确的，这在计算机中自然是无法实现的。那么基本实现思想就是迭代足够多次，并且每个分数也足够精确，这样就能够计算出PI的前n位来。上面这个程序计算800位，迭代公式一共迭代2800次。

1	int a=10000,b,c=2800,d,e,f[2801],g;

这句话中的2800就是迭代次数。

由于float或者double的精度远远不够，因此程序中使用整数类型（实际是长整型），分段运算（每次计算4位）。我们可以看到输出语句 printf("%.4d",e+d/a); 其中%.4就是把计算出来的4位输出，我们看到c每次减少14(c=c-14;)，而c的初始大小为2800，因此一共就分了200段运算，并且每次输出4位，所以一共输出了800位。

由于使用整型数运算，因此有必要乘上一个系数，在这个程序中系数为1000，也就是说，公式如下：

1
2
3

          1          2          3                    k
pi = 2 + --- * (2 + --- * (2 + --- * (2 + ...  (2 + --- * (2 + ... ))...)))
          3          5          7                   2k+1

这里的2k表示2000，也就是f[2801]数组初始化以后的数据,a=10000,a/5=2000,所以下面的程序把f中的每个元素都赋值为2000：

1	for(i=0;i<c;i++) f[i]=a/5;

你可能会觉得奇怪，为什么这里要把一个常数储存到数组中去，请继续往下看。
我们先来跟踪一下程序的运行：

while(c!=0) // 假设这是第一次运行，c=2800，为迭代次数
{
  d=0;
  g=c*2;  // 这里的g是用来做k/(2k+1)中的分子
  b=c;    // 这里的b是用来做k/(2k+1)中的分子
  while(1)
  {
    d=d+f[b]*a; // f中的所有的值都为2000，这里在计算时又把系数扩大了a=10000倍。这样做的目的稍候介绍，你可以看到输出的时候是d/a,所以这不影响计算
    g--;
    f[b]=d%g; // 先不管这一行
    d=d/g; // 第一次运行的g为2*2799+1，你可以看到g做了分母
    g--;
    b--;
    if(b==0) break;
    d=d*b; // 这里的b为2799，可以看到b做了分子。
  }
  c=c-14;
  printf("%.4d",e+d/a);
  e=d%a;
}

只需要粗略的看看上面的程序，我们就大概知道它的确是使用的那个迭代公式来计算 PI 的了，不过不知道到现在为止你是否明白了f数组的用处。如果没有明白，请继续阅读。

d=d/g,这一行的目的是除以2k+1，我们知道之所以程序无法精确计算的原因就是这个除法。即使用浮点数，答案也是不够精确的，因此直接用来计算800位的 PI 是不可能的。那么不精确的成分在哪里？很明显：就是那个余数 d%g。

程序用 f 数组把这个误差储存起来，再下次计算的时候使用。现在你也应该知道为什么 d=d+f[b]*a; 中间需要乘上 a 了吧。把分子扩大之后，才好把误差精确的算出来。

d 如果不乘 10000 这个系数，则其值为 2000，那么运行 d=d/g; 则是 2000/(2*2799+1)，这种整数的除法答案为 0，根本无法迭代下去了。

现在我们知道程序就是把余数储存起来，作为下次迭代的时候的参数，那么为什么这么做就可以使得下次迭代出来的结果为接下来的4位数呢？

这实际上和我们在纸上作除法很类似：

　　　0142
　　/———
　7/　1
　　　10
　　　 7
    ——————
　　　 30
　　　 28
    ——————
　　　　20
　　　　14
    ——————
　　　　 6
.....

我们可以发现，在做除法的时候，我们通常把余数扩大之后再来计算，f中既然储存的是余数，而f[b]*a;则正好把这个余数扩大了a倍，然后如此循环下去，可以计算到任意精度。

这里要说明的是，事实上每次计算出来的d并不一定只有4位数，例如第一次计算的时候，d的值为31415926，输出4位时候，把低四位的值储存在e中间，e=d%a，也就是5926。

最后，这个c=c-14不太好理解。事实上没有这条语句，程序计算出来的仍然正确。只是
因为如果迭代2800次，无论分数如何精确，最后Pi的精度只能够达到800。
你可以把程序改为如下形式尝试一下：

for(i=0;i<800;i++)
{
     {
         d=0;
         g=c*2;
         b=c;
         while(1)
            {
                d=d+f[b]*a;
                g--;
                f[b]=d%g;
                d=d/g;
                g--;
                b--;
                if(b==0) break;
                d=d*b;
    }
    // c=c-14; 不要这句话。
    printf("%.4d",e+d/a);
    e=d%a;
}