mirror of
https://github.com/fluencelabs/musl
synced 2025-06-29 14:42:01 +00:00
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.
41 lines
764 B
C
41 lines
764 B
C
#include "libm.h"
|
|
|
|
/* cosh(x) = (exp(x) + 1/exp(x))/2
|
|
* = 1 + 0.5*(exp(x)-1)*(exp(x)-1)/exp(x)
|
|
* = 1 + x*x/2 + o(x^4)
|
|
*/
|
|
double cosh(double x)
|
|
{
|
|
union {double f; uint64_t i;} u = {.f = x};
|
|
uint32_t w;
|
|
double t;
|
|
|
|
/* |x| */
|
|
u.i &= (uint64_t)-1/2;
|
|
x = u.f;
|
|
w = u.i >> 32;
|
|
|
|
/* |x| < log(2) */
|
|
if (w < 0x3fe62e42) {
|
|
if (w < 0x3ff00000 - (26<<20)) {
|
|
/* raise inexact if x!=0 */
|
|
FORCE_EVAL(x + 0x1p120f);
|
|
return 1;
|
|
}
|
|
t = expm1(x);
|
|
return 1 + t*t/(2*(1+t));
|
|
}
|
|
|
|
/* |x| < log(DBL_MAX) */
|
|
if (w < 0x40862e42) {
|
|
t = exp(x);
|
|
/* note: if x>log(0x1p26) then the 1/t is not needed */
|
|
return 0.5*(t + 1/t);
|
|
}
|
|
|
|
/* |x| > log(DBL_MAX) or nan */
|
|
/* note: the result is stored to handle overflow */
|
|
t = __expo2(x);
|
|
return t;
|
|
}
|