Spigot 演算法之一 計算調和級數的和，spigot級數

Spigot 演算法之一 計算調和級數的和，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



著作權聲明：本文為博主原創文章，未經博主允許不得轉載。