This section shows examples that use functions from the math library.
The following example calls a random number generator to generate an array of numbers and uses a timing function to measure the time it takes to compute the EXP of the given numbers:
Example 8 Random Number Generator#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
To compile the preceding example, place the source code in a file with the suffix F (not f) so that the compiler will automatically invoke the preprocessor, and specify either –DSP or –DDP on the command line to select single or double precision.
This example shows how to use the d_addrans function to generate blocks of random data uniformly distributed over a user-specified range:
Example 9 Using the d_addrans Function/*
* 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;
}
The following Fortran example uses some functions recommended by the IEEE standard:
Example 10 IEEE Recommended Functionsc
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
The output from this program is shown in the following example.
Example 11 Output of Example A-5DOUBLE 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
If using the f95 compiler with the -f77 compatibility option, the following additional messages are displayed.
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)
The following C program calls several of the ieee_values(3m) functions:
#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;
}
Remember to specify both -lsunmath and -lm when linking.
On SPARC-based systems, the output looks like the following:
quiet NaN: NaN = 7ff80000 00000000 nextafter(max_subnormal,0) = 2.2250738585072004e-308 = 000fffff fffffffe single precision min subnormal = 1.40129846e-45 = 00000001
Because the x86 architecture is “little-endian”, the output on x86 is slightly different, such that the high and low order words of the hexadecimal representations of the double precision numbers are reversed:
quiet NaN: NaN = ffffffff 7fffffff
nextafter(max_subnormal,0) = 2.2250738585072004e-308
= fffffffe 000fffff
single precision min subnormal = 1.40129846e-45 = 00000001
Fortran programs that use ieee_values functions should take care to declare those functions' types:
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
On SPARC, the output reads as follows:
6.475175119438025110924438958227646E-4966 Inf 3.4028235E+38
The following example demonstrates how to set the rounding mode to round towards zero:
#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;
}
The following output of the previous rounding direction short program shows the effects of rounding towards zero on 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%
The following output of the previous rounding direction short program shows the effects of rounding towards zero on 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%
To set rounding direction towards zero from a Fortran program, use the following example:
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
The output is as follows:
demo% f95 ieee_flags_demo.f demo% a.out Rounding direction is: tozero
If the program is compiled using the f95 compiler with the -f77 compatibility option, the output includes the following additional messages.
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)
The next example illustrates the use of several of the C99 floating-point environment functions. The norm function computes the Euclidean norm of a vector and uses the environment functions to handle underflow and overflow. The main program calls this function with vectors that are scaled to ensure that underflows and overflows occur, as the retrospective diagnostic output shows
Example 12 C99 Floating-Point Environment Functions#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;
}
On SPARC-based systems, compiling and running this program produces output like the following:
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
The following code example shows the effect of the fesetprec function on x86-based systems. This function is not available on SPARC-based systems. The while loops attempt to determine the available precision by finding the largest power of two that rounds off entirely when it is added to one. As the first loop shows, this technique does not always work as expected on architectures like x86-based systems that evaluate all intermediate results in extended precision. Thus, the fesetprec function can be used to guarantee that all results will be rounded to the desired precision, as the second loop shows.
Example 13 fesetprec Function (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;
}
Compiling on x86 systems with cc A8.c –lm –xarch=386 creates
64 significant bits 53 significant bit
Finally, the following example shows one way to use the environment functions in a multithreaded program to propagate floating-point modes from a parent thread to a child thread and recover exception flags raised in the child thread when it joins with the parent. See the Multithreaded Programming Guide for more information on writing multi-threaded programs.
Example 14 Using Environment Functions in a Multithread Program#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);
...
}