From: Eric Rudd <rudd AT cyberoptics DOT com>
Newsgroups: comp.os.msdos.djgpp
Subject: Re: Code to Fix sinh() in libm.a
Date: Wed, 13 May 1998 17:35:55 -0500
Organization: CyberOptics
Lines: 51
Message-ID: <355A204B.3B8A78EF@cyberoptics.com>
References: <42dba229 DOT 3558c70d AT aol DOT com>
Reply-To: rudd AT cyberoptics DOT com
NNTP-Posting-Host: rudd.cyberoptics.com
Mime-Version: 1.0
Content-Type: text/plain; charset=us-ascii
Content-Transfer-Encoding: 7bit
To: djgpp AT delorie DOT com
DJ-Gateway: from newsgroup comp.os.msdos.djgpp
Precedence: bulk

Kbwms wrote:
> 
> The original complaint was with argument on the order of 1.e-7.  When
> I ran tsinh  and my version of sinh() on this argument, here is what I got:
> 
>  Special Value x = 1.e07 - should produce 1.00000000000000167e-7
> 9.999999999999999547E-08  1.000000000000001675E-07

I encountered this problem myself a while ago, which prompted an
investigation. Here is my current assessment:

1. The libc.a routines are terse and efficient, but not always accurate.
For instance, the sinh(x) routine computes 0.5*(exp(x) - 1./exp(x)),
which has severe cancellation for small x.

2. The libm.a routines were apparently written for a system without
transcendental support in the coprocessor. The polynomial and rational
approximations they use are considerably (sometimes by a factor of 3)
more time-consuming than the routines in libc.a, and even some of those
have accuracy problems in certain ranges.

To solve the problems with poor accuracy for special arguments, certain
alternate forms can be used. For instance, the x87 coprocessor has
special machine instructions for exp(x)-1 and ln(x+1) that preserve good
relative precision for small x, and sinh(x) can be stably computed from
an accurate exp(x)-1 routine. The polynomial approximation given earlier
in this thread is apt to be unduly time-consuming, compared to using the
native x87 transcendental functions.

I considered the following routines in libc.a unacceptable, and have
rewritten them: cosh exp tanh asin acos hypot acosh atanh asinh sinh. I
also wrote C-callable routines to compute exp(x)-1 and ln(x+1), since
those functions are generally useful.

I can almost hear the whine already: "Eric, if you had these wonderful
routines, why on earth did you keep them to yourself?" Well, since the
elementary transcendentals get heavily used, they need to be held to
very high standards of performance, both in speed and accuracy. I think
they measure up, but "thinking" isn't good enough. Since I have not had
time to test them extensively, I didn't consider them to be ready for
public use.  I ran ELEFUNT on some of them, but I don't consider that a
sufficient test. (Intel tests their coprocessor routines on literally
millions of arguments, and even then there was a certain notorious
problem with FDIV.)

However, I would be happy to make the code available for beta testing.
(What is the best way to do this? The routines, compressed, are about
7k)

-Eric Rudd
rudd AT cyberoptics DOT com