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