Numerical Computation Guide

Examples

This appendix provides examples of how to accomplish some popular tasks. The examples are written either in Fortran or ANSI C, and many depend on the current versions of `libm` and `libsunmath`. These examples were tested with the current C and Fortran compilers in the Solaris operating environment.

IEEE Arithmetic

The following examples show one way you can examine the hexadecimal representations of floating-point numbers. Note that you can also use the debuggers to look at the hexadecimal representations of stored data.

The following C program prints a double precision approximation to and single precision infinity:

CODE EXAMPLE A-1

 ```#include ``` ```#include ``` ``` ``` ```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; ``` ```} ```

Double Precision Example

On SPARC, the output of the preceding program looks like this:

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

The following Fortran program prints the smallest normal numbers in each format:

CODE EXAMPLE A-2 Print Smallest Normal Numbers in Each Format

 ``` program print_ieee_values ``` ```c ``` ```c the purpose of the implicit statements is to ensure ``` ```c that the f77_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 ```

On SPARC, the corresponding output reads:

 ```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 ```

(The Fortran compiler on x86 does not support the `real*16` type. To run the preceding example on x86, delete the `real*16` declaration and the code that calculates and prints the quadruple precision value.)

The Math Libraries

This section shows examples that use functions from the math library.

Random Number Generator

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:

CODE EXAMPLE A-3 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 `d_addrans` to generate blocks of random data uniformly distributed over a user-specified range:

CODE EXAMPLE A-4 Using `d_addrans`
 ```/* ``` ``` * test SIZE*LOOPS random arguments to sin in the range ``` ``` * [0, threshold] where ``` ``` * threshold = 3E30000000000000 (3.72529029846191406e-09) ``` ``` */ ``` ```#include ``` ```#include ``` ```#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; ``` ```} ```

IEEE Recommended Functions

This Fortran example uses some functions recommended by the IEEE standard:

CODE EXAMPLE A-5 IEEE Recommended Functions
 ```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 ``` ```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) ``` ```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 ```

The output from this program is shown in CODE EXAMPLE A-6.

CODE EXAMPLE A-6 Output of CODE EXAMPLE A-5
```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.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.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 Special Values

This C program calls several of the `ieee_values`(3m) functions:

 ```#include ``` ```#include ``` ```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, the output looks like this:

 ```quiet NaN: NaN = `7fffffff` ffffffff ``` ```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 (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 `000ffff`f ``` ```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 f77_floatingpoint 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.4751751194380251109244389582276466-4966 ``` ``` Infinity ``` ``` 3.40282E+38 ```

(Recall that the Fortran compiler on x86 does not support the `real*16` type. To run the preceding example on x86, delete all references to `real*16` variables and functions.)

`ieee_flags` -- Rounding Direction

The following example demonstrates how to set the rounding mode to round towards zero:

 ```#include ``` ```#include ``` ```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) The output of this short program shows the effects of rounding towards zero:

 ```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) The output of this short program shows the effects of rounding towards zero:

 ```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:

 ```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:

 ``` Rounding direction is: tozero ``` ``` Note: Rounding direction toward zero. ``` ``` See the Numerical Computation Guide, ieee_flags(3M) ```

C99 Floating Point Environment Functions

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

CODE EXAMPLE A-7 C99 Floating Point Environment Functions
 ```#include ``` ```#include ``` ```#include ``` ```#include ``` ```/* ``` ```* 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, compiling and running this program produces the following:

 ```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 ```

The following example shows the effect of the `fesetprec` function on x86. (This function is not available on SPARC.) 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 that evaluate all intermediate results in extended precision. Thus, the `fesetprec` function may be used to guarantee that all results will be rounded to the desired precision, as the second loop shows.

CODE EXAMPLE A-8 `fesetprec` function (x86)
 ```#include ``` ```#include ``` ```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; ``` ```} ```

The output from this program on x86 systems is:

 ```64 significant bits ``` ```53 significant bits ```

Finally, the following code fragment shows one way to use the environment functions in a multi-threaded 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 Solaris Multithreaded Programming Guide for more information on writing multi-threaded programs.)

CODE EXAMPLE A-9 Using Environment Functions in Multi-Thread Program
 ```#include ``` ```#include ``` ```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); ``` ``` ... ``` ```} ```

Exceptions and Exception Handling

`ieee_flags` -- Accrued Exceptions

Generally, a user program examines or clears the accrued exception bits. Here is a C program that examines the accrued exception flags:

CODE EXAMPLE A-10 Examining the Accrued Exception Flags
 ```#include ``` ```#include ``` ```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; ``` ```} ```

The output from running this program:

 ```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) ```

The same can be done from Fortran:

CODE EXAMPLE A-11 Examining the Accrued Exception Flags - 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 ``` ```f77_floatingpoint.h. ``` ```*/ ``` ```#include ``` ``` 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 ```

The output is as follows:

 ``` Highest priority exception is underflow ``` ``` 0 0 0 1 1 ```

While it is unusual for a user program to set exception flags, it can be done. This is demonstrated in the following C example.

 ```#include ``` ```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; ``` ```} ```

On SPARC, the output from the preceding program is:

 ```out is: division , fp exception code is: 2 ```

On x86, the output is:

 ```out is: division , fp exception code is: 4 ```

`ieee_handler` -- Trapping Exceptions

Note – The examples below apply only to the Solaris operating environment.

Here is a Fortran program that installs a signal handler to locate an exception (for SPARC systems only):

CODE EXAMPLE A-12 Trap on Underflow - 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*1.0d300 ``` ``` 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 ``` ```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 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 ```

The output is:

 ``` 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) Here is a more complex C example:

CODE EXAMPLE A-13 Trap on Invalid, Division by 0, Overflow, Underflow, and Inexact - 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 ``` ```#include ``` ```#include ``` ```#include ``` ```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); ``` ```} ```

The output is similar to the following:

 ```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) The following program shows how you can use `ieee_handler` and the include files to modify the default result of certain exceptional situations:

CODE EXAMPLE A-14 Modifying the Default Result of Exceptional Situations
 ```/* ``` ``` * 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 ``` ```#include ``` ```#include ``` ```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; ``` ``` } ``` ```} ```

As expected, the output is:

 ```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` -- Abort on Exceptions

You can use `ieee_handler` to force a program to abort in case of certain floating-point exceptions:

 ```#include ``` ``` 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` Exception Handling Features

The following examples show how to use some of the exception handling features provided by `libm9x.so`. The first example is based on the following task: given a number x and coefficients a0, a1,..., aN, and b0, b1,..., bN-1, evaluate the function f(x) and its first derivative f'(x), where f is the continued fraction

f(x) = a0 + b0/(x + a1 + b1/(x + ... /(x + aN-1 + bN-1/(x + aN))...).

Computing f is straightforward in IEEE arithmetic: even if one of the intermediate divisions overflows or divides by zero, the default value specified by the standard (a correctly signed infinity) turns out to yield the correct result. Computing f', on the other hand, can be more difficult because the simplest form for evaluating it can have removable singularities. If the computation encounters one of these singularities, it will attempt to evaluate one of the indeterminate forms 0/0, 0*infinity, or infinity/infinity, all of which raise invalid operation exceptions. W. Kahan has proposed a method for handling these exceptions via a feature called "presubstitution".

Presubstitution is an extension of the IEEE default response to exceptions that lets the user specify in advance the value to be substituted for the result of an exceptional operation. Using `libm9x.so`, a program can implement presubstitution easily by installing a handler in the `FEX_CUSTOM` exception handling mode. This mode allows the handler to supply any value for the result of an exceptional operation simply by storing that value in the data structure pointed to by the info parameter passed to the handler. Here is a sample program to compute the continued fraction and its derivative using presubstitution implemented with a `FEX_CUSTOM` handler.

CODE EXAMPLE A-15 Computing the Continued Fraction and its Derivative Using the `FEX_CUSTOM` Handler
 ```#include ``` ```#include ``` ```#include ``` ```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; ``` ```} ```

Several comments about the program are in order. On entry, the function `continued_fraction` saves the current exception handling modes for division by zero and all invalid operation exceptions. It then establishes nonstop exception handling for division by zero and a `FEX_CUSTOM` handler for the three indeterminate forms. This handler will substitute infinity for both 0/0 and infinity/infinity, but it will substitute the value of the global variable `p` for 0*infinity. Note that `p` must be recomputed each time through the loop that evaluates the function in order to supply the correct value to substitute for a subsequent 0*infinity invalid operation. Note also that `p` must be declared `volatile` to prevent the compiler from eliminating it, since it is not explicitly mentioned elsewhere in the loop. Finally, to prevent the compiler from moving the assignment to `p` above or below the computation that can incur the exception for which `p` provides the presubstitution value, the result of that computation is also assigned to a `volatile` variable (called `t` in the program). The final call to `fex_setexcepthandler` restores the original handling modes for division by zero and the invalid operations.

The main program enables logging of retrospective diagnostics by calling the `fex_set_log` function. Before it does so, it raises the inexact flag; this has the effect of preventing the logging of inexact exceptions. (Recall that in `FEX_NONSTOP` mode, an exception is not logged if its flag is raised, as explained in the section Retrospective Diagnostics.) The main program also establishes `FEX_ABORT` mode for the common exceptions to ensure that any unusual exceptions not explicitly handled by `continued_fraction` will cause program termination. Finally, the program evaluates a particular continued fraction at several different points. As the following sample output shows, the computation does indeed encounter intermediate exceptions:

 ```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 ```

(The exceptions that occur in the computation of f'(x) at x = 1, 4, and 5 do not result in retrospective diagnostic messages because they occur at the same site in the program as the exceptions that occur when x = -3.)

The preceding program may not represent the most efficient way to handle the exceptions that can occur in the evaluation of a continued fraction and its derivative. One reason is that the presubstitution value must be recomputed in each iteration of the loop regardless of whether or not it is needed. In this case, the computation of the presubstitution value involves a floating point division, and on modern SPARC and x86 processors, floating point division is a relatively slow operation. Moreover, the loop itself already involves two divisions, and because most SPARC and x86 processors cannot overlap the execution of two different division operations, divisions are likely to be a bottleneck in the loop; adding another division would exacerbate the bottleneck.

It is possible to rewrite the loop so that only one division is needed, and in particular, the computation of the presubstitution value need not involve a division. (To rewrite the loop in this way, one must precompute the ratios of adjacent elements of the coefficients in the `b` array.) This would remove the bottleneck of multiple division operations, but it would not eliminate all of the arithmetic operations involved in the computation of the presubstitution value. Furthermore, the need to assign both the presubstitution value and the result of the operation to be presubstituted to `volatile` variables introduces additional memory operations that slow the program. While those assignments are necessary to prevent the compiler from reordering certain key operations, they effectively prevent the compiler from reordering other unrelated operations, too. Thus, handling the exceptions in this example via presubstitution requires additional memory operations and precludes some optimizations that might otherwise be possible. Can these exceptions be handled more efficiently?

In the absence of special hardware support for fast presubstitution, the most efficient way to handle exceptions in this example may be to use flags, as the following version does:

CODE EXAMPLE A-16 Using Flags to Handle Exceptions
 ```#include ``` ```#include ``` ```#include ``` ```/* ``` ```* 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; ``` ```} ```

In this version, the first loop attempts the computation of f(x) and f'(x) in the default nonstop mode. If the invalid flag is raised, the second loop recomputes f(x) and f'(x) explicitly testing for the appearance of a NaN. Usually, no invalid operation exception occurs, so the program only executes the first loop. This loop has no references to `volatile` variables and no extra arithmetic operations, so it will run as fast as the compiler can make it go. The cost of this efficiency is the need to write a second loop nearly identical to the first to handle the case when an exception occurs. This tradeoff is typical of the dilemmas that floating point exception handling can pose.

Using `libm9x.so` with Fortran Programs

`libm9x.so` is primarily intended to be used from C/C++ programs, but by using the Sun Fortran language interoperability features, you can call some `libm9x.so` functions from Fortran programs as well.

Note – For consistent behavior, do not use both the `libm9x.so` exception handling functions and the `ieee_flags` and `ieee_handler` functions in the same program.

The following example shows a Fortran version of the program to evaluate a continued fraction and its derivative using presubstitution (SPARC only):

CODE EXAMPLE A-17 Evaluating a Continued Fraction and its Derivative Using Presubstitution - 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 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 ``` ``` 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 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(x'ff2')) ``` ``` return ``` ``` end ``` ```c ``` ```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/ ``` ``` 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 ```

The output from this program reads:

 ```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) ```

Miscellaneous

`sigfpe` -- Trapping Integer Exceptions

The previous section showed examples of using `ieee_handler`. In general, when there is a choice between using `ieee_handler` or `sigfpe`, the former is recommended.

Note – `sigfpe` is available only in the Solaris operating environment.

(SPARC) There are instances, such as trapping integer arithmetic exceptions, when `sigfpe` is the handler to be used. The following example traps on integer division by zero:

CODE EXAMPLE A-18 Trapping Integer Exceptions
 ``` /* Generate the integer division by zero exception */ ``` ```#include ``` ```#include ``` ```#include ``` ```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); ``` ```/* ``` ``` * Increment the program counter; the operating system does this ``` ``` * automatically for floating-point exceptions but not for ``` ``` * integer division by zero. ``` ``` */ ``` ``` uap->uc_mcontext.gregs[REG_PC] = ``` ``` uap->uc_mcontext.gregs[REG_nPC]; ``` ```} ```

Calling Fortran from C

Here is a simple example of a C driver calling Fortran subroutines. Refer to the appropriate C and Fortran manuals for more information on working with C and Fortran. The following is the C driver (save it in a file named `driver.c)`:

CODE EXAMPLE A-19 Calling Fortran from C
 ```/* ``` ``` * a demo program that shows: ``` ``` * ``` ``` * 1. how to call f77 subroutine from C, passing an array argument ``` ``` * 2. how to call single precision f77 function from C ``` ``` * 3. how to call double precision f77 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; ``` ```} ```

Save the Fortran subroutines in a file named `drivee.f`:

 ``` 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 ```

Then, perform the compilation and linking:

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

The output looks like this:

 ``` 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 ```

Useful Debugging Commands

TABLE A-1 shows examples of debugging commands for the SPARC architecture.

TABLE A-1   Some Debugging Commands (SPARC)
Action `dbx`` ` `adb`
Set breakpoint
at function
at line number

```stop in myfunct stop at 29 ```

```myfunct:b ```
```23a8:b main+0x40:b ```
Run until breakpoint met ```run ``` ```:r ```
Examine source code ```list ``` ```<pc,10?ia ```
Examine a fp register
IEEE single precision
decimal equivalent (Hex)
IEEE double precision
decimal equivalent (Hex)

```print \$f0 print -fx \$f0 print \$f0f1 print -flx \$f0f1 print -flx \$d0 ```

```<f0=X <f0=f <f0=X; <f1=X <f0=F ```
Examine all fp registers ```regs -F ``` ```\$x for f0-f15 \$X for f16-f31 ```
Examine all registers ```regs ``` ```\$r; \$x; \$X ```
Examine fp status register ```print -fx \$fsr ``` ```<fsr=X ```
Put single precision 1.0 in `f0`
Put double prec 1.0 in `f0/f1`
```assign \$f0=1.0 assign \$f0f1=1.0 ``` ```3f800000>f0 3ff00000>f0; 0>f1 ```
Continue execution ```cont ``` ```:c ```
Single step ```step (or next) ``` ```:s ```
Exit the debugger ```quit ``` ```\$q ```

TABLE A-2 shows examples of debugging commands for the x86 architecture

TABLE A-2   Some Debugging Commands (x86)
Action `dbx`` ` `adb `
Set breakpoint
at function
at line number

```stop in myfunct stop at 29 ```

```myfunct:b ```
```23a8:b main+0x40:b ```
Run until breakpoint met ```run ``` ```:r ```
Examine source code ```list ``` ```<pc,10?ia ```
Examine fp registers ```print \$st0 ... print \$st7 ``` ```\$x ```
Examine all registers ```examine &\$gs/19X ``` ```\$r ```
Examine fp status register ```examine &\$fstat/X ``` ```<fstat=X or \$x ```
Continue execution ```cont ``` ```:c ```
Single step ```step (or next) ``` ```:s ```
Exit the debugger ```quit ``` ```\$q ```

.

The following examples show two ways to set a breakpoint at the beginning of the code corresponding to a routine `myfunction` in `adb`. First you can say:

 ```myfunction:b ```

Second, you can determine the absolute address that corresponds to the beginning of the piece of code corresponding to `myfunction`, and then set a break at that absolute address:

 ```myfunction=X ``` ``` 23a8 ``` ```23a8:b ```

The main subroutine in a Fortran program compiled with `f95` is known as `MAIN_` to `adb`. To set a breakpoint at `MAIN_` in `adb`:

```    MAIN_:b
```

When examining the contents of floating-point registers, the hex value shown by the `dbx` command `regs` `-F` is the base-16 representation, not the number's decimal representation. For SPARC, the `adb` commands `\$x` and `\$X` display both the hexadecimal representation, and the decimal value. For x86, the `adb` command `\$x` displays only the decimal value. For SPARC, the double precision values show the decimal value next to the odd-numbered register.

Because the operating system disables the floating-point unit until it is first used by a process, you cannot modify the floating-point registers until they have been accessed by the program being debugged.

(SPARC) When displaying floating point numbers, you should keep in mind that the size of registers is 32 bits, a single precision floating-point number occupies 32 bits (hence it fits in one register), and double precision floating-point numbers occupy 64 bits (therefore two registers are used to hold a double precision number). In the hexadecimal representation 32 bits correspond to 8-digit numbers. In the following snapshot of FPU registers displayed with `adb`, the display is organized as follows:

<name of fpu register> <IEEE hex value> <single precision> <double precision>

(SPARC) The third column holds the single precision decimal interpretation of the hexadecimal pattern shown in the second column. The fourth column interprets pairs of registers. For example, the fourth column of the `f11` line interprets `f10` and `f11` as a 64-bit IEEE double precision number.

(SPARC) Because `f10` and `f11` are used to hold a double precision value, the interpretation (on the `f10` line) of the first 32 bits of that value, `7ff00000`, as `+NaN`, is irrelevant. The interpretation of all 64 bits, `7ff00000 00000000`, as `+Infinity`, happens to be the meaningful translation.

(SPARC) The `adb` command `\$x`, that was used to display the first 16 floating-point data registers, also displayed `fsr` (the floating-point status register):

 ````\$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) The corresponding output on x86 looks like:

 ``````\$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 ```

Note – (x86) `cw` is the control word; `sw` is the status word.