ports ldouble の利用法

概要

拡張倍精度の算術関数を提供する library である ldouble の使用法や注意点を提示します．

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

1. OS は FreeBSD/i386 の 5 系以降 (ldouble は FreeBSD/amd64 には対応しておりません．)
2. CPU は i386
3. gcc は version 3, 4
4. ldouble の version は記載なし． (配布ファイルの中の最新のものの日付は 2003.5.4 である．)

long double の 2 節 にて述べたように，拡張倍精度 (long double) にて使うことができる算術関数は限定されています． しかし，その問題は qlib と同じく Cephes が配布している library である ldouble を用いればある程度解決します． これがサポートしている関数は添付されている文書 ldoubdoc.html に一覧がありますが，少なくとも libm の倍精度用の関数に対応するものはすべてあります．

とはいえ，特に FreeBSD における ldouble の利用に関しては，次のような問題があります．

1. FreeBSD にて利用するための適切な文書類が用意されていない．
2. ports になったのは比較的最近である． (FreeBSD 4, 5, 6 系で，各々 4.11, 5.3, 6.0 から.)
3. ldouble は拡張倍精度が 96 bit であることを前提に作られており，128 bit に対応していない．[1]

これに FreeBSD の側の long double に関わる問題が重なって，面倒な状態になっています． それは次の二点です．

1. FreeBSD 4 系以前では，printf(), scanf() 等が long double に対応していない．
2. FreeBSD 5 系以降は，gcc 3.4.6 以降で 96 bit の long double を用いると精度が落ちるいう bug がある．
   また gcc 3.3.5 前後においては，long double が 96/128 bit のいずれでも精度が落ちる． [2]

要は，拡張倍精度を使う環境が常に入念に整備され点検された状態にあるわけではないということです． これらの問題には次のように対処します．

1. ldouble が ports になっていない場合には，自分で install する．
2. FreeBSD 4 系では，printf(), scanf() 等で拡張倍精度の入出力を行なわない．
3. 96 bit の拡張倍精度に問題がある場合，128 bit を使い，ldouble も自分で 128 bit に対応させる．
   (なお，gcc でなく clang を使うならば，この作業は不要です．)
4. gcc 3.3.5 は使わない． 必要ならば ports から別の version の gcc を導入する．

次の節以降で，その手順を提示します．

[1] -m96bit-long-double option 付きまたは default で利用する拡張倍精度数を "96 bit"，-m128bit-long-double option 付きで利用するそれを "128 bit" と，そのサイズで呼称して区別することにします．
(ちなみに，clang においては 96 bit のみです．)

[2] これは「FreeBSD における gcc の問題」です． 症状や FreeBSD の version との関係等，詳細は long double を参照．

自力で library を作成した場合，同じ directory に点検用のプログラムも作成されます． makefile でいえばこの部分です．

all: mtstl lparanoi lcalc fltestl nantst testvect monotl libml.a

mtstl はいくつかの関数とその逆関数を組み合わせて精度を確認しています． これを実行してみましょう．

% mtstl

FreeBSD 5.3 (96, 128 bit), 6.2 (128 bit) では次のように出力されるはずです．

Consistency test of math functions.
Max and rms errors for 10000 random arguments.
A = absolute error criterion (but relative if >1):
Otherwise, estimate is of relative error
x =   cbrt(   cube(x) ):  max = 7.65E-20   rms = 4.39E-21
x =   atan(    tan(x) ):  max = 2.01E-19   rms = 3.96E-20
x =    sin(   asin(x) ):  max = 2.15E-19   rms = 3.00E-20
x =   sqrt( square(x) ):  max = 0.00E+00   rms = 0.00E+00
x =    log(    exp(x) ):  max = 5.42E-20 A rms = 1.87E-21 A
x =   log2(   exp2(x) ):  max = 1.08E-19 A rms = 3.37E-21 A
x =  log10(  exp10(x) ):  max = 2.71E-20 A rms = 6.76E-22 A
x =  acosh(   cosh(x) ):  max = 3.13E-18 A rms = 3.21E-20 A
x = pow( pow(x,a),1/a ):  max = 1.25E-17   rms = 1.70E-19
x =   tanh(  atanh(x) ):  max = 1.08E-19   rms = 1.16E-20
x =  asinh(   sinh(x) ):  max = 1.03E-19   rms = 2.94E-21
x =    cos(   acos(x) ):  max = 1.63E-19 A rms = 4.37E-20 A
lgam(x) = log(gamma(x)):  max = 2.31E-19 A rms = 5.93E-20 A
x =  ndtri(   ndtr(x) ):  max = 5.07E-17   rms = 7.02E-19
x =   ndtr(  ndtri(x) ):  max = 1.40E-18   rms = 1.57E-19
Legendre  ellpk,  ellpe:  max = 1.19E-18 A rms = 1.74E-19 A
Relative error and only 100 trials:
x = stdtri(stdtr(k,x) ):  max = 2.75E-18   rms = 5.01E-19
Absolute error and only 2000 trials:
Wronksian of   Yn,   Jn:  max = 4.24E-17 A rms = 9.54E-19 A

FreeBSD 6.2 (96 bit) では次のように出力されるはずです． 上記の結果に比べて精度が悪いです．
計算のすべてまたは大部分を FPU に頼っているものは，比較的精度が落ちません． 一方，多項式に基づいて計算するものは，その係数の精度が落とされるがゆえに精度が悪くなります．

Consistency test of math functions.
Max and rms errors for 10000 random arguments.
A = absolute error criterion (but relative if >1):
Otherwise, estimate is of relative error
x =   cbrt(   cube(x) ):  max = 7.65E-20   rms = 4.51E-21
x =   atan(    tan(x) ):  max = 7.60E-17   rms = 3.91E-17
x =    sin(   asin(x) ):  max = 1.07E-16   rms = 4.24E-17
x =   sqrt( square(x) ):  max = 0.00E+00   rms = 0.00E+00
x =    log(    exp(x) ):  max = 5.12E-18 A rms = 1.20E-19 A
x =   log2(   exp2(x) ):  max = 3.54E-17 A rms = 2.05E-18 A
x =  log10(  exp10(x) ):  max = 1.13E-17 A rms = 2.27E-18 A
x =  acosh(   cosh(x) ):  max = 4.02E-17 A rms = 9.47E-19 A
x = pow( pow(x,a),1/a ):  max = 1.51E-14   rms = 1.87E-16
x =   tanh(  atanh(x) ):  max = 1.90E-17   rms = 5.95E-18
x =  asinh(   sinh(x) ):  max = 3.34E-18   rms = 2.59E-19
x =    cos(   acos(x) ):  max = 1.06E-16 A rms = 5.67E-17 A
lgam(x) = log(gamma(x)):  max = 4.85E-17 A rms = 1.86E-17 A
x =  ndtri(   ndtr(x) ):  max = 4.91E-16   rms = 9.64E-17
x =   ndtr(  ndtri(x) ):  max = 8.71E-16   rms = 1.82E-16
Legendre  ellpk,  ellpe:  max = 5.28E-16 A rms = 1.76E-16 A
Relative error and only 100 trials:
x = stdtri(stdtr(k,x) ):  max = 1.26E-15   rms = 2.41E-16
Absolute error and only 2000 trials:
Wronksian of   Yn,   Jn:  max = 1.51E-14 A rms = 3.40E-16 A

なお，FreeBSD の version や拡張倍精度のサイズによらず，fpsetprec(FP_PE) を実行しないとさらに精度が悪くなります． こういうことは，文書になっていないと困りますよね．

(なお，clang においてもfpsetprec(FP_PE) は必要です)

例えば自作のプログラム s1.c は以下のようにして compile します． (諸々の option は省略．)

% gcc s1.c -o s1 -lml

自作のプログラムを library がある directory 以外のところで作るのであれば，libml.a と protos.h をそこに cp するか，または次のようにしてそれらの所在を指示する必要があります． (protos.h は ../include に，libml.a は ../lib にあるとして，)

% gcc -I../include s1.c -o s1 -L../lib -lml

一方，ports の ldouble を利用するのであれば次のようにします． (shared library を使う場合．)

% gcc s1.c -o s1 -lml

また source は次のようにします． (mathl.h と protos.h とは名前と置き場が違うだけで，中身は同一です．)

#include <stdio.h>

#if 1                          /*  1 を 0 に変えることにより使い分ける． */
#include <mathl.h>             /*  ports の ldouble を使う場合．         */
#else
#include "protos.h"            /*  自前で作成した ldouble を使う場合．   */
#endif

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

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;
}