top page > la vie en rose > here

long double の利用法


概要

i386 + gcc の拡張倍精度 (long double) を用いて,より高精度の数値計算を行なう方法を紹介します.

対象は次のようなものです.

  1. OS は FreeBSD の 5 系以降 (/i386 および /pc98) (4 系にも言及します.)
  2. CPU は i386
  3. gcc は version 2 から 4 まで.
  4. clang は version 3.3

題材は次のようなものです.

  1. 拡張倍精度 (long double) の使用法と注意点
  2. 拡張倍精度用の算術関数の使用法
  3. 拡張倍精度用に関する gcc の bug

目次

1. i386 の拡張倍精度の概説
2. 拡張倍精度の利用
3. 補足と注意点
3.1 計算時間
3.2 定数の取り扱いと演算の制御
3.3 FreeBSD 4 系における拡張倍精度の利用
3.4 FreeBSD 用 gcc における拡張倍精度の bug
3.5 library の問題
4. 資料
4.1 浮動小数点数に関する資料の所在
4.2 変数のサイズ
4.3 FreeBSD 用 gcc における bug の点検法

1. i386 の拡張倍精度の概説

i386 用の gcc の浮動小数点数には単精度,倍精度,拡張倍精度,4 倍精度 (float, double, long double, __float128) の 4 種類があります. 各々の表現精度は 10 進数で約 7, 16, 19, 34 桁です.

ほとんどの i386 CPU には,i387 と呼ばれる浮動小数点数の演算機能 (FPU: Floating Point Unit) が備わっています.[1] その利点は,それらを高速に計算することができることです. i387 は,単精度,倍精度に対しては四則演算だけでなく,いくつかの算術関数も処理します.[2] また拡張倍精度に対しては,四則演算は処理しますが,算術関数は扱うものがより限定されます. 一方,4 倍精度を処理する機能は備わっておりません.

ベクトル演算に関してですが,SSE (Streaming SIMD Extensions) では float を,SSE2 では double を扱うことができますが,long double は SSE4 に至っても扱えないようです. (あくまでカタログ情報ですが.)[3]

本題の i386 の拡張倍精度ですが,単精度や倍精度に比べると実装の仕方が不自然です. 次節以降では,プログラムを書くために必要なの拡張倍精度の「癖」を概説します.

それと,精度は倍精度に比して 3 桁しか向上しておらず,拡張の仕方が中途半端とも言えます.[4] 適宜,4 倍精度や qlib と得失を比較すると良いでしょう.

[1] 過去には i486SX のように FPU がない CPU もありました. (さらに昔には,それがない方がむしろ普通で,別売りの i8087 を重宝していました.) そのような CPU を使ってみると,FPU のありがたみが体感できます.

[2] ある関数の計算機能が FPU にあっても,それを gcc や library (libm) が利用しているとは限りません. さらに,その状況は gcclibm の version や compile option により異なります. その点を明確に記した文書が見つからないので,ソースを見て判別するしかありません.

[3] さすがに需要がないのかもしれません.

[4] 元々これは倍精度の演算の精度を落とさないようにした「内部表現」であり,本来表に出てくるものではなかったという事情もあると思います.

2. 拡張倍精度の利用

拡張倍精度 (long double) は,次のようにすれば利用することができます.

  1. 変数を long double で型宣言する
  2. 定数を long double で扱う場合,数値の後に L を付ける. (実行文だけでなく,数値を読み込む場合も.)
    (小文字の l (エル) でも良いのですが,1 (いち) との判別が面倒です.)
  3. 出力の書式を %Lf, %Le, %LE, %Lg, %LG のいずれかで,入力の書式を %Lf, %Le のいずれかで指定する.
    (こちらは小文字の l (エル) ではダメです.)
  4. 拡張倍精度 (変数,定数,関数) の四則演算を行なうならば,その実行文よりも先に fpsetprec(FP_PE) を実行する. これは一度実行すれば,そのプログラムの終了まで有効です.
    (これは gcc でも clang でも必要です.)

long double 変数の四則演算しか行なわないのであれば,次のようにします.

#include <stdio.h>
#ifdef __FreeBSD__
#include <floatingpoint.h>               /*  必ず include する.       */
#endif

int main(void)
{
   long double t;                        /*  型は long double である. */

#ifdef __FreeBSD__
   fpsetprec(FP_PE);                /*  これを演算の実行文の前におく. */
#endif

   t = 0.25L*3.14159265358979323846L;        /*  定数には L を付ける. */

   printf("t = %Le\n",t);            /* 入出力の書式は L を付加する.  */

   return 0;
}

なお補足ですが,

  1. この例の場合には,gcc に対して特別な option は必要ない.
  2. FreeBSD しか使わないのであれば,#ifdef __FreeBSD__ とそれに対応する #endif は必要ない.
  3. FreeBSD 5 系以後ならば include するのは floatingpoint.h でなく ieeefp.h でも良い.
  4. この例の場合には,math.hincludelibm との link (compile の際の -lm option) は必要ないが,一般には必要である.

一方,算術関数ですが,以下の点に注意が必要です.

  1. long double に対応している算術関数は double のそれよりも遥かに少ない.
    FreeBSD 4.9 の時点で sqrtl(), sinl(), cosl() のみ.
    FreeBSD 5.3, 6.2 の時点ではそれに加えて atan2l(), logl(), expl(), atanl(),及び /usr/include/math.h で定義されている諸関数.(#if 0#endif で comment out されているものが多いので要注意.)
  2. FPU (i387) の命令を builtin 関数として使う. すなわち,ライブラリーの関数を呼び出すのではないので nm で見ても出てこない.
  3. 引数の正当性をチェックしない. (man gcc-funsafe-math-optimizations の部分を参照.)

利用できる関数に関する情報は,/usr/src/contrib/gcc/config/i386/i386.md にありますが,解り難いです.

;; FPU special functions.   ← この下に情報があります.

これらは次のような要領で利用します.

#include <stdio.h>
#ifdef __FreeBSD__
#include <floatingpoint.h>

long double sinl(long double);     /*  sinl() を prototype 宣言する.  */
#endif

int main(void)
{
   long double s,t;

#ifdef __FreeBSD__
   fpsetprec(FP_PE);
#endif

   t = 0.25L*3.14159265358979323846L;
   s = sinl(t);

   printf("sin(pi/4) = %Le\n",s);

   return 0;
}

compile ですが,FreeBSD 4 系では次のようにします.

% gcc -ffast-math -mfancy-math-387 s1.c -o s1

また,FreeBSD 5, 6, 7, 8 系では次のようにします.

% gcc -ffast-math s1.c -o s1

option -ffast-math ではルーズになり過ぎて不安ならば,最低限必要な option に限定すれば良いでしょう.

% gcc -funsafe-math-optimizations s1.c -o s1                   ← FreeBSD 5, 6 (gcc 3)
% gcc -funsafe-math-optimizations -fno-math-errno s1.c -o s1   ← FreeBSD 7, 8 (gcc 4)

また,option -ansi を付けるとエラーになります.

なお,これ以外の算術関数を拡張倍精度で使いたい場合には,library である ldouble を利用することができます. (厳密には,上述の算術関数も ldouble のものを使うことになります.)

3. 補足と注意点

3.1  計算時間

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

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

3.2  定数の取り扱いと演算の制御

数値の後に 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 ではこのような「初期設定」は必要なく,より自然に拡張倍精度を使うことができます.

3.3  FreeBSD 4 系における拡張倍精度の利用

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

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

3.4  FreeBSD 用 gcc における拡張倍精度の bug

後の 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 節参照.)

FreeBSD gcc 96 bit 128 bit 備考
 2.2.6 2.2.2.1  この OS 標準の gcc
 3.4 2.7.2.3  この OS 標準の gcc
 4.9 2.95.4  この OS 標準の gcc
 3.0.4  ports/lang/gcc30
 5.2 3.3.3××  この OS 標準の gcc
 5.3 3.2.3  ports/lang/gcc32
 3.3.5××  ports/lang/gcc33
 3.4.2  この OS 標準の gcc
 6.2 3.4.6×  この OS 標準の gcc
 6.3 3.4.6×  この OS 標準の gcc
 4.2.4×  ports/lang/gcc42
 7.0 4.2.1×  この OS 標準の gcc
 7.2 4.2.1×  この OS 標準の gcc
 7.3 4.2.1×  この OS 標準の gcc
 7.4 4.2.1×  この OS 標準の gcc
 8.0 4.2.1×  この OS 標準の gcc
 8.1 4.2.1×  この OS 標準の gcc
 8.2 4.2.1×  この OS 標準の gcc
 4.7.0×  ports/lang/gcc47
 8.3 4.2.2×  この OS 標準の gcc
 4.8.0×  ports/lang/gcc48
 8.4 4.2.1×  この OS 標準の gcc
 4.9.0×  ports/lang/gcc49

表中の記号の意味は次の通り.
〇 : この問題が生じない.
× : この問題が生じる.
− : この環境では実現できない.

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

   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 の環境においては,clanglong double は 96 bit のみです. 上記の bug はありません.

3.5  library の問題

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

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

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

4. 資料

4.1  浮動小数点数に関する資料の所在

i386 用の FreeBSD にて /usr/include/machine/float.h には次のように記されています. ただこれは libcprintf(), 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

4.2  変数のサイズ

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

#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 回で「ちょうど」取り扱えます. 未使用部分はそのための余白です.

4.3  FreeBSD 用 gcc における bug の点検法

次のような簡単な方法で 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 の段階で生じていることが解ります.