From 60d71b68ef6af3e5fb3e045f7ce2a42a5ff23ac1 Mon Sep 17 00:00:00 2001 From: Raymond Toy <toy.raymond@gmail.com> Date: Sat, 2 Aug 2014 15:30:51 -0700 Subject: [PATCH] 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. --- src/lisp/e_acosh.c | 9 ++++++--- src/lisp/e_atanh.c | 15 ++++++++++----- src/lisp/fdlibm.h | 2 ++ src/lisp/s_asinh.c | 11 +++++++---- 4 files changed, 25 insertions(+), 12 deletions(-) diff --git a/src/lisp/e_acosh.c b/src/lisp/e_acosh.c index 683955410..fb248d544 100644 --- a/src/lisp/e_acosh.c +++ b/src/lisp/e_acosh.c @@ -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)); } } diff --git a/src/lisp/e_atanh.c b/src/lisp/e_atanh.c index 302a89f44..358843698 100644 --- a/src/lisp/e_atanh.c +++ b/src/lisp/e_atanh.c @@ -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; } diff --git a/src/lisp/fdlibm.h b/src/lisp/fdlibm.h index 325d7d42f..72f57f823 100644 --- a/src/lisp/fdlibm.h +++ b/src/lisp/fdlibm.h @@ -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 diff --git a/src/lisp/s_asinh.c b/src/lisp/s_asinh.c index c8094673b..6cae81ff2 100644 --- a/src/lisp/s_asinh.c +++ b/src/lisp/s_asinh.c @@ -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; } -- GitLab