1/* origin: FreeBSD /usr/src/lib/msun/src/e_expf.c */ 2/* 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 4 */ 5/* 6 * ==================================================== 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 8 * 9 * Developed at SunPro, a Sun Microsystems, Inc. business. 10 * Permission to use, copy, modify, and distribute this 11 * software is freely granted, provided that this notice 12 * is preserved. 13 * ==================================================== 14 */ 15 16#include "libm.h" 17 18static const float 19half[2] = {0.5,-0.5}, 20ln2hi = 6.9314575195e-1f, /* 0x3f317200 */ 21ln2lo = 1.4286067653e-6f, /* 0x35bfbe8e */ 22invln2 = 1.4426950216e+0f, /* 0x3fb8aa3b */ 23/* 24 * Domain [-0.34568, 0.34568], range ~[-4.278e-9, 4.447e-9]: 25 * |x*(exp(x)+1)/(exp(x)-1) - p(x)| < 2**-27.74 26 */ 27P1 = 1.6666625440e-1f, /* 0xaaaa8f.0p-26 */ 28P2 = -2.7667332906e-3f; /* -0xb55215.0p-32 */ 29 30float expf(float x) 31{ 32 float_t hi, lo, c, xx, y; 33 int k, sign; 34 uint32_t hx; 35 36 GET_FLOAT_WORD(hx, x); 37 sign = hx >> 31; /* sign bit of x */ 38 hx &= 0x7fffffff; /* high word of |x| */ 39 40 /* special cases */ 41 if (hx >= 0x42aeac50) { /* if |x| >= -87.33655f or NaN */ 42 if (hx > 0x7f800000) /* NaN */ 43 return x; 44 if (hx >= 0x42b17218 && !sign) { /* x >= 88.722839f */ 45 /* overflow */ 46 x *= 0x1p127f; 47 return x; 48 } 49 if (sign) { 50 /* underflow */ 51 FORCE_EVAL(-0x1p-149f/x); 52 if (hx >= 0x42cff1b5) /* x <= -103.972084f */ 53 return 0; 54 } 55 } 56 57 /* argument reduction */ 58 if (hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */ 59 if (hx > 0x3f851592) /* if |x| > 1.5 ln2 */ 60 k = invln2*x + half[sign]; 61 else 62 k = 1 - sign - sign; 63 hi = x - k*ln2hi; /* k*ln2hi is exact here */ 64 lo = k*ln2lo; 65 x = hi - lo; 66 } else if (hx > 0x39000000) { /* |x| > 2**-14 */ 67 k = 0; 68 hi = x; 69 lo = 0; 70 } else { 71 /* raise inexact */ 72 FORCE_EVAL(0x1p127f + x); 73 return 1 + x; 74 } 75 76 /* x is now in primary range */ 77 xx = x*x; 78 c = x - xx*(P1+xx*P2); 79 y = 1 + (x*c/(2-c) - lo + hi); 80 if (k == 0) 81 return y; 82 return scalbnf(y, k); 83} 84