演算精度その2

単純に整数の逆数を足し込んでいくだけのライプニッツ級数は収束が遅いことがわかりました。なので、べき乗項の有るオイラー級数だったらどうなのかをさっそくテスト。大して計算式に変化はありませんのでちょいちょいと変更してコンパイル。

#include <stdio.h>
#include <quadmath.h>
#include <math.h>

void main(void)
{
    char    c1[80];
    long    loop;
    __float128  res=0.0Q,r1=1.0Q;
    
    for (loop=1;loop<100000000;loop++) {
        res=res+1.0Q/(r1*r1);
        r1=r1+2.0Q;
    };
    res=sqrtq(res*8.0Q);
    quadmath_snprintf(c1,sizeof(c1),"%.40Qf",res);
    printf("loop=%ld\tresult=%s\n",loop,c1);
}

結果は、
$ gcc -o paicalc2.o paicalc2.c -lquadmath
$ ./paicalc2.o
 loop=10000000 result=3.1415926217588012757267165560912027187844
$ ./paicalc2.o
 loop=100000000 result=3.1415926504066943431811710021935384220221

と言う訳でライプニッツ級数の1億項の計算に比べて1000万項の計算で1桁精度が向上しました。項数を同じにした場合は2桁の向上ですね。でもそんなものか…

実際にはその後いろんな方法が考案されていて、そちらを使えばってのが身も蓋もない言い方。実際ガウス・ルジャンドル法だと3回目には18桁正しい値が得られるし、昨今の多数桁計算にはマーチンの式が使用されている。
とは言え、今から300年前に手計算だけでこれにたどり着いたのだからすごいことだと思うな。