A P P E N D I X A |
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.
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:
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:
On SPARC, the corresponding output reads:
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:
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:
This Fortran example uses some functions recommended by the IEEE standard:
The output from this program is shown in CODE EXAMPLE A-6.
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 |
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; See the Numerical Computation Guide, ieee_flags(3M) |
This C program calls several of the ieee_values(3m) functions:
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 000fffff single precision min subnormal = 1.40129846e-45 = 00000001 |
Fortran programs that use ieee_values functions should take care to declare those functions' types:
On SPARC, the output reads as follows:
6.4751751194380251109244389582276466-4966 Inf 3.40282E+38 |
The following example demonstrates how to set the rounding mode to round towards zero:
(SPARC) The output of this short program shows the effects of rounding towards zero:
(x86) The output of this short program shows the effects of rounding towards zero:
To set rounding direction towards zero from a Fortran program:
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.
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.
On SPARC, compiling and running this program produces the following:
CODE EXAMPLE A-8 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.
The output from this program on x86 systems is:
64 significant bits 53 significant bits |
Finally, CODE EXAMPLE A-9 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.)
Generally, a user program examines or clears the accrued exception bits. CODE EXAMPLE A-10 is a C program that examines the accrued exception flags.
The output from running this program:
The same can be done from Fortran:
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.
On SPARC, the output from the preceding program is:
out is: division , fp exception code is: 2 |
out is: division , fp exception code is: 4 |
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):
(SPARC) Here is a more complex C example:
The output is similar to the following:
(SPARC) CODE EXAMPLE A-14 shows how you can use ieee_handler and the include files to modify the default result of certain exceptional situations:
You can use ieee_handler to force a program to abort in case of certain floating-point exceptions:
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
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.
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 .) 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:
(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:
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 trade-off is typical of the dilemmas that floating point exception handling can pose.
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):
The output from this program reads:
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. CODE EXAMPLE A-18 traps on integer division by zero.
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):
Save the Fortran subroutines in a file named drivee.f:
Then, perform the compilation and linking:
cc -c driver.c f95 -c drivee.f demo_one: demo_two: demo_three: f95 -o driver driver.o drivee.o |
TABLE A-1 shows examples of debugging commands for the SPARC architecture.
TABLE A-2 shows examples of debugging commands for the x86 architecture.
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):
(x86) The corresponding output on x86 looks like:
Copyright © 2003, Sun Microsystems, Inc. All rights reserved.