数値計算ガイド ホーム目次前ページへ次ページへ索引

第 3 章

数学ライブラリ

この章では、Solaris オペレーティング環境と Sun WorkShop Compiler で提供されている数学ライブラリ libm.a、libm.so、libsunmath.a を含む数学関数の実装について説明します。この章では、各ライブラリの内容を示すとともに、Sun WorkShop Compiler に含まれる数学ライブラリがサポートする機能の一部 (IEEE をサポートする関数、乱数発生機構、IEEE と非 IEEE 形式間でデータを変換する関数など) について説明します。

libm と libsunmath ライブラリの内容は、intro(3M) のマニュアルページにも挙げられています。

数学ライブラリ

libm 数学ライブラリには、Solaris オペレーティングシステムが準拠している標準に必要な関数が含まれています。このライブラリは、静的ライブラリ libm.a と共有ライブラリ libm.so の 2 つの形式で Solaris に含まれています。

標準インストールをした場合の libm のデフォルトディレクトリは、次のとおりです。

libm の標準位置

/opt/SUNWspro/<release>/lib/libm.a
/opt/SUNWspro/lib/libm.so
/usr/lib/libm.a
/usr/lib/libm.so

libm のヘッダーファイル

/usr/include/floatingpoint.h
/usr/include/math.h
/usr/include/sys/ieeefp.h

表 3-1 は、ライブラリの関数の一覧です。

表 3-1   libm の内容  
型 関数名
代数関数 
cbrt、hypot、sqrt
基本超越正関数 
asin、acos、atan、atan2、acosh、asinh、atanh、

exp、expm1、pow、

log、log1p、log10、

sin、cos、tan、sinh、cosh、tanh
高級超越正関数 
j0、j1、jn、y0、y1、yn、

erf、erfc、gamma、lgamma、gamma_r、lgamma_r
整数丸め関数 
ceil、floor、rint
IEEE 規格で勧告の関数 
copysign、fmod、ilogb、nextafter、remainder、 
scalbn、fabs
IEEE 分類の関数 
isnan
旧式の浮動小数点関数 
logb、scalb、significand
浮動小数点トラップ
処理 
libm__sigfpe
エラー処理関数 
matherr


関数 gamma_r と lgamma_r は、gamma と lgamma の再入可能バージョンを意味するので注意してください。

上記で述べているように、標準の数学ライブラリは 2 か所にインストールできます。
/usr/lib に置かれるライブラリは、Solaris オペレーティング環境にバンドルされており、すべてのシステムにインストールされます。一方、/opt/SUNWspro に置かれるライブラリは、Sun WorkShop Compiler に付属しており、コンパイラインストールにアクセスできるシステム上でのみ利用できます。この 2 セットのライブラリは同じインタフェースを使用したまったく同じ関数を提供しますが、わずかな違いがあります。まず、コンパイラに付属した標準の数学ライブラリ内の関数は、Solaris オペレーティング環境にバンドルされたライブラリ内の関数よりも高速な場合があります。また、これらのライブラリの正確度は同等ですが、一定の引数について特定の関数が配布する結果は最下位ビットが異なる場合があります。ほとんどのプログラムはこれらのライブラリ間の数値の違いに特に影響は受けませんが、共有 libm.so にリンクされたプログラムをコンパイラに付属の libm.so を使用したシステムで実行した場合と、デフォルトの Solaris libm.so だけを使用したシステムで実行した場合では、この違いのためにわずかな変化が出ます。また、プログラムを実行する Solaris リリースによっても差が出ることがあります。

プログラムが多様な結果を出すのを防ぐ方法はいくつかあります。これは、使用する数学ライブラリによって異なります。

動的リンクと静的リンク、およびプログラムの実行時に読み込む共有オブジェクトを決定するオプションと環境変数の詳細は、ld(1) とコンパイラのマニュアルページを参照してください。

付加価値数学ライブラリ

Sun 数学ライブラリ

libsunmath ライブラリは、サンの言語製品と共に提供されるライブラリの一部です。libsunmath ライブラリには、サンの旧バージョンの libm に組み込まれていた関数が格納されています。

libsunmath がインストールされるデフォルトのディレクトリ

/opt/SUNWspro/<release>/lib/libsunmath.a
/opt/SUNWspro/lib/libsunmath.so

libsunmath のヘッダーファイルがインストールされるデフォルトのディレクトリ

/opt/SUNWspro/<release>/include/cc/sunmath.h
/opt/SUNWspro/<release>/include/f77/f77_floatingpoint.h

次の表 3-2 は、libsunmath に含まれる関数を示しています。この表では、数学関数ごとに、C プログラムから呼び出される場合の、関数の倍精度名だけを示しています。

表 3-2   libsunmath の内容  
型 関数名
表 3-1 からの関数 matherr を除く単精度および拡張 / 4 倍精度の関数
基本超越正関数 exp2、exp10、log2、sincos
三角関数 (度数) asind、acosd、atand、atan2d、
sind、cosd, sincosd、tand
π で基準化 (スケール)
した三角関数 		asinpi、acospi、atanpi、atan2pi、
sinpi、cospi、sincospi、tanpi
倍精度 π 引数還元のある三角関数 asinp、acosp、atanp、
sinp、cosp、sincosp、tanp
財務関数 annuity、compound
整数丸め関数 aint、anint、irint、nint
IEEE 規格で勧告済みの
関数 		signbit
IEEE 分類の関数 fp_class、isinf、isnormal、issubnormal、iszero
有用な IEEE 値を与える
関数 		min_subnormal、 max_subnormal、
min_normal、 max_normal、
infinity、signaling_nan、quiet_nan
加法的乱数の生成 i_addran_、i_addrans_、i_init_addrans_、i_get_addrans_、i_set_addrans_、
r_addran_、r_addrans_、r_init_addrans_、r_get_addrans_、r_set_addrans_、
d_addran_、d_addrans_、d_init_addrans_、d_get_addrans_、d_set_addrans_、
u_addrans_
線形合同乱数の生成 i_lcran_、i_lcrans_、i_init_lcrans_、i_get_lcrans_、i_set_lcrans_、
r_lcran_、r_lcrans_、d_lcran_、d_lcrans_、
u_lcrans_
桁上げ乗算乱数の生成 i_mwcran_、i_mwcrans_、i_init_mwcrans_、i_get_mwcrans_、i_set_mwcrans、
i_lmwcran_、i_lmwcrans_、i_llmwcran_、i_llmwcrans_、
u_mwcran_、u_mwcrans_、u_lmwcran_、u_lmwcrans、u_llmwcran_、u_llmwcrans_、
r_mwcran_、r_mwcrans_、d_mwcran_、d_mwcrans_,smwcran_
乱数シャフル i_shufrans_、r_shufrans_、d_shufrans_、u_shufrans_
データ変換 convert_external
制御丸めモード、および浮動小数点例外フラグ ieee_flags
浮動小数点トラップ処理 ieee_handler、sigfpe
状態表示 ieee_retrospective
標準外演算の
有効化 / 無効化 		standard_arithmetic、nonstandard_arithmetic


最適化ライブラリ

一部の libm ルーチン中の最適化バージョンは、libmopt ライブラリによって提供されます。一部の libc 中のサポートルーチンは、libcopt ライブラリによって提供されます。また、SPARC では、一部の libc サポートルーチンの代替は libcx によって提供されます。

標準インストールをした場合の libmopt、libcopt、libcx のデフォルトディレクトリは以下のとおりです。

/opt/SUNWspro/<release>/lib/<arch>/libmopt.a
/opt/SUNWspro/<release>/lib/<arch>/libcopt.a
/opt/SUNWspro/<release>/lib/<arch>/libcx.a (SPARC のみ)
/opt/SUNWspro/<release>/lib/<arch>/libcx.so.1 (SPARC のみ)

(<arch> はアーキテクチャ固有のライブラリディレクトリを示す。SPARC では、v7、v8、v8a、v8plus、v8plusa、v8plusb、v9、v9a、v9b など。x86 プラットフォームでは、提供されるディレクトリは f80387 のみ。)

-xarch の詳細は、『Fortran ユーザーズガイド』 、『C ユーザーズガイド』または『C++ ユーザーズガイド』を参照してください。

libcopt に含まれるルーチンは、ユーザーが直接呼び出すことはできません。その代わり、libc 内にサポートルーチンを置き換えて、コンパイラが使用できるようにします。

libmopt に含まれるルーチンは、libm 内の対応するルーチンに置き換えられます。libmopt バージョンの処理速度は、一般的に著しく速くなります。これは、libm バージョンは、ANSI/POSIX、SVID、X/Open、または IEEE 形式の例外のケースを処理するために構成されますが、libmopt ルーチンは IEEE 形式の例外のケースの処理のみをサポートするためです (付録 E を参照してください)。

cc を使用して libmopt と libcopt の両方とリンクするには、コマンド行で -lmopt と -lcopt オプションを指定します (-lm の直前に -lmopt を置き、-lcopt を最後に置くともっともよい結果となる)。ほかのコンパイラを使用してこの 2 つのライブラリとリンクする場合は、コマンド行の任意の位置に -xlibmopt フラグを指定します。

SPARC では、ライブラリ libcx に 128 ビット 4 倍精度浮動小数点算術演算サポートルーチンよりも速いバージョンを含みます。これらのルーチンはユーザーが直接呼び出すことはできず、コンパイラによって呼び出されます。-nocx オプションが指定されていない場合、FORTRAN 77 コンパイラでは自動的に libcx と共にリンクが行われますが、C コンパイラでは行われません。C プログラムで libcx を使用するには、-lcx とともにリンクを行なってください。

libcx の共有バージョン (libcx.so.1) も提供されます。共有バージョンを実行時にあらかじめロードするには、環境変数 LD_PRELOAD を libcx.so.1 ファイルのフルパス名に指定してください。パフォーマンスを最大にするには、使用しているシステムのアーキテクチャに該当する libcx.so.1 を使用します。たとえば UltraSPARC システムでは、ライブラリがデフォルトの位置にインストールされている場合、LD_PRELOAD を次のように設定します。

csh:

setenv LD_PRELOAD /opt/SUNWspro/lib/v8plus/libcx.so.1

sh:

LD_PRELOAD=/opt/SUNWspro/lib/v8plus/libcx.so.1
export LD_PRELOAD

ベクトル数学ライブラリ (SPARC のみ)

SPARC プラットフォームでは、libmvec ライブラリの提供するルーチンは、引数の全部のベクトルに対して共通数学関数を評価します。

標準インストールを行なった場合の libmvec のデフォルトディレクトリは次のとおりです。

/opt/SUNWspro/<release>/lib/<arch>/libmvec.a
/opt/SUNWspro/<release>/lib/<arch>/libmvec_mt.a

表 3-3 は、libmvec 中の関数を一覧しています。

表 3-3   libmvec の内容
種類 関数名
代数関数 vhypot_、vhypotf_、vc_abs_、vz_abs_
指数関数および
関連した関数 		vexp_、vexpf_、vlog_、vlogf_、vpow_、vpowf_、vc_exp_、vz_exp_、vc_log_、vz_log_、vc_pow_、vz_pow_
三角関数 vatan_、vatanf_、vatan2_、vatan2f_、vcos_、
vcosf_、vsin_、vsinf_、vsincos_、vsincosf_


libmvec_mt.a はマルチプロセッサの並列化にもとづいて、ベクトル関数の並列化バージョンを提供します。libmvec_mt.a を使用するには、-xparallel と共にリンクしなければなりません。

詳細は、libmvec(3m) および clibmvec(3m) のマニュアルページを参照してください。

libm9X 数学ライブラリ

libm9x 数学ライブラリには、C99 で規定されている数学および浮動小数点に関した関数の一部が含まれています。Sun WorkShop 6 Compiler リリースでは、改良された浮動小数点例外処理をサポートするために、<fenv.h>
(Floating-Point Environment、浮動小数点環境) 関数および各種の拡張プログラムがこのライブラリに含まれています。

libm9x の標準インストールの場合、デフォルトディレクトリは次のとおりです。

/opt/SUNWspro/lib/libm9x.so

libm9x のヘッダーファイルの標準インストールの場合、デフォルトディレクトリは次のとおりです。

/opt/SUNWspro/<release>/include/cc/fenv.h
/opt/SUNWspro/<release>/include/cc/fenv96.h

次の表 3-4 は、libm9x に含まれる関数を示しています (精度制御関数 fegetprec と fesetprec は、x86 プラットフォームのみに含まれています)。

表 3-4   libm9x の内容
種類 関数名
C99 規格浮動小数点環境関数 feclearexcept、fegetenv、fegetexceptflag、
fegetround、feholdexcept、feraiseexcept、fesetenv、fesetexceptflag、fesetround、fetestexcept、
feupdateenv
精度制御 (x86) fegetprec、fesetprec
例外処理と遡及診断 fex_get_handling、fex_get_log、
fex_getexcepthandler、fex_log_entry、
fex_merge_flags、fex_set_handling、fex_set_log、fex_setexcepthandler


libm9x は、共有ライブラリとしてのみ提供されています。cc コンパイラは、リンク時に共有ライブラリのインストールディレクトリ内でライブラリを自動的に検索することはありません。そのため、cc を使用して libm9x とリンクする場合は、静的リンカーと実行時リンカーの両方を有効にして、ライブラリを検出する必要があります。libm9x を検出するために静的リンカーを有効にするには、次の 3 つの方法のいずれかを使用します。

libm9x を検出するために実行時リンカーを有効にするには、次の 3 つの方法のいずれかを使用します。

ほかの Sun WorkShop Compiler はすべて、自動的に共有ライブラリのインストールディレクトリを検索します。これらのコンパイラのどれかを使用して libm9x とリンクするには、コマンド行で -lm9x と指定してください。libm9x は主に C および C++ プログラムでの使用を想定したものですが、FORTRAN プログラムでも使用できます。例は、付録 A を参照してください。

単精度、倍精度、4 倍精度

多くの数値関数は単精度、倍精度、および 4 倍精度で使用できます。他言語から異なった精度のバージョンを呼び出す例を、表 3-5 に示します。

表 3-5   単精度、倍精度、および 4 倍精度 libm 関数の呼び出し 
言語 単精度 倍精度 4 倍精度
C、 C++ #include <sunmath.h> float x,y,z; x = sinf(y); x = fmodf(y,z); x = max_normalf(); x = r_addran_(); #include <math.h> double x,y,z; x = sin(y); x = fmod(y,z);
#include <sunmath.h> double x, y, z; x = max_normal(); x = d_addran_();
#include <sunmath.h> long double x,y,z; x = sinl(y); x = fmodl(y,z); x = max_normall();
Fortran REAL x,y,z x = sin(y) x = r_fmod(y,z) x = r_max_normal() x = r_addran() REAL*8 x,y,z x = sin(y) x = d_fmod(y,z) x = d_max_normal() x = d_addran() REAL*16 x,y,z x = sin(y) x = q_fmod(y,z) x = q_max_normal()


C では、単精度関数の名前は倍精度名に f を付加し、4 倍精度関数の名前は倍精度名に l を付加することで作成されます。FORTRAN 呼び出し規約は異なっているので、libsunmath は単精度関数、倍精度関数、および 4 倍精度関数用にそれぞれ r_... 、d_... 、および q_... バージョンを提供しています。FORTRAN 組み込み関数は、3 種類の精度のすべてについて総称で呼び出すことができます。

すべての関数が q_... バージョンを持っているわけではありません。libm と libsunmath 関数の名前と定義については、<math.h> と <sunmath.h> を参照してください。

FORTRAN プログラムの中で r_... 関数を real、 d_... 関数を倍精度、q_... 関数を real*16 として宣言することを忘れないでください。そうしないと、結果として型が不整合になることがあります。

注 - x86 バージョンの FORTRAN では単精度と倍精度だけがサポートされ、REAL*16 はサポートされません。ただし、x86 バージョンの C では 4 倍精度がサポートされます。

IEEE サポート関数

この節では、IEEE 推奨関数、有益な値を与える関数、ieee_flags、ieee_retrospective、および standard_arithmetic と nonstandard_arithmetic について説明します。関数 ieee_handler および ieee_flags ついての詳細は、第 4 章を参照してください。

ieee_functions(3m) と ieee_sun(3m)

ieee_functions(3m) と ieee_sun(3m) によって記述される関数は、IEEE 規格で要求される機能または IEEE 規格の付録で推奨される機能を備えています。これらはビットマスク演算として効率的に実装されています。

表 3-6   ieee_functions(3m)  
関数 戻り値 
math.h
ヘッダーファイル 
copysign(x,y)
y の符号ビットを持つ x 
fabs(x)
x の絶対値 
fmod(x, y)
y に関する x の剰余 
ilogb(x)
整数形式である x の基数 2 非バイアス型指数 
nextafter(x,y)
y の方向で x の次に表現可能な数 
remainder(x,y)
y に関する x の剰余 
scalbn(x,n)
x × 2n


　　　　

表 3-7   ieee_sun(3m) 
関数 戻り値
sunmath.h ヘッダーファイル
fp_class(x) 分類関数
isinf(x) 分類関数
isnormal(x) 分類関数
issubnormal(x) 分類関数
iszero(x) 分類関数
signbit(x) 分類関数
nonstandard_arithmetic(void) トグルハードウェア
standard_arithmetic(void) トグルハードウェア
ieee_retrospective(*f)


remainder(x,y) は、IEEE 規格 754-1985 で規定された演算です。remainder(x,y) と fmod(x,y) の相違は、fmod(x,y) が常に x と一致する符号を持つ結果を返すのに対し、remainder(x,y) が返す結果の符号は x または y の符号のどちらとも一致しないことがあるという点です。両関数とも明確な結果を返し、不明確な例外は発生しません。

表 3-8   FORTRAN からの ieee_functions の呼び出し 
IEEE 関数 単精度 倍精度 4 倍精度
copysign(x,y) t=r_copysign(x,y) z=d_copysign(x,y) z=q_copysign(x,y)
ilogb(x) i=ir_ilogb(x) i=id_ilogb(x) i=iq_ilogb(x)
nextafter(x,y) t=r_nextafter(x,y) z=d_nextafter(x,y) z=q_nextafter(x,y)
scalbn(x,n) t=r_scalbn(x,n) z=d_scalbn(x,n) z=q_scalbn(x,n)
signbit(x) i=ir_signbit(x) i=id_signbit(x) i=iq_signbit(x)


表 3-9   FORTRAN からの ieee_sun の呼び出し
IEEE 関数 単精度 倍精度 4 倍精度
signbit(x) i=ir_signbit(x) i=id_signbit(x) i=iq_signbit(x)


注 - ユーザーは d_<関数> および q_<関数> を使用する FORTRAN プログラムの中で d_<関数> を倍精度として宣言し、q_<関数> を REAL*16 として宣言する必要があります。

ieee_values(3m)

ieee_values(3m) のマニュアルページに記述されている特殊関数によって、無限大、NaN、および最大/最小浮動小数点数のような IEEE 値が得られます。
ieee_values(3m) 関数によって記述される 10 進値と 16 進の IEEE 表現を、
表 3-10、表 3-11、表 3-12、表 3-13 に示します。

表 3-10   IEEE 値: 単精度 
IEEE 値 10 進値および IEEE 表現 C、C++
FORTRAN
最大正規数 3.40282347e+38 7f7fffff r = max_normalf();
r = r_max_normal()
最小正規数 1.17549435e-38 00800000 r = min_normalf();
r = r_min_normal()
最大非正規数 1.17549421e-38 007fffff r = max_subnormalf();
r = r_max_subnormal()
最小非正規数 1.40129846e-45 00000001 r = min_subnormalf();
r = r_min_subnormal()
無限大 7f800000 r = infinityf();
r = r_infinity()
シグナルを発生
しない NaN 
NaN

7fffffff
r = quiet_nanf(0);
r = r_quiet_nan(0)
シグナルを発生する NaN 
NaN

7f800001
r = signaling_nanf(0);
r = r_signaling_nan(0)


　　　　

表 3-11   IEEE 値: 倍精度  
IEEE 値 10 進値および IEEE 表現 C、C++
FORTRAN
最大正規数 1.7976931348623157e+308 
7fefffff ffffffff
d = max_normal();
d = d_max_normal()
最小正規数 2.2250738585072014e-308 
00100000 00000000
d = min_normal();
d = d_min_normal()
最大
非正規数 		2.2250738585072009e-308 
000fffff ffffffff
d = max_subnormal();
d = d_max_subnormal()
最小
非正規数 		4.9406564584124654e-324 
00000000 00000001
d = min_subnormal();
d = d_min_subnormal()
無限大 
7ff00000 00000000
d = infinity();
d = d_infinity()
シグナルを発生
しない NaN 
NaN

7fffffff ffffffff
d = quiet_nan(0);
d = d_quiet_nan(0)
シグナルを発生する NaN 
NaN

7ff00000 00000001
d =signaling_nan(0);
d = d_signaling_nan(0)


　　　　　

表 3-12   IEEE 値: 4 倍精度 (SPARC)  
IEEE 値 10 進値および IEEE 表現 C、C++
FORTRAN
最大正規数 1.1897314953572317650857593266280070e+4932 7ffeffff ffffffff ffffffff ffffffff q = max_normall();
q = q_max_normal()
最小正規数 3.3621031431120935062626778173217526e-4932 00010000 00000000 00000000 00000000 q = min_normall();
q = q_min_normal()
最大非正規数 3.3621031431120935062626778173217520e-4932 0000ffff ffffffff ffffffff ffffffff q = max_subnormall();
q = q_max_subnormal()
最小非正規数 6.4751751194380251109244389582276466e-4966 00000000 00000000 00000000 00000001 q = min_subnormall();
q = q_min_subnormal()
無限大 7fff0000 00000000 00000000 00000000 q = infinityl();
q = q_infinity()
シグナルを発生しない NaN NaN 7fff0000 00000000 00000000 00000001 q = quiet_nanl(0);
q = q_quiet_nan(0);
シグナルを発生する NaN NaN 7fff8000 00000000 00000000 00000001 q = signaling_nanl(0);
q = q_signaling_nan(0)


　

表 3-13   IEEE 値 : 拡張倍精度 (x86)  
IEEE 値 10 進値および IEEE 表現
(80 ビット) 		C、C++
最大正規数 1.18973149535723176509e+4932 7ffe ffffffff ffffffff x = max_normall();
最小正規数 3.36210314311209350626e-4932 0001 80000000 00000000 x = min_normall();
最大
非正規数 		3.3621031431120935062626e-4932 0000 7fffffff ffffffff x = max_subnormall();
最小
非正規数 		1.82259976594123730126e-4951 0000 00000000 00000001 x = min_subnormall();
無限大 7fff 00000000 00000000 x = infinityl(); 
シグナルを発
生しない  NaN
NaN 7fff c0000000 00000000 x = q; 
シグナルを発
生する  NaN
NaN 7fff 80000000 00000001 x = signaling_nanl(0);


ieee_flags(3m)

ieee_flags(3m) は、以下の機能に対する SUN インタフェースを提供しています。

ieee_flags(3m) 呼び出しの構文は次の通りです。

i = ieee_flags (action, mode, in, out);

パラメータに対して設定できる値の ASCII 文字列を表 3-14 に示します。

表 3-14   ieee_flags のパラメータ値  
パラメータ C または C++ の型 取り得る値
action char * get、set、clear、clearall
mode char * direction、precision、exception
in char * nearest、tozero、negative、positive、extended、double、single、inexact、division、underflow、overflow、invalid、all、common
out char ** nearest、tozero、negative、positive、extended、double、single、inexact、division、underflow、overflow、invalid、all、common


パラメータについての詳細は、ieee_flags(3m) のマニュアルページを参照してください。

以降に、ieee_flags を使用すると変更可能な算術演算機能を簡単に説明します。ieee_flags および IEEE 例外フラグについての詳細は、第 4 章を参照してください。

mode が direction の時、指定された動作は現在の丸め方向に適用されます。丸め方向は、表現可能な最も近い数に丸める、ゼロ方向に丸める、+· 方向に丸める、-· 方向に丸めるのいずれかに設定できます。IEEE のデフォルトの丸め方向は、表現可能な最も近い数に丸めるよう設定されています。算術演算結果が 隣接する2 つの表現可能な数の間にあるとき、算術結果に近い値が結果となります。算術演算結果が表現可能なもっとも近い 2 つの数のちょうど中間にある場合には、最下位ビットがゼロである値が結果となります。このことを強調するように、表現可能なもっとも近い数への丸めは、もっとも近い偶数値への丸めとも呼ばれます。

ゼロ方向に丸める方法は、IEEE 規格以前に多くのコンピュータが採用していた方法です。これは、数学的には演算結果の切り捨てを意味します。たとえば、2/3 を 6 桁の 10 進数に丸める場合、表現可能な最も近い数に丸めると .666667 になりますが、ゼロ方向に丸めると .666666 になります。

丸め方向の確認、クリアーまたは設定に ieee_flags を使用するときに、4 つの入力パラメータが取り得る値を、表 3-15 に示します。

表 3-15   ieee_flags による丸め方向の入力値
パラメータ 取り得る値
action get、set、clear、clearall
in nearest、tozero、negative、positive
out nearest、tozero、negative、positive


mode が precision の時、指定した動作は現在の丸め精度に適用されます。x86 プラットフォームでは、丸め精度は、単精度、倍精度、拡張精度に設定できます。デフォルトの丸め精度は拡張精度です。このモードでは、浮動小数点レジスタに渡される算術演算の結果は、完全な 64 ビット精度の拡張倍精度レジスタの形式に丸められます。丸め精度が単精度または倍精度の場合は、浮動小数点レジスタに渡される算術演算の結果は、それぞれ上位の 24 ビットまたは 53 ビットに丸められます。拡張丸め精度を使用してもほとんどのプログラムでは少なくとも正確な結果が得られますが、IEEE 算術演算の論理に厳密に準拠している必要があるプログラムの中には、拡張丸め精度モードで正常に動作しないものがあります。このようなプログラムは単精度または倍精度に設定された丸め精度で実行する必要があります。

SPARC アーキテクチャを使用しているシステムでは、丸め精度は設定できません。丸め精度に関連する ieee_flags の呼び出しは、現在の実装では演算結果に影響を与えません。

mode が exception の時、指定された動作は現在の IEEE 例外フラグに適用されます。ieee_flags を使用して IEEE 例外フラグを調べ制御する方法についての詳細は、第 4 章を参照してください。

ieee_retrospective(3m)

libsunmath 関数 ieee_retrospective は、不要な例外および非標準 IEEE モードに関する情報を stderr に出力します。この関数は以下の事項を判断します。

これらの情報はハードウェア浮動小数点状態レジスタに格納されます。

ieee_retrospective は、例外フラグ、および対応するトラップが有効になっている例外についての情報を出力します。これら 2 つの異なった情報を混同しないように気を付けてください。例外フラグが発生している場合、それはプログラム実行中のある時点で発生して無視された例外です。例外に対応するトラップが有効であれば、例外は実際には発生していないことがあります。例外が発生している場合は、SIGFPE シグナルが送信されます。ieee_retrospective メッセージは、調査する必要のある例外についてユーザーに警告する (例外フラグが発生している場合) か、または例外がシグナルハンドラによって処理されていることをユーザーに知らせます (例外がトラップされた場合)。例外がシグナルハンドラをトラップすると、その例外はもう発生しません。例外、シグナル、トラップ、通知された例外の原因調査方法については、第 4 章を参照してください。

ieee_retrospective はいつでも呼び出すことができますが、通常出口の前で呼び出します。FORTRAN 77 でコンパイルされたプログラムでは、デフォルトで ieee_retrospective を出口で呼び出します。

この関数を呼び出すための構文を以下に示します。

C および C++ ieee_retrospective_(fp);
FORTRAN call ieee_retrospective()

C の関数の場合、引数 fp には出力ファイルを指定します。FORTRAN 関数の場合、常に標準エラー出力に出力されます。

この例は、6 種類ある ieee_retrospective 警告メッセージの中の 4 種類を示しています。

 Note: IEEE floating-point exception flags raised:  

    Inexact; Underflow; 

 Rounding direction toward zero 

 IEEE floating-point exception traps enabled: 

    overflow; 

 See the Numerical Computation Guide, ieee_flags(3M), 

    ieee_handler(3M), ieee_sun(3m)


 

[日本語訳]


 

注: 以下の IEEE 浮動小数点例外が発生しました:


 

	不正確、アンダーフロー


 

ゼロ方向への丸め


 

以下の IEEE 浮動小数点例外のトラップが有効です:


 

	オーバーフロー


 

詳細は、『数値計算ガイド』の ieee_flags(3M), ieee_handler(3M),

ieee_sun(3m)に関する説明を参照してください。



  














   警告メッセージは、トラップが有効になっているか、例外が発生した場合にのみ表示されます。






   FORTRAN プログラムからの ieee_retrospective メッセージを抑制するには 3 つの方法があります。1 つめの方法では、処理されていない例外をすべてクリアーし、トラップを無効にし、プログラムが終了する前の、近似値に丸める、拡張精度、標準モードを復元します。これを行うには、次のように ieee_flags、ieee_handler および standard_arithmetic を呼び出します。



  
  

     

character*8 out 


 

i = ieee_flags('clearall', '', '', out) 


 

call ieee_handler('clear', 'all', 0)


 

call standard_arithmetic()



  














   注 - 原因を調査せずに、処理されていない例外をクリアーすることはお勧めしません。






   ieee_retrospective メッセージが表示されないようにするもう 1 つの方法は、標準エラー出力 (stderr) をファイルにリダイレクトする方法です。プログラムが ieee_retrospective 以外のメッセージを stderr に出力する場合は、この方法を使用しないでください。






    3 つめの方法は、次のようにダミーの ieee_retrospective 関数をプログラムに組み込む方法です。



  
  

     

subroutine ieee_retrospective


 

return


 

end



  














   nonstandard_arithmetic(3m) 






   第 2 章で説明したように、IEEE 演算では、段階的アンダーフローを使用して、アンダーフロー例外を処理します。一部の SPARC システム上では、演算のソフトウェアエミュレーションを通じて段階的アンダーフローが実装される場合があります。そのような場合に数多くの計算でアンダーフローが発生すると、性能を低下させる原因となることがあります。






   特定のアプリケーションで上記のような状況が発生しているかどうかを調べる場合は、ieee_retrospective または ieee_flags を使用して、アンダーフロー例外の発生有無を確認し、プログラムで使用されているシステム時間をチェックします。プログラムが多大な時間をシステム呼び出しに費やし、アンダーフロー例外が発生している場合は、段階的アンダーフローが性能低下の原因と考えられます。このような場合には、IEEE 以外の演算機能を使用すると、プログラムの実行速度が上がります。






   nonstandard_arithmetic 関数を呼び出した場合、ゼロへのフラッシュが高速に行われるようなハードウェアモードの SPARC 実装上では、アンダーフローした結果がゼロにフラッシュされます。ただし、段階的アンダーフローの利点が失われるため、高速になる代わりに正確性が低下します。






   standard_arithmetic 関数を呼び出すと、ハードウェアがリセットされ、デフォルトの IEEE 演算機能が復元されます。これら 2 つの関数は、IEEE 754 形式のデフォルト演算機能しか利用できない実装 (SuperSPARCTMなど) では効果がありません。






   C99 浮動小数点環境関数






   この節では、C99 で規定されている <fenv.h> 浮動小数点環境関数について説明します。Sun WorkShop 6 Compiler リリースでは、libm9x.so ライブラリにこれらの関数があります。これらの関数には ieee_flags 関数と同じ機能が多数ありますが、より自然な C インタフェースが使用されており、C99 で定義されているため、将来は移植性がさらに向上する可能性があります。






   注 - 一貫した動作を保つため、libm9x.so に含まれる C99 浮動小数点環境関数と例外処理の拡張機能、および libsunmath 内の ieee_flags と ieee_handler 関数の両方を同じプログラム内で使用することは避けてください。






   例外フラグ関数






   fenv.h ファイルは、5 つの IEEE 浮動小数点例外フラグ (FE_INEXACT、FE_UNDERFLOW、FE_OVERFLOW、FE_DIVBYZERO、および FE_INVALID) のそれぞれに対してマクロを定義しています。また、5 つのフラグマクロすべてのビット単位の論理和となるように、マクロ FE_ALL_EXCEPT を定義します。以下の説明では、excepts パラメータは 5 つのフラグマクロのいずれかのビット単位の論理和であるか、値 FE_ALL_EXCEPT です。fegetexceptflag 関数と fesetexceptflag 関数の場合は、flagp パラメータが型 fexcept_t (この型は fenv.h で定義されている) のオブジェクトを指すポインタでなければなりません。






   C99 は、次に示す例外フラグ関数を定義しています。



    

   表 3-16 	  C99 規格例外フラグ関数

  
 関数

     処理

  

  

    
 feclearexcept (excepts)

    
 指定されたフラグをクリアーする。

  

  

    
 fetestexcept (excepts)

    
 指定されたフラグの設定を返す。

  

  

    
 feraiseexcept (excepts)

    
 指定された例外を発生させる。

  

  

    
 fegetexceptflag (flagp, excepts)

    
 指定された例外を *flagp に保存する。

  

  

    
 fesetexceptflag (flagp, excpts)

    
 指定された例外を *flagp から復元する。

  















   feclearexcept 関数は、指定されたフラグをクリアーします。fetestexcept 関数は、設定されている excepts 引数が指定するフラグのサブセットに対応するマクロ値のビット単位の論理和を返します。たとえば、現在設定されているフラグだけが不正確な場合、アンダーフローが発生し、ゼロ除算が行われ、次の記述が i を FE_DIVBYZERO にセットします。




  

     
  i = fetestexcept(FE_INVALID | FE_DIVBYZRO);
  





   feraiseexcept 関数は、指定された例外のトラップのいずれかが有効であれば、トラップを発生させます (例外トラップの詳細は、第 4 章を参照)。有効なトラップがない場合は、対応するフラグをセットするだけです。






   fegetexceptflag 関数と fesetexceptflag 関数は、特定のフラグの状態を一時的に保存しておき、あとでその状態を復元するのに便利です。fesetexceptflag 関数は、トラップを発生させません。この関数は、指定されたフラグの値を復元するだけです。






   丸め制御






   fenv.h ファイルは、IEEE の 4 つの丸め方向モードである FE_TONEAREST、FE_UPWARD (正の無限大の方向)、FE_DOWNWARD (負の無限大の方向)、および FE_TOZERO のそれぞれについてマクロを定義しています。C99 は、丸め方向モードを制御する 2 つの関数を定義しています。fesetround は、現在の丸め方向をその引数 (これは上記の 4 つのマクロの 1 つでなければならない) で指定された方向に設定します。fegetround は、現在の丸め方向に対応するマクロの値を返します。






   x86 プラットフォームでは、fenv.h ファイルは 3 つの丸め精度モード、FE_FLTPREC (単精度)、FE_DBLPREC (倍精度)、および FE_LDBLPREC (拡張倍精度) のそれぞれについてマクロを定義しています。これらは C99 には含まれませんが、x86 上の libm9x.so は、丸め精度モードを制御する 2 つの関数を提供します。fesetprec は、現在の丸め精度をその引数 (これは上記の 3 つのマクロの 1 つでなければならない) で指定された精度に設定します。fegetprec は、現在の丸め精度に対応するマクロの値を返します。






   環境関数






   fenv.h ファイルは、例外フラグ、丸め制御モード、例外処理モード、標準外モード (標準外モードは SPARC の場合のみ) など、浮動小数点環境全体を表現するデータ型 fenv_t を定義しています。以下の説明では、envp パラメータは型 fenv_t のオブジェクトのポインタでなければなりません。






   C99 は、浮動小数点環境を操作する 4 つの関数を定義しています。libm9x.so は、マルチスレッド対応のプログラムに便利なその他の関数を提供しています。これらの関数の概要を次の表に示します。



    

   表 3-17 	  libm9x.so 浮動小数点環境関数

  
 関数

     処理

  

  

    
 fegetenv(envp)

    
 *envp に環境を保存する。

  

  

    
 fesetenv(envp)

    
 *envp から環境を復元する。

  

  

    
 feholdexcept(envp)

    
 *envp に環境を保存し、連続モードを確立する。

  

  

    
 feupdateenv(envp)

    
 *envp から環境を復元し、例外を発生させる。

  

  

    
 fex_merge_flags(envp)

    
 *envp から例外フラグの論理和をとる。

  















   fegetenv 関数は、浮動小数点環境を保存します。fesetenv 関数は、浮動小数点環境を復元します。fesetenv に対する引数は、fegetenv または feholdexcept の呼び出しによって先に保存されている環境のポインタか、fenv.h で定義されている定数 FE_DFL_ENV のどちらかです。FE_DFL_ENV は、すべての例外フラグのクリアー、最近似への丸め (x86 では拡張倍精度への丸めも含む)、連続した例外処理モード (トラップが無効になる)、および無効にされた標準外モード (SPARC の場合) によるデフォルトの環境を表現します。






   feholdexcept 関数は、現在の環境を保存したあとすべての例外をクリアーし、すべての例外に対して連続した例外処理モードを確立します。feupdateenv 関数は、保存されている環境 (fegetenv または feholdexcept の呼び出しによって保存されているもの、または定数 FE_DFL_ENV) を復元し、そのあとで以前の環境でフラグがセットされている例外を発生させます。それらの例外のいずれかに対して有効になったトラップが、復元される環境に含まれる場合は、トラップが発生します。それ以外の場合は、フラグがセットされます。次のコード例に示すように、これらの 2 つの関数はあわせて使用し、例外を発生させそうなサブルーチンを呼び出すことができます。






   



  
  

     

#include <fenv.h>


 


 

void myfunc(...) {

	fenv_t env;


 


 

	/* 環境を保存し、フラグをクリアーし、トラップを無効にします */

	feholdexcept(&env);

	/* 例外を発生させる可能性がある計算を行います */


 

	...


 

	/* 疑似例外を調べます */

	if (fetestexcept(...)) {

			/* 例外を適切に処理し、それらのフラグをクリアーします */

			...

			feclearexcept(...);

	}

	/* 環境を復元し、関連する例外を発生させます */

	feupdateenv(&env);

{



  














   fex_merge_flags 関数は、トラップを発生させることなく、保存済みの環境の例外フラグの論理和をとり現在の環境に入れるだけです。この関数をマルチスレッド対応のプログラムで使用すると、子スレッドでの計算で発生したフラグに関する情報を親スレッドに保持できます。fex_merge_flags の使用方法を示した例は、付録 Aを参照してください。






   libm と libsunmath の実装上の特徴






   この節では、使用可能な libm と libsunmath の実装上の特徴について説明します。






   アルゴリズム






   SPARC システム上の libm および libsunmath 中の基本関数は、テーブル駆動型および多項式有理数の近似値のアルゴリズムによって実装されています。x86 プラットフォーム上の libm および libsunmath 中の基本関数には、x86 命令セット中に提供されている基本関数カーネル命令を使用して実装されているものと、SPARC システムで使用されているのと同じテーブル駆動型および多項式有理数の近似値アルゴリズムによって実装されているものがあります。






   libm 中の共通の基本演算および libsunmath 中の単精度基本演算に対する、テーブル駆動型および多項式有理数の近似値のアルゴリズムは、最後の place 内の 1 単位 
(ulp = unit-in-the-last-place ) 中だけで正確な結果を算出します。SPARC システムでは、libsunmath 中の共通の 4 倍精度基本演算は、1 ulp 中で正確な結果を算出します (expmll および loglpl 演算を除く)。expmll および loglpl 演算は、2 ulp 中で正確な結果を算出します (共通の関数には、指数関数、対数関数、累乗関数、ラジアン引数の循環三角関数などがあります。双曲線三角関数や高度な超越関数のようなその他の関数は、正確度が低くなります)。これらのエラー境界は、アルゴリズムを直接分析して調べることができます。BeEF (Berkley Elementary Function) 検査プログラムを使用して、これらのルーチンの正確度を検査することもできます。BeEF は、Z.Alex Liu によって開発され、ucbtest パッケージの netlib から入手できます (http://www.netlib.org/fp/ucbtest.tgz)。






   三角関数の引数還元






   範囲 [-π/4, π/4] の外側のラジアン引数に対する三角関数は、通常 π/2 の整数倍を引いて、引数を指定された範囲内に還元して計算します。






   π はマシンで表現可能な数ではないので、何らかの方法で近似する必要があります。最終的に計算された三角関数の誤差は、引数還元の丸め誤差 (近似 π と丸めによる) と、還元された引数の三角関数計算の誤差で決まります。相当に大きい引数の場合でも、引数還元に起因する誤差は他の誤差より大きくありませんが、これに反して相当に小さい引数の場合は、最終結果の相対誤差は引数還元誤差によって左右されます。






   マシンで表現可能な大きい数は、π より大きい範囲で分割されているという理由から、大きな引数の三角関数は本来不正確で、小さい引数はすべて比較的正確であるという誤解が広まっています。






   しかし、計算した 三角関数の値が急に大きくなるという固有の境界はなく、不正確な関数値が役に立たないということもありません。引数還元が一様に行われるならば、基本的同一性と関連性が大きい引数に対しても小さい引数と同様に保護されるので、引数還元が π への近似で行われたという事実はほとんど検知できません。






   libm と libsunmath 三角関数は、引数還元のため「無限に」正確な π を使用します。値 2/π は 16 進数 916 桁まで計算され、参照用テーブルに格納されて引数還元中に使用されます。






   関数 sinpi、cospi、および tanpi のグループ (表 3-2 を参照) は、引数の一定範囲への縮小によって不正確になることを避けるため π で入力引数を基準化 (スケール) します。






   データ変換ルーチン






   libm と libsunmath には、IEEE 形式とそれ以外の形式間で 2 進浮動小数点データの変換に使用されるデータ変換ルーチン convert_external があります。






   サポートする形式には、SPARC (IEEE)、IBM PC、VAX、IBM S/370、および Cray が使用する形式が含まれます。






   Cray で生成されたデータを関数 convert_external を使用して、SPARC システムが希望する IEEE 形式に変換する例については、convert_external(3m) のマニュアルページを参照してください。






   乱数発生機構






   32 ビット整数、単精度浮動小数点、および倍精度浮動小数点の各形式で、一様疑似乱数を生成する 3 つの方法を次に示します。






   また、shufrans(3m) のマニュアルページに述べられた関数をこれらの任意のジェネレータとあわせて使用すると、疑似乱数の配列を動かし、アプリケーションにより大きなランダム性を与えることができます (64 ビット整数の配列を動かす方法はありません)。






   各乱数機能には、一度に 1 つ (関数呼び出しごとに 1 つ) の乱数を生成するルーチンと、一度の呼び出しで乱数の配列を生成するルーチンが含まれています。一度に 1 つの乱数を生成する関数は、次の表 3-18 に示す範囲の数値を配布します。



    

   表 3-18 	  単一値乱数の値の発生範囲

  
 関数

     下限

     上限

  

  

    
 i_addran_

    
 -2147483648

    
 2147483647

  

  

    
 r_addran_

    
 0

    
 0.9999999403953552246

  

  

    
 d_addran_

    
 0

    
 0.9999999999999998890

  

  

    
 i_lcran_

    
 1

    
 21477483646

  

  

    
 r_lcran_

    
 4.656612873077392578E-10

    
 1

  

  

    
 d_lcran_

    
 4.656612875245796923E-10

    
 0.9999999995343387127

  

  

    
 i_mwcran_

    
 0

    
 2147483647

  

  

    
 u_mwcran_

    
 0

    
 4294967295

  

  

    
 i_llmwcan_

    
 0

    
 9223372036854775807

  

  

    
 u_llmwcan_

    
 0

    
 18446744073709551615

  

  

    
 r_mwcan_

    
 0

    
 0.9999999403953552246

  

  

    
 d_mwcran_

    
 0

    
 0.9999999999999998890

  















   一度の呼び出しで乱数の配列全体を生成する関数を使用すると、生成される数値の範囲を指定できます。付録 Aでは、異なる間隔で一様に分布される乱数の配列を生成する方法を示したいくつかの例が挙げられています。






   addrans および mwcrans ジェネレータは、lcrans ジェネレータよりも一般に効率的ですが、これらのジェネレータの理論は精緻ではありません。線形合同アルゴリズムの理論的な特性については、"Communications of the ACM" の 1988 年 10 月号に掲載された S. Park と K. Miller による "Random Number Generators: Good Ones Are Hard To Find" で説明されています。また、加法的乱数ジェネレータについては、Knuth の『The Art of Computer Programming』の第 2 版で説明されています。






   

















サン・マイクロシステムズ株式会社

Copyright information. All rights reserved.






ホーム
  |  
目次
  |  
前ページへ
  |  
次ページへ
  |  
索引