Spigot 算法之一 计算调和级数的和



   我是首先在[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<P;j++)
  {
    x=0;
    for (i=N;i;i--)
    {
	x+=(a[i]*R)/i; a[i]=(a[i]*R)%i;
    }
    x+=a[i]*R;
    printf(FMT_STR,x/R); a[0]=x/R;
    if ( j==0)
     printf(".");
  }
}

 几点说明：
  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<P+1;j++)
  {
    x=0;
    for (i=N;i;i--)
    {
      x+=(a[i]*R)/i; a[i]=(a[i]*R)%i;
    }
    a[0]= a[0]*R + x;
    if ( 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 ,http://numbers.computation.free.fr/Constants/constants.html

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

3. https://en.wikipedia.org/wiki/Spigot_algorithm