Implementing complex elementary functions using exception handling.
Algorithms are developed for reliable and accurate evaluations of the complex elementary functions required in FORTRAN 77 and FORTRAN 9, namely, cabs, csqrt, cexp, clog, csin, and ccos. The algorithms are presented in a pseudocode that has a convenient exception-handling facility. A tight error bound is derived for each algorithm. Corresponding FORTRAN programs for an IEEE environment have also been developed to illustrate the practicality of the algorithms, and these programs have been tested very carefully to help confirm the correctness of the algorithms and their error bounds. The results of these tests are included in the paper, but the FORTRAN programs are not.
INTRODUCTION
Our purpose is to develop algorithms, along with error bounds, for reliable and accurate
evaluations of the complex elementary functions required in Fortran 77 [1] and
Fortran 90 [4], namely cabs, csqrt, cexp, clog, csin, and ccos. These (seemingly oxy-
moronic) complex elementary functions can be expressed in terms of formulas involving
only real arithmetic and real elementary functions. Complex arithmetic is not needed.
If care is taken these formulas can usually be arranged so that serious numerical cancellation
will not occur during their evaluation. (If this cannot be arranged, higher
precision may be necessary at such critical points in the calculations.)
This work was supported by the Natural Sciences and Engineering Research Council of Canada and
the Information Technology Research Centre of Ontario, as well as by the Applied Mathematical
Sciences subprogram of the Office of Energy Research, U. S. Department of Energy, under Contract
W-31-109-Eng-38, and by the STARS Program Office, under AF Order RESD-632.
Authors' addresses: T. E. Hull and T. F. Fairgrieve, Department of Computer Science, University of
Toronto, Toronto, Ontario, Canada M5S 1A4 (e-mail: ftehull,tffg@cs.utoronto.ca); Ping Tak Peter
Tang, Mathematics and Computer Science Division, Argonne National Laboratory, 9700 South Cass
Ave., Argonne, IL U.S.A. 60439-4801 (e-mail: tang@antares.mcs.anl.gov)
The main difficulty that remains in such evaluations is the possibility that overflow
or underflow might occur at some intermediate stages in the calculation. Such exceptions
are often "spurious", in the sense that the final mathematical result is within
the range of machine representable numbers. In these circumstances, which normally
occur only rarely, the algorithms must provide for alternative calculations, which may
be more lengthy, but which are able to circumvent the spurious exceptional situations.
All of this suggests designing the algorithms with the help of exception handling
facilities. Each algorithm would then begin with a direct evaluation of the original
carefully arranged formulas. This would be efficient, and almost always successful. But
if an exception should occur during such a calculation, it would cause control of the
calculations to be transferred to an exception handler that would do what was needed
to circumvent the difficulty, if possible, or, if not possible, to cause an appropriate
exception to be returned by the function itself. To capture this essential feature of our
algorithms, we present them in a pseudo-code which possesses a convenient exception
handling facility.
Section 2 provides basic information on error bounds that we assume about the
real arithmetic and real elementary functions used in evaluating the complex func-
tions. Special attention is given to the real sine and cosine functions. Three functions
needed later are introduced, and our conventions about exceptions returned by the
complex functions are also specified precisely. Section 3 introduces the exception handling
construct used in our pseudo-code. The pseudo-code algorithms are presented in
Section 4, along with error analyses. Special implementations for testing are described
in Section 5, and formulas for the error bounds derived in Section 4 are tabulated in
Section 5, but extensive testing is done only for IEEE binary arithmetic [3] on Sun
systems [6]. Issues to be considered in production implementations are discussed in
Section 6, and concluding remarks are given in Section 7.
We should emphasize that we assume throughout that the function arguments are
exact. If an argument z is only slightly in error, the induced error in the corresponding
value of the function f(z) is generally small. But in special situations, the error can
also be very large. The relative error in f(z); \Deltaf (z)=f(z); is approximately zf 0 (z)=f(z)
times the relative error in z; \Deltaz=z: The relative error magnification factor zf 0 (z)=f(z)
can be very large. For example, in one extreme case, when f(z) is log(z); the factor is
1= log(z) which can be arbitrarily large when z is near 1:
Complex exponentiation is not included in this paper; it is difficult to develop an algorithm
that does much better than simply evaluating cexp(w   clog(z)) to approximate
z   \Lambdaw: We hope to consider this problem in a separate paper.
An earlier version of this paper appeared as an Argonne Preprint [2].
2. BASIC NUMERICAL OPERATIONS
The complex elementary functions described in Section 4 depend on real arithmetic
operations and on real elementary functions. The main purpose of this section is
to introduce a notation for the errors and error bounds associated with these real
operations which can be used later to derive error bounds for the complex functions.
bounds associated with the real sine and cosine functions are discussed further
in a separate subsection, and the implementation of log1p(x)
considered separately, since log1p is not always available as one of the standard real
functions. Two special functions for manipulating exponents are also specified. Possible
exceptions for complex functions are specified in a final subsection.
2.1 Errors and Error Bounds
We assume that input values for the complex elementary functions are normalized complex
floating-point numbers. With this assumption, it can be arranged that all internal
operations which might generate errors (such as rounding errors) will operate only on
normalized real floating-point numbers. The associated error analyses are mostly relatively
straightforward, and produce results which are simple and easy to use. (It would
also be a straightforward matter to make provision for the input or output of other
values, such as complex numbers whose components could be denormalized, but then
the error analyses would be very much more complicated. We also do not distinguish
between signed zeros; for more about the latter in this context, see Kahan [5].)
Except for cabs, whose output is to be a normalized real floating-point number (if
overflow is not returned), the output in each example is to be a normalized complex
floating-point number (if no exception is returned).
The main assumption we make about the real arithmetic is that, if x and y are
normalized real floating-point numbers, and op is one of the four basic arithmetic
operations, then there is a relative roundoff error bound E such that
f l (x op
for some ffl where j ffl j - E, provided no exception occurs. Here f l (x op y) is the rounded
floating-point approximation to x op y produced by the machine. And we assume that
f l (x op
We also assume that corresponding error bounds are known for the real elementary
functions. For example, we assume there is a bound E sqrt such that
f l (
for some ffl sqrt where
Similarly, we assume relative error bounds, namely
for the other real elementary functions used in Section 4. Each such bound
should be only a small multiple of E; but the bound for sin and cos needs special
attention and will be considered in more detail in the next subsection. (We also assume
that there can be no underflow with sin or cos.)
With these assumptions we can derive bounds on the errors in evaluating expressions
that arise in Section 4. For example, we can conclude that the program expression
x   y   sqrt(x), or f l (xy p
x); has the value
xy
and hence is correct to within a relative error bounded by 2E+E sqrt , if we neglect terms
which are small multiples of E 2 , and provided no exception is returned. (The notation
is convenient - note that different occurrences of ffl are not necessarily the same.)
We use results like this to derive relative error bounds for the approximations we
obtain to the complex functions in Section 4. In general, if f is a complex function
with f r and f i its real and imaginary parts respectively, and f c ; f c
r
and f c
are the
corresponding calculated approximations, then the magnitude of the relative error in
f c is
f
r
f r
r
r
r
are the relative errors in f c
r
and f c
respectively. (We are assuming
here that f r and f i are not zero. If either one is zero, such cases have to be treated
separately.)
Let us denote bounds for
f
r
r
r
This is bounded by max(E r which is a useful bound when E r and E i are not very
different, which turns out to be the case with all our examples in Section 4 except for
clog.
In two of the examples in Section 4 (csqrt and cexp) we examine in more detail
the expression from which this bound is obtained and we are able to determine bounds
which are somewhat smaller than max(E r In the case of clog, we are able to obtain
a bound which is dramatically smaller than max(E r
2.2 Errors in Sine and Cosine Approximations
As stated in the previous subsection, we assume a relative error bound E sincos for the
sine and cosine functions. This means that we assume a bound E sincos such that
f l
This assumption is of no use unless we
restrict the range of values of x for which it is to hold. In Section 5 we will give the
result of measuring this value for on the system we use for testing, but
in the remainder of this subsection we will examine in more detail how such a result
depends on the range for x; the accuracy with which the range reduction is done, the
radix, and the precisions used in calculating the approximations.
We consider in detail only one relatively common case, namely the case where
sin(x) and cos(x) are approximated by first finding x 0 , where x 0 is approximately x mod
-=2 and j x using x 0 in some standard way to determine the final
approximation to sin(x) or cos(x):
Suppose that the approximation to -=2 that is used in the range reduction is P Iby2;
where and we are expecting to choose the bound E -=2 for ffl -=2
to be very much smaller than E: We first determine n so that
will be in relatively high precision. (For a good
discussion of how this might be accomplished, see Payne and Hanek [7].) The value of
would now be rounded, to x 00 say, where x is at least - E;
and then x 00 would be used as the argument for a sine or a cosine approximation valid
over (\Gamma-=4; -=4): (It is common practice to use (possibly simulated) high precision
arithmetic in these computations, in which case ffl 00 would likely be bounded by a small
multiple of E 2 :)
The final approximation to sin(x) is then (\Gamma1) n=2 sin(x 00
even, or (\Gamma1) (n\Gamma1)=2 cos(x 00 where the relative errors
introduced at this stage are the errors committed in approximating the sine and cosine
functions over the interval (\Gamma-=4; -=4) and should be bounded by only small multiples
of E. Similar expressions can be obtained to approximate cos(x):
If n is even, the calculated approximation to sin(x) is therefore
f l
where
We first require P Iby2 to be accurate enough so that
can be neglected, and j A j is bounded by a small multiple of E; so that
we can replace Thus, the relative error in
f l
neglecting only small multiples of In this expression
it can be shown that j u Thus the relative error in
f l (sin(x)) is bounded by
neglecting small multiples of E 2 ; and where max is taken over all machine representable
values of x in the appropriate range. (Of course E sincos can be expected to be somewhat
smaller, but this expression shows clearly what factors need to be controlled to keep
reasonably small.)
This result is applicable when n is even. When n is odd, cot(x) must be replaced by
tan(x) and E sin(\Gamma-=4;-=4) by E cos(\Gamma-=4;-=4) . Similar results can be derived for approximations
to cos(x): All of these results can be combined into one bound which is applicable
in approximations to either sin(x) or cos(x); namely
neglecting small multiples of E 2 ; and where max is taken over all machine representable
values of x in the appropriate range.
We now impose an even more stringent requirement on the accuracy of P Iby2: It
must be chosen so that the error contributions involving the tangent and cotangent
terms in the above expression are bounded by a small multiple of E: The choice will
depend on the radix and precision of the number system, as well as on the appropriate
range of values of x in that number system.
To illustrate, let us consider just one example. If x is a 24-bit binary number, and
we have determined that max (j (n-=2) tan(x) almost
stored as accurately as possible in 66 bits (which is the default
case in Sun systems [6, p. 53]), it turns out that so that
(j (n-=2) tan(x)
Under these circumstances, the above bound for the relative errors in the sine and
cosine functions is which indicates that storing P Iby2 in
bits is accurate enough when j x
If we allow the larger range j x
is almost 586E; which suggests that a considerably more accurate value for P Iby2
should be used in the range reduction. (As an option, the Sun system provides access
to a very high precision approximation.)
We will not consider any further detail at this point. The above is enough to
illustrate what sort of considerations should be taken into account in practice. Some
of our bounds in Section 4 depend on It would be helpful if this quantity were
to be carefully documented for the individual systems on which our algorithms might
be implemented.
2.3 The LOG1P Function
In Subsection 4.4 we need to evaluate a real elementary log function in a situation
where the argument of the log function must itself be evaluated and so is known only
approximately. This introduces a serious problem when the value of that argument is
near say 1+x where x is small. It turns out that we can often calculate x to working
precision, so that a real elementary function log1p(x) that approximates log(1
without explicitly adding 1 to x; is ideally suited to our purposes.
Since log1p is not always available, we present here a simple way to specify a reasonably
accurate log1p function using the log function. (For a more accurate implementa-
7function
possible exceptions (domainerror )
real y
real E / appropriately initialized
real oneoverE / appropriately initialized
return domainerror
return log(x)
return x
else
return log(y)
endif
Fig. 1. An implementation of log1p using log. After looking after three special cases, the program
takes care to obtain an accurate approximation to log(1
tion that does not rely on the log function, see Tang [8].) The program in Figure 1 first
looks after the exceptional case when x - \Gamma1: Then, when x is so large that log1p(x)
can be replaced by log(x); the program makes the replacement; it is sufficient to have
x ? 1=E; and this also ensures that the replacement is made in the one case where it
is necessary to do so, which is when x is so large that x
x is the largest machine representable number and the rounding mode is "round-up".
Then, when j x j is so small that log1p(x) can be replaced by x; the program makes
the replacement; it is sufficient to have but is also necessary to make the
replacement if the arithmetic happens to be truncated. (Otherwise, you could have, on
a binary machine for example, so that y would be 1 \Gamma E=2 instead of 1; and
then log(y) would be \GammaE=2 plus possibly a small multiple of depending on how
accurate the log function is - and the returned value would be 0 or a small multiple of
instead of \GammaE 3 ; in either case the relative error is enormous.)
Having looked after all these special cases, the program assigns the rounded value
x to y, and then recovers the error with x. The (approximate) relative
error in y is then the rounded value of the quotient ((y x)=y, which we denote
here by relerr. Then we have
equal x to high accuracy.
to high accuracy.
We have tested the program only in the standard rounding mode of IEEE [3] arithmetic
(where the two elsif clauses are not needed), but we believe the program is
also valid for other arithmetic systems with reasonable rounding conventions, including
truncation.
2.4 The logb and scalb Functions
We need two functions for manipulating exponents. The first returns an integer value.
It is
(j x j)c ;
if x 6= 0. We will not use this function when
The second function is, for x real and n an integer,
unless this overflows or underflows, in which case the appropriate exception is raised.
Both functions are exact, provided of course that there is no overflow or underflow
with the second one.
2.5 Exceptions
In the examples of Section 4, the only exceptions that can be returned are overflow,
underflow and domainerror.
We have adopted the convention that overflow will be returned whenever either one
of the components of the result overflows, or if both overflow. We do not try to provide
the component values themselves, but only minor modifications would be needed to
provide values of \Sigma1, where appropriate, if that were considered to be desirable. An
appropriate normalized value for a component that did not overflow or underflow could
also be provided.
If only one component underflows, and the value of the other component is so much
larger that setting the underflowed component to zero does not alter the error bound
by more than a small multiple of we set that underflowed component to zero and
do not return underflow. (In fact, this is the situation that occurs with csqrt, clog, csin
and ccos.) If the underflowed component is smaller than the normal component by a
factor of at most E (this is quite arbitrary), we still set the underflowed component to
zero and do not return underflow, but we increase the bound as required - otherwise,
if the criterion for setting the underflowed component to zero is not met, or if both
components underflow, we return underflow. (What is described in the last sentence
can only occur with cexp. The setting of the underflowed component to zero has been
accounted for in our error analysis, and it increases the bound by only about 5% in the
system we use for testing.) If underflow is returned, we do not try to provide component
values themselves, but, once again, only minor changes would be needed to provide for
special values, such as denormalized numbers.
3. PSEUDO-CODE
Our notation for the needed exception handling construct is shown in Figure 2. The
calculations in the enable block will normally produce the required result, but, if an
enable
handle
Fig. 2. The exception handling construct. The enable block would normally produce the required
result, but the handler takes over if an exception occurs in the enable block.
exception occurs during the course of these calculations, control is transferred to the
handle block, or handler, where action is taken to circumvent or otherwise cope with
the exceptional situation.
There are various ways in which such constructs can be implemented, depending
on the programming language used. We will discuss some of these ways briefly in
Sections 5 and 6. But for our present purposes it does not matter how such constructs
are implemented. In particular, it does not matter whether the transfer of control
takes place as soon as the first exception occurs, or whether, as is possible in an IEEE
environment [3], the calculation continues to the end of the enable block where a test
is made to determine whether or not an exception occurred. We do not make use of
any intermediate results that might have been obtained in the enable block.
We assume that any indication that exceptions occurred in the enable block disappears
on leaving the handler and also that exception handling constructs can be nested
within handlers of other such constructs. However, we do not assume that any indication
of the type of exception (overflow, or underflow, etc.) which caused the transfer
of control is available in the handler. This seems appropriate considering the "impre-
ciseness" of interrupts in pipelined machines, or in the presence of any parallelism.
We do not allow exits from exception handling constructs, except for possible returns
from within handlers - whether they are returns of values or returns of exceptions.
Otherwise our pseudo-code is reasonably self-explanatory. It is intended to provide
an easy-to-understand description of algorithms for calculating good approximations to
the complex elementary functions. Implementing the programs can be more convenient
in some languages than in others, a matter to which we return in Sections 5 and 6.
4. PSEUDO-CODE ALGORITHMS
In this section we present algorithms in the form of pseudo-code programs for each
of the six complex elementary functions required in Fortran 77 and Fortran 90. The
error bounds derived in this section are repeated in tabular form in Section 5. The
term precision is used in the programs to denote the number of significant digits in
the machine representations of real numbers - for example, in the IEEE [3] binary
representation, single precision is 24.
4.1 Complex Absolute Value CABS
We first consider the absolute value function j z j; where z y: The program in

Figure

3 for calculating an approximation to this function is based on the formula
The result of any such calculation is a representable real value for any value of z, except
that it can overflow in extreme cases when x and y are both very large.
The main difficulties in developing a program for this function are in dealing with
the possibility of spurious overflows or underflows. Since these will occur only rarely in
most calculations, a good strategy is to attempt to calculate immediately the required
approximation directly on the basis of the formula, since this is both efficient and
reasonably accurate - and it works most of the time - and then, if that approach fails,
the program can take all the time it needs to look after the exceptional cases.
The handler first looks after the special cases where x or y is zero. Otherwise
logb(x) and logb(y) both exist and can be used to determine whether or not x and
y differ greatly in magnitude. If they do differ by enough, the smaller of j x j and
can be neglected. Otherwise they are close enough that they can be scaled so
that the corresponding scaled result can be calculated without any spurious overflows
or underflows. Finally, the scaled value of the result can be unscaled to provide the
required result, but it might happen that this unscaling will itself produce a result
which overflows, in which case overflow must be returned.
For the error analysis we first consider the case when no exception occurs. The
analysis proceeds as follows:
f l
and, since these are of the same sign,
f l
so that
f l
neglected. Then
f l
so that
f l
if small multiples of ffl 2 are neglected.
Thus the relative error in the final result is bounded by E +E sqrt if small multiples
of E 2 are neglected. It is a straightforward matter to check that the error cannot exceed
this bound in any of the paths through the handler, unless of course overflow occurs in
the final unscaling. The bound is therefore valid for all values of the input argument,
as long as overflow is not returned.
function
possible exceptions (overflow )
real x;
integer logbx; logby
integer precision / appropriately initialized
x := z:realpart y := z:imagpart
enable - try the simplest formula - it will work most of the time
answer := sqrt(x   \Lambda2
handle - overflow or underflow has occurred
answer
else
logbx := logb(x) logby := logb(y)
- exponents are so different that one of x and y
- can be ignored
answer := max(abs(x); abs(y))
else - scale - we scale so that abs(x) is near 1
scaledx := scalb(x; \Gammalogbx)
scaledy := scalb(y; \Gammalogbx)
scaledanswer := sqrt(scaledx   \Lambda2
enable - now unscale if possible - this might overflow
answer := scalb(scaledanswer; logbx)
handle - must be overflow in scalb
return overflow
endif
endif
return answer
Fig. 3. This program for the absolute value function first attempts to approximate
directly,
which is efficient and reasonably accurate, and is almost always successful. But, if overflow or underflow
occurs in this attempt, the handler takes over and manages to avoid any spurious overflows
or underflows, usually by scaling; however, the final result can still overflow in very exceptional cases
when the scaling is undone - then overflow is returned.
function
possible exceptions (overflow )
real x;
x := z:realpart y := z:imagpart
enable
answer := sqrt(x   \Lambda2
handle - must be overflow or underflow
maxmag := max(abs(x); abs(y))
enable
handle - must be underflow
enable
answer := maxmag   temp
handle - must be overflow
return overflow
return answer
Fig. 4. An alternative to what is in Figure 3. With this approach the error bound is somewhat larger.
To conclude this section we will mention another, very different approach which
leads to the alternative program shown in Figure 4. It should be noted that the error
bound in this approach is somewhat larger, namely 2:25E
4.2 Complex Square Root CSQRT
We now consider the square root function
re i' ,
we have
r e i'=2
r sin('=2)
r
r
where we will adopt the conventions that the real part is non-negative, and, if y is
zero, the imaginary part is also non-negative but otherwise its sign is the same as the
sign of y: Note that neither component can overflow, so the function cannot overflow.
Furthermore, if one of the components underflows, this component can be set to zero
without altering the error bound by more than a small multiple of according
to the criterion described in Subsection 2.5, the function cannot underflow.
To avoid loss of accuracy due to cancellation, we rewrite
finally obtain
according to the sign of y. In this
form the mathematical specification satisfies the sign conventions stated above and it
also avoids any possibility of loss of accuracy due to cancellation. By writing
pin place of
avoids the possibility of underflow in the latter form.
There are two remaining difficulties which are taken into account in the program
in

Figure

5. One is to avoid any spurious overflows or underflows in evaluating an
approximation to t: This is done in a way analogous to what was done for cabs: if such
an exception occurs, the handler first looks after the cases when x or y is zero; then,
if their exponents differ by enough (depending on which has the larger exponent), the
smaller of j x j and j y j can be ignored, and the corresponding expressions for t are
quite simple; finally, in all other cases, scaling can be used to avoid the exceptions -
but the scaling must be done in terms of an even exponent so that the final unscaling is
in terms of an integer exponent. Because two square root operations are required, the
final approximation for t cannot overflow, unlike the case for the final approximation
in cabs.
function
possible exceptions (none )
real x;
real sqrt2 / appropriately initialized
integer logbx; logby; evennearlogbx
integer precision / appropriately initialized
complex answer
x := z:realpart y := z:imagpart
enable
answer:realpart := t=2
answer:imagpart := y=t
answer:realpart := abs(y)=t
answer:imagpart := sign(y)   t=2
else
answer:realpart := temp
answer:imagpart := sign(y)   temp
endif
handle - overflow or underflow has occurred
answer:realpart := temp
answer:imagpart := sign(y)   temp
answer:realpart := sqrt(x)
answer:imagpart := 0
else
answer:realpart := 0
answer:imagpart := sqrt(\Gammax)
endif
else - determine t
logbx := logb(x) logby := logb(y)
- x can be ignored
- y can be ignored
else - scale and unscale - we scale so that
- abs(x) is near 1, with even exponent
evennearlogbx
scaledx := scalb(x; \Gammaevennearlogbx)
scaledy := scalb(y; \Gammaevennearlogbx)
scaledt := sqrt(2   (sqrt(scaledx   \Lambda2
+scaledy
endif
answer:realpart := t=2
enable
answer:imagpart := y=t
handle - underflow has occurred
answer:imagpart := 0
enable
answer:realpart := abs(y)=t
handle - underflow has occurred
answer:realpart := 0
answer:imagpart := sign(y)   t=2
endif
endif
return answer
Fig. 5. This program for the square root function first attempts to approximate
directly, and, if successful, the final approximation to p
z is obtained according to whether x ?
or overflow or underflow occurs in this attempt, the handler takes over and manages to
avoid any such overflow or underflow. There are two places in the handler where underflow might still
occur - but if it does, the component involved is set to 0.
Once the approximation to t has been obtained, the approximations to the real
and imaginary parts of p
z are easily obtained, depending on whether x is greater
than or less than zero. The only remaining difficulty is that, in the handler only, the
approximation to j y j=t might underflow. If this should occur, the program in Figure 5
replaces the underflowed value with zero. It is safe to do this and still preserve the error
bound obtained below, since in any such case, the other part of the approximation to
z (real part if x ? 0 and imaginary part if x ! 0) is always well above the level
required in Subsection 2.5 for such situations.
For the error analysis, we first consider the case when no exception occurs. We can
begin with the final result for cabs, namely,
f l
Then
f l
so that
f l
on binary machines (another ffl is needed on non-binary machines). Then
r
f l
r
so that
f l
r
on binary machines (another 0:5ffl is needed on non-binary machines). Of course, all of
the above is assuming that small multiples of ffl 2 can be neglected.
On binary machines there is a further division by t in one of the components of
the final result when x 6= 0, which means that the relative error in that component is
The error in each component when assuming the
constant sqrt2 is initialized to be within 1+ ffl of
2: It is also easily seen that the errors
in any path through the exception handler cannot be greater. Thus we can conclude
that, in the notation of Subsection 2.1, one of E r and E i can be while the
other is 2E so that their maximum,
is a bound for the relative error in the approximation to p
neglecting small multiples
of binary machines.
This bound can be tightened somewhat by examining the error formula in more
detail. From Subsection 2.1, we have
f
r
r
r
which, for x ? 0; gives
f
r
r
The right hand side of this inequality reaches its maximum value when is as
large as possible. However, for x ? that
f
which is less than 2E The result is the same for x ! 0: The bound is even
smaller when 0: Thus we conclude that the relative error in the approximation to
z is bounded by q
neglecting small multiples of binary machines.
On non-binary machines there is an extra 0:5ffl in the bound for t: There is also
the same further division by t in one component, as well as a division by 2 in the
other component, when x 6= 0; so that the error bounds for both components of the
result are equal. The final relative error bound for non-binary machines is therefore
neglecting small multiples of
4.3 Complex Exponential CEXP
The complex exponential function of z easily expressed in terms of real
elementary functions of x and y. The relationship is
and this form leads immediately to the enable block in the program of Figure 6.
In the handler, exp(x) is first tested for overflow. If it does, then overflow is returned.
Admittedly, this ignores a narrow "fringe" of possible values of x and y where the
overflow could be avoided, namely those values of x and y such that exp(x) overflows
by such a small amount that multiplication by cos(y) and sin(y) would produce values
that are just below the overflow level.
We have chosen to ignore these fringe values as not being worth the trouble to
detect and reinstate. However, if one wishes to include them, one possibility would be
to develop a special procedure to return separately the significand s and exponent e of
exp(x), and to use the results of this procedure to determine scalb(s   cos(y); e) and
scalb(s   sin(y); e). Then overflow can occur only if one or both components do actually
overflow.
Otherwise the handler must cope with underflow. If exp(x) underflows, or even if
only exp(x)   max(j cos(y) underflow must be returned.
Otherwise both components do not underflow, so only one does, and our agreed upon
function
possible exceptions (overflow , underflow )
real x;
real E / appropriately initialized
complex answer
x := z:realpart y := z:imagpart
enable
expx := exp(x)
answer:realpart := expx   cos(y)
answer:imagpart := expx   sin(y)
handle - overflow or underflow has occurred
enable
expx := exp(x)
handle
return overf low
else
return underf low
endif
cosy := cos(y) siny := sin(y)
enable
handle - both components underflow
return underf low
- only one component underflows
- the underflowed component can be set to zero
answer:realpart := expx   cosy
answer:imagpart := 0
else
answer:realpart := 0
answer:imagpart := expx   siny
endif
else
return underf low
endif
return answer
Fig. 6. In this program for the complex exponential function, the handler first deals with overflow or
underflow of exp(x): Then, if only one component underflows, it determines whether or not that one
component can be safely set to zero.
test, which in this case is whether
is used to determine whether the result after setting the underflowed component to 0
can be accepted. Otherwise underflow is returned.
The error analysis is straightforward. In the absence of any exceptions, the relative
error in the real part is ffl and in the imaginary part it is ffl
which gives an overall bound of E neglecting small multiples of
The bound is somewhat larger in the case where one underflowed part is set to zero.
The relative error in such an underflowed part is 1; and the relative error in the overall
approximation is then bounded by
s
s
neglecting small multiples of
Hence, in all cases, the overall relative error in the approximation of e z is bounded
by q
neglecting small multiples of and provided that overflow or underflow is not re-
turned. This bound turns out to be about 5% more than on the
system we use for testing.
4.4 Complex Natural Logarithm CLOG
The complex natural logarithm of z re i' can be expressed in terms of its
components as follows:
where \Gamma-: (The arg function can be approximated by Fortran's ATAN2(y;x):)
To evaluate the components of this function, it is convenient to first introduce M
and m, the maximum and minimum of j x j and j y j, respectively. Then, if
an exception must be returned, which is designated as domainerror in the program in

Figure

7.
The imaginary part of the logarithm function is easily calculated. The arg function
might underflow, but then the real part is simply log(M) and it can be shown that the
underflowed part can be set to zero without any significant increase in the error bound.
The main difficulties with the complex log function arise in evaluating the real part.
The program first deals with the case when and the real part is again simply
log(M): Then the most serious difficulty occurs when x what is the same
near 1, for then an accurate approximation to the log function
would require an accurate approximation to M 2 +m 1, and this cannot be obtained
directly because of possible serious cancellation.
We will postpone dealing with this difficulty, for the time being, and first consider
only those cases where M is not in the interval (1=2;
2). (This ensures that M 2 +m 2
is not in the interval (1=2; 2).) Then the remaining difficulty is that spurious overflow
or underflow may occur in the evaluation of But we can deal with such
exceptions in a way analogous to what we did for cabs and csqrt. As indicated in the
program, m can be ignored if the exponents of M and m are sufficiently different -
and the resulting real part is simply log(M ). Otherwise, scaling can be used to avoid
overflow or underflow. If scaled by a factor radix \Gammascale , the logarithm
must be corrected by adding scale \Theta log(radix) to the logarithm of the scaled value of
It is important that the scaling be chosen so that there is no possibility of
cancellation in this addition.
Now let us consider the case where
2. In this case 1=4 ! M 2 +m
and we want to calculate an accurate approximation to
The best way to do this is to evaluate the expression (M in double the working
precision and then round the result back to working precision. If the evaluation is done
in the working precision, the error analysis becomes more complicated, as we shall
see, and it turns out that the error in the real part can be enormous, although the
overall error bound turns out to be less than double the quite modest bound we obtain
for the case when doubled precision is used. (If doubled precision is not available, it
could be simulated, but, since the simulation is likely to be extremely slow, it may be
worthwhile avoiding the simulation in those cases where there is no cancellation - i.e.,
when M - 1; or where the cancellation is not very serious - i.e., when
are not very close in magnitude.)
We turn now to an error analysis for the program in Figure 7, and at first we will
assume that the arithmetic is binary. The error bound for the imaginary part is of
course unless the imaginary part underflows. If this underflow does occur, the
imaginary part can safely be set to zero, as we have already indicated, and, since the
real part is simply log(M); with an error bound of E log ; an overall relative error bound
is
Otherwise, for the real part we consider first the case where M - 1=2 or M -
and no exception occurs. Here we have
f l
if we neglect small multiples of
know that j log(M 2 +m 2 log 2; so that the relative error bound for the real part
function CLOG
possible exceptions (domainerror)
real x;
integer scale
complex answer
real sqrt2; logradix / appropriately initialized
integer precision / appropriately initialized
x := z:realpart y := z:imagpart
return domainerror
endif
enable
answer:imagpart := arg(z)
handle - must be underflow of arg
answer:imagpart
return answer
now determine the real part
answer:realpart
enable
answer:realpart := 0:5   log(M
handle - must be overflow or underflow in M
- m can be ignored
answer:realpart
else - scale so that exceptions are avoided, and the two
- terms in answer:realpart are of the same sign
must have been overflow
scale
else - must have been underflow
scale := logb(M
endif
scaledM := scalb(M; \Gammascale) scaledm := scalb(m; \Gammascale)
scaledr := scaledM   \Lambda2
answer:realpart := scale   logradix
endif
else
enable - use doubled precision if possible in evaluating the argument of log1p
answer:realpart
handle - must be underflow in m   \Lambda2
answer:realpart
endif
return answer
Fig. 7. This program for the complex logarithm function first looks after three special cases - when
the argument is zero, when the imaginary part underflows, and when 0: The real part is then
calculated in a way depending on whether max(j x outside the interval (1=2;
or not; in
the first case scaling may be needed to cope with spurious overflow or underflow; in the second case the
accuracy of the final result can be very sensitive to how accurately the argument of log1p is calculated.
is 2E= log 2 +E log ; which is bounded by 2:886E +E log ; if we neglect small multiples of
When an exception does occur this bound is obviously still valid when no scaling
is required. When scaling is required this bound is valid for the second term in the
expression for the real part, while the first term, scale   logradix; is bounded by 2E; if
logradix is stored as accurately as possible; since the two terms are of the same sign,
the error in their sum is bounded by 2:886E+E log : Thus, for all cases when M - 1=2 or
2; an overall relative error bound is max(2:886E neglecting small
multiples of
We now consider the case where
2: If an exception occurs in this case,
the real part is log(M) and the error bound is E log : Otherwise we need to determine
how errors in the evaluation of affect log1p of this expression. If
this expression can be evaluated to within a factor of 1 then we have
f l
log1p
if we neglect small multiples of j 2 and fflj:
We know that and we can show that
is an error bound for the real part, where H is a
bound for j j j:
If doubled precision is used, are evaluated exactly, and so are
The sum of these two expressions will suffer a rounding error
of at most ffl 2 ; so that the final rounding to working precision presents log1p with an
argument in which ffl: The final relative error in the real part when
pis therefore bounded by 2:165E neglect small multiples of
The overall relative error, for all M; is therefore bounded by
if we neglect small multiples of E 2 and provided of course that domainerror is not
returned. This completes the error analysis of the program in Figure 7 when the
argument of log1p near the end of the program is calculated to within a factor of
which is the case when doubled precision is used to evaluate that argument.
(The 2:886 and 2:165 in this bound must be replaced by 3:886 and 3:165; respectively,
for non-binary machines because of the multiplication by 0:5 in the expressions being
evaluated.)
When doubled precision is not used, the situation is quite a bit more complicated
since serious cancellation may occur in the evaluation of
there is obviously no cancellation, and the argument for log1p is simply
be accurate to within a factor of 1 as was the case with doubled
precision. If M ? 1; there is also no cancellation, but the argument is now accurate
only to within a factor of 1
leads us to a bound which replaces the 2:165 above with 6:495: If M ! 1; there can be
cancellation, but, if 4m 2 the cancellation is not very serious: it can be
shown that in this case the argument is accurate to within a factor of 1
and the bound is therefore the same as it was for
This leaves us with a situation where more serious cancellation can take place. For
this, we derive the following:
f l
log1p
hi
if we neglect small multiples of From the second term in this expression we see that
the relative error in the real part could be enormous (and this is confirmed by our
tests in the next section). But we continue towards finding a relative error bound in
the overall result. From the above we see that the absolute error in the real part is
bounded by
and, since j M in turn is bounded by
This expression is a bound for E r which can be substituted into the overall
relative error bound developed in Subsection 2.1, namely
f
r
r
r
We can proceed as follows:
f
r
so that we are considering the first part of the bound for E r but the
second part relative to j f r j: This will lead to an overall error bound that is acceptably
small.
The second term is obviously For a bound on the first term, we consider two
cases, one where j y and the other where j y
The imaginary part is arg(z) and, in the first case, j arg(z) j -
so that the
contribution of the first part to the relative error is bounded by
since
In the second case, j arg(z) so the contribution of the first part
to the relative error is
which also turns out to be - 4:457E:
Thus, we end up with an overall relative error bound of 4:457E +E log1p
the case when serious cancellation can take place.
Collecting together all the results we have obtained for the program in Figure 7
when only working precision is used throughout, we have established the following
overall relative error bound for the clog function:
if we neglect small multiples of and provided of course that domainerror is not
returned. (Here each of 2:886; 6:495 and 4:457 must be increased by 1 for non-binary
machines, because of the multiplication by 0:5 in the expressions being evaluated.)
4.5 Complex Sine CSIN
The sine function of z = x+iy can be represented in terms of real elementary functions
as follows:
and the program in Figure 8 is based on this formula.
If overflow occurs in the evaluation of the real and imaginary parts of this function,
cosh(y) or sinh(y) is too large in magnitude (probably both), and the handler returns
overflow. As in the case of cexp some "fringe" values of z (for which the real and
function
possible exceptions (overflow )
real x;
complex answer
x := z:realpart y := z:imagpart
enable
answer:realpart := sin(x)   cosh(y)
answer:imagpart := cos(x)   sinh(y)
handle
enable
coshy := cosh(y)
sinhy := sinh(y)
handle - must be overflow
return overflow
answer:realpart := sin(x)   coshy
enable
answer:imagpart := cos(x)   sinhy
handle - must be underflow
answer:imagpart := 0
return answer
Fig. 8. This program for the complex sine function is straightforward, although it should be acknowledged
that overflow can be returned for some "fringe" values of z whose corresponding real and
imaginary parts are actually slightly below the overflow threshold.
imaginary parts of sin(z) are slightly below the overflow threshold, even though at least
one of cosh(y) and sinh(y) alone do overflow) are neglected here. But, as with cexp, an
auxiliary procedure that computes the exponent and fraction part of exp(y) separately
could be used to avoid this situation, since sinh(y) and cosh(y) both effectively equal
exp(y)=2 in magnitude when these functions overflow.
If underflow occurs it can only occur in the multiplication associated with the imaginary
part. But when that happens the real part is so much larger that the underflowed
part can be set to zero without any significant increase in the error bound.
An upper bound for the overall relative error in this function is
neglecting small multiples of of course, that overflow is not returned.
4.6 Complex Cosine CCOS
The cosine function of z can be represented in terms of real elementary
functions as follows:
function
possible exceptions (overflow )
real x;
complex answer
x := z:realpart y := z:imagpart
enable
answer:realpart := cos(x)   cosh(y)
answer:imagpart := \Gammasin(x)   sinh(y)
handle
enable
coshy := cosh(y)
sinhy := sinh(y)
handle - must be overflow
return overflow
answer:realpart := cos(x)   coshy
enable
answer:imagpart := \Gammasin(x)   sinhy
handle - must be underflow
answer:imagpart := 0
return answer
Fig. 9. This program for the complex cosine function is straightforward, although it should be
acknowledged that overflow can be returned for some "fringe" values of z whose corresponding real
and imaginary parts are slightly below the overflow threshold.
and the program in Figure 9 is based on this formula.
If overflow occurs cosh(y) or sinh(y) is too large in magnitude (probably both), and,
as with csin in the previous section, the handler returns overflow . But also as with
csin, some "fringe" values of z are neglected here, although they could be included with
the help of an auxiliary procedure that returns the exponent and fraction part of exp(y)
separately.
If underflow occurs it can only occur in the multiplication associated with the imaginary
part. But when that happens the real part is so much larger that the underflowed
part can be set to zero without any significant increase in the error bound.
An upper bound for the overall relative error in this function is
neglecting small multiples of of course, that overflow is not returned.
5. SPECIAL IMPLEMENTATIONS AND TESTING
We have implemented the algorithms presented in Section 4 on a Sun 4/40 in Fortran 77
(compiler version 1.4), in order to test their correctness, especially the correctness of
their error bounds. These implementations are special in the sense that they are as

Table

I. Observed error bounds for the single precision real elementary functions in the Sun library
(version 1.4), in units of E; the relative error bound for single precision real arithmetic. The result
given for E sincos is for
1:000 1:152 1:102 2:326 1:382 1:000 1:000
close as possible to the pseudo-code descriptions of Section 4, and in particular have not
been modified in any way to improve their efficiency. (Efficiency and other production
implementation issues will be discussed in the next section.)
Except for a portion of one version of clog, which was an alternative suggested in
Subsection 4.4, the floating point operations in these implementations are in single
precision. (Some care had to be taken to make sure this was the case! We examined
the generated assembly language instructions to ensure that no extended precision was
used in any intermediate calculations.) This enabled us to use the corresponding double
precision results from the Sun system as the "true" results for test purposes. We assume
that these "true" results are correct to single precision accuracy.
The exception handling construct was implemented by allowing the enable block to
be executed and then testing for which exception flags had been raised, as is possible
in an IEEE environment.
It would have been natural, with our interpretation, to use the "ieee handler"
trap-handling facility provided by the Sun system [6, p. 67] but it turned out to be
both inefficient and somewhat difficult to use. Testing the floating point exception flags
using the "ieee flags" subroutine [6, p. 64] is also inefficient, so we instead accessed
the flags by using the math library routine "swapEX ".
To test the correctness of the error bounds for the examples in Section 4, we must
first determine what those bounds turn out to be for the system we are using, and this
requires the determination of E etc. It is convenient to present these in units
of E; and the results are given in Table I. These relative error bounds were determined
by examining all relevant single precision arguments, except that in the case of E sincos
we considered only values of the argument which were less than 10 6 in absolute value,
and in the case of E arg we determined the bound by adding E to E atan ; where the extra
allows for the additional error induced by using f l (y=x) as the argument of atan in
place of y=x: The boundaries of the arg function, where at least one of x or y is 0 or
were also considered. (The bounds for log1p and sinhcosh are 1:000E because
the single precision versions of these functions use correctly rounded results from their
double precision implementations. The bound for our version of log1p in Figure 1,
which uses only single precision, is 2:198E:)
The results in Table I are used to determine the theoretical bounds in column 3 of

Table

II. The observed bounds in column 4 of Table II were obtained by comparing the
results from our implementations with the "true" results provided by the Sun system's
double precision functions for a large number of mostly random input arguments. The
IEEE single precision arguments at which these observed bounds occurred are given in
hexadecimal form in column 5.
The random arguments were constructed from random real parts and random imag-
Table

II. A comparison of theoretical and observed relative error bounds for Sun Fortran (version 1.4)
implementations of the complex elementary function programs in Figures 3-9. Columns 3 and 4 are in
units of E; the relative error bound for real arithmetic. SP and PDP stand for IEEE single precision
and partial IEEE double precision, respectively. (For cexp, real for csin and ccos,
imag
Function Theoretical Bound \Gamma! Based on

Table

I
Observed
Bound
at this
argument
cabs2 2:25E +E sqrt 3.250 2.495 e0723701
csqrt
sqrt
cexp
clog
inary parts, by generating random exponents and random significands within appropriate
ranges. Ten thousand such random arguments were generated for each of the
4 semi-axes, and 100; 000; 000 were generated for each of the 4 quadrants. The origin
was also tested. The regions were restricted where necessary so that overflow and
underflow would be avoided most of the time, while at the same time ensuring good
coverage of the proper domain of each function. Many special cases were also tested,
including many near the boundaries of the regions that separated points which would
probably lead to exceptions being returned and points which would probably not lead
to exceptions being returned. The most important special cases in terms of trying to
observe large errors were those arguments where real and imaginary parts were chosen
to maximize the errors in the relevant real elementary functions; in fact, these special
cases produced most of the observed maximums - so much for random testing!
As can be seen from Table II, all non-exceptional results were in error by no more
than what would be expected on the basis of the theoretical error bounds. In fact,
the theoretical bounds are not much larger than the observed bounds. The sort of
discrepancy shown in the table is not surprising, considering the kinds of reasoning
used in determining theoretical error bounds, especially in the case of the two clog
functions.
(It happens that most of the relatively large discrepancy in the case of cabs2 can be
explained. The IEEE arithmetic in our system satisfies more than what we assumed
in Subsection 2.1. It is accurate to within half a unit in the last place. This can make
a significant difference to the error analysis in cabs2 because of the special form of
an expression such It can be argued that the contribution from such an
expression to the final relative error bound in cabs2 is a maximum when its value is just
greater than 1:5: We have followed this argument through to obtain finally a relative
error bound of only 2:650 for cabs2, in place of 3:250:)
Special mention should be made of the results for clog. The version that makes use
of double precision in the computation of the argument for log1p is not very much more
accurate in terms of its observed overall relative error bound than the version that uses
only single precision. However the real part of the former is much more accurate than
the real part of the latter - the observed bound is 3:604E in the former but approximately
in the latter! - this enormous error occurred at 3f7ffffc + i3a3504f3:
(The observed error bounds in their imaginary parts were equal, only 2:292E:)
Of course we do not claim that these experimental results actually prove the correctness
of the theoretical error bounds, but we do believe the evidence is very convincing.
However, correctness of the programs involves more than correctness of the error bounds
for values of the argument which do not lead to exceptions being returned. It must also
be true that exceptions are returned only when it is reasonable to do so. (Otherwise
a program could be considered correct if it always returned an exception, no matter
what the input argument!)
According to Subsection 2.5, overflow is to be returned when either component
overflows. If one of our function programs does return overflow, for a particular value
of the input argument, our test program considers this to be a correct return if at least
one of the components of the "true" result is within the relative error bound of a true
overflow, specifically if
where HUGE is the largest machine representable number, and ulpup is a unit in the
last place of HUGE in the direction of 1:
The test for an underflow return to be correct is that either
(1) both components of the "true" result are within the relative error
bound of a true underflow, specifically if
or
(2) one component underflows, but is nevertheless within an error bound
of being greater than a value that is within an error bound of the other
component, specifically if
min (j f r
and
min (j f r
where T INY is the smallest positive machine representable number, and ulpdown is a
unit in the last place of T INY in the direction of 0:
These criteria had to be modified slightly, in an obvious way, to allow for the fact
that "fringe" areas were neglected in the algorithm for cexp, csin and ccos.
In each test case for which overflow or underflow was returned, our test program
determined that the appropriate criterion was satisfied.
The exceptional return of domainerror in clog is a special case which was easy to
check separately.
6. PRODUCTION IMPLEMENTATIONS
The exception handling construct in Ada is reasonably close to what we have used in
this paper, except for the unfortunate fact that Ada does not recognize underflow as an
exception. Our construct is simpler than Ada's in that we do not expect the handler to
be told, in effect, what exceptions have occurred. A proposal somewhat like Ada's, but
which did recognize underflow, was proposed for Fortran 90 - but was rejected in favor
of having no exception handling facility at all! Our construct could be implemented in
PL/1, but not very elegantly.
Extensions to existing languages often provide facilities for implementing our con-
struct, for example by providing access to trapping and exception flags in IEEE environments
- as we have already indicated in the preceding section.
In the absence of exception handling facilities, pretesting can be used. For example,
in the case of cabs, j x j and j y j can be tested in advance to make sure that no overflow
or underflow will occur in the evaluation of x if that is the case, evaluate
the expression
carry out the calculations in the handler, using
another pretest to determine whether or not the final unscaling will cause overflow. In
such an implementation, the original pretesting will determine if conditions are satisfied
which are sufficient to ensure that no spurious exceptions will occur. They may not be
necessary, and in such cases this will cause the program to execute the handler more
often than is necessary.
Of course, if higher precision is available, it may be possible to avoid some pretest-
ing, as long as the exponent range in that higher precision is sufficiently broader. For
example, again for cabs, the entire calculation of
can be done in higher preci-
sion, and a test needs to be performed only when the result of this calculation is about
to be coerced to the original precision. But this idea is of little use if the working
precision is already the highest available.
Apart from trying to implement exception handling in a more efficient way, our
implementations can be made slightly more efficient, mostly by writing the programs
so that some of the intermediate "store" operations are avoided.
In addition, the floating point status register would have to be saved on entry to the
function subprograms, and later restored and updated to include the exception flag, if
any, returned by the function.
7. CONCLUDING REMARKS
We have presented algorithms for reliable and accurate evaluations of the complex
elementary functions required in Fortran 77 and Fortran 90. It was convenient to
describe these algorithms with the help of an exception handling construct.
Implementations in Fortran for Sun systems have been tested extensively. The
observed error bounds were between 64% and 98% of the theoretical bounds - which
is not only convincing evidence of the correctness of the theoretical bounds, but also
indicates that the theoretical bounds are quite tight. (It was interesting to discover that
choosing arguments near where we thought the largest error might occur usually led
to observed bounds that were larger than those found even by very extensive random
testing - in one case by over 13%!) In the tests it was also found that exceptions were
returned when, and only when, it was reasonable to do so.

ACKNOWLEDGMENTS

Much of this work was inspired by discussions with members of the Ada Numerics
Working Group under the chairmanship of Gil Myers. Jim Cody was particularly
helpful during early stages of our investigation. We also wish to thank the referees for
helpful suggestions.



--R

American National Standard Programming Language FORTRAN.
Implementing complex elementary functions using exception handling.
IEEE Standard for Binary Floating-Point Arithmetic
Branch cuts for complex elementary functions
Part
Radian reduction for trigonometric functions.

--TR
Table-driven implementation of the logarithm function in IEEE floating-point arithmetic

--CTR
Technical report for floating-point exception handling, ACM SIGPLAN Fortran Forum, v.15 n.3, p.1-28, Dec. 1996
T. E. Hull , Thomas F. Fairgrieve , Ping Tak Peter Tang, Implementing the complex arcsine and arccosine functions using exception handling, ACM Transactions on Mathematical Software (TOMS), v.23 n.3, p.299-335, Sept. 1997
T. E. Hull , Thomas F. Fairgrieve , Ping Tak Peter Tang, Implementing the complex arcsine and arccosine functions using exception handling, ACM Transactions on Mathematical Software (TOMS), v.23 n.3, p.299-335, Sept. 1997
David M. Smith, Algorithm 786: multiple-precision complex arithmetic and functions, ACM Transactions on Mathematical Software (TOMS), v.24 n.4, p.359-367, Dec. 1998
David Bindel , James Demmel , William Kahan , Osni Marques, On computing givens rotations reliably and efficiently, ACM Transactions on Mathematical Software (TOMS), v.28 n.2, p.206-238, June 2002
Milo D. Ercegovac , Jean-Michel Muller, Complex Square Root with Operand Prescaling, Journal of VLSI Signal Processing Systems, v.49 n.1, p.19-30, October   2007
John R. Hauser, Handling floating-point exceptions in numeric programs, ACM Transactions on Programming Languages and Systems (TOPLAS), v.18 n.2, p.139-174, March 1996
implementation;complex elementary functions