Date: Thu, 06 Mar 2003 11:05:45 -0600
From: Eric Rudd <rudd AT cyberoptics DOT com>
Subject: Re: elefunt results
In-reply-to: <200302281831.h1SIVHT22372@speedy.ludd.luth.se>
To: djgpp-workers AT delorie DOT com
Message-id: <3E677FE9.4070703@cyberoptics.com>
Organization: CyberOptics
MIME-version: 1.0
Content-type: text/plain; charset=ISO-8859-1; format=flowed
Content-transfer-encoding: 7bit
X-Accept-Language: en,pdf
User-Agent: Mozilla/5.0 (Windows; U; Win98; en-US; rv:1.2.1) Gecko/20021130
References: <200302281831 DOT h1SIVHT22372 AT speedy DOT ludd DOT luth DOT se>
Reply-To: djgpp-workers AT delorie DOT com
Errors-To: nobody AT delorie DOT com
X-Mailing-List: djgpp-workers AT delorie DOT com
X-Unsubscribes-To: listserv AT delorie DOT com
Precedence: bulk

ams AT ludd DOT luth DOT se wrote:

>Alas there seems to be some failures. I'm not planning to try to
>understand or correct them, but I'm reporting it in case somebody else
>likes to examine it.
>
I will comment on these results, since I wrote the libc elementary math 
functions.

>How current are these test anyway? All capitals is a little weird;
>makes me think of the 80ties.
>
As I understand, these tests *are* old; they were translated from 
FORTRAN.  However, mathematics is not invalid just because it is old.

These ELEFUNT errors occurred because of conscientious design choices on 
my part, rather than unexpected bugs.  Informed and responsible 
numerical analysts don't always agree on the handling of exceptional 
cases, and these are some of the cases where I disagreed with the 
philosophy of ELEFUNT.  Below are my explanations of why these error 
messages are triggered:

>ttan.djl, error message missing:
>TAN WILL BE CALLED WITH THE ARGUMENT     9.0072e+15
> THIS SHOULD TRIGGER AN ERROR MESSAGE
>
>
> TAN RETURNED THE VALUE     1.6063e+00
>
>tsin.djl, error message missing:
> SIN WILL BE CALLED WITH THE ARGUMENT     9.0072e+15
> THIS SHOULD TRIGGER AN ERROR MESSAGE
>
>
> SIN RETURNED THE VALUE    -8.4893e-01
>  
>
Accurate results should not be expected for trig functions of large 
arguments, since tiny perturbations of the arguments give rise to large 
perturbations of the results.  However, these large but finite arguments 
are not outside the mathematical domain of definition of the functions, 
so I felt that a domain error (EDOM) was not called for in this case. 
 It is actually possible to return accurate results for large arguments, 
if multiple precision is used to reduce the arguments.  libm does this, 
but one pays a price in execution time, which I thought was not 
warranted, especially if those arguments are contaminated by roundoff 
error from previous computation.  Those who want trig routines that can 
perform such feats can use libm, but we needed an alternative for the 
user who wanted speed.  (I had to give up very little in accuracy to get 
that speed, thanks to the excellent coprocessors on Intel machines.)  My 
choice was to use the built-in coprocessor value for pi to reduce the 
arguments, and not flag an error.  The trig functions of large arguments 
are thus very inaccurate, but the argument reduction is identical for 
all the trig functions, so they at least obey the usual identities, e.g. 
sin^2 + cos^2 = 1.

>tpower.djl, error message:
> ( 2.0000000e+00) ** (-1.0740000e+03) WILL BE COMPUTED.
> THIS SHOULD NOT TRIGGER AN ERROR MESSAGE
>
>pow(): Output of function out of range (ERANGE)
> THE VALUE RETURNED IS     4.9407e-324
>
I chose to flag an error here because the result is denormalized.  I 
could have chosen to raise ERANGE only when a finite result gets 
assigned a value of zero of infinity.  However, when the results get so 
small that they are denormalized, full precision is no longer attained. 
 In the case that ELEFUNT evaluated above, the result has only one bit 
of significance left in it, and I felt that one ought to know about 
that, since the answer is otherwise well defined in terms of the 
arguments.  The error is reported by setting errno, rather than 
generating an exception, so one can easily ignore it.

The difference between this case and the previous one is subtle: in the 
former case, the results are what one calls "backward-stable," that is, 
the error in the result (large though it may be) can be explained by a 
small perturbation to the argument.  In the latter, the result is not 
backward-stable when it begins to underflow.

-Eric Rudd
rudd AT cyberoptics DOT com