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

Side by Side Diff: fusl/src/math/lgammaf_r.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/e_lgammaf_r.c */ 1 /* origin: FreeBSD /usr/src/lib/msun/src/e_lgammaf_r.c */
2 /* 2 /*
3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. 3 * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
4 */ 4 */
5 /* 5 /*
6 * ==================================================== 6 * ====================================================
7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 7 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 * 8 *
9 * Developed at SunPro, a Sun Microsystems, Inc. business. 9 * Developed at SunPro, a Sun Microsystems, Inc. business.
10 * Permission to use, copy, modify, and distribute this 10 * Permission to use, copy, modify, and distribute this
11 * software is freely granted, provided that this notice 11 * software is freely granted, provided that this notice
12 * is preserved. 12 * is preserved.
13 * ==================================================== 13 * ====================================================
14 */ 14 */
15 15
16 #include "libm.h" 16 #include "libm.h"
17 #include "libc.h" 17 #include "libc.h"
18 18
19 static const float 19 static const float pi = 3.1415927410e+00, /* 0x40490fdb */
20 pi = 3.1415927410e+00, /* 0x40490fdb */ 20 a0 = 7.7215664089e-02, /* 0x3d9e233f */
21 a0 = 7.7215664089e-02, /* 0x3d9e233f */ 21 a1 = 3.2246702909e-01, /* 0x3ea51a66 */
22 a1 = 3.2246702909e-01, /* 0x3ea51a66 */ 22 a2 = 6.7352302372e-02, /* 0x3d89f001 */
23 a2 = 6.7352302372e-02, /* 0x3d89f001 */ 23 a3 = 2.0580807701e-02, /* 0x3ca89915 */
24 a3 = 2.0580807701e-02, /* 0x3ca89915 */ 24 a4 = 7.3855509982e-03, /* 0x3bf2027e */
25 a4 = 7.3855509982e-03, /* 0x3bf2027e */ 25 a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */
26 a5 = 2.8905137442e-03, /* 0x3b3d6ec6 */ 26 a6 = 1.1927076848e-03, /* 0x3a9c54a1 */
27 a6 = 1.1927076848e-03, /* 0x3a9c54a1 */ 27 a7 = 5.1006977446e-04, /* 0x3a05b634 */
28 a7 = 5.1006977446e-04, /* 0x3a05b634 */ 28 a8 = 2.2086278477e-04, /* 0x39679767 */
29 a8 = 2.2086278477e-04, /* 0x39679767 */ 29 a9 = 1.0801156895e-04, /* 0x38e28445 */
30 a9 = 1.0801156895e-04, /* 0x38e28445 */ 30 a10 = 2.5214456400e-05, /* 0x37d383a2 */
31 a10 = 2.5214456400e-05, /* 0x37d383a2 */ 31 a11 = 4.4864096708e-05, /* 0x383c2c75 */
32 a11 = 4.4864096708e-05, /* 0x383c2c75 */ 32 tc = 1.4616321325e+00, /* 0x3fbb16c3 */
33 tc = 1.4616321325e+00, /* 0x3fbb16c3 */ 33 tf = -1.2148628384e-01, /* 0xbdf8cdcd */
34 tf = -1.2148628384e-01, /* 0xbdf8cdcd */ 34 /* tt = -(tail of tf) */
35 /* tt = -(tail of tf) */ 35 tt = 6.6971006518e-09, /* 0x31e61c52 */
36 tt = 6.6971006518e-09, /* 0x31e61c52 */ 36 t0 = 4.8383611441e-01, /* 0x3ef7b95e */
37 t0 = 4.8383611441e-01, /* 0x3ef7b95e */ 37 t1 = -1.4758771658e-01, /* 0xbe17213c */
38 t1 = -1.4758771658e-01, /* 0xbe17213c */ 38 t2 = 6.4624942839e-02, /* 0x3d845a15 */
39 t2 = 6.4624942839e-02, /* 0x3d845a15 */ 39 t3 = -3.2788541168e-02, /* 0xbd064d47 */
40 t3 = -3.2788541168e-02, /* 0xbd064d47 */ 40 t4 = 1.7970675603e-02, /* 0x3c93373d */
41 t4 = 1.7970675603e-02, /* 0x3c93373d */ 41 t5 = -1.0314224288e-02, /* 0xbc28fcfe */
42 t5 = -1.0314224288e-02, /* 0xbc28fcfe */ 42 t6 = 6.1005386524e-03, /* 0x3bc7e707 */
43 t6 = 6.1005386524e-03, /* 0x3bc7e707 */ 43 t7 = -3.6845202558e-03, /* 0xbb7177fe */
44 t7 = -3.6845202558e-03, /* 0xbb7177fe */ 44 t8 = 2.2596477065e-03, /* 0x3b141699 */
45 t8 = 2.2596477065e-03, /* 0x3b141699 */ 45 t9 = -1.4034647029e-03, /* 0xbab7f476 */
46 t9 = -1.4034647029e-03, /* 0xbab7f476 */ 46 t10 = 8.8108185446e-04, /* 0x3a66f867 */
47 t10 = 8.8108185446e-04, /* 0x3a66f867 */ 47 t11 = -5.3859531181e-04, /* 0xba0d3085 */
48 t11 = -5.3859531181e-04, /* 0xba0d3085 */ 48 t12 = 3.1563205994e-04, /* 0x39a57b6b */
49 t12 = 3.1563205994e-04, /* 0x39a57b6b */ 49 t13 = -3.1275415677e-04, /* 0xb9a3f927 */
50 t13 = -3.1275415677e-04, /* 0xb9a3f927 */ 50 t14 = 3.3552918467e-04, /* 0x39afe9f7 */
51 t14 = 3.3552918467e-04, /* 0x39afe9f7 */ 51 u0 = -7.7215664089e-02, /* 0xbd9e233f */
52 u0 = -7.7215664089e-02, /* 0xbd9e233f */ 52 u1 = 6.3282704353e-01, /* 0x3f2200f4 */
53 u1 = 6.3282704353e-01, /* 0x3f2200f4 */ 53 u2 = 1.4549225569e+00, /* 0x3fba3ae7 */
54 u2 = 1.4549225569e+00, /* 0x3fba3ae7 */ 54 u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */
55 u3 = 9.7771751881e-01, /* 0x3f7a4bb2 */ 55 u4 = 2.2896373272e-01, /* 0x3e6a7578 */
56 u4 = 2.2896373272e-01, /* 0x3e6a7578 */ 56 u5 = 1.3381091878e-02, /* 0x3c5b3c5e */
57 u5 = 1.3381091878e-02, /* 0x3c5b3c5e */ 57 v1 = 2.4559779167e+00, /* 0x401d2ebe */
58 v1 = 2.4559779167e+00, /* 0x401d2ebe */ 58 v2 = 2.1284897327e+00, /* 0x4008392d */
59 v2 = 2.1284897327e+00, /* 0x4008392d */ 59 v3 = 7.6928514242e-01, /* 0x3f44efdf */
60 v3 = 7.6928514242e-01, /* 0x3f44efdf */ 60 v4 = 1.0422264785e-01, /* 0x3dd572af */
61 v4 = 1.0422264785e-01, /* 0x3dd572af */ 61 v5 = 3.2170924824e-03, /* 0x3b52d5db */
62 v5 = 3.2170924824e-03, /* 0x3b52d5db */ 62 s0 = -7.7215664089e-02, /* 0xbd9e233f */
63 s0 = -7.7215664089e-02, /* 0xbd9e233f */ 63 s1 = 2.1498242021e-01, /* 0x3e5c245a */
64 s1 = 2.1498242021e-01, /* 0x3e5c245a */ 64 s2 = 3.2577878237e-01, /* 0x3ea6cc7a */
65 s2 = 3.2577878237e-01, /* 0x3ea6cc7a */ 65 s3 = 1.4635047317e-01, /* 0x3e15dce6 */
66 s3 = 1.4635047317e-01, /* 0x3e15dce6 */ 66 s4 = 2.6642270386e-02, /* 0x3cda40e4 */
67 s4 = 2.6642270386e-02, /* 0x3cda40e4 */ 67 s5 = 1.8402845599e-03, /* 0x3af135b4 */
68 s5 = 1.8402845599e-03, /* 0x3af135b4 */ 68 s6 = 3.1947532989e-05, /* 0x3805ff67 */
69 s6 = 3.1947532989e-05, /* 0x3805ff67 */ 69 r1 = 1.3920053244e+00, /* 0x3fb22d3b */
70 r1 = 1.3920053244e+00, /* 0x3fb22d3b */ 70 r2 = 7.2193557024e-01, /* 0x3f38d0c5 */
71 r2 = 7.2193557024e-01, /* 0x3f38d0c5 */ 71 r3 = 1.7193385959e-01, /* 0x3e300f6e */
72 r3 = 1.7193385959e-01, /* 0x3e300f6e */ 72 r4 = 1.8645919859e-02, /* 0x3c98bf54 */
73 r4 = 1.8645919859e-02, /* 0x3c98bf54 */ 73 r5 = 7.7794247773e-04, /* 0x3a4beed6 */
74 r5 = 7.7794247773e-04, /* 0x3a4beed6 */ 74 r6 = 7.3266842264e-06, /* 0x36f5d7bd */
75 r6 = 7.3266842264e-06, /* 0x36f5d7bd */ 75 w0 = 4.1893854737e-01, /* 0x3ed67f1d */
76 w0 = 4.1893854737e-01, /* 0x3ed67f1d */ 76 w1 = 8.3333335817e-02, /* 0x3daaaaab */
77 w1 = 8.3333335817e-02, /* 0x3daaaaab */ 77 w2 = -2.7777778450e-03, /* 0xbb360b61 */
78 w2 = -2.7777778450e-03, /* 0xbb360b61 */ 78 w3 = 7.9365057172e-04, /* 0x3a500cfd */
79 w3 = 7.9365057172e-04, /* 0x3a500cfd */ 79 w4 = -5.9518753551e-04, /* 0xba1c065c */
80 w4 = -5.9518753551e-04, /* 0xba1c065c */ 80 w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */
81 w5 = 8.3633989561e-04, /* 0x3a5b3dd2 */ 81 w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
82 w6 = -1.6309292987e-03; /* 0xbad5c4e8 */
83 82
84 /* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */ 83 /* sin(pi*x) assuming x > 2^-100, if sin(pi*x)==0 the sign is arbitrary */
85 static float sin_pi(float x) 84 static float sin_pi(float x) {
86 { 85 double_t y;
87 » double_t y; 86 int n;
88 » int n; 87
89 88 /* spurious inexact if odd int */
90 » /* spurious inexact if odd int */ 89 x = 2 * (x * 0.5f - floorf(x * 0.5f)); /* x mod 2.0 */
91 » x = 2*(x*0.5f - floorf(x*0.5f)); /* x mod 2.0 */ 90
92 91 n = (int)(x * 4);
93 » n = (int)(x*4); 92 n = (n + 1) / 2;
94 » n = (n+1)/2; 93 y = x - n * 0.5f;
95 » y = x - n*0.5f; 94 y *= 3.14159265358979323846;
96 » y *= 3.14159265358979323846; 95 switch (n) {
97 » switch (n) { 96 default: /* case 4: */
98 » default: /* case 4: */ 97 case 0:
99 » case 0: return __sindf(y); 98 return __sindf(y);
100 » case 1: return __cosdf(y); 99 case 1:
101 » case 2: return __sindf(-y); 100 return __cosdf(y);
102 » case 3: return -__cosdf(y); 101 case 2:
103 » } 102 return __sindf(-y);
103 case 3:
104 return -__cosdf(y);
105 }
104 } 106 }
105 107
106 float __lgammaf_r(float x, int *signgamp) 108 float __lgammaf_r(float x, int* signgamp) {
107 { 109 union {
108 » union {float f; uint32_t i;} u = {x}; 110 float f;
109 » float t,y,z,nadj,p,p1,p2,p3,q,r,w; 111 uint32_t i;
110 » uint32_t ix; 112 } u = {x};
111 » int i,sign; 113 float t, y, z, nadj, p, p1, p2, p3, q, r, w;
112 114 uint32_t ix;
113 » /* purge off +-inf, NaN, +-0, tiny and negative arguments */ 115 int i, sign;
114 » *signgamp = 1; 116
115 » sign = u.i>>31; 117 /* purge off +-inf, NaN, +-0, tiny and negative arguments */
116 » ix = u.i & 0x7fffffff; 118 *signgamp = 1;
117 » if (ix >= 0x7f800000) 119 sign = u.i >> 31;
118 » » return x*x; 120 ix = u.i & 0x7fffffff;
119 » if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */ 121 if (ix >= 0x7f800000)
120 » » if (sign) { 122 return x * x;
121 » » » *signgamp = -1; 123 if (ix < 0x35000000) { /* |x| < 2**-21, return -log(|x|) */
122 » » » x = -x; 124 if (sign) {
123 » » } 125 *signgamp = -1;
124 » » return -logf(x); 126 x = -x;
125 » } 127 }
126 » if (sign) { 128 return -logf(x);
127 » » x = -x; 129 }
128 » » t = sin_pi(x); 130 if (sign) {
129 » » if (t == 0.0f) /* -integer */ 131 x = -x;
130 » » » return 1.0f/(x-x); 132 t = sin_pi(x);
131 » » if (t > 0.0f) 133 if (t == 0.0f) /* -integer */
132 » » » *signgamp = -1; 134 return 1.0f / (x - x);
133 » » else 135 if (t > 0.0f)
134 » » » t = -t; 136 *signgamp = -1;
135 » » nadj = logf(pi/(t*x)); 137 else
136 » } 138 t = -t;
137 139 nadj = logf(pi / (t * x));
138 » /* purge off 1 and 2 */ 140 }
139 » if (ix == 0x3f800000 || ix == 0x40000000) 141
140 » » r = 0; 142 /* purge off 1 and 2 */
141 » /* for x < 2.0 */ 143 if (ix == 0x3f800000 || ix == 0x40000000)
142 » else if (ix < 0x40000000) { 144 r = 0;
143 » » if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */ 145 /* for x < 2.0 */
144 » » » r = -logf(x); 146 else if (ix < 0x40000000) {
145 » » » if (ix >= 0x3f3b4a20) { 147 if (ix <= 0x3f666666) { /* lgamma(x) = lgamma(x+1)-log(x) */
146 » » » » y = 1.0f - x; 148 r = -logf(x);
147 » » » » i = 0; 149 if (ix >= 0x3f3b4a20) {
148 » » » } else if (ix >= 0x3e6d3308) { 150 y = 1.0f - x;
149 » » » » y = x - (tc-1.0f); 151 i = 0;
150 » » » » i = 1; 152 } else if (ix >= 0x3e6d3308) {
151 » » » } else { 153 y = x - (tc - 1.0f);
152 » » » » y = x; 154 i = 1;
153 » » » » i = 2; 155 } else {
154 » » » } 156 y = x;
155 » » } else { 157 i = 2;
156 » » » r = 0.0f; 158 }
157 » » » if (ix >= 0x3fdda618) { /* [1.7316,2] */ 159 } else {
158 » » » » y = 2.0f - x; 160 r = 0.0f;
159 » » » » i = 0; 161 if (ix >= 0x3fdda618) { /* [1.7316,2] */
160 » » » } else if (ix >= 0x3F9da620) { /* [1.23,1.73] */ 162 y = 2.0f - x;
161 » » » » y = x - tc; 163 i = 0;
162 » » » » i = 1; 164 } else if (ix >= 0x3F9da620) { /* [1.23,1.73] */
163 » » » } else { 165 y = x - tc;
164 » » » » y = x - 1.0f; 166 i = 1;
165 » » » » i = 2; 167 } else {
166 » » » } 168 y = x - 1.0f;
167 » » } 169 i = 2;
168 » » switch(i) { 170 }
169 » » case 0: 171 }
170 » » » z = y*y; 172 switch (i) {
171 » » » p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); 173 case 0:
172 » » » p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); 174 z = y * y;
173 » » » p = y*p1+p2; 175 p1 = a0 + z * (a2 + z * (a4 + z * (a6 + z * (a8 + z * a10))));
174 » » » r += p - 0.5f*y; 176 p2 = z * (a1 + z * (a3 + z * (a5 + z * (a7 + z * (a9 + z * a11)))));
175 » » » break; 177 p = y * p1 + p2;
176 » » case 1: 178 r += p - 0.5f * y;
177 » » » z = y*y; 179 break;
178 » » » w = z*y; 180 case 1:
179 » » » p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ 181 z = y * y;
180 » » » p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); 182 w = z * y;
181 » » » p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); 183 p1 = t0 + w * (t3 + w * (t6 + w * (t9 + w * t12))); /* parallel comp */
182 » » » p = z*p1-(tt-w*(p2+y*p3)); 184 p2 = t1 + w * (t4 + w * (t7 + w * (t10 + w * t13)));
183 » » » r += (tf + p); 185 p3 = t2 + w * (t5 + w * (t8 + w * (t11 + w * t14)));
184 » » » break; 186 p = z * p1 - (tt - w * (p2 + y * p3));
185 » » case 2: 187 r += (tf + p);
186 » » » p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); 188 break;
187 » » » p2 = 1.0f+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); 189 case 2:
188 » » » r += -0.5f*y + p1/p2; 190 p1 = y * (u0 + y * (u1 + y * (u2 + y * (u3 + y * (u4 + y * u5)))));
189 » » } 191 p2 = 1.0f + y * (v1 + y * (v2 + y * (v3 + y * (v4 + y * v5))));
190 » } else if (ix < 0x41000000) { /* x < 8.0 */ 192 r += -0.5f * y + p1 / p2;
191 » » i = (int)x; 193 }
192 » » y = x - (float)i; 194 } else if (ix < 0x41000000) { /* x < 8.0 */
193 » » p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); 195 i = (int)x;
194 » » q = 1.0f+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); 196 y = x - (float)i;
195 » » r = 0.5f*y+p/q; 197 p = y *
196 » » z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */ 198 (s0 + y * (s1 + y * (s2 + y * (s3 + y * (s4 + y * (s5 + y * s6))))));
197 » » switch (i) { 199 q = 1.0f + y * (r1 + y * (r2 + y * (r3 + y * (r4 + y * (r5 + y * r6)))));
198 » » case 7: z *= y + 6.0f; /* FALLTHRU */ 200 r = 0.5f * y + p / q;
199 » » case 6: z *= y + 5.0f; /* FALLTHRU */ 201 z = 1.0f; /* lgamma(1+s) = log(s) + lgamma(s) */
200 » » case 5: z *= y + 4.0f; /* FALLTHRU */ 202 switch (i) {
201 » » case 4: z *= y + 3.0f; /* FALLTHRU */ 203 case 7:
202 » » case 3: z *= y + 2.0f; /* FALLTHRU */ 204 z *= y + 6.0f; /* FALLTHRU */
203 » » » r += logf(z); 205 case 6:
204 » » » break; 206 z *= y + 5.0f; /* FALLTHRU */
205 » » } 207 case 5:
206 » } else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */ 208 z *= y + 4.0f; /* FALLTHRU */
207 » » t = logf(x); 209 case 4:
208 » » z = 1.0f/x; 210 z *= y + 3.0f; /* FALLTHRU */
209 » » y = z*z; 211 case 3:
210 » » w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); 212 z *= y + 2.0f; /* FALLTHRU */
211 » » r = (x-0.5f)*(t-1.0f)+w; 213 r += logf(z);
212 » } else /* 2**58 <= x <= inf */ 214 break;
213 » » r = x*(logf(x)-1.0f); 215 }
214 » if (sign) 216 } else if (ix < 0x5c800000) { /* 8.0 <= x < 2**58 */
215 » » r = nadj - r; 217 t = logf(x);
216 » return r; 218 z = 1.0f / x;
219 y = z * z;
220 w = w0 + z * (w1 + y * (w2 + y * (w3 + y * (w4 + y * (w5 + y * w6)))));
221 r = (x - 0.5f) * (t - 1.0f) + w;
222 } else /* 2**58 <= x <= inf */
223 r = x * (logf(x) - 1.0f);
224 if (sign)
225 r = nadj - r;
226 return r;
217 } 227 }
218 228
219 weak_alias(__lgammaf_r, lgammaf_r); 229 weak_alias(__lgammaf_r, lgammaf_r);
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698