演算精度の話

大先輩のブログでライプニッツ級数で円周率の計算って話が出ていました。もともとはプログラム電卓でどこまでやれるのかということが主眼なので、でかいPCを持ち出すのは野暮の極みといったところですが、ついつい悪乗り。

gccでもVer.4.8以降だと4倍精度小数点つまり128bit floatが使用できるとのことでしたので、さっそくテスト。
ちなみにソース(paicalc.c)はこんな感じ。面倒なので定義通りの記述。なんかねぇ…

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

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

コンパイルして実行した結果は、
$ gcc -o paicalc.o paicalc.c -lquadmath
$ ./paicalc.o
 loop=100000000 result=3.1415926485897931884626429145303967199477

円周率の頭の方の分かっている数値は3.141592653589793…なので4倍精度だからといって7桁程度しか出てませんが、式から言えば当たり前ですね。
それとこのquadmathライブラリは当たり前ですが、ソフト処理なのでしょうから通常の64bit演算に比べると大幅に時間がかかります。倍精度の場合は0.76秒で済みましたが、4倍精度の場合は19.4秒と25倍以上も時間がかかっています。

オイラー式だとべき乗級数なのでもうちょっと収束が早いかな。ただ結果が円周率の2乗なので、開平演算の誤差が乗ってきそう… まぁ明日試してみますか。