本节显示使用数学库中函数的示例。
以下示例调用随机数字生成器以生成数字数组,并使用计时函数测量计算指定数字的 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
如果使用 f95 编译器和 -f77 兼容性选项,将显示以下附加消息。
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 上的输出略有不同,其中将保留双精度数值的十六进制表示形式的高位字和低位字:
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
如果使用 f95 编译器和 -f77 兼容性选项编译程序,则输出中将包括以下附加消息。
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
以下代码示例显示 fesetprec 函数在基于 x86 的系统上的效果。此函数在基于 SPARC 的系统上不可用。while 循环通过查找在加 1 时完全舍入的 2 的最大次幂,尝试确定可用精度。显示第一个循环时,在使用扩展精度计算所有中间结果的类似于基于 x86 的系统的体系结构上,此技术并不总是按照预期工作。因此,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; }
在 x86 系统上使用 cc A8.c –lm –xarch=386 编译时将生成
64 significant bits 53 significant bit
最后,以下示例显示一种方法,在多线程程序中使用环境函数将浮点模式从父线程传播到子线程,并恢复在子线程与父线程连接时在其中引发的异常标志。有关编写多线程程序的更多信息,请参见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); ... }