no substitute for getting the right answer
David Hough
sun!Eng!David.Hough
Tue Oct 25 10:42:33 PDT 1994
Elementary transcendental functions are used in many unpredictable ways,
which is why the implementer of the standard library should do as good a
job as possible. Persons creating specific applications for which they
can prove that they can use faster, less accurate transcendental functions
can code their own far easier than persons creating specific applications
that require full precision and correct exception handling can code their
own when the library is defective. fdlibm incorporates "infinite-precision"
trigonometric argument reduction.
The whole issue of large arguments to trig functions boils down to just this:
should there be an exception raised on loss of significance, or cancellation,
which is the situation when an EXACT subtraction magnifies any relative error
that may exist in the arguments. Whether this is a problem or not depends
on whether there is any relative error in the arguments, and what is done
with the result, neither of which is known to the floating-point subtraction
hardware. Note that a conscientious implementation of significance loss
detection in software would have to detect trig cases of rather modest size;
as discussed in my article in the March 1980 IEEE Computer.
In fact skillful exploitation of cancellation is part of the
mechanism for building higher-precision floating point out of lower precision;
dealing with an exception in that case would get in the way, which is probably
why significance exceptions have seldom been implemented in hardware
and when implemented, as on IBM 360 I think,
have invariably been disabled in software. The "damage" of cancellation is
always done earlier,
when error contaminated the operands; hence IEEE 754 specifies
an exception on inexact rather than loss of significance. That almost
all floating-point programs generate inexact exceptions reinforces the
observation that if you don't know the error (or its bound) then you don't
know the answer, regardless of what the computer prints; that argues for
manual or automatic error bounding (interval arithmetic).
Kulisch's papers usually incorporate
drastic examples of cancellation of inexact intermediate results
to demonstrate how even increasing precision may fail to yield any correct
significant figures in conventional floating-point arithmetic. The lack
of significance in the results is immediately demonstrated by any kind of
error analysis.
The following recent USENET exchange illustrates these points; one release
of HP libraries evidently didn't even enforce the identity sin*sin+cos*cos = 1,
with unfortunate results little alleviated by an error message:
=============================================================================
Article: 5643 of comp.lang.fortran
From: oliveriaadown.engin.umich.edu (Roque Donizete de Oliveira)
Newsgroups: comp.lang.fortran,comp.sys.hp.hpux
Subject: Re: how to detect TOTAL LOSS OF PRECISION with HP's f77 compiler
Date: 24 Oct 1994 07:51:12 GMT
Organization: University of Michigan Engineering, Ann Arbor
To followup on my original question, after 2 days of painful debugging.
I have a complex nonlinear integral equation I want to solve for C
0 = f(C) + integral_{Pi/4}^{Pi/2} { g(Z) dtheta }
where the argument of g(Z) looks like
Z = Z(C,theta) = C / DCOS(theta)
At the place where I got the exception
C = (-.2721096508319146D+01 , -.2400237696120182D-05 )
This value of C is pretty close to exact solution, so the nonlinear
solver is doing ok.
To avoid the singularity at theta=Pi/2 we have used an upper limit
a tiny bit (10^-13) smaller than Pi/2. This is ok (we noticed the
integral converges, and the main contribution comes from near
theta=Pi/2, and the integrand, as a whole, is bound at theta=Pi/2).
The problem arises because near theta=Pi/2 (these points are picked
by the numerical integrator) Z becomes very large, like
Z = ( -.2721605324800001D+14 , -.2400686515371547D+08 )
However, the g(Z) function, looks like (in this range of Z):
g(Z) = some_other_terms + CDEXP(- Z**2 * 0.5)
and therefore the argument of the complex exponential became very
large (like (-3.703567771986978E+26 ,-6.533721203410762E+20)).
Since a exponential of a complex argument can be written as
CDEXP(x + i*y) = DEXP(x) * ( DCOS(y) + i*DSIN(y) )
we are left with trigonometric functions of very large arguments.
Unfortunately, HP's trigonometric functions break for arguments
larger than about 10**8. They return 0 for the function value
and the "warning" "TOTAL LOSS OF PRECISION".
Incidentally, IBM's and Sun's fortran libraries happily accept
trigonometric functions for arguments as large as the
largest double precision number (about 10**308), even though
the result is correct only for (|y| < 2**50 Pi =3.53 10**15
and x <= 174.673), as stated in the IBM documentation.
I think what I'll end up doing is to modify g(Z) to avoid
computing CDEXP(- Z**2 * 0.5) if the real part of the
argument is very large and negative. Or, if someone has
an implemtation of CDEXP that handles these cases quietly
I would appreciate to hear about it.
And to those who send me "me too" replies on this subject,
I suggest using an error handler to catch these errors.
For example, create a C program called traps.c containing
#include <math.h>
#include <stdio.h>
#include <signal.h>
int matherr(x)
register struct exception *x;
{
switch (x->type) {
case DOMAIN:
/* change sqrt to return sqrt(-arg1), not 0 */
if (!strcmp(x->name, "sqrt")) {
/*x->retval = sqrt(-x->arg1);*/
return (0); /* print message and set errno */
}
else if (!strcmp(x->name, "sqrtf")) {
/*x->retval = sqrtf(-x->arg1);*/
return (0); /* print message and set errno */
}
case SING:
/* all other domain or sing errors, print message and abort */
fprintf(stderr, "domain error in %s\n", x->name);
abort( );
case PLOSS:
/* print detailed error message */
fprintf(stderr, "loss of significance in %s(%g) = %g\n", x->name, x->arg1, x->retval);
return (1); /* take no other action */
case TLOSS:
/* print detailed error message */
fprintf(stderr, "total loss of precision in %s(%g) = %g\n",x->name, x->arg1, x->retval);
/*return (1); /* take no other action */
abort( );
}
return (0); /* all other errors, execute default procedure */
}
void traps()
{
/* enable trapping of floating point exceptions */
signal(SIGFPE,matherr);
}
If you add a line
CALL TRAPS()
at the beginning of your fortran program, and link it with traps.o
then you'll catch these errors (it will crash and create a core file).
You can then use xdb and the trace command to find out where the
problem happened.
Roque
=============================================================================
Article: 5650 of comp.lang.fortran
From: shenkinastill3.chem.columbia.edu (Peter Shenkin)
Newsgroups: comp.lang.fortran,comp.sys.hp.hpux
Subject: Re: how to detect TOTAL LOSS OF PRECISION with HP's f77 compiler
Date: 24 Oct 1994 16:13:45 GMT
Organization: MacroModel Development Group, Chemistry, Columbia U., NY, NY
In article <38fp1g$8tasrvr1.engin.umich.edu>,
Roque Donizete de Oliveira <oliveriaadown.engin.umich.edu> wrote:
>Unfortunately, HP's trigonometric functions break for arguments
>larger than about 10**8. They return 0 for the function value
>and the "warning" "TOTAL LOSS OF PRECISION".
Well, now the nature of the problem is becoming clear. If the
value of the argument is so large that a change by one in the last
significant digit amounts to more than 360 degrees, then it's
in some sense stupid to try to determine a cosine, since it
is likely that numerical rounding has already removed all information
necessary to determine what the cosine is.
There has been discussion in many newsgroups about what the
"correct" behavior is; one view is that the sytem should point
out the problem, as HP does. Another view is that the trig function
should just return the value that it would have if its argument
had its given value, assumed known to high precision.
In either case you need to worry whether what you're doing here
is sensible. The fact that other platforms don't give you this
error might simply mean that those platforms are keeping you
unaware of a real numerical problem you need to deal with.
-P.
--
*********** World music: What bluegrass is to a Bulgarian. **********
*Peter S. Shenkin, Box 768 Havemeyer Hall, Chemistry, Columbia Univ.,*
* New York, NY 10027; shenkinacolumbia.edu; (212) 854-5143 *
************ Alphabetical order makes strange bedfellows. ************
=============================================================================
Article: 5661 of comp.lang.fortran
From: oliveriaadown.engin.umich.edu (Roque Donizete de Oliveira)
Newsgroups: comp.lang.fortran,comp.sys.hp.hpux
Subject: Re: how to detect TOTAL LOSS OF PRECISION with HP's f77 compiler
Date: 25 Oct 1994 03:06:22 GMT
Organization: University of Michigan Engineering, Ann Arbor
C This program illustrates the "very small" range of arguments
C of trigonometric functions that are handled by HP fortran
C libraries under HPUX 9.01. For arguments above about 10**8 the
C trigonometric functions give the error "TOTAL LOSS OF PRECISION"
C and return 0. In HPUX 9.05 it seems this limit was raised to
C about 10**15 and it doesn't give you a warning anymore when
C the argument is too large (the returned results are bogus though).
C I don't think I like the behaviour in HPUX 9.05 (I do appreciate
C the increased range of arguments though).
implicit none
double precision a,c,s,t
1 write(*,*) 'enter a'
read(*,*) a
write(6,*) 'calling DCOS'
c = DCOS(a)
write(6,*) 'calling DSIN'
s = DSIN(a)
write(6,*) 'calling DTAN'
t = DTAN(a)
write(*,*) 'a=',a,' Tan=',t,' Sin=',s,' Cos=',c
go to 1
stop
end
C down% uname -a
C HP-UX down A.09.01 A 9000/710 2007005523 two-user license
C f77 junk.f
C a.out
C enter a
C 1.d9
C *** MATH LIBRARY ERROR 32: DCOS(X): TOTAL LOSS OF PRECISION
C *** MATH LIBRARY ERROR 31: DSIN(X): TOTAL LOSS OF PRECISION
C *** MATH LIBRARY ERROR 30: DTAN(X): TOTAL LOSS OF PRECISION
C a= 1000000000.0 Tan= .0 Sin= .0 Cos= .0
C
C biscayne% uname -a
C HP-UX biscayne A.09.05 E 9000/710 2009185315 8-user license
C
C a.out
C enter a
C 1.d9 <--- correct result
C calling DCOS
C calling DSIN
C calling DTAN
C a= 1000000000.0 Tan= .6514522021451413 Sin= .5458434494486995 Cos=
C .8378871813639024
C enter a
C 3.d15 <--- sensible result but wrong
C calling DCOS
C calling DSIN
C calling DTAN
C a= 3000000000000000. Tan= .1266142734353373 Sin= .1256114272401405
C Cos= .9920795170403197
C enter a
C 1.d17 <--- totally bogus result
C calling DCOS
C calling DSIN
C calling DTAN
C a= 1.000000000000000E+17 Tan= -1.31357109754878 Sin= 1548.316605399928
C Cos= -1178.70788135427
C enter a
C I hit Control-C at this point
C
C Some comments:
C Even though the results at a=3.d15 may seem reasonable (you don't
C get magnitude of sin and cos greater than 1), according to
C Mathematica they are wrong. Further tests revealed that for
C arguments up to 10**15 the results agree pretty well with Mathematica.
C I'm not saying a fortran implementation of the trig functions
C should match the results of a multiprecision system.
C In my problem I was not evaluating SIN and COS of large arguments
C directly, I was using CDEXP(x + I*y). Even though under HPUX 9.01
C I would get tons of warnings due to COS(y) and SIN(y), and the
C the returned result set to 0, the real exponential DEXP(x) turned
C out to be practically 0 (since x was very large and negative) and
C therefore the "TOTAL LOSS OF PRECISION" problem didn't affect the
C final result.
C It seems the behaviour under HPUX 9.05 is more in line with that
C of other vendors (the IBM RS/6000 performed a bit better in this
C case, returning reliable results up to about 10**20).
C
C I prefer the warnings (so the user knows something is happening)
C but I also prefer a compile-time option to turn off the warnings
C (once the user has indeed verified the warnings are harmless).
C Quietly returning 0 for very large arguments (like a=1.d30) may
C just cause problems later in the user program.
C
More information about the Numeric-interest
mailing list