Date: Sun, 10 May 1998 12:27:41 +0300 (IDT)
From: Eli Zaretskii <eliz AT is DOT elta DOT co DOT il>
To: Kbwms <Kbwms AT aol DOT com>
cc: Jens Bischoff <j DOT bischoff AT airbus DOT dasa DOT de>, DJ Delorie <dj AT delorie DOT com>,
        djgpp-workers AT delorie DOT com
Subject: Re: Code to Fix sinh() in libm.a
In-Reply-To: <5bc654c6.3551fb66@aol.com>
Message-ID: <Pine.SUN.3.91.980510122651.16252N-100000@is>
MIME-Version: 1.0
Content-Type: TEXT/PLAIN; charset=US-ASCII
Precedence: bulk


On Thu, 7 May 1998, Kbwms wrote:

> Appended is the code, with driver attached, for a solution to the
> sinh() problem that I reported recently.  The current version of
> the function in libm.a calculates inaccurate values for small
> arguments.

Thanks.  However, I think this needs some more work.

For starters, it is less robust than the library version for extreme
values of the argument.  For example, if the argument causes `exp' to
overflow, so will this version of `sinh', for no good reason; the
current library version is much more defensive about large argument
values.  (Personally, I would just add a separate approximation region
to the original code in libm.a sources, rather than rewriting the
whole thing from scratch.)

Next, please run the new version through tsinh test from the ELEFUNT
suite (it's in v2/djtst201.zip).  The number of argument values in the
small test program that you used is too small; ELEFUNT uses thousands
of random values, and some specific extreme values on top of that.

Also, please note that the way to submit patches is to prepare a file
with context diffs, like this:

	diff -c old-version new-version > diffs.txt

then send the diffs.  (A port of GNU `diff' program is available in
v2gnu/dif271b.zip.)  For the best results, run `diff' from the root of
the DJGPP installation, and use the relative path names of the source
files, like this:

  diff -c src/libm/src/e_sinh.old src/libm/src/e_sinh.c > diffs.txt

(it doesn't really matter how you call the old version--it can be
e.g. e_sinh.bak or anything else, provided that the new version has
the right name).  This makes applying the diffs a snap.

Finally, patches to library sources should be sent to
djgpp-workers AT delorie DOT com, so that people who work on DJGPP
development could get a chance to critique them (I have already
redirected this thread there).

> ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
> /* ------ */
> /* sinh.c */
> /* ------ */
> #include <math.h>
> /* ------------------------------------------ */
> /* sinh - Returns hyperbolic sine of argument */
> /* ------------------------------------------ */
> double
> sinh(double x)
> {
>     /* --------------------------------------- */
>     /* Approxination no. 1906 from Hart, et al */
>     /* --------------------------------------- */
>     long double C[] = {
> 	    .10000000000000000000428964e1,
> 	    .16666666666666666026158977e0,
> 	    .83333333333336099945076e-2,
> 	    .19841269840742957068490e-3,
> 	    .27557319739041248550000e-5,
> 	    .25051838782995770000000e-7,
> 	    .16130886652601000000000e-9
>     };
> 
>     long double RetVal;
> 
>     if (fabs(x) <= 0.5)
>     {
> 	long double Xsq = x*x;
> 
> 	RetVal = (Xsq*(C[1] +
> 		  Xsq*(C[2] +
> 		  Xsq*(C[3] +
> 		  Xsq*(C[4] +
> 		  Xsq*(C[5] +
> 		  Xsq*(C[6])))))) + C[0]) * x;
>     }
>     else
>     {
> 	long double D = exp(x);
> 
> 	RetVal = 0.5L * D - 0.5L / D;
>     }
> 
>     return RetVal;
> }
> # if defined TEST
> #include <stdio.h>
> int
> main()
> {
>     double  Ans, Arg;
>     Arg = 1.e-18;
>     Ans = sinh(Arg);
>     printf("sinh(%+.15g) = %+.17e\n", Arg, Ans);
>     Ans = sinh(-Arg);
>     printf("sinh(%+.15g) = %+.17e\n", -Arg, Ans);
>     printf("Verify: sinh(1e-18) = 1.00000000000000000000E-18\n");
>     Arg = 1.0e-7;
>     Ans = sinh(Arg);
>     printf("sinh(%+.15g) = %+.17e\n", Arg, Ans);
>     Ans = sinh(-Arg);
>     printf("sinh(%+.15g) = %+.17e\n", -Arg, Ans);
>     printf("Verify: sinh(1e-7)  = 1.00000000000000166667E-7\n");
>     Arg = .5;
>     Ans = sinh(Arg);
>     printf("sinh(%+.15g) = %+.17e\n", Arg, Ans);
>     Ans = sinh(-Arg);
>     printf("sinh(%+.15g) = %+.17e\n", -Arg, Ans);
>     printf("Verify: sinh(.5)    = 5.21095305493747361622E-1\n");
>     Arg = .6;
>     Ans = sinh(Arg);
>     printf("sinh(%+.15g) = %+.17e\n", Arg, Ans);
>     Ans = sinh(-Arg);
>     printf("sinh(%+.15g) = %+.17e\n", -Arg, Ans);
>     printf("Verify: sinh(.6)    = 6.36653582148241271123E-1\n");
>     Arg = 1;
>     Ans = sinh(Arg);
>     printf("sinh(%+.15g) = %+.17e\n", Arg, Ans);
>     Ans = sinh(-Arg);
>     printf("sinh(%+.15g) = %+.17e\n", -Arg, Ans);
>     printf("Verify: sinh(1)     = 1.17520119364380145688E0\n");
>     printf("Verify: sinh(10)    = 1.10132328747033933772E4\n");
>     printf("Verify: sinh(23)    = 4.87240172312445129997E9\n");
>     Arg = 10;
>     Ans = sinh(Arg);
>     printf("sinh(%+.15g) = %+.17e\n", Arg, Ans);
>     Ans = sinh(-Arg);
>     printf("sinh(%+.15g) = %+.17e\n", -Arg, Ans);
>     printf("Verify: sinh(10)    = 1.10132328747033933772E4\n");
>     Arg = 23;
>     Ans = sinh(Arg);
>     printf("sinh(%+.15g) = %+.17e\n", Arg, Ans);
>     Ans = sinh(-Arg);
>     printf("sinh(%+.15g) = %+.17e\n", -Arg, Ans);
>     printf("Verify: sinh(23)    = 4.87240172312445129997E9\n");
> 
>     return 0;
> }
> # endif
>