From: Kbwms <Kbwms AT aol DOT com>
Message-ID: <886e92e9.3550fe0c@aol.com>
Date: Wed, 6 May 1998 20:19:23 EDT
To: j DOT bischoff AT airbus DOT dasa DOT de
Cc: djgpp AT delorie DOT com
Mime-Version: 1.0
Subject: Re: Problem with sinh() in libm.a
Content-type: text/plain; charset=US-ASCII
Content-transfer-encoding: 7bit
Precedence: bulk

Dear Jens Bischoff,

On 05-06-98 at 06:08:09 EST you wrote:
>
> Yes there is a problem. The math library uses the following algorithm:
>
> For x <= 22:  sinh(x) ~= (E + E/(E+1))/2, with E = exp(x) - 1
>
> This is just another notation for sinh(x) = 0.5 * (exp(x) - exp(-x)).
>

According to Hart, et al, *Computer Approximations*, page 99, the first
expression is the preferred one.  On the other hand, it helps not a bit
for tiny arguments.  Steve Moshier, *Methods and Programs for Mathematical
Functions*, does what you describe below for arguments less than 1.

> Now it seems to me that the cumulated relative error of this
> expression is higher than masch.eps (~2.2*10^-15) for small x,
> and therefore greater than it should be. Simply spoken: The algorithm
> has to be improved.
>
> The FPU shows a better result since it internally uses all 80 bits,
> so mach.eps ~= 1.08*10^19 (about 4 digits exacter).
>
> To fix this problem you can either patch the math library so that
> it uses the FPU directly, or you can implement a better algorithm
> (for example: for |x| < arsinh(1) calculate sinh(x) as
> sinh(x) = c0 + c1 * x^2 + c2 * x^4 + c3 * x^6 + ... + c6 * x^12.
> Assume c0 = 1, and obtain the other coefficients ci by the minmax
> approximation of sinh(x)/x as a function of x^2.  This algorithm
> is used by the IBM VS Fortran R2 compiler).

Great. Build a fix that works and send it in.

Thanks for your response.


Sincerely,

K.B. Williams