ISO/IEC JTC1/SC22/WG5 - N1149 To: WG5 From: John Reid Subject: Floating point exceptions The development body has commenced work and has written a first draft technical report based on the enable construct. The development body decided that it must develop a full proposal at least in outline in order to be sure that the limited proposal is consistent with it. The full proposal is based on the Edinburgh proposal, but with an important addition: as well as intrinsic conditions, conditions may be declared by Fortran code. This was a wish expressed strongly by WG5 in Edinburgh, but no way was seen at that time to satisfy this requirement without significantly enlarging the proposal. However, the suggestion of Richard Hanson to use a derived type in an intrinsic module, provides this and is actually an overall simplification. Because it is a simplification, we have retained it in the subset. Both the full and subset proposals are explained in a recent article in Fortran Forum (Sept 1995), but more work is needed and I hope that members of the development body can get together in San Diego. An alternative to the enable construct is to use a collection of procedures and the development body was asked by WG5 to pursue this approach too. I have constructed a draft technical report based on this approach and it is attached (without the detailed edits to the standard). I have given the development body a preview, but it has not had long enough yet to consider it carefully. The procedures approach is easy to understand and is quite close to what is available in C. Its main defect is that, without the enable and handle blocks, there is no way to scope statically extra testing (such as for integer overflow) nor is there is a branch target for when things go seriously wrong. Nevertheless, I would be far more confident of meeting our deadlines if the technical report is based on this approach. I think we should seriously consider putting enable on one side for later development and concentrate on the procedures approach. ................................................................. DRAFT PROCEDURES TECHNICAL REPORT 0. NOTES This is a draft Technical Report using procedures in a module for exception handling. It is based on Keith Bierman's message of 7 September. For IEEE hardware, it supports exceptions and other aspects of IEEE arithmetic. Non-IEEE processors are required to support only OVERFLOW and DIVIDE_BY_ZERO, and may support more. We could reduce this proposal to one that supports exceptions only. This would involve removal of the data type IEEE_RND_VALUE and the reduction of the set of procedures to IEEE_GET_FLAG(FLAG), IEEE_CLEAR_FLAGS, IEEE_SET_FLAGS, IEEE_GET_FSR, IEEE_HALT, and IEEE_SET_FSR. More work is needed yet on the interaction of these features with PURE and FORALL. 1. RATIONALE Exception handling is required for the development of robust and efficient numerical software. In particular, it is necessary in order to be able to write portable scientific libraries. In numerical Fortran programming, current practice is to employ whatever exception handling mechanisms are provided by the system/vendor. This clearly inhibits the production of fully portable numerical libraries and programs. It is particularly frustrating now that IEEE arithmetic is so widely used, since built into it are the five conditions: overflow, invalid, divide_by_zero, underflow, and inexact. Our aim to to provide support for these conditions. This proposal involves a standard module, IEEE_ARITHMETIC, containing three derived types, some named constants of these types, and a collection of simple procedures. They allow the flags to be tested, cleared, set, saved, or restored. To facilitate maximum performance, each of the proposed functions does very little processing of arguments. In many cases, a processor may generate only a few inline machine code instructions rather than library calls. In order to allow for the maximum number of processors to provide the maximum value to users, we do NOT require complete IEEE conformance. What a vendor must do is to provide the module, and support the overflow and divide-by-zero flags and the basic inquiry functions. A user must utilize the inquiry functions to determine if they can count on specific features of the IEEE standard, before using other inquiry or value setting procedures. Note that a processor ought not implement these as "macros", as IEEE conformance is often controlled by compiler switches. A processor which offers a switch to turn off a facility should adjust the values returned for these inquiries. For example, a processor which allows gradual underflow to be turned off (replaced with flush to zero) should return .FALSE. for IEEE_SUPPORT_DENORMAL(X) when a source file is processed with that option on. Naturally it should return .TRUE. when that option is not in effect. The most important use of a floating-point exception handling facility is to make possible the development of much more efficient software than is otherwise possible. The following "hypotenuse" function, sqrt(x**2 + y**2), illustrates the use of the facility in developing efficient software. REAL FUNCTION HYPOT(X, Y) ! In rare circumstances this may lead to the setting of the OVERFLOW flag USE IEEE_ARITHMETIC REAL X, Y REAL SCALED_X, SCALED_Y, SCALED_RESULT LOGICAL OLD_OVERFLOW, OLD_UNDERFLOW INTRINSIC SQRT, ABS, EXPONENT, MAX, DIGITS, SCALE ! Store the old flags and clear them OLD_OVERFLOW = IEEE_GET_FLAG(OVERFLOW) OLD_UNDERFLOW = IEEE_GET_FLAG(UNDERFLOW) CALL IEEE_CLEAR_FLAGS(UNDERFLOW, OVERFLOW) ! Try a fast algorithm first HYPOT = SQRT( X**2 + Y**2 ) IF (IEEE_GET_FLAG(UNDERFLOW) .OR. IEEE_GET_FLAG(OVERFLOW) ) THEN CALL IEEE_CLEAR_FLAGS(UNDERFLOW, OVERFLOW) IF ( X==0.0 .OR. Y==0.0 ) THEN HYPOT = ABS(X) + ABS(Y) ELSE IF ( 2*ABS(EXPONENT(X)-EXPONENT(Y)) > DIGITS(X)+1 ) THEN HYPOT = MAX( ABS(X), ABS(Y) )! one of X and Y can be ignored ELSE ! scale so that ABS(X) is near 1 SCALED_X = SCALE( X, -EXPONENT(X) ) SCALED_Y = SCALE( Y, -EXPONENT(X) ) SCALED_RESULT = SQRT( SCALED_X**2 + SCALED_Y**2 ) HYPOT = SCALE( SCALED_RESULT, EXPONENT(X) ) ! may cause overflow END IF END IF IF(OLD_OVERFLOW) CALL IEEE_SET_FLAGS(OVERFLOW) IF(OLD_UNDERFLOW) CALL IEEE_SET_FLAGS(UNDERFLOW) END FUNCTION HYPOT An attempt is made to evaluate this function directly in the fastest possible way. (Note that with hardware support, exception checking is very efficient; without language facilities, reliable code would require programming checks that slow the computation significantly.) The fast algorithm will work almost every time, but if an exception occurs during this fast computation, a safe but slower way evaluates the function. This slower evaluation may involve scaling and unscaling, and in (very rare) extreme cases this unscaling can cause overflow (after all, the true result might overflow if X and Y are both near the overflow limit). If the overflow or underflow flag is set on entry, it is reset on return, so that earlier exceptions are not lost. Can all this be accomplished without the help of an exception handling facility? Yes, it can - in fact, the alternative code can do the job, but of course it is much less efficient. That's the point. The HYPOT function is special, in this respect, in that the normal and alternative codes try to accomplish the same task. This is not always the case. In fact, it very often happens that the alternative code concentrates on handling the exceptional cases and is not able to handle all of the non-exceptional cases. When this happens, a program which cannot take advantage of hardware flags could have a structure like the following: if ( in the first exceptional region ) then handle this case else if ( in the second exceptional region ) then handle this case : else execute the normal code end But this is not only inefficient, it also "inverts" the logic of the computation. For other examples, see Hull, Fairgrieve and Tang (1994) and Demmel and Li (1994). The HYPOT algorithm is used in example 4 of section 15.10. The code for the HYPOT function can be generalized in an obvious way to compute the Euclidean Norm, sqrt( x(1)**2 + x(2)**2 + ... + x(n)**2 ), of an n-vector; the generalization of the alternative code is not so obvious (though straightforward) and will be much slower relative to the normal code than is the case with the HYPOT function. Of the five IEEE conditions, underflow and inexact are obviously milder than the other three (inexact is particularly mild and signals almost all the time in typical codes). We therefore group the other three together as USUAL. There is a need for further intrinsic conditions in connection with reliable computation. Examples are a. INSUFFICIENT_STORAGE for when the processor is unable to find sufficient storage to continue execution. b. INTEGER_OVERFLOW and INTEGER_DIVIDE_BY_ZERO for when an intrinsic integer operation has a very large result or has a zero denominator. c. INTRINSIC for when an intrinsic procedure has been unsuccessful. d. SYSTEM_ERROR for when a system error occurs. This proposal has been designed to allow such enhancements in the future. References Demmel, J.W. and Li, X. (1994). Faster Numerical Algorithms via Exception Handling. IEEE Transactions on Computers, 43, no. 8, 983-992. Hull, T.E., Fairgrieve, T.F., and Tang, T.P.T. (1994). Implementing complex elementary functions using exception handling. ACM Trans. Math. Software 20, 215-244. 2. TECHNICAL SPECIFICATION This proposal involves a standard module, IEEE_ARITHMETIC, containing three derived types: 1. IEEE_FLAG, with the constants INVALID, OVERFLOW, DIVIDE_BY_ZERO, USUAL, UNDERFLOW, INEXACT, and ALL. 2. IEEE_RND_VALUE, with the constants NEAREST, TO_ZERO, UP, DOWN, and OTHER 3. IEEE_FSR_VALUE, for saving the current floating point status. and a collection of simple procedures. There must be flags that are initially clear and that are set when an exception occurs. The value of a flag is determined by the function IEEE_GET_FLAG(flag) Flags may be cleared by the subroutine IEEE_CLEAR_FLAGS(flag-list) and may be set by the subroutine IEEE_SET_FLAGS(flag-list) The minimum requirement is for the support of OVERFLOW and DIVIDE_BY_ZERO. Whether the other exceptions are supported may be determined by the inquiry function IEEE_SUPPORT_FLAG(FLAG) The extent of further support for the IEEE standard may be determined by the inquiry functions: IEEE_DATATYPE (X) Inquire if the processor supports IEEE arithmetic for reals of the same kind type parameter as the argument. IEEE_SUPPORT_ALL() Inquire if processor supports all the IEEE facilities defined in this standard for all kinds of reals. IEEE_SUPPORT_DENORMAL(X) Inquire if the processor supports IEEE gradual underflow for reals of the same kind type parameter as the argument. IEEE_SUPPORT_FLAG(FLAG) Inquire if the processor supports an exception or set of exceptions. IEEE_SUPPORT_HALTING() Inquire if the processor supports the ability to control during program execution whether to abort or continue execution after an exception. IEEE_SUPPORT_INF(X) Inquire if processor supports the IEEE infinity facility for reals of the same kind type parameter as the argument. IEEE_SUPPORT_NAN(X) Inquire if processor supports the IEEE Not-A-Number facility for reals of the same kind type parameter as the argument. IEEE_SUPPORT_ROUNDING (RND_VALUE) Inquire if processor supports changing to a particular rounding mode for IEEE kinds of reals. IEEE_SUPPORT_SQRT(X) Inquire if the processor supports IEEE square root for reals of the same kind type parameter as the argument. Implementations with the appropriate extent of IEEE support provide: A. Elemental functions: IEEE_COPYSIGN(X,Y) IEEE copysign function. IEEE_INFINITY(X) Generate IEEE infinity with the same sign as X. IEEE_IS_DENORMAL(X) Determine if value is IEEE denormalized. IEEE_IS_INF(X) Determine if value is IEEE Infinity. IEEE_IS_NAN(X) Determine if value is IEEE Not-a-Number. IEEE_IS_NORMAL(X) Whether a value is normal, that is, neither an Infinity, a NaN, nor denormalized. IEEE_LOGB(X) Unbiased exponent in the IEEE floating point format. IEEE_NAN(X) Generate IEEE Not-a-Number. IEEE_NEXTAFTER(X,Y) Returns the next representable neighbor of X in the direction toward Y. IEEE_SCALB (X,I) Returns X * 2**I. IEEE_SQRT(X) Square root as defined in the IEEE standard. IEEE_UNORDERED(X,Y) IEEE unordered function. True if X or Y is a NaN and false otherwise. B. Scalar functions: IEEE_GET_FSR() Save the current values of the set of flags that define the current state of the floating point environment (usually the floating point status register). The result is of type IEEE_FSR_VALUE. IEEE_GET_ROUNDING_MODE() Save the current IEEE rounding mode. The result is of type IEEE_RND_VALUE. C. Subroutines: IEEE_CLEAR_FLAGS (FLAG1 [, FLAG2, ...]) Clear one or more exception flags. IEEE_HALT(FLAG,HALTING) Controls continuation or halting on exceptions. FLAG is of type IEEE_FLAG and value INVALID, OVERFLOW, DIVIDE_BY_ZERO, USUAL, UNDERFLOW, INEXACT, or ALL. IEEE_SET_FLAGS(FLAG1 [, FLAG2, ...]) Set one or more exception flags. IEEE_SET_FSR(FSR_VALUE) Restore the values of the set of flags that define the current state of the floating point environment (usually the floating point status register). FSR_VALUE is of type IEEE_FSR_VALUE and has been set by a call of IEEE_GET_FSR. IEEE_SET_ROUNDING_MODE(RND_VALUE) Set the current IEEE rounding mode. RND_VALUE is of type IEEE_RND_VALUE, with value NEAREST, TO_ZERO, UP, or DOWN. Note: An IEEE processor is required to support single precision, and an extended precision. This is usually REAL and DOUBLE PRECISION but it might not always be. In addition, some vendors support a "REAL*16" which does not follow the "logical extension" to the IEEE 754 standard (e.g. IBM's RS6K "doubled double" implementation of REAL*16). Users should be able to determine if a particular datatype is going to behave as "one would expect". But implementors should not be constrained to supply only datatypes which conform to IEEE 754 (for performance, or compatibility with older processors or eventual improvements in state of the art in floating point implementation). If a processor supports a particular IEEE exception, the associated flag is always clear. 3. EDITS TO THE DRAFT STANDARD (N1122) These are written but to keep this paper short are not included here.