我是首先在[1] 注意到 Spigot-Algorithm的，这个算法公布的相当早。见[2].  [1] 给出几个令人惊异的程序。仅仅用非常少的代码就能够计算e,pi,log(2)等常数。

当中那个4行代码计算圆周率的程序被网友称作外星人写的程序。但我一直没有勇气去分析和学习它，近期最终决定学习这个 Spigot-Algorithm，先看了文献【3】。明确了其基本思想。遂计划尝试编写各种计算常数的代码。并写一个系列博客。从这篇開始，我将讲述假设使用这个算法计算各种常数或者级数的和。

    数列a[n]={1/1,1/2,1/3,1/4... 1/n} 被称作调和数列。调和数列的和f(n)= 1/1 + 1/2 + 1/3 +1/4 + ... 1/n 被称作调和级数。调和级数是最简单的级数之中的一个。

用Spigot 算法来计算这个级数也最为简单。

   以下给出代码。

详细说明以后补上。

#define R 10		//进制。可改为100。1000，10000。#define FMT_STR "%d"	//当R=100,1000,10000时。对应的，须要改为"%02d","%03d","%04d",#define N 10		//计算交错级数的前10项#define P 20		//当R=10^k时，可打印前P*k位有效数字int a[N+1],i,j,x;void main(){  for (i=N;i;a[i--]=1);  for (j=0;j

 几点说明：

  1. 数组的长度和计算的项数有关，而和终于精度无关。

  2. a[1] to a[n]存储第j轮计算时。各项的分子。

  3. 每轮计算中。得到2位10进制数。最高位直接输出，次高位缓存在a[0]

  4. 当N大于10时。计算结果错误。这是由于交错级数收敛非常慢. 当得到次高位时就急于输出最高位仍是冒险的方法。

改进的方法是当得到第5位时再输出最高位，然后将第2到5位存储在a[0],以下是改动后的代码，可正确计算交错级数前100项的值。

#define R 10  //进制，可改为100。1000，10000，#define FMT_STR "%d" //当R=100,1000,10000时，对应的，须要改为"%02d","%03d","%04d",#define N 100  //计算交错级数的前N项#define P 100  //当R=10^k时，可打印前P*k位有效数字int a[N+1],i,j,x;void main(){  for (i=N;i;i--)  a[i]=1;  for (j=0;j
     
      2)    {      printf(FMT_STR,a[0]/(R*R*R*R));      a[0]%=(R*R*R*R);      if ( j==3)        printf(".");    }  }}
     

參考文献：

1. Tiny programs for constants ,

1. M. Abramowitz and I. Stegun, Handbook of Mathematical Functions, Dover, New York, (1964)

3.