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.