posting to sci.math.num-analysis on trigonometric functions

[David G. Hough on validgh]

Should computer arithmetic prevent subtraction of two close but unequal
numbers?   After all, when "cancellation" occurs in this case, any error
inherent in the original operands will be magnified, possibly leaving no
significant figures in the result.   After all, most floating-point operands
are inherently inexact.

So some possible approaches to dealing with this threat are 

1) set the result to zero, and on C implementations set errno to something,
and continue

2) abort the computation

3) set a flag that can be tested later in the program, but by function calls
that aren't portable across different machines

4) use computational methods such as interval arithmetic to maintain error
bounds on all the variables in the computation, including the results of
the cancelling subtraction

5) use analytical methods to determine error bounds for all the variables

===============================================================================

Which of these approaches make sense to you?   It's instructive to consider
this simple problem (z = x-y) before deciding what to do about (z=sin(x)).

#1 represents the System V Validation Suite.   Conforming implementations
have an interesting anomaly in their trigonometric functions near an 
arbitrary internal boundary.   Trigonometric identities fail in a big way.
Many people, however, find zero z = x-y confusing when x != y.

#2 represents a very conservative approach: abort on any exception.
Not many floating-point programs pass this rigorous standard.   In fact,
by aborting on the rounding error that precedes the cancellation, even more
dangerous results can be avoided.

#3 represents an approach somewhat like the current state of IEEE arithmetic.
The reason IEEE arithmetic has an inexact exception flag, rather than a
loss-of-significance exception flag, is that the inexactness is what's
important.   No error is created by the cancellation step.   The portability
problem is being worked on by NCEG.

#4 and #5 are the right answer for dealing with rounding error, whether 
or not revealed by cancellation.   If you don't care about the error, then
you don't care about the answer, unless you have some other method of
verifying it.

===============================================================================

With trigonometric functions, ultimately the problem is that the natural
transcendental functions are... transcendental.    They don't mesh well with
ordinary arithmetic.   You can change the units to be degrees and 10, say,
rather than pi and e, but then you get into  problems with the complex
numbers: exp (i * pi) = -1 represents a fundamental relationship that can't
be defined away conveniently.

So one higher-level solution, redefining the base of logarithms and the
units with which we measure angles, is of some limited value in some cases.
It might have been better, overall, if the first Fortran had defined
	sinpi(x) = sin(pi * x)
as well as sin(x), since many mathematical formulas have explicit pi factors
in the trig arguments.   Then argument reduction of sinpi, at least, could
have been exact.

But there is a growing consensus that there's no good reason to compromise
sin(x).  The computed value should be within a small rounding error
of the correct value of sin(x), for all x, even large ones, even potentially
inexact ones, just as almost all computers now, except Crays, accept their
obligation to compute x-y to within a small rounding error of the correct
value of x-y, even though x and y might be inexact, and there might be 
no significant digits in the result.   There's no way of knowing locally
whether the potential
magnification of error by cancellation is of global significance.
There might be no error in the operands; it may be that what's important
is preserving relationships like trigonometric identities, so that preserving
them approximately is more important that grossly violating them near a
boundary.

Which is all a long-winded way of saying that the subject of computer
approximations to elementary transcendental functions is a subject that 
deserves more careful treatment than a few usenet postings.   
A good place to start is Peter Tang's papers in recent issues of ACM TOMS.