Date: Tue, 18 Jun 2013 09:29:41 +1000 (EST) From: Bruce Evans <brde@optusnet.com.au> To: David Schultz <das@FreeBSD.org> Cc: freebsd-numerics@FreeBSD.org, Steve Kargl <sgk@troutmask.apl.washington.edu> Subject: Re: Implementation for coshl. Message-ID: <20130618091306.A787@besplex.bde.org> In-Reply-To: <20130610055834.GA68643@zim.MIT.EDU> References: <20130610003645.GA16444@troutmask.apl.washington.edu> <20130610110740.V24058@besplex.bde.org> <20130610055834.GA68643@zim.MIT.EDU>
next in thread | previous in thread | raw e-mail | index | archive | help
[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
Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20130618091306.A787>