この付録では、一般的なタスクを行う方法の例を示します。例は Fortran または ANSI C で記述されており、その多くは現在のバージョンの libm および libsunmath に依存しています。これらの例は、Oracle Solaris 10 Update 10 OS 以降のリリースの Oracle Solaris Studio 12.4 でテストされています。C の例は、–lsunmath –lm オプションを使用してコンパイルされています。
次の例は、浮動小数点数の 16 進表現を検査できる 1 つの方法を示しています。格納されているデータの 16 進表現を参照するには、デバッガを使用することもできます。
次の C のプログラムは、倍精度の近似値を π と単精度の無限大に出力しています。
使用例 A-1 倍精度の例#include <math.h>
#include <sunmath.h>
int main() {
union {
float flt;
unsigned un;
} r;
union {
double dbl;
unsigned un[2];
} d;
/* double precision */
d.dbl = M_PI;
(void) printf("DP Approx pi = %08x %08x = %18.17e \n",
d.un[0], d.un[1], d.dbl);
/* single precision */
r.flt = infinityf();
(void) printf("Single Precision %8.7e : %08x \n",
r.flt, r.un);
return 0;
}
–lsunmath を指定してコンパイルされた SPARC ベースのシステムでは、前述の プログラム の出力は次のようになります。
DP Approx pi = 400921fb 54442d18 = 3.14159265358979312e+00 Single Precision Infinity: 7f800000
次の Fortran プログラムは、最小の正規数をそれぞれの形式で出力します。
使用例 A-2 各形式での最小の正規数の出力 (続き)
program print_ieee_values
c
c the purpose of the implicit statements is to ensure
c that the floatingpoint pseudo-intrinsic functions
c are declared with the correct type
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
このセクションでは、数学ライブラリの関数を使用した例を示します。
次の例では、数値の配列を生成する乱数ジェネレータを呼び出し、 指定した数の 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 tarray(2), u1, u2 c c compute EXP on random numbers in [-ln2/2,ln2/2] lb = -0.3465735903 ub = 0.3465735903 c c generate array of random numbers #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 start the clock call dtime(tarray) u1 = tarray(1) c c compute exponentials do 16 i=1,SIZE y = exp(x(i)) 16 continue c c get the elapsed time call dtime(tarray) 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 関数の使用/*
* test SIZE*LOOPS random arguments to sin in the range
* [0, threshold] where
* threshold = 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 {
unsigned u[2];
double d;
} upperbound;
upperbound.u[0] = 0x3e300000;
upperbound.u[1] = 0x00000000;
/* initialize the random number generator */
d_init_addrans_();
/* test (SIZE * LOOPS) arguments to sin */
for (j = 0; j < LOOPS; j++) {
/*
* generate a vector, x, of length SIZE,
* of random numbers to use as
* input to the trig functions.
*/
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]);
/* is sin(x) == x? It ought to, for tiny 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;
}
次の Fortran の例では、IEEE 標準規格が推奨するいくつかの関数を使用しています。
使用例 A-5 IEEE が推奨する関数c
c Demonstrate how to call 5 of the more interesting IEEE
c recommended functions from Fortran. These are implemented
c with "bit-twiddling", and so are as efficient as you could
c hope. The IEEE standard for floating-point arithmetic
c doesn't require these, but recommends that they be
c included in any IEEE programming environment.
c
c For example, to accomplish
c y = x * 2**n,
c since the hardware stores numbers in base 2,
c shift the exponent by n places.
c
c Refer to
c ieee_functions(3m)
c libm_double(3f)
c libm_single(3f)
c
c The 5 functions demonstrated here are:
c
c ilogb(x): returns the base 2 unbiased exponent of x in
c integer format
c signbit(x): returns the sign bit, 0 or 1
c copysign(x,y): returns x with y's sign bit
c nextafter(x,y): next representable number after x, in
c the direction y
c scalbn(x,n): x * 2**n
c
c function double precision single precision
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)
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.5 is 1 Starting from 4.9406564584124654E-324, the next representable number towards -Inf is 0.0000000000000000E+000 Starting from 4.9406564584124654E-324, the next representable number towards 1.0 is 9.8813129168249309E-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.5 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
-f77 互換性オプションを指定して f95 コンパイラを使用した場合は、次の追加のメッセージが表示されます。
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_flags(3M), ieee_handler(3M)
次の 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 = 7ff80000 00000000 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 the purpose of the implicit statements is to insure c that the floating-point pseudo-instrinsic c functions are declared with the correct type 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.475175119438025110924438958227646E-4966 Inf 3.4028235E+38
次の例は、丸めモードをゼロへの丸めに設定する方法を示しています。
#include <math.h>
#include <sunmath.h>
int main()
{
int i;
double x, y;
char *out_1, *out_2, *dummy;
/* get prevailing rounding direction */
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);
/* set rounding direction */
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);
/*
* restore original rounding direction before printf, since
* printf is also affected by the current rounding direction
*/
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 -lsunmath -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 -lsunmath -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
出力は次のようになります。
demo% f95 ieee_flags_demo.f demo% a.out Rounding direction is: tozero
-f77 互換性オプションを指定して f95 コンパイラを使用してプログラムをコンパイルした場合は、次の追加のメッセージが出力に表示されます。
demo% f95 ieee_flags_demo.f -f77
demo% a.out
Note: Rounding direction toward zero
IEEE floating-point exception traps enabled:
overflow; division by zero; invalid operation;
See the Numerical Computation Guide, ieee_flags(3M), ieee_handler(3M)
次の例では、いくつかの C99 浮動小数点環境関数を使用しています。norm 関数は、アンダーフローおよびオーバーフローを処理するため、ベクトルのユークリッドノルムを計算し、この環境関数を使用します。遡及診断出力が示すように、メインプログラムは、アンダーフローとオーバーフローを発生させるように位取りされたベクトルを使用してこの関数を呼び出します。
使用例 A-7 C99 浮動小数点環境関数#include <stdio.h>
#include <math.h>
#include <sunmath.h>
#include <fenv.h>
/*
* Compute the euclidean norm of the vector x avoiding
* premature underflow or overflow
*/
double norm(int n, double *x)
{
fenv_t env;
double s, b, d, t;
int i, f;
/* save the environment, clear flags, and establish nonstop
exception handling */
feholdexcept(&env);
/* attempt to compute the dot product x.x */
d = 1.0; /* scale factor */
s = 0.0;
for (i = 0; i < n; i++)
s += x[i] * x[i];
/* check for underflow or overflow */
f = fetestexcept(FE_UNDERFLOW | FE_OVERFLOW);
if (f & FE_OVERFLOW) {
/* first attempt overflowed, try again scaling down */
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)) {
/* first attempt underflowed, try again scaling up */
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;
}
}
/* hide any underflows that have occurred so far */
feclearexcept(FE_UNDERFLOW);
/* restore the environment, raising any other exceptions
that have occurred */
feupdateenv(&env);
/* take the square root and undo any scaling */
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 -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;
}
cc A8.c –lm –xarch=386 を指定して x86 システムでコンパイルすると、次のように出力されます。
64 significant bits 53 significant bit
最後に、次の例は、マルチスレッドプログラムで環境関数を使用して、親スレッドから子スレッドに浮動小数点モードを伝播し、子スレッドが親スレッドに結合されるときに子スレッド内で発生した例外フラグを回復する 1 つの方法を示しています。マルチスレッドプログラムの記述については、Multithreaded Programming Guide を参照してください。
使用例 A-9 マルチスレッドプログラムでの環境関数の使用#include <thread.h>
#include <fenv.h>
fenv_t env;
void * child(void *p)
{
/* inherit the parent's environment on entry */
fesetenv(&env);
...
/* save the child's environment before exit */
fegetenv(&env);
}
void parent()
{
thread_t tid;
void *arg;
...
/* save the parent's environment before creating the child */
fegetenv(&env);
thr_create(NULL, NULL, child, arg, NULL, &tid);
...
/* join with the child */
thr_join(tid, NULL, &arg);
/* merge exception flags raised in the child into the
parent's environment */
fex_merge_flags(&env);
...
}
通常、累積例外ビットの検査またはクリアはユーザープログラムで行います。次の例は、累積例外フラグを 検査する C プログラムです。
使用例 A-10 累積例外フラグの検査#include <sunmath.h>
#include <sys/ieeefp.h>
int main()
{
int code, inexact, division, underflow, overflow, invalid;
double x;
char *out;
/* cause an underflow exception */
x = max_subnormal() / 2.0;
/* this statement insures that the previous */
/* statement is not optimized away */
printf("x = %g\n",x);
/* find out which exceptions are raised */
code = ieee_flags("get", "exception", "", &out);
/* decode the return value */
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" is the raised exception with the highest priority */
printf(" Highest priority exception is: %s\n", out);
/* The value 1 means the exception is raised, */
/* 0 means it isn't. */
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)
使用例 A-11 累積例外フラグの検査 – Fortran/*
A Fortran example that:
* causes an underflow exception
* uses ieee_flags to determine which exceptions are raised
* decodes the integer value returned by ieee_flags
* clears all outstanding exceptions
Remember to save this program in a file with the suffix .F, so that
the c preprocessor is invoked to bring in the header file
floatingpoint.h.
*/
#include <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 Cause an underflow exception
x = d_max_subnormal() / 2.0
c Find out which exceptions are raised
accrued = ieee_flags('get', 'exception', '', out)
c Decode value returned by ieee_flags using bit-shift intrinsics
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 The exception with the highest priority is returned in "out"
print *, "Highest priority exception is ", out
c The value 1 means the exception is raised; 0 means it is not
print *, inv, over, div, under, inx
c Clear all outstanding exceptions
i = ieee_flags('clear', 'exception', 'all', out)
end
出力は次のようになります。
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 set 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
次の例は、例外を特定するためのシグナルハンドラをインストールする Fortran プログラムです (SPARC ベースシステムの場合のみ)。
使用例 A-12 アンダーフローでのトラップ (SPARC) program demo
c declare signal handler function
external fp_exc_hdl
double precision d_min_normal
double precision x
c set up signal handler
i = ieee_handler('set', 'common', fp_exc_hdl)
if (i.ne.0) print *, 'ieee trapping not supported here'
c cause an underflow exception (it will not be trapped)
x = d_min_normal() / 13.0
print *, 'd_min_normal() / 13.0 = ', x
c cause an overflow exception
c the value printed out is unrelated to the result
x = 1.0d300
x = x * x
print *, '1.0d300*1.0d300 = ', x
end
c
c the floating-point exception handling function
c
integer function fp_exc_hdl(sig, sip, uap)
integer sig, code, addr
character label*16
c
c The structure /siginfo/ is a translation of siginfo_t
c from <sys/siginfo.h>
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 See <sys/machsig.h> for list of FPE codes
c Figure out the name of the 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 Print information about the signal that happened
write (*,77) code, label, addr
77 format ('floating-point exception code ', i2, ',',
* a17, ',', ' at address ', z8 )
end
前のコードを –f77 を指定してコンパイルすると、出力は次のようになります。
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_flags(3M),
ieee_handler(3M)
次の例は、SPARC ベースシステムでの より複雑な C の例です。
使用例 A-13 無効、ゼロ除算、オーバーフロー、アンダーフロー、および不正確でのトラップ (SPARC)/*
* Generate the 5 IEEE exceptions: invalid, division,
* overflow, underflow and inexact.
*
* Trap on any floating point exception, print a message,
* and continue.*
* Note that you could also inquire about raised exceptions by
* i = ieee("get","exception","",&out);* where out contains the name of the highest exception
* raised, and i can be decoded to find out about all the
* exceptions raised.
*/
#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);
}
出力は、次のようになります。
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)
次のコードは、SPARC で ieee_handler およびインクルードファイルを使用して、特定の例外状況のデフォルトの結果を変更する方法を示しています。
使用例 A-14 例外状況のデフォルトの結果の変更/*
* Cause a division by zero exception and use the
* signal handler to substitute MAXDOUBLE (or MAXFLOAT)
* as the result.
*
* compile with the flag -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;
/*
* Use ieee_handler to establish division_handler as the
* signal handler to use for the IEEE exception division.
*/
if (ieee_handler("set","division",division_handler)!=0) {
printf(" IEEE trapping not supported here.\n");
}
/* Cause a division-by-zero exception */
x = 1.0;
y = 0.0;
z = x / y;
/*
* Check to see that the user-supplied value, MAXDOUBLE,
* is indeed substituted in place of the IEEE default
* value, infinity.
*/
printf("double precision division: %g/%g = %g \n",x,y,z);
/* Cause a division-by-zero exception */
r = 1.0;
s = 0.0;
t = r / s;
/*
* Check to see that the user-supplied value, MAXFLOAT,
* is indeed substituted in place of the IEEE default
* value, infinity.
*/
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;
/* Get instruction that caused exception. */
inst = uap->uc_mcontext.fpregs.fpu_q->FQu.fpq.fpq_instr;
/*
* Decode the destination register. Bits 29:25 encode the
* destination register for any SPARC floating point
* instruction.
*/
mask = 0x1f;
rd = (mask & (inst >> 25));
/*
* Is this a single precision or double precision
* instruction? Bits 5:6 encode the precision of the
* opcode; if bit 5 is 1, it's sp, else, dp.
*/
mask = 0x1;
single_prec = (mask & (inst >> 5));
/* put user-defined value into destination register */
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_handler を使用して強制的にプログラムを中止させることができます。
#include <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
次の例は、libm によって提供される例外処理機能のいくつかを使用する方法を示しています。最初の例は、数値 x と係数 a0、a1、...、aN、および b0、b1、...、bN-1 の場合に、関数 f(x) とその一次導関数 f'(x) を評価するタスクに基づいています。ここで f() は次の連分数です。
f(x) =a0 + b0 / (x + a1 + b1 / (x + ... / (x + aN-1 + bN-1 / (x + aN))...))
IEEE 演算では、f() の計算は単純です。中間除算の 1 つがオーバーフローしたり、ゼロ除算を行う場合でも、標準規格によって指定されているデフォルトの値 (正しい符号が付けられた無限大) から正しい結果が算出されます。一方、f'() の計算は、これを評価するもっとも簡潔な形式が、削除可能な特異点を持つことができるため、f の計算よりも難しくなります。計算中にこれらの特異点の 1 つが出現すると、不定形 0/0、0*無限大、または無限大/無限大のいずれか 1 つの評価が試みられます (これらはすべて無効な演算例外を発生します)。W. Kahan は、前置換という機能によってこれらの例外を処理する方法を提唱しています。
前置換は、例外に対する IEEE のデフォルトの応答を拡張したものであり、例外演算の結果を置換する値をユーザーが前もって指定できます。libm の例外処理機能を使用すると、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();
}
/*
* Evaluate the continued fraction given by coefficients a[j] and
* b[j] at the point x; return the function value in *pf and the
* derivative in *pf1
*/
void continued_fraction(int N, double *a, double *b,
double x, double *pf, double *pf1)
{
fex_handler_t oldhdl; /* for saving/restoring handlers */
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;
/* the following assignment to the volatile variable t
is needed to maintain the correct sequencing between
assignments to p and evaluation of f1 */
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;
}
/* For the following coefficients, x = -3, 1, 4, and 5 will all
encounter intermediate exceptions */
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); /* prevent logging of 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 はループ内で明示的に記述されていないため、コンパイラによって削除されないように 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>
/*
* Evaluate the continued fraction given by coefficients a[j] and
* b[j] at the point x; return the function value in *pf and the
* derivative in *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 番目のループを最初のループとほとんど同じように記述する必要があります。このトレードオフは、浮動小数点例外処理において生じることがある典型的なジレンマです。
libm の例外処理機能は主に C/C++ プログラムからの使用を想定していますが、Sun の Fortran 言語の相互運用性機能を使用すると、Fortran プログラムからも一部の libm 関数を呼び出すことができます。
次の例では、前置換を使用して連分数とその導関数を評価するための Fortran バージョンのプログラムを示します (SPARC のみ)。
使用例 A-17 前置換を使用した連分数とその導関数の評価 (SPARC)c
c Presubstitution handler
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; see <fenv.h> for this and other constants
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 Evaluate the continued fraction given by coefficients a(j) and
c b(j) at the point x; return the function value in f and the
c derivative in 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
data ixff2/x'ff2'/
data ix2/x'2'/
data ixb0/x'b0'/
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(ixff2))
c x'2' = FEX_DIVBYZERO, 0 = FEX_NONSTOP
call fex_set_handling(%val(ix2), %val(0), %val(0))
c x'b0' = FEX_INV_ZDZ | FEX_INV_IDI | FEX_INV_ZMI, 3 = FEX_CUSTOM
call fex_set_handling(%val(ixb0), %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 the following assignment to the volatile variable t
c is needed to maintain the correct sequencing between
c assignments to p and evaluation of f1
t = f1
p = b(j-1) * d1 / b(j)
f = a(j) + q
end do
call fex_setexcepthandler(oldhdl, %val(ixff2))
return
end
c Main program
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/
data ixffa/x'ffa'/
external fex_set_handling
c$pragma c(fex_set_handling)
c x'ffa' = FEX_COMMON, 1 = FEX_ABORT
call fex_set_handling(%val(ixffa), %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
–f77 フラグを指定してコンパイルされたこのプログラムの出力は次のようになります。
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; Underflow; 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_handler を使用した例を示しました。通常、ieee_handler または sigfpe のどちらを使用するかを選択する場合は、前者をお勧めします。
整数演算例外をトラップする場合などに、ハンドラとして sigfpe を使用することがあります。Example A–18 では、SPARC ベースのシステムで整数のゼロ除算をトラップしています。
使用例 A-18 整数例外のトラップ/* Generate the integer division by zero exception */
#include <signal.h>
#include <siginfo.h>
#include <ucontext.h>
void int_handler(int sig, siginfo_t *sip, ucontext_t *uap);
int main() {
int a, b, c;
/*
* Use sigfpe(3) to establish "int_handler" as the signal handler
* to use on integer division by zero
*/
/*
* Integer division-by-zero aborts unless a signal
* handler for integer division by zero is set up
*/
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);
/*
* automatically for floating-point exceptions but not for
* integer division by zero.
*/
uap->uc_mcontext.gregs[REG_PC] =
uap->uc_mcontext.gregs[REG_nPC];
}
次に、 Fortran のサブルーチンを呼び出す C ドライバの簡単な例を示します。C および Fortran での動作の詳細は、Oracle Solaris Studio 12.4: C ユーザーガイド およびOracle Solaris Studio 12.4: Fortran ユーザーズガイド を参照してください。C ドライバを次に示します (driver.c というファイルに保存します)。
使用例 A-19 C からの Fortran の呼び出し/*
* a demo program that shows:
* 1. how to call f95 subroutine from C, passing an array argument
* 2. how to call single precision f95 function from C
* 3. how to call double precision f95 function from 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;
/* pass an array to a fortran function (print the array) */
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");
}
/* call a single precision fortran function */
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");
/* call a double precision fortran function */
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 f95 -c drivee.f demo_one: demo_two: demo_three: f95 -o driver driver.o drivee.o
出力は次のようになります。
from the fortran routine:
array[ 1 ][ 1 ] = 0.0E+0
array[ 1 ][ 2 ] = 1.0
array[ 1 ][ 3 ] = 2.0
array[ 2 ][ 1 ] = 2.0
array[ 2 ][ 2 ] = 3.0
array[ 2 ][ 3 ] = 4.0
array[ 3 ][ 1 ] = 4.0
array[ 3 ][ 2 ] = 5.0
array[ 3 ][ 3 ] = 6.0
array[ 4 ][ 1 ] = 6.0
array[ 4 ][ 2 ] = 7.0
array[ 4 ][ 3 ] = 8.0
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
次の表は、SPARC アーキテクチャーのデバッグコマンドの例を示しています。
|
浮動小数点数を表示する場合、レジスタのサイズは 32 ビットであり、1 つの単精度浮動小数点数は 32 ビットを占有すること (つまり、1 つのレジスタに収まります)、倍精度浮動小数点数は 64 ビットを占有する (つまり、1 つの倍精度数を保持するには 2 つのレジスタが使用されます) に留意してください。16 進数表現では、32 ビットは 8 桁の 16 進数に相当します。adb を使用して FPU レジスタを表示したスナップショットでは、表示は次のような形式になります。
<fpu レジスタ名> <IEEE 16 進数の値> <単精度> <倍精度>
3 番目の列には、2 番目の列に表示されている 16 進数を単精度の 10 進数に変換した値が示されます。4 番目の列は、レジスタのペアを解釈しています。たとえば、f11 の行の 4 番目の列は、f10 および f11 を 64 ビットの IEEE 倍精度数として解釈しています。
f10 および f11 は 1 つの倍精度値の保持に使用されているため、その値の最初の 32 ビットの (f10 行の) 7ff00000 が +NaN として解釈されても無意味になります。64 ビット全体 7ff00000 00000000 の解釈である +Infinity は、意味のある解釈になっています。
最初の 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 アーキテクチャーのデバッグコマンドの例を示しています。
|
次の例は、adb のルーチン myfunction に対応するコードの先頭にブレークポイントを設定する 2 つの方法を示しています。最初に次のコマンドを使用できます。
myfunction:b
次に、myfunction に対応するコードの先頭に対応する絶対アドレスを判別して、その絶対アドレスにブレークを設定できます。
myfunction=X 23a8 23a8:b
f95 を指定してコンパイルされた Fortran プログラムのメインサブルーチンは、 adb への MAIN_ と認識されています。adb の MAIN_ にブレークポイントを設定するには、次のコマンドを実行します。
MAIN_:b
浮動小数点レジスタの内容を検査する場合、dbx コマンドの regs –F で表示される 16 進数値は、数値の 10 進数表現ではなく base-16 表現です。SPARC ベースシステムでは、adb コマンドの $x および $X は、16 進数表現および 10 進数値の両方を表示します。x86 ベースシステムでは、adb コマンドの $x は 10 進数値のみを表示します。SPARC ベースシステムでは、倍精度値の場合、奇数レジスタの横に 10 進数値が表示されます。
オペレーティングシステムではプロセスが最初に使用するまで浮動小数点ユニットが無効にされるため、デバッグしているプログラムがアクセスするまで、浮動小数点レジスタを変更できません。
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