long double の利用法

あくまで当方で行なった限りの計算 (上記のような簡単なものでなく，もっと時間の掛るもの) では計算時間は double の場合よりも少し増えるだけです． 私的なプログラムでの経験を普遍化するつもりはありませんが，この程度の話でも目安にはなるでしょう．

当然のことですが，注目したいのは自分のプログラムの実行速度です． それと，計算時間だけでなく結果の確認も必要です． ベンチマーク用のプログラムは自前で用意しておくのが良いでしょう．

数値の後に L を付けない定数は倍精度として処理されます． (厳密には「L, l, F, f を付けない」です．)
従って，拡張倍精度の t に

   t = M_PI;                    /*  定数 M_PI は 21 桁あるが，L が付いていない． */

   t = 3.14159265358979323846;  /*  この数値は M_PI と同じ．  */

のように代入すると，倍精度の表現精度を越える数値が一旦倍精度に丸められ (すなわち精度が落ち)，その後に t に代入されることになります．

また，fpsetprec(FP_PE) を実行しないと，四則演算の精度は double と同じになってしまいます． 例えば，

   x = 0.25L*3.14159265358979323846L;
   y = sinl(x) + 0.5L;

は右辺の演算の精度が倍精度になってしまいます． ただし，これを実行しなくても演算を伴わない代入，例えば

   x = 3.14159265358979323846L;
   y = sinl(x);                 /*  関数の戻り値の代入も「演算を伴わない代入」扱いです．  */

は拡張倍精度が維持されます． また入出力

   static char *d = "1.000000000000000010L";
   sscanf(d,"%Lf",&t);
   printf("t = %21.19Lf\n",t);

も問題なく実行されます． (これはかえって混乱する，迷惑な仕様ですね． おそらく副作用はないので，必ず実行するようにして下さい．)

ちなみに，Linux ではこのような「初期設定」は必要なく，より自然に拡張倍精度を使うことができます．

FreeBSD の 4 系においても，拡張倍精度の計算は 5 系以降と同様に可能です． しかし，困ったことに libc の printf(), scanf() のような入出力が拡張倍精度をサポートしていません． そのため，拡張倍精度数の入出力は倍精度で行なわれてしまいます． この問題は gcc の version を 3 に上げても (すなわち ports/lang/gcc30 等を利用しても) 解決しません． 従って，計算結果の数値そのものを高い精度で入出力したい場合には対応できず，FreeBSD の version を上げるしか方策はありません．

一方，連立一次方程式や交代級数を拡張倍精度で計算するが，最終的に出力する結果は倍精度で良いといった用途では問題なく利用できるでしょう． (ただしその場合にも，行列の要素や級数の係数は倍精度から拡張倍精度に変換するのでなく，最初から拡張倍精度で求めておくべきでしょう．)

後の 4.2 節にて「-m96bit-long-double と -m128bit-long-double とで生成される long double はサイズは違うが精度は違わない」と書きました． 本来そうなるはずですし，昔の FreeBSD では実際にそうなっていました．

しかし，下表に「×」で示すところでは，拡張倍精度の定数の精度が double に落とされてしまうという問題が生じます．
(-m96bit-long-double option 付きまたは default で利用する拡張倍精度数を "96 bit"，-m128bit-long-double option 付きで利用するそれを "128 bit" と，そのサイズで呼称して区別することにします．)

この表以外の version の FreeBSD, gcc をお使いの方は，個々に点検されることをお勧めします． (4.3 節参照．)

再度確認しますが，精度が倍精度に落とされてしまうのは定数です． 入出力および変数や関数の演算においては拡張倍精度が維持されます．

   long double x;

   x = 3.14159265358979323846L;       ← これは精度が落ちる．

   x = sqrtl(1.5L);                   ← これは精度が落ちる．

   static char *d = "3.14159265358979323846L";
   sscanf(d,"%Lf",&x);                ← これは精度が落ちない．(定数ではないので)

比較的新しい FreeBSD をお使いならば，(gcc 3.3.5 およびその前後を除けば) 128 bit ではこの問題は生じないので，この bug が直るまではそちらを使えば良いでしょう．

参考までに，この bug は FreeBSD/{i386,pc98} + intel 特有のもので，FreeBSD/amd64 + AMD, FreeBSD/amd64 + intel, Mac OS X + intel, Linux + intel では生じていないそうです．[1][2][3]

[1] この情報，および上記の表の FreeBSD 6.3, 7.0 についての情報をいただいた方々に御礼申し上げます．
  FreeBSD-users-jp 91681, 91685, 91686 (2008.6)

[2] gcc ですが，FreeBSD としては最後の GPL 2 である 4.2.1 よりも後の version (GPL 3) は採用しないらしいです． それでは困るので，別の c compiler (LLVM C) に移行中です． (ここから XFree86 を想起するのは私だけ？)
  なお，ports からはより新しい version の gcc が利用できます．

[3] FreeBSD/i386 の環境においては，clang の long double は 96 bit のみです． 上記の bug はありません．

FreeBSD 6.2 の時点ですが，次のような問題に遭遇しました． 既に修正されていれば良いのですが．

(1) long double のための関数は意外と用意されていない． 例えば fmod() に対応する fmodl() がない．
  (libm だけでなく ports の ldouble にも．)

(2) long double として 128 bit を用いた場合，printf() で long double の変数の値を一つ出力すると，その後に続く変数の出力がおかしくなる． この問題は，printf() で一つの変数しか出力しないようにすれば避けられる．

i386 用の FreeBSD にて /usr/include/machine/float.h には次のように記されています． ただこれは libc の printf(), scanf() のような入出力にしか効いていないようです．

#define FLT_RADIX       2               /* b */

#define FLT_MANT_DIG    24              /* p */
#define FLT_DIG         6               /* floor((p-1)*log10(b))+(b == 10) */
#define FLT_MAX_EXP     128             /* emax */
#define FLT_MAX_10_EXP  38              /* floor(log10((1-b**(-p))*b**emax)) */

#define DBL_MANT_DIG    53
#define DBL_DIG         15
#define DBL_MAX_EXP     1024
#define DBL_MAX_10_EXP  308

#define LDBL_MANT_DIG   64              ← FreeBSD 4 系においては，LDBL_* は DBL_* と同一です．
#define LDBL_DIG        18
#define LDBL_MAX_EXP    16384
#define LDBL_MAX_10_EXP 4932

• RADIX : 基数．
• MANT_DIG : 仮数の bit 数． (FLT (float), DBL (double) には仮数部の bit 数に 1 が加えられている．)
• MAX_EXP : 指数 e の取り得る最大値．(この値から指数部の bit 数を求めることができる．)
• DIG : 10 進数の精度 (桁)．上の表の値は辛めのように思う．
• MAX_10_EXP : 10 進数で表現可能な桁数．(表現精度ではない！)

関連する man page や info は以下の通り． (jman でなく man でも良いですが．)

% jman math
% jman ieee
% jman fpsetprec
% jman fenv
% jman gcc
% info gcc     ← man gcc の内容は概して古いので，こちらも参照する．

また libm に関連するソースは，以下のところにあります． FreeBSD 4 系は，

/usr/src/lib/libm

FreeBSD 5, 6 系 (おそらくそれ以降もしばらくは) は，

/usr/src/lib/msun

gcc に関連するソースは，以下のところにあります．

/usr/src/contrib/gcc

浮動小数点数の種々の型のサイズは，次のようにして知ることができます．

#include <stdio.h>

int main(void)
{
   printf("float       = %2d byte\n",sizeof(float));
   printf("double      = %2d byte\n",sizeof(double));
   printf("long double = %2d byte\n",sizeof(long double));

   return 0;
}

gcc の 3 系以降には -m96bit-long-double, -m128bit-long-double という option があります． 前者が default です． これは long double のサイズを 96 または 128 bit に規定するものです． ただし，これはあくまで long double のサイズをどうとるかの問題であり，どちらを選んでも実際に使われるのは 80 bit のみで，精度は変わりません． そのことは，次のようなプログラムにより確認することができます． (この例は，3.4 節に記したような bug がない場合です．)

#include <stdio.h>

int main(void)
{
   long double lx;

   lx = 3.14159265358979323846L;

   ana_ldouble(lc);
   printf("  lx = %22.20Lf\n",lc);

   return 0;
}

void ana_ldouble( long double x )
{
   int    i;
   unsigned int s;
   typedef union {  long double     l;
                    unsigned char   c[16]; }  exd;
   exd e;

   s = sizeof(long double);

   for ( i = 0; i < 16; i++ )  e.c[i] = '\0';  /*  clear  */
   e.l = x;

   printf("sizeof(long double) = %ud\n",s);
   for ( i = s-1; i <= 0; i-- )  printf("%02x ",e.c[i]&0xff);
}

-m96bit-long-double 付きまたは default で compile した場合の出力は次のようになります． 青字の 00 は未使用部分です．

sizeof(long double) = 12
00 00 40 00 c9 0f da a2 21 68 c2 35   lx = 3.14159265358979323851

-m128bit-long-double 付きで compile した場合の出力は次のようになります．

sizeof(long double) = 16
00 00 00 00 00 00 40 00 c9 0f da a2 21 68 c2 35   lx = 3.14159265358979323851

96 bit は 32 bit で 3 回，128 bit は 64 bit で 2 回で「ちょうど」取り扱えます． 未使用部分はそのための余白です．

次のような簡単な方法で bug の有無を点検することができます． (問題を生じる代表として FreeBSD 6.2 + gcc 3.4.6 を，生じない代表として FreeBSD 5.3 + gcc 3.4.2 を用い，各々の OS の version だけを書くことにします．)

#include <stdio.h>

int main(void)
{
   double c;
   long double lc;

   c  = 3.14159265358979323846;     /*  こちらは比較のため              */
   lc = 3.14159265358979323846L;    /*  bug があれば，こちらに差が出る  */

   return 0;
}

これを s.c とし，FreeBSD 6.2 にて次のように二通りに compile して assembler source を出力します．

% gcc -m96bit-long-double  -S s.c -o  s1.a    (-m96bit-long-double は default)
% gcc -m128bit-long-double -S s.c -o  s2.a

こうして得られる source の中の，これらの数値に関連した部分は次のようになります．

変数 c
        .long   1413754136         ← 仮数部 (2)
        .long   1074340347         ← 符合，指数部，仮数部 (1) が混在

変数 lc
        movl    $560513024, %eax   ← 仮数部 (2). 128 bit ではここは $560513589 です．
        movl    $-921707870, %edx  ← 仮数部 (1)
        movl    $16384, %ecx       ← 符合と指数部

一方，FreeBSD 5.3 ではこれら二つの assembler source は同じであり，かつ FreeBSD 6.2 の 128 bit のそれとも同じです．

この差異の有無と，差異の内容とから，問題としている bug の有無が判定できます．

参考までに，この数値から double の変数 c の仮数部 (mant) を求めると次のようになります．

double : 1074340347 1413754136   mant = 09 21 fb 54 44 2d 18

これを long double の変数 lc の仮数部と比較するために，左に 3 bit シフトし，implicit bit を持たない表現へ変換するために最上位の 1 bit を立てます． その結果と lc の 96, 128 bit の場合の仮数部を並記すると次のようになります．

double :                        mant = c9 0f da a2 21 68 c0
 96bit : -921707870 560513024   mant = c9 0f da a2 21 68 c0 00   ← double と同じ数値である．
128bit : -921707870 560513589   mant = c9 0f da a2 21 68 c2 35

すなわち，96 bit の拡張倍精度の場合，定数の精度が倍精度に落とされています． またそれは実行時でなく compile の段階で生じていることが解ります．