Szabolcs Nagy
1aec620f93
changed the algorithm: large input is not special cased (when exp(-x) is small compared to exp(x)) and the threshold values are reevaluated (fdlibm code had a log(2)/2 cutoff for which i could not find justification, log(2) seems to be a better threshold and this was verified empirically) the new code is simpler, makes smaller binaries and should be faster for common cases the old comments were removed as they are no longer true for the new algorithm and the fdlibm copyright was dropped as well because there is no common code or idea with the original anymore except for trivial ones.
45 lines
827 B
C
45 lines
827 B
C
#include "libm.h"
|
|
|
|
#if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
|
|
long double coshl(long double x)
|
|
{
|
|
return cosh(x);
|
|
}
|
|
#elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
|
|
long double coshl(long double x)
|
|
{
|
|
union {
|
|
long double f;
|
|
struct{uint64_t m; uint16_t se; uint16_t pad;} i;
|
|
} u = {.f = x};
|
|
unsigned ex = u.i.se & 0x7fff;
|
|
uint32_t w;
|
|
long double t;
|
|
|
|
/* |x| */
|
|
u.i.se = ex;
|
|
x = u.f;
|
|
w = u.i.m >> 32;
|
|
|
|
/* |x| < log(2) */
|
|
if (ex < 0x3fff-1 || (ex == 0x3fff-1 && w < 0xb17217f7)) {
|
|
if (ex < 0x3fff-32) {
|
|
FORCE_EVAL(x + 0x1p120f);
|
|
return 1;
|
|
}
|
|
t = expm1l(x);
|
|
return 1 + t*t/(2*(1+t));
|
|
}
|
|
|
|
/* |x| < log(LDBL_MAX) */
|
|
if (ex < 0x3fff+13 || (ex == 0x3fff+13 && w < 0xb17217f7)) {
|
|
t = expl(x);
|
|
return 0.5*(t + 1/t);
|
|
}
|
|
|
|
/* |x| > log(LDBL_MAX) or nan */
|
|
t = expl(0.5*x);
|
|
return 0.5*t*t;
|
|
}
|
|
#endif
|