fluencelabs/musl
1
0
Fork 0
mirror of https://github.com/fluencelabs/musl synced 2025-06-30 15:11:55 +00:00
Code Issues Projects Releases Wiki Activity

math: finished cosh.c cleanup

Browse Source 
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.
This commit is contained in:
Szabolcs Nagy
2012-12-16 19:23:51 +01:00
parent 58bba42d1b
commit 1aec620f93
3 changed files with 49 additions and 142 deletions
Download Patch File Download Diff File Expand all files Collapse all files

73
src/math/cosh.c
View File

@ -1,71 +1,40 @@
/* origin: FreeBSD /usr/src/lib/msun/src/e_cosh.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.
* ====================================================
*/
/* cosh(x)
* Method :
* mathematically cosh(x) if defined to be (exp(x)+exp(-x))/2
* 1. Replace x by |x| (cosh(x) = cosh(-x)).
* 2.
* [ exp(x) - 1 ]^2
* 0 <= x <= ln2/2 : cosh(x) := 1 + -------------------
* 2*exp(x)
*
* exp(x) + 1/exp(x)
* ln2/2 <= x <= 22 : cosh(x) := -------------------
* 2
* 22 <= x <= lnovft : cosh(x) := exp(x)/2
* lnovft <= x <= ln2ovft: cosh(x) := exp(x/2)/2 * exp(x/2)
* ln2ovft < x : cosh(x) := inf (overflow)
*
* Special cases:
* cosh(x) is |x| if x is +INF, -INF, or NaN.
* only cosh(0)=1 is exact for finite x.
*/
#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 ix;
uint32_t w;
double t;
/* |x| */
u.i &= (uint64_t)-1/2;
x = u.f;
ix = u.i >> 32;
w = u.i >> 32;
/* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */
if (ix < 0x3fd62e43) {
t = expm1(x);
if (ix < 0x3c800000)
/* |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| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|))/2; */
if (ix < 0x40360000) {
/* |x| < log(DBL_MAX) */
if (w < 0x40862e42) {
t = exp(x);
return 0.5*t + 0.5/t;
/* note: if x>log(0x1p26) then the 1/t is not needed */
return 0.5*(t + 1/t);
}
/* |x| in [22, log(maxdouble)] return 0.5*exp(|x|) */
if (ix < 0x40862e42)
return 0.5*exp(x);
/* |x| in [log(maxdouble), overflowthresold] */
if (ix <= 0x408633ce)
return __expo2(x);
/* overflow (or nan) */
x *= 0x1p1023;
return x;
/* |x| > log(DBL_MAX) or nan */
/* note: the result is stored to handle overflow */
t = __expo2(x);
return t;
}