diff options
Diffstat (limited to 'lib/mlibc/options/ansi/musl-generic-math/expm1f.c')
-rw-r--r-- | lib/mlibc/options/ansi/musl-generic-math/expm1f.c | 111 |
1 files changed, 0 insertions, 111 deletions
diff --git a/lib/mlibc/options/ansi/musl-generic-math/expm1f.c b/lib/mlibc/options/ansi/musl-generic-math/expm1f.c deleted file mode 100644 index 297e0b4..0000000 --- a/lib/mlibc/options/ansi/musl-generic-math/expm1f.c +++ /dev/null @@ -1,111 +0,0 @@ -/* origin: FreeBSD /usr/src/lib/msun/src/s_expm1f.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" - -static const float -o_threshold = 8.8721679688e+01, /* 0x42b17180 */ -ln2_hi = 6.9313812256e-01, /* 0x3f317180 */ -ln2_lo = 9.0580006145e-06, /* 0x3717f7d1 */ -invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */ -/* - * Domain [-0.34568, 0.34568], range ~[-6.694e-10, 6.696e-10]: - * |6 / x * (1 + 2 * (1 / (exp(x) - 1) - 1 / x)) - q(x)| < 2**-30.04 - * Scaled coefficients: Qn_here = 2**n * Qn_for_q (see s_expm1.c): - */ -Q1 = -3.3333212137e-2, /* -0x888868.0p-28 */ -Q2 = 1.5807170421e-3; /* 0xcf3010.0p-33 */ - -float expm1f(float x) -{ - float_t y,hi,lo,c,t,e,hxs,hfx,r1,twopk; - union {float f; uint32_t i;} u = {x}; - uint32_t hx = u.i & 0x7fffffff; - int k, sign = u.i >> 31; - - /* filter out huge and non-finite argument */ - if (hx >= 0x4195b844) { /* if |x|>=27*ln2 */ - if (hx > 0x7f800000) /* NaN */ - return x; - if (sign) - return -1; - if (x > o_threshold) { - x *= 0x1p127f; - return x; - } - } - - /* argument reduction */ - if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ - if (hx < 0x3F851592) { /* and |x| < 1.5 ln2 */ - if (!sign) { - hi = x - ln2_hi; - lo = ln2_lo; - k = 1; - } else { - hi = x + ln2_hi; - lo = -ln2_lo; - k = -1; - } - } else { - k = invln2*x + (sign ? -0.5f : 0.5f); - t = k; - hi = x - t*ln2_hi; /* t*ln2_hi is exact here */ - lo = t*ln2_lo; - } - x = hi-lo; - c = (hi-x)-lo; - } else if (hx < 0x33000000) { /* when |x|<2**-25, return x */ - if (hx < 0x00800000) - FORCE_EVAL(x*x); - return x; - } else - k = 0; - - /* x is now in primary range */ - hfx = 0.5f*x; - hxs = x*hfx; - r1 = 1.0f+hxs*(Q1+hxs*Q2); - t = 3.0f - r1*hfx; - e = hxs*((r1-t)/(6.0f - x*t)); - if (k == 0) /* c is 0 */ - return x - (x*e-hxs); - e = x*(e-c) - c; - e -= hxs; - /* exp(x) ~ 2^k (x_reduced - e + 1) */ - if (k == -1) - return 0.5f*(x-e) - 0.5f; - if (k == 1) { - if (x < -0.25f) - return -2.0f*(e-(x+0.5f)); - return 1.0f + 2.0f*(x-e); - } - u.i = (0x7f+k)<<23; /* 2^k */ - twopk = u.f; - if (k < 0 || k > 56) { /* suffice to return exp(x)-1 */ - y = x - e + 1.0f; - if (k == 128) - y = y*2.0f*0x1p127f; - else - y = y*twopk; - return y - 1.0f; - } - u.i = (0x7f-k)<<23; /* 2^-k */ - if (k < 23) - y = (x-e+(1-u.f))*twopk; - else - y = (x-(e+u.f)+1)*twopk; - return y; -} |