Implementation for coshl.

Bruce Evans brde at optusnet.com.au
Mon Jun 17 23:29:50 UTC 2013


[Reply to my latest saved mail in this thread.  Mostly not in response to it.]

On Sun, 9 Jun 2013, David Schultz wrote:

> On Mon, Jun 10, 2013, Bruce Evans wrote:
>> On Sun, 9 Jun 2013, Steve Kargl wrote:

>>> The former would require a refactoring of
>>> s_expl.c into a kernel __kernel_expl(x, hi, lo).  I have no plans on
>>> pursuing this at the this time.
>>
>> But you need this soon for __ldexp_exp() and __ldexp_cexp(), which are
>> needed for hyperbolic functions (the large args already fixed for float
>> and double precision) and for cexp() and trig and hyperbolic complex
>> functions.  It is much easier to implement these using a kernel.  I do
>> this only for float precision.
>
> That's fine for now.  We need a better k_exp, too.  I believe
> Bruce was working on one.

The wrappers turned out to be fairly easy.  Especially __ldexp_exp*(()
and __ldexp_cexp*(), given the main wrapper.

Are you going to finish coshl() and sinhl()?  I'm getting closer to
hacking on them.  I'm now in the middle of doing this for cosh():

diff -u2 e_cosh.c~ e_cosh.c
% --- e_cosh.c~	Thu May 30 18:14:16 2013
% +++ e_cosh.c	Mon Jun 17 23:54:12 2013
% @@ -39,10 +39,23 @@
%  #include "math_private.h"
% 
% -static const double one = 1.0, half=0.5, huge = 1.0e300;
% +static const double one = 1.0, half = 0.5;
% +static const volatile double huge = 1.0e300, tiny = 1.0e-300;

Part of fixing old bugs.

% +static const double
% +/*
% + * Domain [-1, 1], range ~[-2.2746e-18, 2.2748e-18]:
% + * |cosh(x) - c(x)| < 2**-58.6
% + */
% +C2  =  5.0000000000000000e-1,		/*  0x10000000000000.0p-53 */
% +C4  =  4.1666666666665422e-2,		/*  0x155555555554a2.0p-57 */
% +C6  =  1.3888888889057476e-3,		/*  0x16c16c16c29bca.0p-62 */
% +C8  =  2.4801587216389553e-5,		/*  0x1a01a01881edc9.0p-68 */
% +C10 =  2.7557340350974160e-7,		/*  0x127e50a5570530.0p-74 */
% +C12 =  2.0873996906899188e-9,		/*  0x11ee3d8f036513.0p-81 */
% +C14 =  1.1653016291730710e-11;		/*  0x19a010a270ce5a.0p-89 */

New method for small x.

% 
%  double
%  __ieee754_cosh(double x)
%  {
% -	double t,w;
% +	double_t t,x2;
%  	int32_t ix;
%

double -> double_t prepares for using the kernel and avoids discarding
accuracy on i386.

% @@ -54,13 +67,14 @@
%  	if(ix>=0x7ff00000) return x*x; 
% 
% -    /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
% -	if(ix<0x3fd62e43) {
% -	    t = expm1(fabs(x));
% -	    w = one+t;
% -	    if (ix<0x3c800000) return w;	/* cosh(tiny) = 1 */
% -	    return one+(t*t)/(w+w);
% +    /* |x| in [0,1), return 1 or c(x) */

Fix comment to match code.  Update it for new code.

% +	if(ix<0x3ff00000) {
% +	    if (ix<0x3c800000)
% +		return one+tiny;	/* cosh(tiny) = 1(+) with inexact */

Fix minor bugs in setting inexact

% +	    x2 = x*x;
% +	    return x2*x2*(x2*x2*x2*(x2*(x2*C14 + C12) + C10) +
% +		x2*(x2*C8 + C6) + C4) + x2*C2 + 1;

More accuracy for |x| < 1.

%  	}
% 
% -    /* |x| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */
% +    /* |x| in [1,22], return (exp(|x|)+1/exp(|x|)/2; */
%  	if (ix < 0x40360000) {
%  		t = __ieee754_exp(fabs(x));

I just need to change this and later cases to use the kernel and its
wrappers.  In float precision, this is:

% -    /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */
% -	if (ix < 0x41100000) {
% -		t = __ieee754_expf(fabsf(x));
% -		return half*t+half/t;
% +    /* |x| in [1, 12), return accurate exp(|x|)/2+1/exp(|x|)/2 */
% +	if (ix<0x41400000) {
% +	    k_hexpf(fabsf(x), &hi, &lo);
% +	    return lo + 0.25F/(hi + lo) + hi;
%  	}

Bruce


More information about the freebsd-numerics mailing list