Annotation of src/lib/libm/src/e_expf.c, Revision 1.9
1.1 jtc 1: /* e_expf.c -- float version of e_exp.c.
2: * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3: */
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
1.8 simonb 11: * software is freely granted, provided that this notice
1.1 jtc 12: * is preserved.
13: * ====================================================
14: */
15:
1.7 lukem 16: #include <sys/cdefs.h>
1.3 jtc 17: #if defined(LIBM_SCCS) && !defined(lint)
1.9 ! wiz 18: __RCSID("$NetBSD: e_expf.c,v 1.8 1999/07/02 15:37:39 simonb Exp $");
1.1 jtc 19: #endif
20:
21: #include "math.h"
22: #include "math_private.h"
23:
1.6 phil 24: static const float huge = 1.0e+30;
1.4 jtc 25:
1.1 jtc 26: static const float
27: one = 1.0,
28: halF[2] = {0.5,-0.5,},
29: twom100 = 7.8886090522e-31, /* 2**-100=0x0d800000 */
30: o_threshold= 8.8721679688e+01, /* 0x42b17180 */
31: u_threshold= -1.0397208405e+02, /* 0xc2cff1b5 */
32: ln2HI[2] ={ 6.9313812256e-01, /* 0x3f317180 */
33: -6.9313812256e-01,}, /* 0xbf317180 */
34: ln2LO[2] ={ 9.0580006145e-06, /* 0x3717f7d1 */
35: -9.0580006145e-06,}, /* 0xb717f7d1 */
36: invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
37: P1 = 1.6666667163e-01, /* 0x3e2aaaab */
38: P2 = -2.7777778450e-03, /* 0xbb360b61 */
39: P3 = 6.6137559770e-05, /* 0x388ab355 */
40: P4 = -1.6533901999e-06, /* 0xb5ddea0e */
41: P5 = 4.1381369442e-08; /* 0x3331bb4c */
42:
1.9 ! wiz 43: float
! 44: __ieee754_expf(float x) /* default IEEE double exp */
1.1 jtc 45: {
46: float y,hi,lo,c,t;
1.2 jtc 47: int32_t k,xsb;
48: u_int32_t hx;
1.1 jtc 49:
1.7 lukem 50: hi = lo = 0;
51: k = 0;
1.1 jtc 52: GET_FLOAT_WORD(hx,x);
53: xsb = (hx>>31)&1; /* sign bit of x */
54: hx &= 0x7fffffff; /* high word of |x| */
55:
56: /* filter out non-finite argument */
57: if(hx >= 0x42b17218) { /* if |x|>=88.721... */
58: if(hx>0x7f800000)
59: return x+x; /* NaN */
60: if(hx==0x7f800000)
61: return (xsb==0)? x:0.0; /* exp(+-inf)={inf,0} */
62: if(x > o_threshold) return huge*huge; /* overflow */
63: if(x < u_threshold) return twom100*twom100; /* underflow */
64: }
65:
66: /* argument reduction */
1.8 simonb 67: if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
1.1 jtc 68: if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
69: hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
70: } else {
71: k = invln2*x+halF[xsb];
72: t = k;
73: hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
74: lo = t*ln2LO[0];
75: }
76: x = hi - lo;
1.8 simonb 77: }
1.1 jtc 78: else if(hx < 0x31800000) { /* when |x|<2**-28 */
79: if(huge+x>one) return one+x;/* trigger inexact */
80: }
81: else k = 0;
82:
83: /* x is now in primary range */
84: t = x*x;
85: c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
1.8 simonb 86: if(k==0) return one-((x*c)/(c-(float)2.0)-x);
1.1 jtc 87: else y = one-((lo-(x*c)/((float)2.0-c))-hi);
88: if(k >= -125) {
1.2 jtc 89: u_int32_t hy;
1.1 jtc 90: GET_FLOAT_WORD(hy,y);
91: SET_FLOAT_WORD(y,hy+(k<<23)); /* add k to y's exponent */
92: return y;
93: } else {
1.2 jtc 94: u_int32_t hy;
1.1 jtc 95: GET_FLOAT_WORD(hy,y);
96: SET_FLOAT_WORD(y,hy+((k+100)<<23)); /* add k to y's exponent */
97: return y*twom100;
98: }
99: }
CVSweb <webmaster@jp.NetBSD.org>