Chromium Code Reviews
chromiumcodereview-hr@appspot.gserviceaccount.com (chromiumcodereview-hr) | Please choose your nickname with Settings | Help | Chromium Project | Gerrit Changes | Sign out
(84)

Side by Side Diff: fusl/src/math/exp2f.c

Issue 1714623002: [fusl] clang-format fusl (Closed) Base URL: git@github.com:domokit/mojo.git@master
Patch Set: headers too Created 4 years, 10 months ago
Use n/p to move between diff chunks; N/P to move between comments. Draft comments are only viewable by you.
Jump to:
View unified diff | Download patch
OLDNEW
1 /* origin: FreeBSD /usr/src/lib/msun/src/s_exp2f.c */ 1 /* origin: FreeBSD /usr/src/lib/msun/src/s_exp2f.c */
2 /*- 2 /*-
3 * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG> 3 * Copyright (c) 2005 David Schultz <das@FreeBSD.ORG>
4 * All rights reserved. 4 * All rights reserved.
5 * 5 *
6 * Redistribution and use in source and binary forms, with or without 6 * Redistribution and use in source and binary forms, with or without
7 * modification, are permitted provided that the following conditions 7 * modification, are permitted provided that the following conditions
8 * are met: 8 * are met:
9 * 1. Redistributions of source code must retain the above copyright 9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer. 10 * notice, this list of conditions and the following disclaimer.
(...skipping 11 matching lines...) Expand all
22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT 22 * HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY 23 * LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 24 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
25 * SUCH DAMAGE. 25 * SUCH DAMAGE.
26 */ 26 */
27 27
28 #include "libm.h" 28 #include "libm.h"
29 29
30 #define TBLSIZE 16 30 #define TBLSIZE 16
31 31
32 static const float 32 static const float redux = 0x1.8p23f / TBLSIZE, P1 = 0x1.62e430p-1f,
33 redux = 0x1.8p23f / TBLSIZE, 33 P2 = 0x1.ebfbe0p-3f, P3 = 0x1.c6b348p-5f,
34 P1 = 0x1.62e430p-1f, 34 P4 = 0x1.3b2c9cp-7f;
35 P2 = 0x1.ebfbe0p-3f,
36 P3 = 0x1.c6b348p-5f,
37 P4 = 0x1.3b2c9cp-7f;
38 35
39 static const double exp2ft[TBLSIZE] = { 36 static const double exp2ft[TBLSIZE] = {
40 0x1.6a09e667f3bcdp-1, 37 0x1.6a09e667f3bcdp-1, 0x1.7a11473eb0187p-1, 0x1.8ace5422aa0dbp-1,
41 0x1.7a11473eb0187p-1, 38 0x1.9c49182a3f090p-1, 0x1.ae89f995ad3adp-1, 0x1.c199bdd85529cp-1,
42 0x1.8ace5422aa0dbp-1, 39 0x1.d5818dcfba487p-1, 0x1.ea4afa2a490dap-1, 0x1.0000000000000p+0,
43 0x1.9c49182a3f090p-1, 40 0x1.0b5586cf9890fp+0, 0x1.172b83c7d517bp+0, 0x1.2387a6e756238p+0,
44 0x1.ae89f995ad3adp-1, 41 0x1.306fe0a31b715p+0, 0x1.3dea64c123422p+0, 0x1.4bfdad5362a27p+0,
45 0x1.c199bdd85529cp-1, 42 0x1.5ab07dd485429p+0,
46 0x1.d5818dcfba487p-1,
47 0x1.ea4afa2a490dap-1,
48 0x1.0000000000000p+0,
49 0x1.0b5586cf9890fp+0,
50 0x1.172b83c7d517bp+0,
51 0x1.2387a6e756238p+0,
52 0x1.306fe0a31b715p+0,
53 0x1.3dea64c123422p+0,
54 0x1.4bfdad5362a27p+0,
55 0x1.5ab07dd485429p+0,
56 }; 43 };
57 44
58 /* 45 /*
59 * exp2f(x): compute the base 2 exponential of x 46 * exp2f(x): compute the base 2 exponential of x
60 * 47 *
61 * Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927. 48 * Accuracy: Peak error < 0.501 ulp; location of peak: -0.030110927.
62 * 49 *
63 * Method: (equally-spaced tables) 50 * Method: (equally-spaced tables)
64 * 51 *
65 * Reduce x: 52 * Reduce x:
66 * x = k + y, for integer k and |y| <= 1/2. 53 * x = k + y, for integer k and |y| <= 1/2.
67 * Thus we have exp2f(x) = 2**k * exp2(y). 54 * Thus we have exp2f(x) = 2**k * exp2(y).
68 * 55 *
69 * Reduce y: 56 * Reduce y:
70 * y = i/TBLSIZE + z for integer i near y * TBLSIZE. 57 * y = i/TBLSIZE + z for integer i near y * TBLSIZE.
71 * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z), 58 * Thus we have exp2(y) = exp2(i/TBLSIZE) * exp2(z),
72 * with |z| <= 2**-(TBLSIZE+1). 59 * with |z| <= 2**-(TBLSIZE+1).
73 * 60 *
74 * We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a 61 * We compute exp2(i/TBLSIZE) via table lookup and exp2(z) via a
75 * degree-4 minimax polynomial with maximum error under 1.4 * 2**-33. 62 * degree-4 minimax polynomial with maximum error under 1.4 * 2**-33.
76 * Using double precision for everything except the reduction makes 63 * Using double precision for everything except the reduction makes
77 * roundoff error insignificant and simplifies the scaling step. 64 * roundoff error insignificant and simplifies the scaling step.
78 * 65 *
79 * This method is due to Tang, but I do not use his suggested parameters: 66 * This method is due to Tang, but I do not use his suggested parameters:
80 * 67 *
81 * Tang, P. Table-driven Implementation of the Exponential Function 68 * Tang, P. Table-driven Implementation of the Exponential Function
82 * in IEEE Floating-Point Arithmetic. TOMS 15(2), 144-157 (1989). 69 * in IEEE Floating-Point Arithmetic. TOMS 15(2), 144-157 (1989).
83 */ 70 */
84 float exp2f(float x) 71 float exp2f(float x) {
85 { 72 double_t t, r, z;
86 » double_t t, r, z; 73 union {
87 » union {float f; uint32_t i;} u = {x}; 74 float f;
88 » union {double f; uint64_t i;} uk; 75 uint32_t i;
89 » uint32_t ix, i0, k; 76 } u = {x};
77 union {
78 double f;
79 uint64_t i;
80 } uk;
81 uint32_t ix, i0, k;
90 82
91 » /* Filter out exceptional cases. */ 83 /* Filter out exceptional cases. */
92 » ix = u.i & 0x7fffffff; 84 ix = u.i & 0x7fffffff;
93 » if (ix > 0x42fc0000) { /* |x| > 126 */ 85 if (ix > 0x42fc0000) { /* |x| > 126 */
94 » » if (u.i >= 0x43000000 && u.i < 0x80000000) { /* x >= 128 */ 86 if (u.i >= 0x43000000 && u.i < 0x80000000) { /* x >= 128 */
95 » » » x *= 0x1p127f; 87 x *= 0x1p127f;
96 » » » return x; 88 return x;
97 » » } 89 }
98 » » if (u.i >= 0x80000000) { /* x < -126 */ 90 if (u.i >= 0x80000000) { /* x < -126 */
99 » » » if (u.i >= 0xc3160000 || (u.i & 0x0000ffff)) 91 if (u.i >= 0xc3160000 || (u.i & 0x0000ffff))
100 » » » » FORCE_EVAL(-0x1p-149f/x); 92 FORCE_EVAL(-0x1p-149f / x);
101 » » » if (u.i >= 0xc3160000) /* x <= -150 */ 93 if (u.i >= 0xc3160000) /* x <= -150 */
102 » » » » return 0; 94 return 0;
103 » » } 95 }
104 » } else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */ 96 } else if (ix <= 0x33000000) { /* |x| <= 0x1p-25 */
105 » » return 1.0f + x; 97 return 1.0f + x;
106 » } 98 }
107 99
108 » /* Reduce x, computing z, i0, and k. */ 100 /* Reduce x, computing z, i0, and k. */
109 » u.f = x + redux; 101 u.f = x + redux;
110 » i0 = u.i; 102 i0 = u.i;
111 » i0 += TBLSIZE / 2; 103 i0 += TBLSIZE / 2;
112 » k = i0 / TBLSIZE; 104 k = i0 / TBLSIZE;
113 » uk.i = (uint64_t)(0x3ff + k)<<52; 105 uk.i = (uint64_t)(0x3ff + k) << 52;
114 » i0 &= TBLSIZE - 1; 106 i0 &= TBLSIZE - 1;
115 » u.f -= redux; 107 u.f -= redux;
116 » z = x - u.f; 108 z = x - u.f;
117 » /* Compute r = exp2(y) = exp2ft[i0] * p(z). */ 109 /* Compute r = exp2(y) = exp2ft[i0] * p(z). */
118 » r = exp2ft[i0]; 110 r = exp2ft[i0];
119 » t = r * z; 111 t = r * z;
120 » r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4); 112 r = r + t * (P1 + z * P2) + t * (z * z) * (P3 + z * P4);
121 113
122 » /* Scale by 2**k */ 114 /* Scale by 2**k */
123 » return r * uk.f; 115 return r * uk.f;
124 } 116 }
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698