数値計算ガイド

付録 E

規格への準拠

Solaris 環境に対応する Sun WorkShop Compiler 言語製品の中のコンパイラ、ヘッダーファイル、ライブラリは、複数の規格、すなわち System V Interface Definition (SVID) 第 3 版、X/Open、および ANSI C をサポートしています。その結果、数学ライブラリ libm とそれに関連したファイルも、C プログラムが各規格に準拠するように修正されました。この修正の変更点は、主として例外処理関係であるため、ユーザーのプログラムは通常は影響を受けません。

SVID の歴史

SVID に従った例外処理と IEEE 規格が表わしている立場との相違点を理解するためには、両者が発展してきた状況を検討してみる必要があります。SVID にある考え方の源は、その多くが UNIX の誕生間もない時期、つまりコンピュータ本体上に初めて実装された時期に発しています。このような初期の環境に共通しているのは、有理浮動小数点演算 +、-、*、/ が不可分 (atomic) な機械命令であること、また一方では sqrt、浮動小数点形式での整数値への変換、ならびに基本超越整関数が多数の不可分な機械命令から成るサブルーチンであることです。

これらの環境では浮動小数点例外を多様な方法で処理しますが、一様性を持たせようとすると、不可分 (atomic) な各浮動小数点命令の前後にソフトウェア内で引数と結果をチェックしなければ実現できないと考えられます。しかし、これを行うと性能に及ぼす影響が大きくなりすぎると考えられるので、SVID ではゼロによる除算やオーバーフローなど浮動小数点例外の影響は明示していません。

サブルーチンによって実現される演算は、単一の不可分 (atomic) な浮動小数点命令に比べると速度が劣ります。引数と結果について特別なエラーチェックを行なっても性能にはほとんど影響がないので、SVID ではそのようなチェックを必須にしています。例外が検出されると、デフォルト結果が指定され、不適格なオペランドの場合には errno が EDOM に設定されます。またオーバーフロー、あるいはアンダーフローする結果が出る場合には、errno が ERANGE に設定されます。さらに、例外の詳細を含むレコード付きで関数 matherr() が呼び出されます。このことは、UNIX が開発された当初に対象としたマシンにはほとんど負担をかけませんが、基本的な演算 +、-、*、/ における一般的な例外が全く未指定であるため、価値はそれ相応に小さいものとなります。

IEEE 754 の歴史

IEEE 規格は、以前の実装との互換性は目標でなかったと明白に述べています。代わりに、効率とユーザーの要求内容とを念頭に置いて例外処理の機構が開発されました。この機構は、単純な有理演算 ( +、-、*、/ ) と、さらに複雑な剰余、平方根、フォーマット間変換などの演算との両方を通じて一様です。規格では超越関数については規定していませんが、規格の創設者は、準拠システム内の基本超越整関数にも同じ例外処理機構が適用されることを期待していました。

IEEE 例外処理の要素には、あらかじめ要求された場合にのみ、適当なデフォルト結果と計算の中断が含まれます。

SVID の将来の方向

現在の SVID (第 3 版または SVR4) では、将来の発展について一定の方向が明らかにされています。方向の 1 つは IEEE 規格との互換性です。特に、SVID の将来のバージョンにより、大きな有限数用に用意された HUGE への参照が、IEEE システム上での無限大である HUGE_VAL で置換されます。たとえば HUGE_VAL は、浮動小数点オーバーフローの結果として返されます。例外を発生させる入力引数について libm 関数が返す値は、後出の表 E-1 の IEEE 欄に示すものとなります。errno は今後設定される必要がなくなります。

SVID の実装

以下の表に示す libm 関数により、SVID に対応するオペランドチェックまたは結果チェックが行われます。-xlibmil 経由で libm のインライン展開テンプレートを使用する C プログラムから呼び出されたとき、平方根用のハードウェア命令 fsqrt[sd] が関数コールの代わりに使用されるため、sqrt 関数は SVID に準拠しない唯一の関数です。

表 E-1   例外のケースと libm 関数 

関数	errno	エラー メッセージ	SVID	X/Open	IEEE
acos(|x|>1)	EDOM	DOMAIN	0.0	0.0	NaN
acosh(x<1)	EDOM	DOMAIN	NaN	NaN	NaN
asin(|x|>1)	EDOM	DOMAIN	0.0	0.0	NaN
atan2((+-0,+-0)	EDOM	DOMAIN	0.0	0.0	+-0.0,+-pi
atanh(|x|>1)	EDOM	DOMAIN	NaN	NaN	NaN
atanh(+-1)	EDOM/ERANGE	SING	+-HUGE (EDOM)	+-HUGE_VAL (ERANGE)	+-infinity
cosh overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
exp overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
exp underflow	ERANGE	-	0.0	0.0	0.0
fmod(x,0)	EDOM	DOMAIN	x	NaN	NaN
gamma(0 or -integer)	EDOM	SING	HUGE	HUGE_VAL	infinity
gamma overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
hypot overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
j0(|x| > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
j1(|x| > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
jn(|x| > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
lgamma (0 or -integer)	EDOM	SING	HUGE	HUGE_VAL	infinity
lgamma overflow	ERANGE	-	HUGE	HUGE_VAL	infinity
log(0)	EDOM/ERANGE	SING	-HUGE (EDOM)	-HUGE_VAL (ERANGE)	infinity
log(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
log10(0)	EDOM/ERANGE	SING	-HUGE (EDOM)	-HUGE_VAL (ERANGE)	-infinity
log10(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
loglp(-1)	EDOM/ERANGE	SING	-HUGE (EDOM)	-HUGE_VAL (ERANGE)	-infinity
loglp(x<-1)	EDOM	DOMAIN	NaN	NaN	NaN
pow(0,0)	EDOM	DOMAIN	0.0	1.0 (no error)	1.0 (no error)
pow(NaN,0)	EDOM	DOMAIN	NaN	NaN	1.0 (no error)
pow(0,neg)	EDOM	DOMAIN	0.0	-HUGE_VAL	+-infinity
pow (neg, non-integer)	EDOM	DOMAIN	0.0	NaN	NaN
pow overflow	ERANGE	-	+-HUGE	+-HUGE_VAL	+-infinity
pow underflow	ERANGE	-	+-0.0	+-0.0	+-0.0
remainder(x,0)	EDOM	DOMAIN	NaN	NaN	NaN
scalb overflow	ERANGE	-	+-HUGE_VAL	+-HUGE_VAL	+-infinity
scalb underflow	ERANGE	-	+-0.0	+-0.0	+-0.0
sinh overflow	ERANGE	-	+-HUGE	+-HUGE_VAL	+-infinity
sqrt(x<0)	EDOM	DOMAIN	0.0	NaN	NaN
y0(0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	-infinity
y0(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
y0(x > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
y1(0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	-infinity
y1(x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
y1(x > X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer
yn(n,0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	-infinity
yn(n,x<0)	EDOM	DOMAIN	-HUGE	-HUGE_VAL	NaN
yn(n, x> X_TLOSS)	ERANGE	TLOSS	0.0	0.0	correct answer

例外のケースと libm 関数についての一般的注意事項

表 E-1 には、各規格の影響を受ける libm 関数がすべてリストされています。値 X_TLOSS は、<values.h> に定義されています。SVID では、<math.h> に対し HUGE を MAXFLOAT と定義するように要求しています。これはおよそ 3.4e+38 です。HUGE_VAL は libc で無限大と定義されています。errno は、C プログラムおよび C++ プログラムにアクセス可能なグローバル変数です。

<errno.h> では、errno 用として考えられる値を 120 個程度定義しています。数学ライブラリが使用するものが 2 つあります。ドメインエラー用の EDOM と範囲エラー用の ERANGE です。intro(3) と perror(3) を参照してください。

• ANSI C コンパイラは、-Xt、-Xa、-Xc、-Xs を切り換え、特にコンパイラによって実施される規格準拠レベルを管理します。これらのスイッチの詳細については、
  cc(1) を参照してください。
• libm に関する限り、-Xt と -Xa によってそれぞれ SVID と X/Open が動作します。-Xc は厳密な ANSI C 動作に相当します。

  付加スイッチ -libmieee は、指定されると、IEEE 754 に従って値の返却を行います。libm と、libsunmath のデフォルト動作は、Solaris 2.6、Solaris 7 および Solaris 8 上では SVID 準拠になっています。

• 厳密な ANSI C (-Xc) の場合、errno は必ず設定され、matherr() の呼び出しは行われず、X/Open 値が返されます。
• SVID (-Xt または -Xs) の場合、関数 matherr() が例外についての情報付きで呼び出されます。この情報にはデフォルトの SVID 戻り値となる値が含まれています。

  ユーザーが与えた matherr() により、戻り値が変化することもあります。
  matherr(3m) を参照してください。ユーザーが与えた matherr() がなければ、libm は errno の設定を行い、メッセージを標準エラー出力ファイルに出力し、表 E-1 内で SVID 欄に示されている値を返します。

• X/Open (-Xa) の場合、動作は matherr() が起動され、それに合わせて errno が設定される点で SVID の場合と同じです。ただし、標準エラー出力ファイルにはエラーメッセージの出力はなく、また多くの場合、X/Open の戻り値は IEEE の戻り値と同じです。
• libm の例外処理に対して、-Xs は -Xt と同じ動作を行います。すなわち、-Xs でコンパイルされたプログラムは、表 E-1 に示されている libm 関数の SVID 準拠バージョンを使用します。
• sqrt が負の引数に出会った場合、インラインハードウェア浮動小数点付きでコンパイルされたプログラムは、効率上の観点から EDOM を設定したり matherr() を呼び出したりするため必要となる特別なチェックは行いません。そのままでは EDOM が設定される可能性がある場合には、関数値用に NaN が返されます。

  sqrt() 関数コールを fsqrt[sd] 命令で置換する C プログラムは、IEEE 浮動小数点規格には準拠しますが、System V Interface Definition のエラー処理要求には今後準拠しなくなる可能性があります。

libm についての注意

SVID では、PLOSS (Partial Loss of Significance) と TLOSS (Total Loss of Significance) の 2 つの浮動小数点例外を規定しています。sqrt(-1) と異なり、これらには固有の数学的意味がありません。また、exp(+-10000) と異なり、これらは浮動小数点記憶形式の固有の制限を反映しません。

その代わり PLOSS と TLOSS は、fmod 用の特定アルゴリズム、および明確な境界により不意に正確度がおちる三角関数用の特定アルゴリズムの制限を反映します。

libm のアルゴリズムは、ほとんどの IEEE の実装と同様、そのような不意の低下を被ることがなく、PLOSS のシグナルを出すことはありません。SVID 準拠要件を満たすために、ベッセル関数では、正確な結果を安全に計算できるにもかかわらず、大きい入力引数については TLOSS のシグナルを出します。

sin、cos、および tan の実装では、無限大での本質的な特異点を、無限引数について EDOM を設定して NaN を返すことにより、ほかの本質的な特異点と同様に処理します。

同様に SVID では、x/y がオーバーフローすると考えられる場合には、fmod(x,y) は 0 でなくてはならないことを規定しています。しかし、IEEE 剰余関数から派生した fmod の libm の実装では、x/y を明示的に計算せず、必ず厳密な結果をもたらします。

LIA-1 準拠

ここでは、LIA-1 は ISO/IEC 10967-1:1994 Information Technology Language Independent Arithmetic Part 1 のことを差します。言語に依存しない数値演算についての基準です。

C コンパイラ (cc) および FORTRAN 77 コンパイラ (f77) は、Sun WorkShop Compiler 6.0 に含まれており、以下の点で LIA-1 に準拠しています。
行頭のアルファベットは、LIA-1 の第 8 節で使用されているものと対応しています。

a. データ型 (LIA 5.1)

LIA-1 準拠している型には、C の int および FORTRAN の INTEGER があります。この他にも LIA-1 準拠の型はありますが、ここでは取り上げません。その他特定の言語に対する仕様は、言語が LIA-1 に準拠してから決められます。

b. パラメータ (LIA 5.1)

#include <values.h> defines MAXINT #define TRUE 1 #define FALSE 0 #define BOUNDED TRUE #define MODULO TRUE #define MAXINT 2147483647 #define MININT -2147483648 logical bounded, modulo integer maxint, minint parameter (bounded = .TRUE.) parameter (modulo = .TRUE.)

#include <values.h> defines MAXINT

d. DIV/REM/MOD (5.1.3)

C の / および % と FORTRAN の / および mod() によって、DIVtI(x,y) と REMtI(x,y) が提供されます。また、modaI(x,y) は、以下のコードによって使用できます。

int modaI(int x, int y){ int t = x % y; if (y < 0 && t > 0) t -= y; else if (y > 0 && t < 0) t += y; return t; }

または

integer function modaI (x, y) integer x, y, t t = mod(x,y) if (y .lt. 0 .and. t .gt. 0) t = t - y if (y .gt. 0 .and. t .lt. 0) t = t + y modaI = t return end

i. 記数法 (LIA 5.1.3)

LIA-1 の整数演算で認識される記数法を示します。

表 E-2   LIA-1 準拠の記数法  

LIA	C	FORTRAN (C と異なる場合)
addI(x,y)	x+y
subI(x,y)	x-y
mulI(x,y)	x*y
divtI(x,y)	x/y
remtI(x,y)	x%y	mod(x,y)
modaI(x,y)	x%y	mod(x,y)
negI(x)	-x
absI(x)	abs(x) #include <stdlib.h>	abs(x)
signI(x)	#define signI((x) (x>0?1:(x<0?-1:0))	以下を参照
eqI(x,y)	x==y	x.eq.y
neqI(x,y)	x!=y	x.ne.y
lssI(x,y)	x<y	x.lt.y
leqI(x,y)	x<=y	x.le.y
gtrI(x,y)	x>y	x.gt.y
geqI(x,y)	x>=y	x.ge.y

signI(x) の FORTRAN でのコード例を示します。

integer function signi(x) integer x, t if (x .gt. 0) t=1 if (x .lt. 0) t=-1 if (x .eq. 0) t=0 return end

j. 式の評価

デフォルトでは、最適化が指定されていない場合は、式は C の場合は int、FORTRAN の場合は INTEGER の精度で評価されます。括弧も評価されます。a+b+c または a*b*c などの、括弧で囲まれていない結合式の評価順序は、指定されていません。

k. パラメータの受け取り方

ソースコード中にある上記の定義をインクルードします。

n. 通知

整数の例外は x/0、x%0、mod(x,0) です。デフォルトでは、これらの例外が SIGFPE を生成します。シグナルハンドラが指定されていない場合は、プロセスを終了してメモリーダンプを行います。

o. 選択のしくみ

signal(3) および signal(3F) が使用され、ユーザーが SIGFPE に対する例外処理を行うことができるようになります。