Math Library Basics

Math libraries in most computer systems, including HP-UX, are collections of frequently used mathematical functions. The functions take one or more arguments and return one or more results. When an application source file contains a use of a math function, the compiler automatically generates a call to the appropriate routine name in the appropriate math library. For example, suppose your Fortran program contains the following declaration and statement:

       DOUBLE PRECISION A, B, Y
          .
          .
          .
       Y = A**B

The statement raises A to the power B and assigns the result to Y. The Fortran compiler emits a call to FTN_DTOD, which is one of the functions in libcl.a and libcl.sl. libcl is the Fortran and Pascal library.

Conceptually, the definition of a math library function is simple; in this example, FTN_DTOD merely raises a DOUBLE PRECISION value to the power of another DOUBLE PRECISION value. However, there are practical questions about math library functions that the programmer should consider:

How accurate is the result returned?
How efficient is the function (that is, how fast does it run)?
What happens when an error occurs (for example, overflow or invalid arguments)?
What are the calling conventions for the function?
What functions are available?

The answers to questions 1 and 2 are implementation-dependent. The answers to questions 3, 4, and 5 are determined by a variety of programming environment specifications, such as XPG4.2.

Appendix A “The C Math Library” partially answers question 5 by listing the functions in the HP-UX C math library. The Fortran equivalents are the intrinsic functions; see the HP Fortran 90 Programmer's Reference or the HP FORTRAN/9000 Programmer's Reference for a list of these functions. The following section, which describes what happens when a program calls a math library function, provides some answers to questions 1 through 4.

Anatomy of a Math Library Function Call

Figure 4-1 “Anatomy of a Math Library Call” shows a generalized flowchart of a math library function call, applicable to all languages and standards. The figure and the discussion that follows assume that the function takes one argument, but they also apply to functions that take more than one argument.

Figure 4-1 Anatomy of a Math Library Call

Steps A through E of this process are described below.

Step A

The compiler-generated call to a math library function includes code that

Converts the argument to the required format, if necessary
Places the argument in the appropriate place (that is, where the function expects it)
Calls the function

Step B

The function determines whether the argument is valid. All functions check for NaN arguments and make additional checks specific to the function. For example, a logarithmic function checks for negative or zero arguments; an exponential function checks to see if the argument is so large that the result would overflow. If the argument is valid, control passes to Step C, execution of the function. If the argument is not valid, control passes to Step D.

Step C

The execution of the function involves performing the appropriate transformation on the argument(s). For most math functions, the result has the same data type as the argument, which is usually double-precision.

The execution of the function actually consists of one of many possible paths through the function code. The nature of the argument usually determines the path taken. For example, there may be special paths for particular argument ranges, or there may be an argument reduction phase. Math library implementors tune the most frequently executed paths for optimal efficiency. Consequently, arguments that are unusual (for example, very large or very small) can cause noticeable performance degradation.

Step D

If an argument proves to be invalid, control passes to the error-handling function of the library. When this happens, performance slows considerably; usually it is 2 to 3 orders of magnitude slower than for a valid argument. This follows from the basic assumption in math library design that errors are rare events.

If you do not supply an error-handling function in your program, a math library call that encounters an illegal argument does some or all of the following, depending on which programming language you are using and which standards are being enforced:

Supplies a system-defined default result: NaN for invalid operations, a huge value for overflows, and so on
Sets the globally accessible error code variable errno
Sets some state in the hardware floating-point status register
Prints an error message to stderr
Returns to the calling application, returning the default result

Your program may cause any of the above steps to be modified or stopped and may also provide an error-handling subroutine or function to be invoked. “Math Library Error Handling for C” describes how these steps are performed for the C programmer in the XPG4.2 and SVID environments; “Math Library Error Handling for Fortran” describes how they are performed for the Fortran programmer. Users of other languages should refer to the appropriate language manual for details.

Step E

When the function exits, the compiler-generated call to the function passes the result from where the function left it to where it is needed, converting it to the required format if necessary.

Math Library Error Handling for C

The XPG4.2 and SVID specifications specify similar math library error handling, so we describe them together. The XPG4.2 specification is a superset of the ANSI C standard.

To differentiate between XPG and SVID, HP-UX in the past implemented two separate C math libraries, libm and libM. The SVID libraries were libm.a and libm.sl; the XPG libraries were libM.a and libM.sl.

Because the SVID3 and XPG4.2 standards are essentially identical, HP-UX now supplies only one library, libm. The libM library is now obsolete. All versions of libM are provided only as soft links to the corresponding versions of libm (see “Locations of the Math Libraries at Release 10.30” for details).

The flowchart in Figure 4-2 “C Math Library Error Handling for the sqrt Function” summarizes the error handling in the SVID and XPG4.2 math libraries. We use the sqrt function as an example.

Figure 4-2 C Math Library Error Handling for the sqrt Function

If a C library function such as sqrt encounters an invalid argument, it ordinarily returns a default result that indicates a failure—a result that the function could not ordinarily return. The default result is usually a NaN, zero, or HUGE_VAL, depending on the function and the argument value.

In addition, some functions set the global value errno, defined in the header file errno.h. (Many HP-UX system calls also set this value.) When a library call fails, the value of errno may be set to an appropriate code, also defined by errno.h. The math library functions set errno to either EDOM or to ERANGE. EDOM indicates a domain error—that is, an error in the argument. ERANGE indicates a range error—usually an overflow in the result.

To detect errors, an application should check both the returned value and errno. It may set errno to 0 (that is, no error) at the beginning of a section of code and then examine it after one or more library calls to see if it is nonzero. A nonzero value indicates that some sort of library error has occurred.

For example, you could use the following C code fragment to print an error message when a call to pow fails:

errno = 0;
x = z * pow(y, w);
if (isnan(x) || errno)
    fprintf(stderr, "error in pow function: errno = %d\n", errno);

If y is negative and w is not an integer value, this fragment would print out

error in pow function: errno = 33

If you look up 33 in the system include file /usr/include/sys/errno.h, you find that it indicates a domain error (EDOM).




	NOTE: C math library functions formerly used a function called `matherr`, which was required by the SVID2 specification but is not specified by ANSI C, SVID3, or XPG4.2. At HP-UX Release 10.30, this function is no longer provided.

Math Library Error Handling for Fortran

If a Fortran intrinsic function encounters an invalid argument, it returns a default result, just as a C function does. The default result depends on both the function and the nature of the argument. For example, a negative argument to the DLOG function causes it to return a NaN value. The most generally useful method of detecting errors is for the application to check for an anomalous result and then to take appropriate action.

What happens after the error depends on the following factors:

Whether an exception trap for the error is enabled
Whether the program contains an ON EXTERNAL ERROR statement (HP FORTRAN/9000 only)

We discuss the possible sequences of events in the following sections. The flowchart in Figure 4-3 “Fortran 77 Math Library Error Handling” summarizes HP FORTRAN/9000 math library error handling. The flowchart in Figure 4-4 “Fortran 90 Math Library Error Handling” summarizes HP Fortran 90 math library error handling.

Figure 4-3 Fortran 77 Math Library Error Handling

Figure 4-4 Fortran 90 Math Library Error Handling

Enabling Exception Traps for Invalid Operations

You can enable an exception trap for any of the five floating-point exception conditions described in “Exception Conditions”. The exception condition that indicates an invalid argument to a math library function is the invalid operation condition. You can enable a trap for this condition in any of several ways:

By compiling with the option +FPV, which enables the invalid operation trap
(HP Fortran 90) By compiling with the option +fp_exception, which enables traps for invalid operation, overflow, underflow, and division by zero.
(HP FORTRAN/9000) By compiling with the option +T, which enables traps for invalid operation, overflow, underflow, and division by zero. (If your program contains an ON EXTERNAL ERROR statement, you must use this option.)
By calling the fesettrapenable routine with the argument FE_INVALID, which enables the invalid operation trap

See Chapter 5 “Manipulating the Floating-Point Status Register” and Chapter 6 “Floating-Point Trap Handling” for information about all of these methods.

If you do not enable a trap for invalid operations, all that happens in the case of an invalid argument is that the invalid operation exception flag in the status register is set and the default result is returned. You can retrieve the value of the exception flags by calling the fetestexcept routine, described in “Exception Bits”.

The ON EXTERNAL ERROR Statement (HP FORTRAN/9000 only)




	NOTE: HP Fortran 90 does not support the `ON EXTERNAL ERROR` statement, although it does support the `ON` statement for other kinds of error handling.

If you use +T to enable a trap for invalid operations, what happens next depends on whether your program contains an ON EXTERNAL ERROR statement. If it does not, the function merely returns the default result.

If your program contains an ON statement, you must compile with the +T option in order to enable trap handling. If your program contains an ON statement and you do not specify +T, you get a compile-time warning.

If your program does contain an ON statement, what happens depends on the action you specify in the statement. You can specify any of the following:

ABORT. If you specify ABORT, the program exits.
IGNORE (usually not a good idea). If you specify IGNORE, the default result is returned.
CALL sub (call a subroutine). If you specify CALL sub, the user-defined trap-handling subroutine sub is called. The subroutine must have three or four arguments (four if the function takes two arguments); see the HP FORTRAN/9000 Programmer's Guide for details. If the subroutine sets the first argument to 0, the program exits. Otherwise, the default result is returned.




	NOTE: A subroutine that handles a math library error takes a different number of arguments from a subroutine that handles an IEEE exception (as described in “Using the ON Statement (Fortran only)”). See the HP FORTRAN/9000 Programmer's Guide for details.

A Program Example

The following program illustrates Fortran 77 math library error handling. It calls DLOG with a negative argument, which is invalid.

See “Run-Time Mode Control: The fenv(5) Suite” for information about the fetestexcept and fegettrapenable routines.

Example 4-1 Sample Program: liberr77.f

C  HP FORTRAN/9000 only
C  Compile
C    1) As is;
C    2) With comment removed from ON EXTERNAL ERROR ABORT
C       statement, and with +T
C    3) With comment removed from ON EXTERNAL ERROR CALL MYSUB
C       statement, and with +T
 
      PROGRAM LIBERR77
$ALIAS FETESTEXCEPT = 'fetestexcept' (%val)
$ALIAS FEGETTRAPENABLE = 'fegettrapenable'
 
      INTEGER FE_INEXACT, FE_UNDERFLOW, FE_OVERFLOW
      INTEGER FE_DIVBYZERO, FE_INVALID, FE_ALL_EXCEPT
      PARAMETER (FE_INEXACT = Z'08000000')
      PARAMETER (FE_UNDERFLOW = Z'10000000')
      PARAMETER (FE_OVERFLOW = Z'20000000')
      PARAMETER (FE_DIVBYZERO = Z'40000000')
      PARAMETER (FE_INVALID = Z'80000000')
      PARAMETER (FE_ALL_EXCEPT = Z'f8000000')
      EXTERNAL FETESTEXCEPT, FEGETTRAPENABLE
      INTEGER  FETESTEXCEPT, FEGETTRAPENABLE
      DOUBLE PRECISION X, Y
      LOGICAL TEST
      INTEGER FLAGS, TRAPS
 
C  Abort if a function results in an error
C      ON EXTERNAL ERROR ABORT
 
C  Call mysub if a function results in an error
C  Other possible action is IGNORE
C  Setting I to 0 aborts the program
C      ON EXTERNAL ERROR CALL MYSUB
      X = 1.2345D0
      X = X*1.1D0
      Y = DLOG(0.0D0-X)
      FLAGS = FETESTEXCEPT(FE_ALL_EXCEPT)
      TRAPS = FEGETTRAPENABLE()
      WRITE(*,10) Y, FLAGS, TRAPS
      TEST = FLAGS .AND. FE_INVALID
      IF (TEST) THEN
        PRINT *, 'invalid operation occurred'
      ENDIF
 10   FORMAT(G, Z10.8, Z10.8)
      END
 
      SUBROUTINE MYSUB(I, A, X)
      INTEGER I
      DOUBLE PRECISION A, X
      PRINT *, 'error no. is ', I
      PRINT *, 'result is ', A
      PRINT *, 'arg is ', X
      I = 0
      RETURN
      END

If you compile this program with the ON statements commented out and without the +T option, no traps are enabled. Therefore, the math library error sets the appropriate exception flag. The output of the program is as follows:

$ f77 liberr77.f -lm
liberr77.f:
   MAIN liberr:
   mysub:
$ a.out
NaN                        88000000  00000000
 invalid operation occurred

The exception flag value indicates that both an invalid operation and an inexact result condition were generated. (Table 5-1 “IEEE Exception Bits” shows the bit values for each of these conditions.).

If you compile with +T, the exception flags that correspond to the traps set by +T are not set after an error, and it is not useful to check them.

If you remove the comment from the ON EXTERNAL ERROR ABORT statement and compile with +T, the invalid operation trap is enabled. When you run the program, it exits:

$ f77 +T liberr77.f -lm
liberr77.f:
   MAIN liberr:
   mysub:
$ a.out
$

Finally, if you remove the comment from the ON EXTERNAL ERROR CALL MYSUB statement and compile with +T, the program calls the error-handling subroutine MYSUB before exiting:

$ f77 +T liberr77.f -lm
liberr77.f:
   MAIN liberr:
   mysub:
$ a.out
 error no. is  8
 result is  NaN
 arg is  -1.35795
$

HP does not support the use of the ON statement to handle math library errors in Fortran 90. The following is a Fortran 90 version of the program.

Example 4-2 Sample Program: liberr90.f

C  Compile
C    1) As is;
C    2) With the +fp_exception option
 
      PROGRAM LIBERR90
!$HP$ ALIAS FETESTEXCEPT = 'fetestexcept' (%val)
!$HP$ ALIAS FEGETTRAPENABLE = 'fegettrapenable'
 
      PARAMETER (FE_INEXACT = Z'08000000')
      PARAMETER (FE_UNDERFLOW = Z'10000000')
      PARAMETER (FE_OVERFLOW = Z'20000000')
      PARAMETER (FE_DIVBYZERO = Z'40000000')
      PARAMETER (FE_INVALID = Z'80000000')
      PARAMETER (FE_ALL_EXCEPT = Z'f8000000')
      EXTERNAL FETESTEXCEPT, FEGETTRAPENABLE
      INTEGER  FETESTEXCEPT, FEGETTRAPENABLE
      DOUBLE PRECISION X, Y
      INTEGER FLAGS, TRAPS 
 
      X = 1.2345D0
      X = X*1.1D0
      Y = DLOG(0.0D0-X)
      FLAGS = FETESTEXCEPT(FE_ALL_EXCEPT)
      TRAPS = FEGETTRAPENABLE()
      WRITE(*,10) Y, FLAGS, TRAPS
 10   FORMAT(G, Z10.8, Z10.8)
      END

If you compile the program as we show it, it runs as follows. Because no traps are enabled, the function returns a NaN and sets the appropriate exception flag.

$ f90 liberr90.f -lm
liberr90.f
   program LIBERR90
 
29 Lines Compiled
$ ./a.out
NaN                        88000000  00000000

If you compile the program with the +fp_exception option, the invalid operation trap is enabled, and the program aborts with an error message:

$ f90 +fp_exception liberr90.f -lm
liberr90.f
program LIBERR90
 
15 Lines Compiled
$ ./a.out
PROGRAM ABORTED :  IEEE invalid operation
 
PROCEDURE TRACEBACK:
( 0)  0x00009a0c   _start + 0x64  [./a.out]