Mercurial > pmdwin

annotate fmgen/e_expf.c @ 3:8ad174416431

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Convert the last EUC-JP file to UTF-8.
author Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
date Tue, 21 May 2013 10:45:37 +0200
parents c55ea9478c80
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
1 /* e_expf.c -- float version of e_exp.c.
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
2 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
3 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
4
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
5 /*
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
6 * ====================================================
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
8 *
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
9 * Developed at SunPro, a Sun Microsystems, Inc. business.
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
10 * Permission to use, copy, modify, and distribute this
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
11 * software is freely granted, provided that this notice
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
12 * is preserved.
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
13 * ====================================================
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
14 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
15
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
16 #include <stdint.h>
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
17
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
18 static const float
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
19 one = 1.0,
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
20 halF[2] = {0.5,-0.5,},
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
21 ln2HI[2] ={ 6.9313812256e-01, /* 0x3f317180 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
22 -6.9313812256e-01,}, /* 0xbf317180 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
23 ln2LO[2] ={ 9.0580006145e-06, /* 0x3717f7d1 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
24 -9.0580006145e-06,}, /* 0xb717f7d1 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
25 invln2 = 1.4426950216e+00, /* 0x3fb8aa3b */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
26 P1 = 1.6666667163e-01, /* 0x3e2aaaab */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
27 P2 = -2.7777778450e-03, /* 0xbb360b61 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
28 P3 = 6.6137559770e-05, /* 0x388ab355 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
29 P4 = -1.6533901999e-06, /* 0xb5ddea0e */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
30 P5 = 4.1381369442e-08; /* 0x3331bb4c */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
31
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
32 typedef union {
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
33 float value;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
34 uint32_t word;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
35 } _f32;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
36
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
37 /* Get a 32 bit int from a float. */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
38 #define GET_FLOAT_WORD(i,d) \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
39 do { \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
40 _f32 gf_u; \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
41 gf_u.value = (d); \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
42 (i) = gf_u.word; \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
43 } while (0)
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
44
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
45 /* Set a float from a 32 bit int. */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
46 #define SET_FLOAT_WORD(d,i) \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
47 do { \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
48 _f32 sf_u; \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
49 sf_u.word = (i); \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
50 (d) = sf_u.value; \
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
51 } while (0)
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
52
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
53 float expf(float x) /* default IEEE double exp */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
54 {
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
55 float y,hi=0,lo=0,c,t;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
56 int32_t k,xsb;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
57 uint32_t hx;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
58
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
59 GET_FLOAT_WORD(hx,x);
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
60 xsb = (hx>>31)&1; /* sign bit of x */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
61 hx &= 0x7fffffff; /* high word of |x| */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
62
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
63 /* filter out non-finite argument */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
64 if(hx >= 0x42b17218) { /* if |x|>=88.721... */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
65 if(hx>0x7f800000)
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
66 return x+x; /* NaN */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
67 if(hx==0x7f800000)
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
68 return (xsb==0)? x:0.0f; /* exp(+-inf)={inf,0} */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
69 }
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
70
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
71 /* argument reduction */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
72 if(hx > 0x3eb17218) { /* if |x| > 0.5 ln2 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
73 if(hx < 0x3F851592) { /* and |x| < 1.5 ln2 */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
74 hi = x-ln2HI[xsb]; lo=ln2LO[xsb]; k = 1-xsb-xsb;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
75 } else {
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
76 k = invln2*x+halF[xsb];
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
77 t = k;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
78 hi = x - t*ln2HI[0]; /* t*ln2HI is exact here */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
79 lo = t*ln2LO[0];
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
80 }
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
81 x = hi - lo;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
82 }
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
83 else k = 0;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
84
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
85 /* x is now in primary range */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
86 t = x*x;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
87 c = x - t*(P1+t*(P2+t*(P3+t*(P4+t*P5))));
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
88 if(k==0) return one-((x*c)/(c-(float)2.0)-x);
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
89 else y = one-((lo-(x*c)/((float)2.0-c))-hi);
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
90 {
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
91 uint32_t hy;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
92 GET_FLOAT_WORD(hy,y);
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
93 SET_FLOAT_WORD(y,hy+(k<<23)); /* add k to y's exponent */
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
94 return y;
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
95 }
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
96 }
c55ea9478c80 Hello Gensokyo!
Emmanuel Gil Peyrot <linkmauve@linkmauve.fr>
parents:
diff changeset
97