gccの4倍精度ライブラリー

gcc 4.6 から4倍精度計算ができる様になりました．昔は4倍精度で計算できるマシンがあったみたいだけどね．

今回は、c, (c++はやってない), fortranの言語で実装するとしてどの言語が最も実装しやすいか比較する．
gcc のversionは4.8 を用いる．結論としてはfortranが一番楽だと思う．
cおよびc++は
http://homepage3.nifty.com/maikaze/float128.html
で解説されているようにc(, c++)では幾分実装が不自然な所もあるので，一応比較してみる．
計算する内容は1/3, 4*arctan(1)=pi,sqrt(2)である．

まずfortranから

program fort_quadmath
 implicit none
 real (kind(0D0)) dx, dy, dpi,dsqrt2
 real (kind(0Q0)) qx, qy, qpi,qsqrt2
 dx = 1.0D0;
 dy = 3.0D0;
 dsqrt2 = sqrt(2.0D0)
 dpi = 4.0D0*atan(1.0D0)
 qx = 1.0Q0;
 qy = 3.0Q0;
 qsqrt2 = sqrt(2.0Q0)
 qpi = 4.0Q0*atan(1.0Q0)
 print *, "---- double precision ----"
 print *, "1/3     =", dx/dy
 print *, "pi      =", dpi
 print *, "sqrt(2) =", dsqrt2
 print *, "----  quad  precision ----" 
 print *, "1/3     =", qx/qy
 print *, "pi      =", qpi
 print *, "sqrt(2) =", qsqrt2
end program fort_quadmath 

コンパイルは

 gfortran -o fort_quadmath fort_quadmath.f03 

であり，実行結果は

$ ./fort_quadmath 
 ---- double precision ----
 1/3     =  0.33333333333333331     
 pi      =   3.1415926535897931     
 sqrt(2) =   1.4142135623730951     
 ----  quad  precision ----
 1/3     =  0.333333333333333333333333333333333317      
 pi      =   3.14159265358979323846264338327950280      
 sqrt(2) =   1.41421356237309504880168872420969798   

 pi           3.1415926535897932384626433832795028841971693993751
 sqrt(2)      1.4142135623730950488016887242096980785696718753769

厳密解と比較して32桁くらいまでは遜色無く一致している．
実装に関して変数の型の宣言と、リテラル(0Q0)さえ気をつけていれば，関数の型を意識する必要がない(というか関数が代表名で与えられているので，変数の型によって関数が勝手に陰的に変わるんだと理解している．)

一方c言語に関して

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

int main(void)
{
 char c1[80],c2[80],c3[80],c4[80],c5[80];
 double dx,dy,dpi;
 long double lx,ly,lpi;
 __float128 x,y, pi1,pi2,sqrt1,sqrt2;
 dx = 1.0;
 dy = 3.0;
 dpi = 4.0*atan(1.0);
 lx = 1.0L;
 ly = 3.0L;
 lpi = 4.0L*atan(1.0L);
 x = 1.0Q;
 y = 3.0Q;
 pi1 = 4.0Q*atan(1.0Q);
 pi2 = 4.0Q*atanq(1.0Q); 
 sqrt1 = sqrt(2.0Q);
 sqrt2 = sqrtq(2.0Q);  
 
 printf("---- double precision ----\n");
 printf("1/3 = %.30lf\n", dx/dy);
 printf("pi  = %lf\n", dpi); 
 printf("---- long double prescision ----\n");
 printf("1/3 = %.30Lf\n", lx/ly);
 printf("pi  = %.30Lf\n", lpi);
 printf("---- quad precision ----\n");
 quadmath_snprintf(c1,80,"%.40Qf", x/y);
 quadmath_snprintf(c2,80,"%.40Qf", pi1); 
 quadmath_snprintf(c3,80,"%.40Qf", pi2);  
 quadmath_snprintf(c4,80,"%.40Qf", sqrt1);  
 quadmath_snprintf(c5,80,"%.40Qf", sqrt2);    
 printf("1/3 = %s\n", c1);
 printf("pi1  = %s\n", c2);
 printf("pi2  = %s\n",c3);
 printf("sqrt1  = %s\n",c4);
 printf("sqrt2  = %s\n",c5);  
 return 0;
}

のコンパイルは

gcc -o c_quadmath c_quadmath.c -lquadmath

で行う．実行結果は

$ ./c_quadmath 
---- double precision ----
1/3 = 0.333333333333333314829616256247
pi  = 3.141593
---- long double prescision ----
1/3 = 0.333333333333333333342368351437
pi  = 3.141592653589793115997963468544
---- quad precision ----
1/3 = 0.3333333333333333333333333333333333172839
pi1  = 3.1415926535897931159979634685441851615906
pi2  = 3.1415926535897932384626433832795027974791
sqrt1  = 1.4142135623730951454746218587388284504414
sqrt2  = 1.4142135623730950488016887242096981769402

 
pi     3.1415926535897932384626433832795028841971693993751
sqrt(2)  1.4142135623730950488016887242096980785696718753769

の厳密解と比較して pi1とsqrt1がdouble程度の精度しかない． これは代入する変数の型と関数の返り値の型が一致していないため，関数の方の型で丸められてしまったせいだと考えられる．(c++でも？)
既存のプログラムの精度を倍精度から4倍精度に書き直すにはc言語では関数も書き換えなければならない． (gccがversion upしたらなんとかなるんか？) 精度の仕様変更に関しては圧倒的にfortranの方がいじりやすいと個人的に思う．
 (fortran はあんまり使ってないけどね・・・)