diff options
Diffstat (limited to 'src/math')
| -rw-r--r-- | src/math/cosh.c | 73 | ||||
| -rw-r--r-- | src/math/coshf.c | 49 | ||||
| -rw-r--r-- | src/math/coshl.c | 69 |
3 files changed, 49 insertions, 142 deletions
diff --git a/src/math/cosh.c b/src/math/cosh.c index 2fcdea8f..100f8231 100644 --- a/src/math/cosh.c +++ b/src/math/cosh.c @@ -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; } diff --git a/src/math/coshf.c b/src/math/coshf.c index 2bed1258..b09f2ee5 100644 --- a/src/math/coshf.c +++ b/src/math/coshf.c @@ -1,54 +1,33 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/e_coshf.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" float coshf(float x) { union {float f; uint32_t i;} u = {.f = x}; - uint32_t ix; + uint32_t w; float t; /* |x| */ u.i &= 0x7fffffff; x = u.f; - ix = u.i; + w = u.i; - /* |x| in [0,0.5*ln2], return 1+expm1(|x|)^2/(2*exp(|x|)) */ - if (ix < 0x3eb17218) { - t = expm1f(x); - if (ix < 0x39800000) + /* |x| < log(2) */ + if (w < 0x3f317217) { + if (w < 0x3f800000 - (12<<23)) { + FORCE_EVAL(x + 0x1p120f); return 1; + } + t = expm1f(x); return 1 + t*t/(2*(1+t)); } - /* |x| in [0.5*ln2,9], return (exp(|x|)+1/exp(|x|))/2; */ - if (ix < 0x41100000) { + /* |x| < log(FLT_MAX) */ + if (w < 0x42b17217) { t = expf(x); - return 0.5f*t + 0.5f/t; + return 0.5f*(t + 1/t); } - /* |x| in [9, log(maxfloat)] return 0.5f*exp(|x|) */ - if (ix < 0x42b17217) - return 0.5f*expf(x); - - /* |x| in [log(maxfloat), overflowthresold] */ - if (ix <= 0x42b2d4fc) - return __expo2f(x); - - /* |x| > overflowthresold or nan */ - x *= 0x1p127f; - return x; + /* |x| > log(FLT_MAX) or nan */ + t = __expo2f(x); + return t; } diff --git a/src/math/coshl.c b/src/math/coshl.c index 3ea56e00..d09070bb 100644 --- a/src/math/coshl.c +++ b/src/math/coshl.c @@ -1,35 +1,3 @@ -/* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_coshl.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. - * ==================================================== - */ -/* coshl(x) - * Method : - * mathematically coshl(x) if defined to be (exp(x)+exp(-x))/2 - * 1. Replace x by |x| (coshl(x) = coshl(-x)). - * 2. - * [ exp(x) - 1 ]^2 - * 0 <= x <= ln2/2 : coshl(x) := 1 + ------------------- - * 2*exp(x) - * - * exp(x) + 1/exp(x) - * ln2/2 <= x <= 22 : coshl(x) := ------------------- - * 2 - * 22 <= x <= lnovft : coshl(x) := expl(x)/2 - * lnovft <= x <= ln2ovft: coshl(x) := expl(x/2)/2 * expl(x/2) - * ln2ovft < x : coshl(x) := inf (overflow) - * - * Special cases: - * coshl(x) is |x| if x is +INF, -INF, or NaN. - * only coshl(0)=1 is exact for finite x. - */ - #include "libm.h" #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 @@ -45,41 +13,32 @@ long double coshl(long double x) 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; - uint32_t mx,lx; /* |x| */ u.i.se = ex; x = u.f; - mx = u.i.m >> 32; - lx = u.i.m; + w = u.i.m >> 32; - /* |x| in [0,0.5*ln2], return 1+expm1l(|x|)^2/(2*expl(|x|)) */ - if (ex < 0x3fff-2 || (ex == 0x3fff-2 && mx < 0xb17217f7)) { - t = expm1l(x); - if (ex < 0x3fff-64) + /* |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| in [0.5*ln2,22], return (exp(|x|)+1/exp(|x|)/2; */ - if (ex < 0x3fff+4 || (ex == 0x3fff+4 && mx < 0xb0000000)) { + /* |x| < log(LDBL_MAX) */ + if (ex < 0x3fff+13 || (ex == 0x3fff+13 && w < 0xb17217f7)) { t = expl(x); - return 0.5*t + 0.5/t; - } - - /* |x| in [22, ln(maxdouble)] return 0.5*exp(|x|) */ - if (ex < 0x3fff+13 || (ex == 0x3fff+13 && mx < 0xb1700000)) - return 0.5*expl(x); - - /* |x| in [log(maxdouble), log(2*maxdouble)) */ - if (ex == 0x3fff+13 && (mx < 0xb174ddc0 || - (mx == 0xb174ddc0 && lx < 0x31aec0eb))) { - t = expl(0.5*x); - return 0.5*t*t; + return 0.5*(t + 1/t); } - /* |x| >= log(2*maxdouble) or nan */ - return x*0x1p16383L; + /* |x| > log(LDBL_MAX) or nan */ + t = expl(0.5*x); + return 0.5*t*t; } #endif |