我是首先在[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
版权声明：本文为博主原创文章，未经博主允许不得转载。