From 482ccd2f7497a79ca83e998f54e823e7cedaaa6e Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Tue, 11 Dec 2012 23:06:20 +0100 Subject: [PATCH] math: rewrite inverse hyperbolic functions to be simpler/smaller modifications: * avoid unsigned->signed integer conversion * do not handle special cases when they work correctly anyway * more strict threshold values (0x1p26 instead of 0x1p28 etc) * smaller code, cleaner branching logic * same precision as the old code: acosh(x) has up to 2ulp error in [1,1.125] asinh(x) has up to 1.6ulp error in [0.125,0.5], [-0.5,-0.125] atanh(x) has up to 1.7ulp error in [0.125,0.5], [-0.5,-0.125] --- src/math/acosh.c | 61 +++++++++------------------------------ src/math/acoshf.c | 47 ++++++++----------------------- src/math/acoshl.c | 58 ++++++++------------------------------ src/math/asinh.c | 69 ++++++++++++++------------------------------- src/math/asinhf.c | 62 ++++++++++++++-------------------------- src/math/asinhl.c | 72 ++++++++++++++++------------------------------- src/math/atanh.c | 67 ++++++++++--------------------------------- src/math/atanhf.c | 50 +++++++++----------------------- src/math/atanhl.c | 69 ++++++++++++--------------------------------- 9 files changed, 149 insertions(+), 406 deletions(-) diff --git a/src/math/acosh.c b/src/math/acosh.c index 15f51c6e..4ce9b3d1 100644 --- a/src/math/acosh.c +++ b/src/math/acosh.c @@ -1,54 +1,19 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/e_acosh.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - */ -/* acosh(x) - * Method : - * Based on - * acosh(x) = log [ x + sqrt(x*x-1) ] - * we have - * acosh(x) := log(x)+ln2, if x is large; else - * acosh(x) := log(2x-1/(sqrt(x*x-1)+x)) if x>2; else - * acosh(x) := log1p(t+sqrt(2.0*t+t*t)); where t=x-1. - * - * Special cases: - * acosh(x) is NaN with signal if x<1. - * acosh(NaN) is NaN without signal. - */ - #include "libm.h" -static const double -ln2 = 6.93147180559945286227e-01; /* 0x3FE62E42, 0xFEFA39EF */ - +/* acosh(x) = log(x + sqrt(x*x-1)) */ double acosh(double x) { - double t; - int32_t hx; - uint32_t lx; + union {double f; uint64_t i;} u = {.f = x}; + unsigned e = u.i >> 52 & 0x7ff; + + /* x < 1 domain error is handled in the called functions */ - EXTRACT_WORDS(hx, lx, x); - if (hx < 0x3ff00000) { /* x < 1 */ - return (x-x)/(x-x); - } else if (hx >= 0x41b00000) { /* x > 2**28 */ - if (hx >= 0x7ff00000) /* x is inf of NaN */ - return x+x; - return log(x) + ln2; /* acosh(huge) = log(2x) */ - } else if ((hx-0x3ff00000 | lx) == 0) { - return 0.0; /* acosh(1) = 0 */ - } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ - t = x*x; - return log(2.0*x - 1.0/(x+sqrt(t-1.0))); - } else { /* 1 < x < 2 */ - t = x-1.0; - return log1p(t + sqrt(2.0*t+t*t)); - } + if (e < 0x3ff + 1) + /* |x| < 2, up to 2ulp error in [1,1.125] */ + return log1p(x-1 + sqrt((x-1)*(x-1)+2*(x-1))); + if (e < 0x3ff + 26) + /* |x| < 0x1p26 */ + return log(2*x - 1/(x+sqrt(x*x-1))); + /* |x| >= 0x1p26 or nan */ + return log(x) + 0.693147180559945309417232121458176568; } diff --git a/src/math/acoshf.c b/src/math/acoshf.c index 0f7aae2a..4596085e 100644 --- a/src/math/acoshf.c +++ b/src/math/acoshf.c @@ -1,42 +1,17 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/e_acoshf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" -static const float -ln2 = 6.9314718246e-01; /* 0x3f317218 */ - +/* acosh(x) = log(x + sqrt(x*x-1)) */ float acoshf(float x) { - float t; - int32_t hx; + union {float f; int32_t i;} u = {.f = x}; - GET_FLOAT_WORD(hx, x); - if (hx < 0x3f800000) { /* x < 1 */ - return (x-x)/(x-x); - } else if (hx >= 0x4d800000) { /* x > 2**28 */ - if (hx >= 0x7f800000) /* x is inf of NaN */ - return x + x; - return logf(x) + ln2; /* acosh(huge)=log(2x) */ - } else if (hx == 0x3f800000) { - return 0.0f; /* acosh(1) = 0 */ - } else if (hx > 0x40000000) { /* 2**28 > x > 2 */ - t = x*x; - return logf(2.0f*x - 1.0f/(x+sqrtf(t-1.0f))); - } else { /* 1 < x < 2 */ - t = x-1.0f; - return log1pf(t + sqrtf(2.0f*t+t*t)); - } + if (u.i < 0x3f800000+(1<<23)) + /* x < 2, invalid if x < 1 or nan */ + /* up to 2ulp error in [1,1.125] */ + return log1pf(x-1 + sqrtf((x-1)*(x-1)+2*(x-1))); + if (u.i < 0x3f800000+(12<<23)) + /* x < 0x1p12 */ + return logf(2*x - 1/(x+sqrtf(x*x-1))); + /* x >= 0x1p12 */ + return logf(x) + 0.693147180559945309417232121458176568f; } diff --git a/src/math/acoshl.c b/src/math/acoshl.c index a4024516..472c71cb 100644 --- a/src/math/acoshl.c +++ b/src/math/acoshl.c @@ -1,28 +1,3 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_acoshl.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ -/* acoshl(x) - * Method : - * Based on - * acoshl(x) = logl [ x + sqrtl(x*x-1) ] - * we have - * acoshl(x) := logl(x)+ln2, if x is large; else - * acoshl(x) := logl(2x-1/(sqrtl(x*x-1)+x)) if x>2; else - * acoshl(x) := log1pl(t+sqrtl(2.0*t+t*t)); where t=x-1. - * - * Special cases: - * acoshl(x) is NaN with signal if x<1. - * acoshl(NaN) is NaN without signal. - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -31,29 +6,20 @@ long double acoshl(long double x) return acosh(x); } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -static const long double -ln2 = 6.931471805599453094287e-01L; /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ - +/* acosh(x) = log(x + sqrt(x*x-1)) */ long double acoshl(long double x) { - long double t; - uint32_t se,i0,i1; + union { + long double f; + struct{uint64_t m; int16_t se; uint16_t pad;} i; + } u = {.f = x}; - GET_LDOUBLE_WORDS(se, i0, i1, x); - if (se < 0x3fff || se & 0x8000) { /* x < 1 */ - return (x-x)/(x-x); - } else if (se >= 0x401d) { /* x > 2**30 */ - if (se >= 0x7fff) /* x is inf or NaN */ - return x+x; - return logl(x) + ln2; /* acoshl(huge) = logl(2x) */ - } else if (((se-0x3fff)|i0|i1) == 0) { - return 0.0; /* acosh(1) = 0 */ - } else if (se > 0x4000) { /* x > 2 */ - t = x*x; - return logl(2.0*x - 1.0/(x + sqrtl(t - 1.0))); - } - /* 1 < x <= 2 */ - t = x - 1.0; - return log1pl(t + sqrtl(2.0*t + t*t)); + if (u.i.se < 0x3fff + 1) + /* x < 2, invalid if x < 1 or nan */ + return log1pl(x-1 + sqrtl((x-1)*(x-1)+2*(x-1))); + if (u.i.se < 0x3fff + 32) + /* x < 0x1p32 */ + return logl(2*x - 1/(x+sqrtl(x*x-1))); + return logl(x) + 0.693147180559945309417232121458176568L; } #endif diff --git a/src/math/asinh.c b/src/math/asinh.c index 11bbd71a..4152dc39 100644 --- a/src/math/asinh.c +++ b/src/math/asinh.c @@ -1,55 +1,28 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_asinh.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ -/* asinh(x) - * Method : - * Based on - * asinh(x) = sign(x) * log [ |x| + sqrt(x*x+1) ] - * we have - * asinh(x) := x if 1+x*x=1, - * := sign(x)*(log(x)+ln2)) for large |x|, else - * := sign(x)*log(2|x|+1/(|x|+sqrt(x*x+1))) if|x|>2, else - * := sign(x)*log1p(|x| + x^2/(1 + sqrt(1+x^2))) - */ - #include "libm.h" -static const double -ln2 = 6.93147180559945286227e-01, /* 0x3FE62E42, 0xFEFA39EF */ -huge= 1.00000000000000000000e+300; - +/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */ double asinh(double x) { - double t,w; - int32_t hx,ix; + union {double f; uint64_t i;} u = {.f = x}; + unsigned e = u.i >> 52 & 0x7ff; + unsigned s = u.i >> 63; - GET_HIGH_WORD(hx, x); - ix = hx & 0x7fffffff; - if (ix >= 0x7ff00000) /* x is inf or NaN */ - return x+x; - if (ix < 0x3e300000) { /* |x| < 2**-28 */ - /* return x inexact except 0 */ - if (huge+x > 1.0) - return x; - } - if (ix > 0x41b00000) { /* |x| > 2**28 */ - w = log(fabs(x)) + ln2; - } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */ - t = fabs(x); - w = log(2.0*t + 1.0/(sqrt(x*x+1.0)+t)); - } else { /* 2.0 > |x| > 2**-28 */ - t = x*x; - w =log1p(fabs(x) + t/(1.0+sqrt(1.0+t))); + /* |x| */ + u.i &= (uint64_t)-1/2; + x = u.f; + + if (e >= 0x3ff + 26) { + /* |x| >= 0x1p26 or inf or nan */ + x = log(x) + 0.693147180559945309417232121458176568; + } else if (e >= 0x3ff + 1) { + /* |x| >= 2 */ + x = log(2*x + 1/(sqrt(x*x+1)+x)); + } else if (e >= 0x3ff - 26) { + /* |x| >= 0x1p-26, up to 1.6ulp error in [0.125,0.5] */ + x = log1p(x + x*x/(sqrt(x*x+1)+1)); + } else { + /* |x| < 0x1p-26, raise inexact if x != 0 */ + FORCE_EVAL(x + 0x1p1000); } - if (hx > 0) - return w; - return -w; + return s ? -x : x; } diff --git a/src/math/asinhf.c b/src/math/asinhf.c index efe3af94..fc9f0911 100644 --- a/src/math/asinhf.c +++ b/src/math/asinhf.c @@ -1,48 +1,28 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_asinhf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" -static const float -ln2 = 6.9314718246e-01, /* 0x3f317218 */ -huge= 1.0000000000e+30; - +/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */ float asinhf(float x) { - float t,w; - int32_t hx,ix; + union {float f; uint32_t i;} u = {.f = x}; + uint32_t i = u.i & 0x7fffffff; + unsigned s = u.i >> 31; - GET_FLOAT_WORD(hx, x); - ix = hx & 0x7fffffff; - if (ix >= 0x7f800000) /* x is inf or NaN */ - return x+x; - if (ix < 0x31800000) { /* |x| < 2**-28 */ - /* return x inexact except 0 */ - if (huge+x > 1.0f) - return x; - } - if (ix > 0x4d800000) { /* |x| > 2**28 */ - w = logf(fabsf(x)) + ln2; - } else if (ix > 0x40000000) { /* 2**28 > |x| > 2.0 */ - t = fabsf(x); - w = logf(2.0f*t + 1.0f/(sqrtf(x*x+1.0f)+t)); - } else { /* 2.0 > |x| > 2**-28 */ - t = x*x; - w =log1pf(fabsf(x) + t/(1.0f+sqrtf(1.0f+t))); + /* |x| */ + u.i = i; + x = u.f; + + if (i >= 0x3f800000 + (12<<23)) { + /* |x| >= 0x1p12 or inf or nan */ + x = logf(x) + 0.693147180559945309417232121458176568f; + } else if (i >= 0x3f800000 + (1<<23)) { + /* |x| >= 2 */ + x = logf(2*x + 1/(sqrtf(x*x+1)+x)); + } else if (i >= 0x3f800000 - (12<<23)) { + /* |x| >= 0x1p-12, up to 1.6ulp error in [0.125,0.5] */ + x = log1pf(x + x*x/(sqrtf(x*x+1)+1)); + } else { + /* |x| < 0x1p-12, raise inexact if x!=0 */ + FORCE_EVAL(x + 0x1p120f); } - if (hx > 0) - return w; - return -w; + return s ? -x : x; } diff --git a/src/math/asinhl.c b/src/math/asinhl.c index dc5dd71f..db966246 100644 --- a/src/math/asinhl.c +++ b/src/math/asinhl.c @@ -1,25 +1,3 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/ld80/s_asinhl.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ -/* asinhl(x) - * Method : - * Based on - * asinhl(x) = signl(x) * logl [ |x| + sqrtl(x*x+1) ] - * we have - * asinhl(x) := x if 1+x*x=1, - * := signl(x)*(logl(x)+ln2)) for large |x|, else - * := signl(x)*logl(2|x|+1/(|x|+sqrtl(x*x+1))) if|x|>2, else - * := signl(x)*log1pl(|x| + x^2/(1 + sqrtl(1+x^2))) - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -28,35 +6,33 @@ long double asinhl(long double x) return asinh(x); } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -static const long double -ln2 = 6.931471805599453094287e-01L, /* 0x3FFE, 0xB17217F7, 0xD1CF79AC */ -huge = 1.000000000000000000e+4900L; - +/* asinh(x) = sign(x)*log(|x|+sqrt(x*x+1)) ~= x - x^3/6 + o(x^5) */ long double asinhl(long double x) { - long double t,w; - int32_t hx,ix; + union { + long double f; + struct{uint64_t m; uint16_t se; uint16_t pad;} i; + } u = {.f = x}; + unsigned e = u.i.se & 0x7fff; + unsigned s = u.i.se >> 15; - GET_LDOUBLE_EXP(hx, x); - ix = hx & 0x7fff; - if (ix == 0x7fff) - return x + x; /* x is inf or NaN */ - if (ix < 0x3fde) { /* |x| < 2**-34 */ - /* return x, raise inexact if x != 0 */ - if (huge+x > 1.0) - return x; - } - if (ix > 0x4020) { /* |x| > 2**34 */ - w = logl(fabsl(x)) + ln2; - } else if (ix > 0x4000) { /* 2**34 > |x| > 2.0 */ - t = fabsl(x); - w = logl(2.0*t + 1.0/(sqrtl(x*x + 1.0) + t)); - } else { /* 2.0 > |x| > 2**-28 */ - t = x*x; - w =log1pl(fabsl(x) + t/(1.0 + sqrtl(1.0 + t))); + /* |x| */ + u.i.se = e; + x = u.f; + + if (e >= 0x3fff + 32) { + /* |x| >= 0x1p32 or inf or nan */ + x = logl(x) + 0.693147180559945309417232121458176568L; + } else if (e >= 0x3fff + 1) { + /* |x| >= 2 */ + x = logl(2*x + 1/(sqrtl(x*x+1)+x)); + } else if (e >= 0x3fff - 32) { + /* |x| >= 0x1p-32 */ + x = log1pl(x + x*x/(sqrtl(x*x+1)+1)); + } else { + /* |x| < 0x1p-32, raise inexact if x!=0 */ + FORCE_EVAL(x + 0x1p1000); } - if (hx & 0x8000) - return -w; - return w; + return s ? -x : x; } #endif diff --git a/src/math/atanh.c b/src/math/atanh.c index dbe241d1..84a84c69 100644 --- a/src/math/atanh.c +++ b/src/math/atanh.c @@ -1,58 +1,21 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/e_atanh.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunSoft, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - * - */ -/* atanh(x) - * Method : - * 1.Reduced x to positive by atanh(-x) = -atanh(x) - * 2.For x>=0.5 - * 1 2x x - * atanh(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x - * - * For x<0.5 - * atanh(x) = 0.5*log1p(2x+2x*x/(1-x)) - * - * Special cases: - * atanh(x) is NaN if |x| > 1 with signal; - * atanh(NaN) is that NaN with no signal; - * atanh(+-1) is +-INF with signal. - * - */ - #include "libm.h" -static const double huge = 1e300; - +/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ double atanh(double x) { - double t; - int32_t hx,ix; - uint32_t lx; + union {double f; uint64_t i;} u = {.f = x}; + unsigned e = u.i >> 52 & 0x7ff; + unsigned s = u.i >> 63; + + /* |x| */ + u.i &= (uint64_t)-1/2; + x = u.f; - EXTRACT_WORDS(hx, lx, x); - ix = hx & 0x7fffffff; - if ((ix | ((lx|-lx)>>31)) > 0x3ff00000) /* |x| > 1 */ - return (x-x)/(x-x); - if (ix == 0x3ff00000) - return x/0.0; - if (ix < 0x3e300000 && (huge+x) > 0.0) /* x < 2**-28 */ - return x; - SET_HIGH_WORD(x, ix); - if (ix < 0x3fe00000) { /* x < 0.5 */ - t = x+x; - t = 0.5*log1p(t + t*x/(1.0-x)); - } else - t = 0.5*log1p((x+x)/(1.0-x)); - if (hx >= 0) - return t; - return -t; + if (e < 0x3ff - 1) { + /* |x| < 0.5, up to 1.7ulp error */ + x = 0.5*log1p(2*x + 2*x*x/(1-x)); + } else { + x = 0.5*log1p(2*x/(1-x)); + } + return s ? -x : x; } diff --git a/src/math/atanhf.c b/src/math/atanhf.c index 2be780bb..ca106efc 100644 --- a/src/math/atanhf.c +++ b/src/math/atanhf.c @@ -1,42 +1,20 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/e_atanhf.c */ -/* - * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. - */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ - #include "libm.h" -static const float huge = 1e30; - +/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ float atanhf(float x) { - float t; - int32_t hx,ix; + union {float f; uint32_t i;} u = {.f = x}; + unsigned s = u.i >> 31; + + /* |x| */ + u.i &= 0x7fffffff; + x = u.f; - GET_FLOAT_WORD(hx, x); - ix = hx & 0x7fffffff; - if (ix > 0x3f800000) /* |x| > 1 */ - return (x-x)/(x-x); - if (ix == 0x3f800000) - return x/0.0f; - if (ix < 0x31800000 && huge+x > 0.0f) /* x < 2**-28 */ - return x; - SET_FLOAT_WORD(x, ix); - if (ix < 0x3f000000) { /* x < 0.5 */ - t = x+x; - t = 0.5f*log1pf(t + t*x/(1.0f-x)); - } else - t = 0.5f*log1pf((x+x)/(1.0f-x)); - if (hx >= 0) - return t; - return -t; + if (u.i < 0x3f800000 - (1<<23)) { + /* |x| < 0.5, up to 1.7ulp error */ + x = 0.5f*log1pf(2*x + 2*x*x/(1-x)); + } else { + x = 0.5f*log1pf(2*x/(1-x)); + } + return s ? -x : x; } diff --git a/src/math/atanhl.c b/src/math/atanhl.c index 931bae32..b4c5e58b 100644 --- a/src/math/atanhl.c +++ b/src/math/atanhl.c @@ -1,31 +1,3 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_atanh.c */ -/* - * ==================================================== - * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. - * - * Developed at SunPro, a Sun Microsystems, Inc. business. - * Permission to use, copy, modify, and distribute this - * software is freely granted, provided that this notice - * is preserved. - * ==================================================== - */ -/* atanhl(x) - * Method : - * 1.Reduced x to positive by atanh(-x) = -atanh(x) - * 2.For x>=0.5 - * 1 2x x - * atanhl(x) = --- * log(1 + -------) = 0.5 * log1p(2 * --------) - * 2 1 - x 1 - x - * - * For x<0.5 - * atanhl(x) = 0.5*log1pl(2x+2x*x/(1-x)) - * - * Special cases: - * atanhl(x) is NaN if |x| > 1 with signal; - * atanhl(NaN) is that NaN with no signal; - * atanhl(+-1) is +-INF with signal. - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -34,31 +6,26 @@ long double atanhl(long double x) return atanh(x); } #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 -static const long double huge = 1e4900L; - +/* atanh(x) = log((1+x)/(1-x))/2 = log1p(2x/(1-x))/2 ~= x + x^3/3 + o(x^5) */ long double atanhl(long double x) { - long double t; - int32_t ix; - uint32_t se,i0,i1; + union { + long double f; + struct{uint64_t m; uint16_t se; uint16_t pad;} i; + } u = {.f = x}; + unsigned e = u.i.se & 0x7fff; + unsigned s = u.i.se >> 15; + + /* |x| */ + u.i.se = e; + x = u.f; - GET_LDOUBLE_WORDS(se, i0, i1, x); - ix = se & 0x7fff; - if ((ix+((((i0&0x7fffffff)|i1)|(-((i0&0x7fffffff)|i1)))>>31)) > 0x3fff) - /* |x| > 1 */ - return (x-x)/(x-x); - if (ix == 0x3fff) - return x/0.0; - if (ix < 0x3fe3 && huge+x > 0.0) /* x < 2**-28 */ - return x; - SET_LDOUBLE_EXP(x, ix); - if (ix < 0x3ffe) { /* x < 0.5 */ - t = x + x; - t = 0.5*log1pl(t + t*x/(1.0 - x)); - } else - t = 0.5*log1pl((x + x)/(1.0 - x)); - if (se <= 0x7fff) - return t; - return -t; + if (e < 0x3fff - 1) { + /* |x| < 0.5, up to 1.7ulp error */ + x = 0.5*log1pl(2*x + 2*x*x/(1-x)); + } else { + x = 0.5*log1pl(2*x/(1-x)); + } + return s ? -x : x; } #endif -- GitLab