このセクションでは、数学ライブラリの関数を使用した例を示します。
次の例では、数値の配列を生成する乱数ジェネレータを呼び出し、 指定した数の EXP を計算するためにかかる時間を測定する関数を使用します。
使用例 8 乱数ジェネレータ#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 関数を使用して、ユーザーが指定した範囲に均一に分散するランダムデータのブロックを生成する方法を示します。
使用例 9 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 標準規格が推奨するいくつかの関数を使用しています。
使用例 10 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
このプログラムからの出力は次の例のようになります。
使用例 11 例 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 関数は、アンダーフローおよびオーバーフローを処理するため、ベクトルのユークリッドノルムを計算し、この環境関数を使用します。遡及診断出力が示すように、メインプログラムは、アンダーフローとオーバーフローを発生させるように位取りされたベクトルを使用してこの関数を呼び出します。
使用例 12 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 関数を使用すると、すべての結果を必要な精度に丸めることができます。
使用例 13 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を参照してください。
使用例 14 マルチスレッドプログラムでの環境関数の使用#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);
...
}