数値計算ガイド

付録 A

例

この付録では、よく行われる作業の例を示します。例は FORTRAN、および ANSI C で書かれており、その多くは現バージョンの libm と libsunmath に依存しています。これらの例は、Solaris オペレーティングシステム上の C および FORTRAN コンパイラでテストされています。

IEEE の演算機能

浮動小数点数の 16 進表現の検証方法を例で示します。格納されたデータの 16 進表現を見るには、デバッガを使用することもできます。

以下の C のプログラムでは、倍精度の近似値を π と単精度の無限大に出力します。

コード例 A-1 倍精度の例

#include <math.h> #include <sunmath.h> main() { union { float flt; unsigned un; } r; union { double dbl; unsigned un[2]; } d; /* 倍精度 */ d.dbl = M_PI; (void) printf("DP Approx pi = %08x %08x = %18.17e \n", d.un[0], d.un[1], d.dbl); /* 単精度 */ r.flt = infinityf(); (void) printf("Single Precision %8.7e : %08x \n", r.flt, r.un); exit(0); }

SPARC では、上記のプログラムの出力は次のようになります。

DP Approx pi = 400921fb 54442d18 = 3.14159265358979312e+00 Single Precision Infinity : 7f800000

次の FORTRAN プログラムは、最小の正規数をそれぞれの形式で出力します。

コード例 A-2 最小の正規数をそれぞれの形式で出力

program print_ieee_values c c implicit 文は、f77_floatingpoint 疑似組み込み関数が、 c 正しい型で宣言されているかを確認します。 c implicit real*16 (q) implicit double precision (d) implicit real (r) real*16 z double precision x real r c z = q_min_normal() write(*,7) z, z 7 format('min normal, quad: ',1pe47.37e4,/,' in hex ',z32.32) c x = d_min_normal() write(*,14) x, x 14 format('min normal, double: ',1pe23.16,' in hex ',z16.16) c r = r_min_normal() write(*,27) r, r 27 format('min normal, single: ',1pe14.7,' in hex ',z8.8) c end

SPARC では、対応する出力が次のようになります。

min normal, quad: 3.3621031431120935062626778173217526026E- 4932 in hex 00010000000000000000000000000000 min normal, double:2.2250738585072014-308 in hex 0010000000000000 min normal, single:1.1754944E-38 in hex 00800000

(x86 では、FORTRAN コンパイラは real*16 型をサポートしません。前記の例を x86 で実行する場合は、real*16 宣言と、4 倍精度の値の計算と出力を行うコードを削除してください。)

数学ライブラリ

この節では、数学ライブラリの関数を使用した例を示します。

乱数生成

次の例では、数の配列を生成する乱数関数を呼び出し、与えられた数の EXP を計算するのにかかる時間を測定する関数を使用します。

コード例 A-3 乱数生成 

#ifdef DP #define GENERIC double precision #else #define GENERIC real #endif #define SIZE 400000 program example c implicit GENERIC (a-h,o-z) GENERIC x(SIZE), y, lb, ub real time(2), u1, u2 c c [-ln2/2,ln2/2] における乱数での EXP を計算 lb = -0.3465735903 ub = 0.3465735903 c c 乱数の配列を生成 #ifdef DP call d_init_addrans() call d_addrans(x,SIZE,lb,ub) #else call r_init_addrans() call r_addrans(x,SIZE,lb,ub) #endif c c クロック開始 call dtime(tarray) u1 = tarray(1) c c 指数を計算 do 16 i=1,SIZE y = exp(x(i)) 16 continue c c 経過時間の取得 call dtime(time) u2 = tarray(1) print *,'time used by EXP is ',u2-u1 print *,'last values for x and exp(x) are ',x(SIZE),y c call flush(6) end

前記の例をコンパイルするには、コンパイラが自動的にプリプロセッサを呼び出すようにソースコードを接尾辞 F (f ではない) のついたファイルに入れ、コマンド行で -DSP または -DDP を指定して単精度または倍精度を選択します。

以下の例では、d_addrans を使用して、ユーザーが指定した範囲で均一に分散する乱数データのブロックを発生させる方法を示します。

コード例 A-4 d_addrans を使用する 

/* * sin への SIZE*LOOPS 乱数引数をテストします。 * 範囲は [0、しきい値 ] とし、 * しきい値 = 3E30000000000000 (3.72529029846191406e-09) とします。 */ #include <math.h> #include <sunmath.h> #define SIZE 10000 #define LOOPS 100 int main() { double x[SIZE], y[SIZE]; int i, j, n; double lb, ub; union { unsignedu[2]; doubled; } upperbound; upperbound.u[0] = 0x3e300000; upperbound.u[1] = 0x00000000; /* 乱数関数を初期化 */ d_init_addrans_(); /* sin について (SIZE * LOOPS) 個の引数をテスト */ for (j = 0; j < LOOPS; j++) { /* * 長さ SIZE の乱数のベクトル x を生成し、 * 三角関数の関数への入力として使用 */ n = SIZE; ub = upperbound.d; lb = 0.0; d_addrans_(x, &n, &lb, &ub); for (i = 0; i < n; i++) y[i] = sin(x[i]); /*sin(x) == x ？　 x が小さい場合になる */ for (i = 0; i < n; i++) if (x[i] != y[i]) printf( " OOPS: %d sin(%18.17e)=%18.17e \n", i, x[i], y[i]); } printf(" comparison ended; no differences\n"); ieee_retrospective_(); return(0); }

IEEE が推奨する関数

次の FORTRAN の例では、IEEE 標準が推奨する関数を使用しています。

コード例 A-5 IEEE が推奨する関数 

c

c このプログラムでは、IEEE 推奨の関数のうちの 5 つを FORTRAN から呼び出す

c 方法を示します。

c たとえば、y=x*2**n を行うには、ハードウェアでは数値が基数 2 で記憶され

c るため、指数を n だけシフトします。

c 浮動小数点の演算機能に関する IEEE 規格では、これらの関数を必須にはしてい

c ませんが、IEEE 浮動小数点環境には含めることを勧めています。

c

c 以下の項目も参照してください。

c

c ieee_functions(3m)

c libm_double(3f)

c libm_single(3f)

c

c ここでは次の 5 つの関数について例を示します：

c

c ilogb(x): x の基数 2 でのバイアスされていない指数を整数形式で返す

c signbit(x): 0 または 1 の符号ビットを返す

c copysign(x,y): y の符号ビットを付けた x を返す

c nextafter(x,y): x から y 方向に次に出現する表現可能な数値を返す

c scalbn(x,n): x * 2**n

c

c 関数 倍精度 単精度

c --------------------------------------------------------

c ilogb(x) i = id_ilogb(x) i = ir_ilogb(r)

c signbit(x) i = id_signbit(x) i = ir_signbit(r)

c copysign(x,y) x = d_copysign(x,y) r = r_copysign(r,s)

c nextafter(x,y) z = d_nextafter(x,y) r = r_nextafter(r,s)

c scalbn(x,n) x = d_scalbn(x,n) r = r_scalbn(r,n)

c

program ieee_functions_demo

implicit double precision (d)

implicit real (r)

double precision x, y, z, direction

real r, s, t, r_direction

integer i, scale

print *

print *, 'DOUBLE PRECISION EXAMPLES:'

print *

x = 32.0d0

i = id_ilogb(x)

write(*,1) x, i

1 format(' The base 2 exponent of ', F4.1, ' is ', I2)

x = -5.5d0

y = 12.4d0

z = d_copysign(x,y)

write(*,2) x, y, z

2 format(F5.1, ' was given the sign of ', F4.1,

* ' and is now ', F4.1)

x = -5.5d0

i = id_signbit(x)

print *, 'The sign bit of ', x, ' is ', i

x = d_min_subnormal()

direction = -d_infinity()

y = d_nextafter(x, direction)

write(*,3) x

3 format(' Starting from ', 1PE23.16E3,

- ', the next representable number ')

write(*,4) direction, y

4 format(' towards ', F4.1, ' is ', 1PE23.16E3)

x = d_min_subnormal()

direction = 1.0d0

y = d_nextafter(x, direction)

write(*,3) x

write(*,4) direction, y

x = 2.0d0

scale = 3

y = d_scalbn(x, scale)

write (*,5) x, scale, y

5 format(' Scaling ', F4.1, ' by 2**', I1, ' is ', F4.1)

print *

print *, 'SINGLE PRECISION EXAMPLES:'

print *

r = 32.0

i = ir_ilogb(r)

write (*,1) r, i

r = -5.5

i = ir_signbit(r)

print *, 'The sign bit of ', r, ' is ', i

r = -5.5

s = 12.4

t = r_copysign(r,s)

write (*,2) r, s, t

r = r_min_subnormal()

r_direction = -r_infinity()

s = r_nextafter(r, r_direction)

write(*,3) r

write(*,4) r_direction, s

r = r_min_subnormal()

r_direction = 1.0e0

s = r_nextafter(r, r_direction)

write(*,3) r

write(*,4) r_direction, s

r = 2.0

scale = 3

s = r_scalbn(r, scale)

write (*,5) r, scale, y

print *

end

このプログラムからの出力は次のようになります。

コード例 A-6 コード例 A-5 のプログラムからの出力

DOUBLE PRECISION EXAMPLES: The base 2 exponent of 32.0 is 5 -5.5 was given the sign of 12.4 and is now 5.5 The sign bit of -5.5000000000000 is 1 Starting from 4.9406564584124654-324, the next representable number towards -Inf is 0.0000000000000000E+000 Starting from 4.9406564584124654-324, the next representable number towards 1.0 is 9.8813129168249309-324 Scaling 2.0 by 2**3 is 16.0 SINGLE PRECISION EXAMPLES: The base 2 exponent of 32.0 is 5 The sign bit of -5.50000 is 1 -5.5 was given the sign of 12.4 and is now 5.5 Starting from 1.4012984643248171E-045, the next representable number towards -Inf is 0.0000000000000000E+000 Starting from 1.4012984643248171E-045, the next representable number towards 1.0 is 2.8025969286496341E-045 Scaling 2.0 by 2**3 is 16.0 Note:IEEE floating-point exception flags raised: Inexact; Underflow; See the Numerical Computation Guide, ieee_flags(3M) [日本語訳] 注: 以下の IEEE 浮動小数点例外が発生しました: 不正確、アンダーフロー 詳細は、『数値計算ガイド』の ieee_flags(3M)に関する説明を参照してくださ い。

IEEE の特殊な値

以下の C プログラムでは、ieee_values(3m) の関数を呼び出しています。

#include <math.h> #include <sunmath.h> int main() { double x; float r; x = quiet_nan(0); printf("quiet NaN: %.16e = %08x %08x \n", x, ((int *) &x)[0], ((int *) &x)[1]); x = nextafter(max_subnormal(), 0.0); printf("nextafter(max_subnormal,0) = %.16e\n",x); printf(" = %08x %08x\n", ((int *) &x)[0], ((int *) &x)[1]); r = min_subnormalf(); printf("single precision min subnormal = %.8e = %08x\n", r, ((int *) &r)[0]); return 0; }

リンクするときには、-lsunmath と -lm の両方を指定してください。

SPARC では、出力は次のようになります。

quiet NaN: NaN = 7fffffff ffffffff nextafter(max_subnormal,0) = 2.2250738585072004e-308 = 000fffff fffffffe single precision min subnormal = 1.40129846e-45 = 00000001

x86 アーキテクチャは「リトルエンディアン」であるため、x86 での出力は少々異なります (倍精度数の 16 進表現で上位と下位のワードが逆になります)。

quiet NaN: NaN = ffffffff 7fffffff nextafter(max_subnormal,0) = 2.2250738585072004e-308 = fffffffe 000fffff single precision min subnormal = 1.40129846e-45 = 00000001

ieee_values 関数を使用する FORTRAN プログラムでは、関数の型を宣言するように注意しなければなりません。

program print_ieee_values c c implicit 文は、f77_floatingpoint の疑似組み込み関数が、 c 正しい型で宣言されているかどうかを確認します。 c implicit real*16 (q) implicit double precision (d) implicit real (r) real*16 z, zero, one double precision x real r c zero = 0.0 one = 1.0 z = q_nextafter(zero, one) x = d_infinity() r = r_max_normal() c print *, z print *, x print *, r c end

SPARC では、出力は次のようになります。

6.4751751194380251109244389582276466-4966 Infinity 3.40282E+38

x86 の FORTRAN コンパイラは、real*16 型をサポートしないことに注意してください。前記の例を x86 で実行するには、real*16 の変数および関数に対する参照をすべて削除する必要があります。

ieee_flags - 丸め方向

以下の例は、丸めモードをゼロ切り捨てモードに設定する方法を示しています。

#include <math.h> #include <sunmath.h> int main() { int i; double x, y; char *out_1, *out_2, *dummy; /* 一般的な丸めの方向を取得する */ i = ieee_flags("get", "direction", "", &out_1); x = sqrt(.5); printf("With rounding direction %s, \n", out_1); printf("sqrt(.5) = 0x%08x 0x%08x = %16.15e\n", ((int *) &x)[0], ((int *) &x)[1], x); /* 丸め方向の設定 */ if (ieee_flags("set", "direction", "tozero", &dummy) != 0) printf("Not able to change rounding direction!\n"); i = ieee_flags("get", "direction", "", &out_2); x = sqrt(.5); /* * printf も現在の丸め方向に影響を受けるため、 * printf の前に元の丸め方向を復元する */ if (ieee_flags("set", "direction", out_1, &dummy) != 0) printf("Not able to change rounding direction!\n"); printf("\nWith rounding direction %s,\n", out_2); printf("sqrt(.5) = 0x%08x 0x%08x = %16.15e\n", ((int *) &x)[0], ((int *) &x)[1], x); return 0; }

SPARC では、この短いプログラムの出力は、ゼロ切り捨てモードの効果を示します。

demo% cc rounding_direction.c -lm demo% a.out With rounding direction nearest, sqrt(.5) = 0x3fe6a09e 0x667f3bcd = 7.071067811865476e-01 With rounding direction tozero, sqrt(.5) = 0x3fe6a09e 0x667f3bcc = 7.071067811865475e-01 demo%

x86 では、この短いプログラムの出力は、ゼロ切り捨てモードの効果を示します。

demo% cc rounding_direction.c -lm demo% a.out With rounding direction nearest, sqrt(.5) = 0x667f3bcd 0x3fe6a09e = 7.071067811865476e-01 With rounding direction tozero, sqrt(.5) = 0x667f3bcc 0x3fe6a09e = 7.071067811865475e-01 demo%

FORTRAN プログラムでゼロ切り捨てモードを設定します。

program ieee_flags_demo character*16 out i = ieee_flags("set", "direction", "tozero", out) if (i.ne.0) print *, 'not able to set rounding direction' i = ieee_flags("get", "direction", "", out) print *, "Rounding direction is: ", out end

出力は以下のようになります。

Rounding direction is: tozero Note: Rounding direction toward zero. See the Numerical Computation Guide, ieee_flags(3M) [日本語訳] 注: ゼロ方向への丸め 詳細は、『数値計算ガイド』の ieee_flags(3M)に関する説明を参照してくださ い。

C99 浮動小数点環境関数

次に、C99 浮動小数点環境関数をいくつか使用した例を示します。norm 関数は、アンダーフローとオーバーフローを処理するため、ベクトルのユークリッド正規形を計算し、C99 浮動小数点環境関数を使用します。メインプログラムは、遡及診断出力が示すように、アンダーフローとオーバーフローを発生させるようにスケールされたベクトルを使用してこの関数を呼び出します。

コード例 A-7 C99 浮動小数点環境関数 

#include <stdio.h>

#include <math.h>

#include <sunmath.h>

#include <fenv.h>

/*

* 早すぎるアンダーフローまたはオーバーフローを避けるベクトル x の

* ユークリッド正規形を計算します

*/

double norm(int n, double *x)

{

fenv_t env;

double s, b, d, t;

int i, f;

/* 環境を保存してフラグをクリアし、連続した例外処理を確立します */

feholdexcept(&env);

/* ドット積 x.x の計算を試みます */

d = 1.0; /* 桁移動子 */

s = 0.0;

for (i = 0; i < n; i++)

s += x[i] * x[i];

/* アンダーフローまたはオーバーフローを調べます */

f = fetestexcept(FE_UNDERFLOW | FE_OVERFLOW);

if (f & FE_OVERFLOW) {

/* 最初の試行がオーバーフローしたら、スケーリングを 小さくしてもう一度実行してください */

feclearexcept(FE_OVERFLOW);

b = scalbn(1.0, -640);

d = 1.0 / b;

s = 0.0;

for (i = 0; i < n; i++) {

t = b * x[i];

s += t * t;

}

}

else if (f & FE_UNDERFLOW && s < scalbn(1.0, -970)) {

/* 最初の試行がアンダーフローしたら、スケーリングを 大きくしてもう一度実行してください */

b = scalbn(1.0, 1022);

d = 1.0 / b;

s = 0.0;

for (i = 0; i < n; i++) {

t = b * x[i];

s += t * t;

}

}

/* これまでに発生したアンダーフローをすべて隠します */

feclearexcept(FE_UNDERFLOW);

/* 環境を復元し、これまでに起きたほかの例外を発生させます */

feupdateenv(&env);

/* 平方根を取得し、スケーリングを解除します */

return d * sqrt(s);

}

int main()

{

double x[100], l, u;

int n = 100;

fex_set_log(stdout);

l = 0.0;

u = min_normal();

d_lcrans_(x, &n, &l, &u);

printf("norm: %g\n", norm(n, x));

l = sqrt(max_normal());

u = l * 2.0;

d_lcrans_(x, &n, &l, &u);

printf("norm: %g\n", norm(n, x));

return 0;

}

SPARC でこのプログラムをコンパイルして実行すると、次のように出力されます。

demo% cc norm.c -R/opt/SUNWspro/lib -L/opt/SUNWspro/lib -lm9x -lsunmath -lm demo% a.out Floating point underflow at 0x000153a8 __d_lcrans_, nonstop mode 0x000153b4 __d_lcrans_ 0x00011594 main Floating point underflow at 0x00011244 norm, nonstop mode 0x00011248 norm 0x000115b4 main norm: 1.32533e-307 Floating point overflow at 0x00011244 norm, nonstop mode 0x00011248 norm 0x00011660 main norm: 2.02548e+155

次の例は、x86 での fesetprec 関数の結果を示しています (この関数は SPARC では使用できません)。while ループは、1 に加えられるときに完全に丸められる、2 の最大の累乗を探すことにより、使用できる精度の決定を試みています。最初のループが示すように、すべての中間結果を拡張倍精度で評価する x86 のようなアーキテクチャでは、この手法が期待どおりに動作しない場合があります。そのため、2 番目のループが示すように、fesetprec 関数を使用して、すべての結果が希望する精度に丸められるように指定できます。

コード例 A-8 fesetprec 関数 (x86)  

#include <math.h>

#include <fenv.h>

int main()

{

double x;

x = 1.0;

while (1.0 + x != 1.0)

x *= 0.5;

printf("%d significant bits\n", -ilogb(x));

fesetprec(FE_DBLPREC);

x = 1.0;

while (1.0 + x != 1.0)

x *= 0.5;

printf("%d significant bits\n", -ilogb(x));

return 0;

}

x86 システムでは、このプログラムの出力は次のようになります。

64 significant bits 53 significant bits

次のコードフラグメントは、マルチスレッド対応のプログラムで環境関数を使用して、親スレッドから子スレッドに浮動小数点モードを伝え、子スレッドが親スレッドに結合されるときに子スレッド内で発生する例外フラグを回復する方法を示しています (マルチスレッド対応のプログラムの記述については、『Solaris でのマルチスレッドのプログラミング』を参照)。

コード例 A-9 マルチスレッド対応のプログラムで環境関数を使用する  

#include <thread.h>

#include <fenv.h>

fenv_t env;

void child(void *p)

{

/* 入口で親の環境を継承します */

fesetenv(&env);

...

/* 終了前に子の環境を保存します */

fegetenv(&env);

}

void parent()

{

thread_t tid;

void *arg;

...

/* 子を作成する前に親の環境を保存します */

fegetenv(&env);

thr_create(NULL, NULL, child, arg, NULL, &tid);

...

/* 子と結合します */

thr_join(tid, NULL, &arg);

/* 子で発生した例外フラグを親の環境にマージします */

fex_merge_flags(&env);

...

}

例外と例外処理

ieee_flags - 例外

一般的には、例外ビットの検証やクリアはユーザープログラムで行います。自然発生した例外フラグを検証する C プログラムを示します。

コード例 A-10 例外ビットの検証  

#include <sunmath.h> #include <sys/ieeefp.h> int main() { int code, inexact, division, underflow, overflow, invalid; double x; char *out; /* アンダーフロー例外を発生させる */ x = max_subnormal() / 2.0; /* この文は、前の文が最適化されていないことを示す */ printf("x = %g\n",x); /* どの例外が発生したかを判定する */ code = ieee_flags("get", "exception", "", &out); /* 戻り値を解読する */ inexact = (code >> fp_inexact) & 0x1; underflow = (code >> fp_underflow) & 0x1; division = (code >> fp_division) & 0x1; overflow = (code >> fp_overflow) & 0x1; invalid = (code >> fp_invalid) & 0x1; /* "out" は発生した例外のうちで最も優先度が高い */ printf(" Highest priority exception is: %s\n", out); /* 1 は例外が発生したことを示し、 */ /* 0 は例外が発生していないことを示す */ printf("%d %d %d %d %d\n", invalid, overflow, division, underflow, inexact); ieee_retrospective_(); return(0); }

このプログラムを実行した結果の出力です。

demo% a.out x = 1.11254e-308 Highest priority exception is: underflow 0 0 0 1 1 Note:IEEE floating-point exception flags raised: Inexact; Underflow; See the Numerical Computation Guide, ieee_flags(3M) [日本語訳] 注: 以下の IEEE 浮動小数点例外が発生しました: 不正確、アンダーフロー 詳細は、『数値計算ガイド』の ieee_flags(3M)に関する説明を参照してくださ い。

FORTRAN でも同じことが行えます。

コード例 A-11 例外ビットの検証 (FORTRAN) 

/* これは、次のような FORTRAN プログラム例 です： * アンダーフロー例外を発生させる * ieee_flags を使用して、どの例外が発生したかを判定する * ieee_flags が返す整数値を解読する * 未処理の例外をすべてクリアにする C のプリプロセッサが起動され、ヘッダーファイル f77_floatingpoint.h が取 り込まれるようにするために、このプログラムは接尾辞 .F のファイルに保存してく ださい。 */ #include <f77_floatingpoint.h> program decode_accrued_exceptions double precision x integer accrued, inx, div, under, over, inv character*16 out double precision d_max_subnormal c アンダーフロー例外を発生させる x = d_max_subnormal() / 2.0 c どの例外が発生したかを判定する accrued = ieee_flags('get', 'exception', '', out) c ビットシフトの組み込みを使用して、ieee_flags が返した値を解読する inx = and(rshift(accrued, fp_inexact) , 1) under = and(rshift(accrued, fp_underflow), 1) div = and(rshift(accrued, fp_division) , 1) over = and(rshift(accrued, fp_overflow) , 1) inv = and(rshift(accrued, fp_invalid) , 1) c 最高の優先度をもつ例外が "out" に返される print *, "Highest priority exception is ", out c 1 は例外が発生したことを示し、0 は例外が発生していないことを示す print *, inv, over, div, under, inx c 未処理の例外をすべてクリアにする i = ieee_flags('clear', 'exception', 'all', out) end

この FORTRAN プログラムの出力です。

Highest priority exception is underflow 0 0 0 1 1

通常、ユーザープログラムで例外フラグを設定することはありませんが、できないことではありません。以下の C の例でそれを示します。

#include <sunmath.h> int main() { int code; char *out; if (ieee_flags("clear", "exception", "all", &out) != 0) printf("could not clear exceptions\n"); if (ieee_flags("set", "exception", "division", &out) != 0) printf("could not get exception\n"); code = ieee_flags("get", "exception", "", &out); printf("out is: %s , fp exception code is: %X \n", out, code); return(0); }

SPARC では、前記のプログラムの出力は次のようになります。

out is: division , fp exception code is: 2

x86 では、出力は次のようになります。

out is: division , fp exception code is: 4

ieee_handler - 例外のトラップ

---

注 - 以下の例は、Solaris オペレーティング環境だけに適用できます。

---

SPARC では、例外を特定するためのシグナルハンドラを設定する FORTRAN および C プログラムの例を示します。

FORTRAN プログラム例

コード例 A-12 アンダーフローでのトラップ (SPARC) 

program demo

c シグナルハンドラ関数の宣言

external fp_exc_hdl

double precision d_min_normal

double precision x

c シグナルハンドラの設定

i = ieee_handler(`set', `common', fp_exc_hdl)

if (i.ne.0) print *, `ieee trapping not supported here'

c アンダーフロー例外を発生させる (トラップはされない)

x = d_min_normal() / 13.0

print *, `d_min_normal() / 13.0 = `, x

c オーバーフロー例外を発生させる

c 出力された値は結果とは関係ない

x = 1.0d300*1.0d300

print *, `1.0d300*1.0d300 = `, x

end

c

c 浮動小数点例外処理関数

c

integer function fp_exc_hdl(sig, sip, uap)

integer sig, code, addr

character label*16

c

c 構造体 /siginfo/ は <sys/siginfo.h> による siginfo_t の解釈

c

structure /fault/

integer address

end structure

structure /siginfo/

integer si_signo

integer si_code

integer si_errno

record /fault/ fault

end structure

record /siginfo/ sip

c FPE コードの一覧は <sys/machsig.h> を参照

c SIGFPE 名を検出

code = sip.si_code

if (code.eq.3) label = 'division'

if (code.eq.4) label = 'overflow'

if (code.eq.5) label = 'underflow'

if (code.eq.6) label = 'inexact'

if (code.eq.7) label = 'invalid'

addr = sip.fault.address

c 発生したシグナルに関する情報を出力

write (*,77) code, label, addr

77 format ('floating-point exception code ', i2, ',',

* a17, ',', ' at address ', z8 )

end

出力は次のとおりです。

d_min_normal() / 13.0 = 1.7115952757748-309 floating-point exception code 4, overflow , at address 1131C 1.0d300*1.0d300 = 1.0000000000000+300 Note: IEEE floating-point exception flags raised: Inexact; Underflow; IEEE floating-point exception traps enabled: overflow; division by zero; invalid operation; See the Numerical Computation Guide, ieee_handler(3M), ieee_flags(3M) [日本語訳] 注: 以下の IEEE 浮動小数点例外が発生しました: 不正確、アンダーフロー 以下の IEEE 浮動小数点例外のトラップが有効です: オーバーフロー、ゼロによる除算、無効な演算 詳細は、『数値計算ガイド』の ieee_flags(3M), ieee_handler(3M)に関す る説明を参照してください。

C プログラム例

コード例 A-13 無効な演算、ゼロ除算、オーバーフロー、アンダーフロー、不正確結果でのトラップ(SPARC) 

/*

* 5 つの IEEE 例外 (無効な演算、ゼロ除算、オーバーフロー、

* アンダーフロー、不正確結果) を発生させます。

* すべての浮動小数点例外をトラップし、 * メッセージを出力してから処理を続けます。 * i = ieee("get","exception","",&out); * 上記の式によって発生した例外について問い合わせをすることも可能です。 * 発生した最高の例外名が含まれ、i はすべての例外を検出するように * 解釈されます。

*/

#include <sunmath.h>

#include <signal.h>

#include <siginfo.h>

#include <ucontext.h>

extern void trap_all_fp_exc(int sig, siginfo_t *sip,

ucontext_t *uap);

int main()

{

double x, y, z;

char *out;

/*

* Use ieee_handler to establish "trap_all_fp_exc"

* as the signal handler to use whenever any floating

* point exception occurs.

*/

if (ieee_handler("set", "all", trap_all_fp_exc) != 0)

printf(" IEEE trapping not supported here.\n");

/* disable trapping (uninteresting) inexact exceptions */

if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)

printf("Trap handler for inexact not cleared.\n");

/* raise invalid */

if (ieee_flags("clear", "exception", "all", &out) != 0)

printf(" could not clear exceptions\n");

printf("1. Invalid: signaling_nan(0) * 2.5\n");

x = signaling_nan(0);

y = 2.5;

z = x * y;

/* raise division */

if (ieee_flags("clear", "exception", "all", &out) != 0)

printf(" could not clear exceptions\n");

printf("2. Div0: 1.0 / 0.0\n");

x = 1.0;

y = 0.0;

z = x / y;

/* raise overflow */

if (ieee_flags("clear", "exception", "all", &out) != 0)

printf(" could not clear exceptions\n");

printf("3. Overflow: -max_normal() - 1.0e294\n");

x = -max_normal();

y = -1.0e294;

z = x + y;

/* raise underflow */

if (ieee_flags("clear", "exception", "all", &out) != 0)

printf(" could not clear exceptions\n");

printf("4. Underflow: min_normal() * min_normal()\n");

x = min_normal();

y = x;

z = x * y;

/* enable trapping on inexact exception */

if (ieee_handler("set", "inexact", trap_all_fp_exc) != 0)

printf("Could not set trap handler for inexact.\n");

/* raise inexact */

if (ieee_flags("clear", "exception", "all", &out) != 0)

printf(" could not clear exceptions\n");

printf("5. Inexact: 2.0 / 3.0\n");

x = 2.0;

y = 3.0;

z = x / y;

/* don't trap on inexact */

if (ieee_handler("set", "inexact", SIGFPE_IGNORE) != 0)

printf(" could not reset inexact trap\n");

/* check that we're not trapping on inexact anymore */

if (ieee_flags("clear", "exception", "all", &out) != 0)

printf(" could not clear exceptions\n");

printf("6. Inexact trapping disabled; 2.0 / 3.0\n");

x = 2.0;

y = 3.0;

z = x / y;

/* find out if there are any outstanding exceptions */

ieee_retrospective_();

/* exit gracefully */

return 0;

}

void trap_all_fp_exc(int sig, siginfo_t *sip, ucontext_t *uap) {

char *label = "undefined";

/* see /usr/include/sys/machsig.h for SIGFPE codes */

switch (sip->si_code) {

case FPE_FLTRES:

label = "inexact";

break;

case FPE_FLTDIV:

label = "division";

break;

case FPE_FLTUND:

label = "underflow";

break;

case FPE_FLTINV:

label = "invalid";

break;

case FPE_FLTOVF:

label = "overflow";

break;

}

printf(

" signal %d, sigfpe code %d: %s exception at address %x\n",

sig, sip->si_code, label, sip->_data._fault._addr);

}

この C プログラムの結果は次のようになります。

1. Invalid: signaling_nan(0) * 2.5 signal 8, sigfpe code 7: invalid exception at address 10da8 2. Div0: 1.0 / 0.0 signal 8, sigfpe code 3: division exception at address 10e44 3. Overflow: -max_normal() - 1.0e294 signal 8, sigfpe code 4: overflow exception at address 10ee8 4. Underflow: min_normal() * min_normal() signal 8, sigfpe code 5: underflow exception at address 10f80 5. Inexact: 2.0 / 3.0 signal 8, sigfpe code 6: inexact exception at address 1106c 6. Inexact trapping disabled; 2.0 / 3.0 Note: IEEE floating-point exception traps enabled: underflow; overflow; division by zero; invalid operation; See the Numerical Computation Guide, ieee_handler(3M) 注: 以下の IEEE 浮動小数点例外が発生しました: 不正確、アンダーフロー 以下の IEEE 浮動小数点例外のトラップが有効です: オーバーフロー、ゼロによる除算、無効な演算 詳細は、『数値計算ガイド』の ieee_flags(3M), ieee_handler(3M)に関す る説明を参照してください。

SPARC では、次のプログラムは、ある例外の状況から起きたデフォルト結果を修正するための ieee_handler とインクルードファイルの使用方法を示しています。

コード例 A-14 例外状況から起きたデフォルト結果を修正する  

/*

* ゼロ除算の例外を発生させ、シグナルハンドラを使用して結果を MAXDOUBLE

* (または MAXFLOAT) で置き換えます。

*

* フラグ -Xa でコンパイルします。

*/

#include <values.h>

#include <siginfo.h>

#include <ucontext.h>

void division_handler(int sig, siginfo_t *sip, ucontext_t *uap);

int main() {

double x, y, z;

float r, s, t;

char *out;

/*

* ieee_handler を使用し、division_handler を IEEE 例外

* "division" が発生したときに使用するシグナルハンドラに設定する

*/

if (ieee_handler("set","division",division_handler)!=0) {

printf(" IEEE trapping not supported here.\n");

}

/* ゼロ除算の例外を発生させる */

x = 1.0;

y = 0.0;

z = x / y;

/*

* ユーザーが提供した値 MAXDOUBLE が IEEE のデフォルト値である

* 無限大の代わりに置換されているかどうかを確認する

*/

printf("double precision division: %g/%g = %g \n",x,y,z);

/* ゼロ除算の例外を発生させる */

r = 1.0;

s = 0.0;

t = r / s;

/*

* ユーザーが提供した値 MAXFLOAT が IEEE のデフォルト値である

* 無限大の代わりに置換されているかどうかを確認する

*/

printf("single precision division: %g/%g = %g \n",r,s,t);

ieee_retrospective_();

return 0;

}

void division_handler(int sig, siginfo_t *sip, ucontext_t *uap) {

int inst;

unsigned rd, mask, single_prec=0;

float f_val = MAXFLOAT;

double d_val = MAXDOUBLE;

long *f_val_p = (long *) &f_val;

/* 例外を発生させる命令 */

inst = uap->uc_mcontext.fpregs.fpu_q->FQu.fpq.fpq_instr;

/*

* 移動先のレジスタを復号化する。

* ビット 29:25 は、SPARC 浮動小数点命令に対して

* 移動先レジスタを符号化する

*/

mask = 0x1f;

rd = (mask & (inst >> 25));

/*

* これは単精度命令か倍精度命令か?

* ビット 5:6 は opcode の精度を符号化する。

* ビット 5 が 1 であれば sp となり、それ以外は dp となる

*/

mask = 0x1;

single_prec = (mask & (inst >> 5));

/* ユーザーが定義した値を移動先レジスタに入れる */

if (single_prec) {

uap->uc_mcontext.fpregs.fpu_fr.fpu_regs[rd] =

f_val_p[0];

} else {

uap->uc_mcontext.fpregs.fpu_fr.fpu_dregs[rd/2] = d_val;

}

}

期待される出力は次のとおりです。

double precision division: 1/0 = 1.79769e+308 single precision division: 1/0 = 3.40282e+38 Note: IEEE floating-point exception traps enabled: division by zero; See the Numerical Computation Guide, ieee_handler(3M) [日本語訳] 注: 以下の IEEE 浮動小数点例外のトラップが有効です: ゼロによる除算 詳細は、『数値計算ガイド』の ieee_handler(3M)に関する説明を参照してくだ さい。

ieee_handler - 例外での異常終了

特定の浮動小数点例外の場合、ieee_handler を使用して強制的にプログラムを異常終了させることができます。

#include <f77_floatingpoint.h> program abort c ieeer = ieee_handler('set', 'division', SIGFPE_ABORT) if (ieeer .ne. 0) print *, ' ieee trapping not supported' r = 14.2 s = 0.0 r = r/s c print *, 'you should not see this; system should abort' c end

libm9x.so の例外処理機能

この節では、libm9x.so が提供する例外処理機能の一部の使用方法を示した例を取りあげます。最初の例は、数値 x と係数 a₀、a₁、...a_N、および b₀、b₁、....b_N-1 の場合に、関数 f(x) とその最初の導関数 f'(x) を評価するタスクに基づいています。この場合 f は連分数 f(x) = a₀ + b₀/(x + a₁ + b₁/(x + ... /(x + a_N-1 + b_N-1/(x + a_N))...))です。IEEE 演算では、f の計算は簡単です。中間除算の 1 つがオーバーフローしたりゼロ除算を行う場合でも、標準によって指定されたデフォルトの値 (正しく符号が付いた無限大) が正しい結果を出します。一方、f' の計算は、これを評価するもっとも簡潔なフォームが、削除可能な特異点を持てるため、f の計算よりも困難です。計算中にこれらの特異点の 1 つが出現すると、不定形 0/0、0*無限大、無限大/無限大の 1 つの評価が試みられます (これらの不定形はすべて無効な演算例外を発生します)。W. Kahan は、「前置換」という機能によりこれらの例外を処理する方法を提唱しています。

前置換は、例外に対する IEEE のデフォルトの応答 (例外演算の結果を置換する値をユーザーが前もって指定することを許可する) を拡張したものです。libm9x.so を使用すると、FEX_CUSTOM 例外処理モードでハンドラをインストールし、プログラムで簡単に前置換を実装できます。このモードでは、ハンドラに渡された info パラメータが示すデータ構造に値を格納するだけで、ハンドラは例外演算の結果に対して任意の値を供給できます。次に、FEX_CUSTOM ハンドラによって実装した前置換を使用して連分数とその導関数を計算するプログラム例を示します。

コード例 A-15 FEX_CUSTOM ハンドラを使用して連分数とその導関数を計算する 

#include <stdio.h>

#include <sunmath.h>

#include <fenv.h>

volatile double p;

void handler(int ex, fex_info_t *info)

{

info->res.type = fex_double;

if (ex == FEX_INV_ZMI)

info->res.val.d = p;

else

info->res.val.d = infinity();

}

/*

* x の位置で係数 a[j] と b[j] によって指定される連分数を

* 評価し、関数値を *pf、導関数の値を *pf1 で返します

*/

void continued_fraction(int N, double *a, double *b,

double x, double *pf, double *pf1)

{

fex_handler_t oldhdl; /* ハンドラの保存または復元のため */

volatile double t;

double f, f1, d, d1, q;

int j;

fex_getexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);

fex_set_handling(FEX_DIVBYZERO, FEX_NONSTOP, NULL);

fex_set_handling(FEX_INV_ZDZ | FEX_INV_IDI | FEX_INV_ZMI,

FEX_CUSTOM, handler);

f1 = 0.0;

f = a[N];

for (j = N - 1; j >= 0; j--) {

d = x + f;

d1 = 1.0 + f1;

q = b[j] / d;

/* 変数 p に対する代入と f1 の評価の間で正しい順序づけを 維持するため、volatile 変数 t に対する次の代入が必要です */

t = f1 = (-d1 / d) * q;

p = b[j-1] * d1 / b[j];

f = a[j] + q;

}

fex_setexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);

*pf = f;

*pf1 = f1;

}

/* 次の係数、x = -3、1、4、および 5 はすべて中間例外に出会います */

double a[] = { -1.0, 2.0, -3.0, 4.0, -5.0 };

double b[] = { 2.0, 4.0, 6.0, 8.0 };

int main()

{

double x, f, f1;

int i;

feraiseexcept(FE_INEXACT); /* 不正確な演算のログ記録を防止します */

fex_set_log(stdout);

fex_set_handling(FEX_COMMON, FEX_ABORT, NULL);

for (i = -5; i <= 5; i++) {

x = i;

continued_fraction(4, a, b, x, &f, &f1);

printf("f(% g) = %12g, f'(% g) = %12g\n", x, f, x, f1);

}

return 0;

}

プログラムについてのコメントが順を追って付けられています。入口で、関数 continued_fraction は、ゼロ除算およびすべての無効な演算例外に対する現在の例外処理モードを保存します。続いて、ゼロ除算に対する連続した例外処理と、3 つの不定形に対する FEX_CUSTOM ハンドラを確立します。このハンドラは 0/0 と 無限大/無限大の両方を無限大で置換しますが、0*無限大は大域変数 p の値で置換します。正しい値を供給して後続の 0*無限大無効演算を置換するように、p は関数を評価するループを通して毎回計算し直す必要があります。また、p はループ内で明示的に記述されていないため、コンパイラが p を削除しないように、p を volatile と宣言する必要があります。最後に、コンパイラが p に対する代入を、例外を発生させる可能性のある演算 (これに対し p が前置換の値を提供する) より上または下に移動させないように、その演算の結果も volatile 変数 (プログラム内では t と呼ばれる) に代入します。fex_setexcepthandler の最後の呼び出しが、ゼロ除算と無効演算に対する本来の処理モードを復元します。

メインプログラムは、fex_set_log 関数を呼び出して、遡及診断のログ記録を有効にします。この呼び出しを行う前に、メインプログラムは不正確フラグを発生させます。これにより、不正確演算の記録を防ぎます (「遡及診断」で説明しているように、FEX_NONSTOP モードでは例外のフラグが発生すると例外がログに記録されません)。メインプログラムはまた、共通例外に対して FEX_ABORT モードを確立し、continued_fraction によって明示的に処理されない異常な例外がプログラムを停止することを防ぎます。プログラムは、最後に数か所で特定の連分数を評価します。次の出力例が示すように、演算は中間例外に出会います。

f(-5) = -1.59649, f'(-5) = -0.1818 f(-4) = -1.87302, f'(-4) = -0.428193 Floating point division by zero at 0x08048dbe continued_fraction, nonstop mode 0x08048dc1 continued_fraction 0x08048eda main Floating point invalid operation (inf/inf) at 0x08048dcf continued_fraction, handler: handler 0x08048dd2 continued_fraction 0x08048eda main Floating point invalid operation (0*inf) at 0x08048dd2 continued_fraction, handler: handler 0x08048dd8 continued_fraction 0x08048eda main f(-3) = -3, f'(-3) = -3.16667 f(-2) = -4.44089e-16, f'(-2) = -3.41667 f(-1) = -1.22222, f'(-1) = -0.444444 f( 0) = -1.33333, f'( 0) = 0.203704 f( 1) = -1, f'( 1) = 0.333333 f( 2) = -0.777778, f'( 2) = 0.12037 f( 3) = -0.714286, f'( 3) = 0.0272109 f( 4) = -0.666667, f'( 4) = 0.203704 f( 5) = -0.777778, f'( 5) = 0.0185185

(x = 1、4、および 5 の場合の計算 f'(x) で発生する例外は、プログラム内で x = -3 のときに起きる例外と同じサイトで発生するため、遡及診断メッセージに出力されません。)

前記のプログラムは、連分数とその導関数の評価で起きる例外を処理する方法として、もっとも効率がよい方法を示しているとは言えません。その 1 つの理由は、必要の有無にかかわらず、ループが繰り返されるごとに前置換の値を計算し直す必要があることです。この場合、前置換の値の計算には浮動小数点除算が含まれますが、最新の SPARC および x86 プロセッサでは浮動小数点除算は比較的遅い演算です。また、ループ自体にすでに 2 つの除算が含まれており、ほとんどの SPARC および x86 プロセッサは 2 つの除算演算の例外をオーバーラップできないため、除算がループ内のボトルネックとなりやすいのです。別の除算を追加するとボトルネックはさらに悪化します。

1 つの除算だけですむようにループを書き直すことができます。実際、前置換値の計算は除算を必要としません (このようにループを書き直すには、b 配列内の係数の隣接要素比を前もって計算する必要があります)。これにより、複数の除算を含む演算のボトルネックは排除されますが、前置換値の計算に含まれるすべての算術演算が除かれるわけではありません。また、前置換値と演算結果の両方が volatile 変数に前置換されるように割り当てる必要があるため、プログラムの速度を低下させるメモリー演算が増えます。これらの代入は、コンパイラがある種のキー操作を再要求することを防止するために欠かせませんが、コンパイラがほかの無関係な演算を再要求することも防止してしまいます。この例のように前置換を介して例外を処理すると、メモリー演算が増え、通常では可能な最適化が行われなくなります。これらの例外を、更に効率よく処理することは可能でしょうか。

高速な前置換に対する特別なハードウェアサポートがない場合、この例における例外をもっとも効率よく処理するには、次の例に示すようにフラグを使用することでしょう。

コード例 A-16 例外処理にフラグを使用する 

#include <stdio.h>

#include <math.h>

#include <fenv.h>

/*

* x の位置で係数 a[j] と b[j] によって指定される連分数を * 評価し、関数値を *pf、導関数の値を *pf1 で返します

*/

void continued_fraction(int N, double *a, double *b,

double x, double *pf, double *pf1)

{

fex_handler_t oldhdl;

fexcept_t oldinvflag;

double f, f1, d, d1, pd1, q;

int j;

fex_getexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);

fegetexceptflag(&oldinvflag, FE_INVALID);

fex_set_handling(FEX_DIVBYZERO | FEX_INV_ZDZ | FEX_INV_IDI |

FEX_INV_ZMI, FEX_NONSTOP, NULL);

feclearexcept(FE_INVALID);

f1 = 0.0;

f = a[N];

for (j = N - 1; j >= 0; j--) {

d = x + f;

d1 = 1.0 + f1;

q = b[j] / d;

f1 = (-d1 / d) * q;

f = a[j] + q;

}

if (fetestexcept(FE_INVALID)) {

/* recompute and test for NaN */

f1 = pd1 = 0.0;

f = a[N];

for (j = N - 1; j >= 0; j--) {

d = x + f;

d1 = 1.0 + f1;

q = b[j] / d;

f1 = (-d1 / d) * q;

if (isnan(f1))

f1 = b[j] * pd1 / b[j+1];

pd1 = d1;

f = a[j] + q;

}

}

fesetexceptflag(&oldinvflag, FE_INVALID);

fex_setexcepthandler(&oldhdl, FEX_DIVBYZERO | FEX_INVALID);

*pf = f;

*pf1 = f1;

}

この例では、最初のループはデフォルトの連続モードで f(x) と f'(x) の計算を試みています。無効フラグが発生すると、2 番目のループは NaN の外形をテストして f(x) と f'(x) を明示的に再計算します。通常、無効演算例外は発生しないため、プログラムは最初のループだけを実行します。このループには volatile 変数に対する参照も余分な算術演算も含まれないため、ループはコンパイラが可能とするかぎりの速度で実行されます。この効率を得るためには、例外が発生するケースを処理するために、2 番目のループを最初のループとほとんど同じように記述しなければなりません。このトレードオフは、浮動小数点例外処理が直面する典型的なジレンマです。

FORTRAN プログラムでの libm9x.so の使用

libm9x.so は主に C および C++ プログラムでの使用を想定したものですが、
Sun の FORTRAN 言語における相互運用性の機能を活用すれば、FORTRAN プログラムからも一部の libm9x.so 関数を呼び出すことができます。

---

注 - 一貫した動作を保つため、libm9x.so の例外処理関数と、ieee_flags および ieee_handler 関数の両方を同じプログラム内で使用することは避けてください。

---

次に、前置換を使用して連分数とその導関数を評価する、FORTRAN によるプログラム例を示します (SPARC のみ)。

コード例 A-17 前置換を使用して連分数とその導関数を評価する (SPARC) 

c

c 前置換ハンドラ

c

subroutine handler(ex, info)

structure /fex_numeric_t/

integer type

union

map

integer i

end map

map

integer*8 l

end map

map

real f

end map

map

real*8 d

end map

map

real*16 q

end map

end union

end structure

structure /fex_info_t/

integer op, flags

record /fex_numeric_t/ op1, op2, res

end structure

integer ex

record /fex_info_t/ info

common /presub/ p

double precision p, d_infinity

volatile p

c 4 = fex_double; 定数については <fenv.h> を参照してください。

info.res.type = 4

c x'80' = FEX_INV_ZMI

if (loc(ex) .eq. x'80') then

info.res.d = p

else

info.res.d = d_infinity()

endif

return

end

c

c x の位置で係数 a[j] と b[j] によって指定される連分数を c 評価し、関数値を f、導関数の値を f1 で返します

c

subroutine continued_fraction(n, a, b, x, f, f1)

integer n

double precision a(*), b(*), x, f, f1

common /presub/ p

integer j, oldhdl

dimension oldhdl(24)

double precision d, d1, q, p, t

volatile p, t

external fex_getexcepthandler, fex_setexcepthandler

external fex_set_handling, handler

c$pragma c(fex_getexcepthandler, fex_setexcepthandler)

c$pragma c(fex_set_handling)

c x'ff2' = FEX_DIVBYZERO | FEX_INVALID

call fex_getexcepthandler(oldhdl, %val(x'ff2'))

c x'2' = FEX_DIVBYZERO, 0 = FEX_NONSTOP

call fex_set_handling(%val(x'2'), %val(0), %val(0))

c x'b0' = FEX_INV_ZDZ | FEX_INV_IDI | FEX_INV_ZMI, 3 = FEX_CUSTOM

call fex_set_handling(%val(x'b0'), %val(3), handler)

f1 = 0.0d0

f = a(n+1)

do j = n, 1, -1

d = x + f

d1 = 1.0d0 + f1

q = b(j) / d

f1 = (-d1 / d) * q

c

c 変数 p に対する代入と f1 の評価の間で正しい順序づけ c を維持するため、volatile 変数 t に対する次の代入が c 必要です

t = f1

p = b(j-1) * d1 / b(j)

f = a(j) + q

end do

call fex_setexcepthandler(oldhdl, %val(x'ff2'))

return

end

c

c メインプログラム

c

program cf

integer i

double precision a, b, x, f, f1

dimension a(5), b(4)

data a /-1.0d0, 2.0d0, -3.0d0, 4.0d0, -5.0d0/

data b /2.0d0, 4.0d0, 6.0d0, 8.0d0/

external fex_set_handling

c$pragma c(fex_set_handling)

c x'ffa' = FEX_COMMON, 1 = FEX_ABORT

call fex_set_handling(%val(x'ffa'), %val(1), %val(0))

do i = -5, 5

x = dble(i)

call continued_fraction(4, a, b, x, f, f1)

write (*, 1) i, f, i, f1

end do

1 format('f(', I2, ') = ', G12.6, ', f''(', I2, ') = ', G12.6)

end

このプログラムの出力は次のようになります。

f(-5) = -1.59649 , f'(-5) = -.181800 f(-4) = -1.87302 , f'(-4) = -.428193 f(-3) = -3.00000 , f'(-3) = -3.16667 f(-2) = -.444089E-15, f'(-2) = -3.41667 f(-1) = -1.22222 , f'(-1) = -.444444 f( 0) = -1.33333 , f'( 0) = 0.203704 f( 1) = -1.00000 , f'( 1) = 0.333333 f( 2) = -.777778 , f'( 2) = 0.120370 f( 3) = -.714286 , f'( 3) = 0.272109E-01 f( 4) = -.666667 , f'( 4) = 0.203704 f( 5) = -.777778 , f'( 5) = 0.185185E-01 Note: IEEE floating-point exception flags raised: Inexact; Division by Zero; Invalid Operation; IEEE floating-point exception traps enabled: overflow; division by zero; invalid operation; See the Numerical Computation Guide, ieee_flags(3M), ieee_handler(3M) [日本語訳] 注: 以下の IEEE 浮動小数点例外が発生しました: 不正確、ゼロによる除算、無効な演算 以下の IEEE 浮動小数点例外のトラップが有効です: オーバーフロー、ゼロによる除算、無効な演算 詳細は、『数値計算ガイド』の ieee_flags(3M), ieee_handler(3M)に関す る説明を参照してください。

その他

sigfpe - 整数例外のトラップ

前節では、ieee_handler を使用した例を示しました。一般的に、ieee_handler と sigfpe のどちらかの使用を選択する場合は、前者をお勧めします。

---

注 - sigfpe は、Solaris オペレーティング環境のみで使用可能です。

---

SPARCでは、たとえば、整数演算例外をトラップする際に使用するハンドラが sigfpe だとします。次の例は整数のゼロ除算でトラップし、結果をユーザーの指定した値に置き換えています。

コード例 A-18 整数例外のトラップ 

/* 整数のゼロ除算例外を生成する */ #include <siginfo.h> #include <ucontext.h> #include <signal.h> void int_handler(int sig, siginfo_t *sip, ucontext_t *uap); int main() { int a, b, c; /* * sigfpe(3) を使用して、整数ゼロ除算で使用するシグナルハンドラとして * int_handler を設定する */ /* * 整数のゼロ除算に対するシグナルハンドラが設定されていないと、 * 整数のゼロ除算は異常終了する */ sigfpe(FPE_INTDIV, int_handler); a = 4; b = 0; c = a / b; printf("%d / %d = %d\n\n", a, b, c); return(0); } void int_handler(int sig, siginfo_t *sip, ucontext_t *uap) { printf("Signal %d, code %d, at addr %x\n", sig, sip->si_code, sip->_data._fault._addr); /* * プログラムカウンタを増やす。オペレーティングシステムは、 * これを浮動小数点例外に対して自動的に行う。 * 整数のゼロ除算に対しては行わない。 */ uap->uc_mcontext.gregs[REG_PC] = uap->uc_mcontext.gregs[REG_nPC]; }

FORTRAN を呼び出す C

FORTRAN のサブルーチンを呼び出す C ドライバの例を示します。C と FORTRAN での動作に関する詳細は、適当な C および FORTRAN のマニュアルを参照してください。以下は C ドライバです (ファイル driver.c に保存します)。

コード例 A-19 FORTRAN を呼び出す C 

/* * 以下の方法を示すデモプログラムです。 * * 1. 配列引数を渡して、f77 サブルーチンを C から呼び出す方法 * 2. 単精度 f77 関数を C から呼び出す方法 * 3. 倍精度 f77 関数を C から呼び出す方法 */ extern int demo_one_(double *); extern float demo_two_(float *); extern double demo_three_(double *); int main() { double array[3][4]; float f, g; double x, y; int i, j; for (i = 0; i < 3; i++) for (j = 0; j < 4; j++) array[i][j] = i + 2*j; g = 1.5; y = g; /* 配列を fortran 関数に渡す (配列の出力) */ demo_one_(&array[0][0]); printf(" from the driver\n"); for (i = 0; i < 3; i++) { for (j = 0; j < 4; j++) printf(" array[%d][%d] = %e\n", i, j, array[i][j]); printf("\n"); } /* 単精度 fortran 関数を呼び出す */ f = demo_two_(&g); printf( " f = sin(g) from a single precision fortran function\n"); printf(" f, g: %8.7e, %8.7e\n", f, g); printf("\n"); /* 倍精度 fortran 関数を呼び出す */ x = demo_three_(&y); printf( " x = sin(y) from a double precision fortran function\n"); printf(" x, y: %18.17e, %18.17e\n", x, y); ieee_retrospective_(); return(0); }

ファイル drivee.f に FORTRAN のサブルーチンを保存します。

subroutine demo_one(array) double precision array(4,3) print *, 'from the fortran routine:' do 10 i =1,4 do 20 j = 1,3 print *, ' array[', i, '][', j, '] = ', array(i,j) 20 continue print * 10 continue return end real function demo_two(number) real number demo_two = sin(number) return end double precision function demo_three(number) double precision number demo_three = sin(number) return end

それから、コンパイルとリンクを行います。

cc -c driver.c f77 -c drivee.f demo_one: demo_two: demo_three: f77 -o driver driver.o drivee.o

出力は次のようになります。

from the fortran routine: array[ 1][ 1] = 0. array[ 1][ 2] = 1.0000000000000 array[ 1][ 3] = 2.0000000000000 array[ 2][ 1] = 2.0000000000000 array[ 2][ 2] = 3.0000000000000 array[ 2][ 3] = 4.0000000000000 array[ 3][ 1] = 4.0000000000000 array[ 3][ 2] = 5.0000000000000 array[ 3][ 3] = 6.0000000000000 array[ 4][ 1] = 6.0000000000000 array[ 4][ 2] = 7.0000000000000 array[ 4][ 3] = 8.0000000000000 from the driver array[0][0] = 0.000000e+00 array[0][1] = 2.000000e+00 array[0][2] = 4.000000e+00 array[0][3] = 6.000000e+00 array[1][0] = 1.000000e+00 array[1][1] = 3.000000e+00 array[1][2] = 5.000000e+00 array[1][3] = 7.000000e+00 array[2][0] = 2.000000e+00 array[2][1] = 4.000000e+00 array[2][2] = 6.000000e+00 array[2][3] = 8.000000e+00 f = sin(g) from a single precision fortran function f, g: 9.9749500e-01, 1.5000000e+00 x = sin(y) from a double precision fortran function x, y: 9.97494986604054446e-01, 1.50000000000000000e+00

効果的なデバッグコマンド

表 A-1 は、SPARC アーキテクチャのデバッグコマンドの例です。

表 A-1   デバッグコマンド(SPARC) 

処理	dbx	adb
ブレークポイントの設定
関数	stop in myfunct	myfunct:b
行番号	stop at 29
絶対アドレス		23a8:b
相対アドレス		main+0x40:b
ブレークポイントに来るまで実行	run	:r
ソースコードの表示	list	<pc,10?ia
浮動小数点レジスタの表示
IEEE 単精度	print $f0	<f0=X
等価の 10 進数	print -fx $f0	<f0=f
IEEE 倍精度	print $f0f1	<f0=X; <f1=X
等価の 10 進数	print -flx $f0f1 print -flx $d0	<f0=F
全浮動小数点レジスタの表示	regs -F	$x for f0-f15 $X for f16-f31
全レジスタの表示	regs	$r; $x; $X
浮動小数点状態レジスタの表示	print -fx $fsr	<fsr=X
f0 に単精度 1.0 を設定	assign $f0=1.0	3f800000>f0
f0/f1 に倍精度 1.0 を設定	assign $f0f1=1.0	3ff00000>f0; 0>f1
実行を継続	cont	:c
シングルステップ	step (or next)	:s
デバッガを終了	quit	$q

表 A-2 は、x86 アーキテクチャのデバッグコマンドの例です。

表 A-2   デバッグコマンド (x86)

処理	dbx	adb
ブレークポイントの設定
関数	stop in myfunct	myfunct:b
行番号	stop at 29
絶対アドレス		23a8:b
相対アドレス		main+0x40:b
ブレークポイントに来るまで実行	run	:r
ソースコードの表示	list	<pc,10?ia
浮動小数点レジスタの表示	print $st0 ... print $st7	$x
全レジスタの表示	examine &$gs/19X	$r
浮動小数点状態レジスタの表示	examine &$fstat/X	<fstat=X または $x
実行を継続	cont	:c
シングルステップ	step (or next)	:s
デバッガを終了	quit	$q

adb のルーチン myfunction に対応するコードの先頭に、ブレークポイントを設定する 2 つの方法を示します。最初の方法では、次のように記述するだけですみます。

myfunction:b

2 番目の方法では、myfunction に対応しているコード部の先頭に相当する絶対アドレスを調べた後、その絶対アドレスにブレークを設定します。

myfunction=X 23a8 23a8:b

f95 でコンパイルされる FORTRAN プログラム内のメインサブルーチンは、adb への MAIN_ として知られています。adb の MAIN_ でブレークポイントを設定するには、次のように指定します。

MAIN_:b

浮動小数点レジスタの内容を調べる場合、dbx コマンドの &$f0/X によって表示される 16 進数値は IEEE 表現であり、その数字の 10 進表現ではありません。SPARC では adb コマンドの $x と $X は、16 進数表現と 10 進数値の両方を表示します。x86 では、adb コマンドの $x は 10 進数のみを表示します。SPARC では倍精度の値については、奇数番号のレジスタの横に 10 進数値が表示されます。

プロセスが浮動小数点ユニットを使用するまで、それはオペレーティングシステムによって使用不可になっているため、デバッグしているプログラムがアクセスするまで、浮動小数点ユニットを変更することはできません。

SPARC では、浮動小数点数を表示する場合、レジスタのサイズは 32 ビットであり、1 つの単精度浮動小数点数は 32 ビットを占め (したがって 1 つのレジスタに収まる)、倍精度浮動小数点数は 64 ビットを占める (したがって倍精度数 1 つを保持するには 2 つのレジスタが使用される) ことを覚えておいてください。16 進数表現では、32 ビットは 8 桁の数字に相当します。adb で表示した浮動小数点ユニットでは、表示は次のような形で行われます。

< 浮動小数点レジスタ名> <IEEE 16 進数の値> <単精度> <倍精度>

SPARC では、3 番目の欄には、2 番目の欄に示された IEEE 16 進数を単精度の 10 進数に変換した値が保持されています。4 番目の欄では、レジスタの対を解釈しています。

SPARC では、f10 と f11 を 64 ビットの IEEE 倍精度数として解釈しています。1 つの倍精度値を保持するのに f10 と f11 を使うので、その値の最初の 32 ビットの (f10 行での) 7ff00000 は +NaN とは解釈されません。64 ビット 7ff00000 00000000 全体の解釈である +無限大は、意味のある解釈になります。

SPARC では、最初の 16 の浮動小数点データレジスタを表示するために使用されている adb コマンド $x は、fsr (浮動小数点ステータスレジスタ) も表示しています。

$x fsr 40020 f0 400921fb +2.1426990e+00 f1 54442d18 +3.3702806e+12 +3.1415926535897931e+00 f2 2 +2.8025969e-45 f3 0 +0.0000000e+00 +4.2439915819305446e-314 f4 40000000 +2.0000000e+00 f5 0 +0.0000000e+00 +2.0000000000000000e+00 f6 3de0b460 +1.0971904e-01 f7 0 +0.0000000e+00 +1.2154188766544394e-10 f8 3de0b460 +1.0971904e-01 f9 0 +0.0000000e+00 +1.2154188766544394e-10 f10 7ff00000 +NaN f11 0 +0.0000000e+00 +Infinity f12 ffffffff -NaN f13 ffffffff -NaN -NaN f14 ffffffff -NaN f15 ffffffff -NaN -NaN

対応する x86 での出力は、次のとおりです。

$x 80387 chip is present. cw 0x137f sw 0x3920 cssel 0x17 ipoff 0x2d93 datasel 0x1f dataoff 0x5740 st[0] +3.24999988079071044921875 e-1 VALID st[1] +5.6539133243479549034419688 e73 EMPTY st[2] +2.0000000000000008881784197 EMPTY st[3] +1.8073218308070440556016047 e-1 EMPTY st[4] +7.9180300235748291015625 e-1 EMPTY st[5] +4.201639036693904927233234 e-13 EMPTY st[6] +4.201639036693904927233234 e-13 EMPTY st[7] +2.7224999213218694649185636 EMPTY

---

注 - (x86 のみ) cw は制御ワード、sw はステータスワードです。

---