diff options
author | Szabolcs Nagy <nsz@port70.net> | 2012-12-11 23:06:20 +0100 |
---|---|---|
committer | Szabolcs Nagy <nsz@port70.net> | 2012-12-11 23:06:20 +0100 |
commit | 482ccd2f7497a79ca83e998f54e823e7cedaaa6e (patch) | |
tree | 0e2e1fd4db0a277ece85e32a23ecc68997404514 /src/math/acoshl.c | |
parent | 64623cd59a5e72c6322548bca3827a75d5d11918 (diff) | |
download | musl-482ccd2f7497a79ca83e998f54e823e7cedaaa6e.tar.gz musl-482ccd2f7497a79ca83e998f54e823e7cedaaa6e.tar.bz2 musl-482ccd2f7497a79ca83e998f54e823e7cedaaa6e.tar.xz musl-482ccd2f7497a79ca83e998f54e823e7cedaaa6e.zip |
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]
Diffstat (limited to 'src/math/acoshl.c')
-rw-r--r-- | src/math/acoshl.c | 58 |
1 files changed, 12 insertions, 46 deletions
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 |