Commit 60d71b68 authored by Raymond Toy's avatar Raymond Toy
Browse files

Fix up the inverse hyperbolics.

 * fdlibm.h:
   * Declare fdlibm_log1p and __ieee754_log
 * e_acosh.c:
   * Use unions to get the parts of a double.
   * Use fdlibm_log1p instead of log1p.
 * e_atanh.c:
   * Use unions to get the parts of a double.
   * Use fdlibm_log1p instead of log1p.
 * s_asinh.c:
   * Rename from asinh to fdlibm_asinh.
   * Use unions to get the parts of a double.
   * Use fdlibm_log1p instead of log1p.
parent 1b876280
......@@ -45,7 +45,10 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
{
double t;
int hx;
hx = __HI(x);
union { int i[2]; double d; } ux;
ux.d = x;
hx = ux.i[HIWORD];
if(hx<0x3ff00000) { /* x < 1 */
return (x-x)/(x-x);
} else if(hx >=0x41b00000) { /* x > 2**28 */
......@@ -53,13 +56,13 @@ ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */
return x+x;
} else
return __ieee754_log(x)+ln2; /* acosh(huge)=log(2x) */
} else if(((hx-0x3ff00000)|__LO(x))==0) {
} else if(((hx-0x3ff00000)|ux.i[LOWORD])==0) {
return 0.0; /* acosh(1) = 0 */
} else if (hx > 0x40000000) { /* 2**28 > x > 2 */
t=x*x;
return __ieee754_log(2.0*x-one/(x+sqrt(t-one)));
} else { /* 1<x<2 */
t = x-one;
return log1p(t+sqrt(2.0*t+t*t));
return fdlibm_log1p(t+sqrt(2.0*t+t*t));
}
}
......@@ -50,19 +50,24 @@ static double zero = 0.0;
double t;
int hx,ix;
unsigned lx;
hx = __HI(x); /* high word */
lx = __LO(x); /* low word */
union { int i[2]; double d; } ux;
ux.d = x;
hx = ux.i[HIWORD]; /* high word */
lx = ux.i[LOWORD]; /* low word */
ix = hx&0x7fffffff;
if ((ix|((lx|(-lx))>>31))>0x3ff00000) /* |x|>1 */
return (x-x)/(x-x);
if(ix==0x3ff00000)
return x/zero;
if(ix<0x3e300000&&(huge+x)>zero) return x; /* x<2**-28 */
__HI(x) = ix; /* x <- |x| */
ux.d = x;
ux.i[HIWORD] = ix; /* x <- |x| */
x = ux.d;
if(ix<0x3fe00000) { /* x < 0.5 */
t = x+x;
t = 0.5*log1p(t+t*x/(one-x));
t = 0.5*fdlibm_log1p(t+t*x/(one-x));
} else
t = 0.5*log1p((x+x)/(one-x));
t = 0.5*fdlibm_log1p((x+x)/(one-x));
if(hx>=0) return t; else return -t;
}
......@@ -52,7 +52,9 @@ extern double fdlibm_sin(double x);
extern double fdlibm_cos(double x);
extern double fdlibm_tan(double x);
extern double fdlibm_expm1(double x);
extern double fdlibm_log1p(double x);
extern double __ieee754_exp(double x);
extern double __ieee754_log(double x);
#endif
......@@ -34,15 +34,18 @@ ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */
huge= 1.00000000000000000000e+300;
#ifdef __STDC__
double asinh(double x)
double fdlibm_asinh(double x)
#else
double asinh(x)
double fdlibm_asinh(x)
double x;
#endif
{
double t,w;
int hx,ix;
hx = __HI(x);
union { int i[2]; double d; } ux;
ux.d = x;
hx = ux.i[HIWORD];
ix = hx&0x7fffffff;
if(ix>=0x7ff00000) return x+x; /* x is inf or NaN */
if(ix< 0x3e300000) { /* |x|<2**-28 */
......@@ -55,7 +58,7 @@ huge= 1.00000000000000000000e+300;
w = __ieee754_log(2.0*t+one/(sqrt(x*x+one)+t));
} else { /* 2.0 > |x| > 2**-28 */
t = x*x;
w =log1p(fabs(x)+t/(one+sqrt(one+t)));
w =fdlibm_log1p(fabs(x)+t/(one+sqrt(one+t)));
}
if(hx>0) return w; else return -w;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment