From 1aec620f9366c29d761fe42b3e02bd8024685db3 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Sun, 16 Dec 2012 19:23:51 +0100 Subject: math: finished cosh.c cleanup 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. --- src/math/coshf.c | 49 ++++++++++++++----------------------------------- 1 file changed, 14 insertions(+), 35 deletions(-) (limited to 'src/math/coshf.c') 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; } -- cgit v1.2.1