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

Side by Side Diff: fusl/src/math/tgammal.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: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.c */ 1 /* origin: OpenBSD /usr/src/lib/libm/src/ld80/e_tgammal.c */
2 /* 2 /*
3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net> 3 * Copyright (c) 2008 Stephen L. Moshier <steve@moshier.net>
4 * 4 *
5 * Permission to use, copy, modify, and distribute this software for any 5 * Permission to use, copy, modify, and distribute this software for any
6 * purpose with or without fee is hereby granted, provided that the above 6 * purpose with or without fee is hereby granted, provided that the above
7 * copyright notice and this permission notice appear in all copies. 7 * copyright notice and this permission notice appear in all copies.
8 * 8 *
9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES 9 * THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES
10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF 10 * WITH REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF
(...skipping 33 matching lines...) Expand 10 before | Expand all | Expand 10 after
44 * IEEE -40,+40 10000 3.6e-19 7.9e-20 44 * IEEE -40,+40 10000 3.6e-19 7.9e-20
45 * IEEE -1755,+1755 10000 4.8e-18 6.5e-19 45 * IEEE -1755,+1755 10000 4.8e-18 6.5e-19
46 * 46 *
47 * Accuracy for large arguments is dominated by error in powl(). 47 * Accuracy for large arguments is dominated by error in powl().
48 * 48 *
49 */ 49 */
50 50
51 #include "libm.h" 51 #include "libm.h"
52 52
53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 53 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
54 long double tgammal(long double x) 54 long double tgammal(long double x) {
55 { 55 return tgamma(x);
56 » return tgamma(x);
57 } 56 }
58 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384 57 #elif LDBL_MANT_DIG == 64 && LDBL_MAX_EXP == 16384
59 /* 58 /*
60 tgamma(x+2) = tgamma(x+2) P(x)/Q(x) 59 tgamma(x+2) = tgamma(x+2) P(x)/Q(x)
61 0 <= x <= 1 60 0 <= x <= 1
62 Relative error 61 Relative error
63 n=7, d=8 62 n=7, d=8
64 Peak error = 1.83e-20 63 Peak error = 1.83e-20
65 Relative error spread = 8.4e-23 64 Relative error spread = 8.4e-23
66 */ 65 */
67 static const long double P[8] = { 66 static const long double P[8] = {
68 4.212760487471622013093E-5L, 67 4.212760487471622013093E-5L, 4.542931960608009155600E-4L,
69 4.542931960608009155600E-4L, 68 4.092666828394035500949E-3L, 2.385363243461108252554E-2L,
70 4.092666828394035500949E-3L, 69 1.113062816019361559013E-1L, 3.629515436640239168939E-1L,
71 2.385363243461108252554E-2L, 70 8.378004301573126728826E-1L, 1.000000000000000000009E0L,
72 1.113062816019361559013E-1L,
73 3.629515436640239168939E-1L,
74 8.378004301573126728826E-1L,
75 1.000000000000000000009E0L,
76 }; 71 };
77 static const long double Q[9] = { 72 static const long double Q[9] = {
78 -1.397148517476170440917E-5L, 73 -1.397148517476170440917E-5L, 2.346584059160635244282E-4L,
79 2.346584059160635244282E-4L, 74 -1.237799246653152231188E-3L, -7.955933682494738320586E-4L,
80 -1.237799246653152231188E-3L, 75 2.773706565840072979165E-2L, -4.633887671244534213831E-2L,
81 -7.955933682494738320586E-4L, 76 -2.243510905670329164562E-1L, 4.150160950588455434583E-1L,
82 2.773706565840072979165E-2L, 77 9.999999999999999999908E-1L,
83 -4.633887671244534213831E-2L,
84 -2.243510905670329164562E-1L,
85 4.150160950588455434583E-1L,
86 9.999999999999999999908E-1L,
87 }; 78 };
88 79
89 /* 80 /*
90 static const long double P[] = { 81 static const long double P[] = {
91 -3.01525602666895735709e0L, 82 -3.01525602666895735709e0L,
92 -3.25157411956062339893e1L, 83 -3.25157411956062339893e1L,
93 -2.92929976820724030353e2L, 84 -2.92929976820724030353e2L,
94 -1.70730828800510297666e3L, 85 -1.70730828800510297666e3L,
95 -7.96667499622741999770e3L, 86 -7.96667499622741999770e3L,
96 -2.59780216007146401957e4L, 87 -2.59780216007146401957e4L,
(...skipping 18 matching lines...) Expand all
115 /* Stirling's formula for the gamma function 106 /* Stirling's formula for the gamma function
116 tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x)) 107 tgamma(x) = sqrt(2 pi) x^(x-.5) exp(-x) (1 + 1/x P(1/x))
117 z(x) = x 108 z(x) = x
118 13 <= x <= 1024 109 13 <= x <= 1024
119 Relative error 110 Relative error
120 n=8, d=0 111 n=8, d=0
121 Peak error = 9.44e-21 112 Peak error = 9.44e-21
122 Relative error spread = 8.8e-4 113 Relative error spread = 8.8e-4
123 */ 114 */
124 static const long double STIR[9] = { 115 static const long double STIR[9] = {
125 7.147391378143610789273E-4L, 116 7.147391378143610789273E-4L, -2.363848809501759061727E-5L,
126 -2.363848809501759061727E-5L, 117 -5.950237554056330156018E-4L, 6.989332260623193171870E-5L,
127 -5.950237554056330156018E-4L, 118 7.840334842744753003862E-4L, -2.294719747873185405699E-4L,
128 6.989332260623193171870E-5L, 119 -2.681327161876304418288E-3L, 3.472222222230075327854E-3L,
129 7.840334842744753003862E-4L, 120 8.333333333333331800504E-2L,
130 -2.294719747873185405699E-4L,
131 -2.681327161876304418288E-3L,
132 3.472222222230075327854E-3L,
133 8.333333333333331800504E-2L,
134 }; 121 };
135 122
136 #define MAXSTIR 1024.0L 123 #define MAXSTIR 1024.0L
137 static const long double SQTPI = 2.50662827463100050242E0L; 124 static const long double SQTPI = 2.50662827463100050242E0L;
138 125
139 /* 1/tgamma(x) = z P(z) 126 /* 1/tgamma(x) = z P(z)
140 * z(x) = 1/x 127 * z(x) = 1/x
141 * 0 < x < 0.03125 128 * 0 < x < 0.03125
142 * Peak relative error 4.2e-23 129 * Peak relative error 4.2e-23
143 */ 130 */
144 static const long double S[9] = { 131 static const long double S[9] = {
145 -1.193945051381510095614E-3L, 132 -1.193945051381510095614E-3L, 7.220599478036909672331E-3L,
146 7.220599478036909672331E-3L, 133 -9.622023360406271645744E-3L, -4.219773360705915470089E-2L,
147 -9.622023360406271645744E-3L, 134 1.665386113720805206758E-1L, -4.200263503403344054473E-2L,
148 -4.219773360705915470089E-2L, 135 -6.558780715202540684668E-1L, 5.772156649015328608253E-1L,
149 1.665386113720805206758E-1L, 136 1.000000000000000000000E0L,
150 -4.200263503403344054473E-2L,
151 -6.558780715202540684668E-1L,
152 5.772156649015328608253E-1L,
153 1.000000000000000000000E0L,
154 }; 137 };
155 138
156 /* 1/tgamma(-x) = z P(z) 139 /* 1/tgamma(-x) = z P(z)
157 * z(x) = 1/x 140 * z(x) = 1/x
158 * 0 < x < 0.03125 141 * 0 < x < 0.03125
159 * Peak relative error 5.16e-23 142 * Peak relative error 5.16e-23
160 * Relative error spread = 2.5e-24 143 * Relative error spread = 2.5e-24
161 */ 144 */
162 static const long double SN[9] = { 145 static const long double SN[9] = {
163 1.133374167243894382010E-3L, 146 1.133374167243894382010E-3L, 7.220837261893170325704E-3L,
164 7.220837261893170325704E-3L, 147 9.621911155035976733706E-3L, -4.219773343731191721664E-2L,
165 9.621911155035976733706E-3L, 148 -1.665386113944413519335E-1L, -4.200263503402112910504E-2L,
166 -4.219773343731191721664E-2L, 149 6.558780715202536547116E-1L, 5.772156649015328608727E-1L,
167 -1.665386113944413519335E-1L, 150 -1.000000000000000000000E0L,
168 -4.200263503402112910504E-2L,
169 6.558780715202536547116E-1L,
170 5.772156649015328608727E-1L,
171 -1.000000000000000000000E0L,
172 }; 151 };
173 152
174 static const long double PIL = 3.1415926535897932384626L; 153 static const long double PIL = 3.1415926535897932384626L;
175 154
176 /* Gamma function computed by Stirling's formula. 155 /* Gamma function computed by Stirling's formula.
177 */ 156 */
178 static long double stirf(long double x) 157 static long double stirf(long double x) {
179 { 158 long double y, w, v;
180 » long double y, w, v;
181 159
182 » w = 1.0/x; 160 w = 1.0 / x;
183 » /* For large x, use rational coefficients from the analytical expansion. */ 161 /* For large x, use rational coefficients from the analytical expansion. */
184 » if (x > 1024.0) 162 if (x > 1024.0)
185 » » w = (((((6.97281375836585777429E-5L * w 163 w = (((((6.97281375836585777429E-5L * w + 7.84039221720066627474E-4L) * w -
186 » » + 7.84039221720066627474E-4L) * w 164 2.29472093621399176955E-4L) *
187 » » - 2.29472093621399176955E-4L) * w 165 w -
188 » » - 2.68132716049382716049E-3L) * w 166 2.68132716049382716049E-3L) *
189 » » + 3.47222222222222222222E-3L) * w 167 w +
190 » » + 8.33333333333333333333E-2L) * w 168 3.47222222222222222222E-3L) *
191 » » + 1.0; 169 w +
192 » else 170 8.33333333333333333333E-2L) *
193 » » w = 1.0 + w * __polevll(w, STIR, 8); 171 w +
194 » y = expl(x); 172 1.0;
195 » if (x > MAXSTIR) { /* Avoid overflow in pow() */ 173 else
196 » » v = powl(x, 0.5L * x - 0.25L); 174 w = 1.0 + w * __polevll(w, STIR, 8);
197 » » y = v * (v / y); 175 y = expl(x);
198 » } else { 176 if (x > MAXSTIR) { /* Avoid overflow in pow() */
199 » » y = powl(x, x - 0.5L) / y; 177 v = powl(x, 0.5L * x - 0.25L);
200 » } 178 y = v * (v / y);
201 » y = SQTPI * y * w; 179 } else {
202 » return y; 180 y = powl(x, x - 0.5L) / y;
181 }
182 y = SQTPI * y * w;
183 return y;
203 } 184 }
204 185
205 long double tgammal(long double x) 186 long double tgammal(long double x) {
206 { 187 long double p, q, z;
207 » long double p, q, z;
208 188
209 » if (!isfinite(x)) 189 if (!isfinite(x))
210 » » return x + INFINITY; 190 return x + INFINITY;
211 191
212 » q = fabsl(x); 192 q = fabsl(x);
213 » if (q > 13.0) { 193 if (q > 13.0) {
214 » » if (x < 0.0) { 194 if (x < 0.0) {
215 » » » p = floorl(q); 195 p = floorl(q);
216 » » » z = q - p; 196 z = q - p;
217 » » » if (z == 0) 197 if (z == 0)
218 » » » » return 0 / z; 198 return 0 / z;
219 » » » if (q > MAXGAML) { 199 if (q > MAXGAML) {
220 » » » » z = 0; 200 z = 0;
221 » » » } else { 201 } else {
222 » » » » if (z > 0.5) { 202 if (z > 0.5) {
223 » » » » » p += 1.0; 203 p += 1.0;
224 » » » » » z = q - p; 204 z = q - p;
225 » » » » } 205 }
226 » » » » z = q * sinl(PIL * z); 206 z = q * sinl(PIL * z);
227 » » » » z = fabsl(z) * stirf(q); 207 z = fabsl(z) * stirf(q);
228 » » » » z = PIL/z; 208 z = PIL / z;
229 » » » } 209 }
230 » » » if (0.5 * p == floorl(q * 0.5)) 210 if (0.5 * p == floorl(q * 0.5))
231 » » » » z = -z; 211 z = -z;
232 » » } else if (x > MAXGAML) { 212 } else if (x > MAXGAML) {
233 » » » z = x * 0x1p16383L; 213 z = x * 0x1p16383L;
234 » » } else { 214 } else {
235 » » » z = stirf(x); 215 z = stirf(x);
236 » » } 216 }
237 » » return z; 217 return z;
238 » } 218 }
239 219
240 » z = 1.0; 220 z = 1.0;
241 » while (x >= 3.0) { 221 while (x >= 3.0) {
242 » » x -= 1.0; 222 x -= 1.0;
243 » » z *= x; 223 z *= x;
244 » } 224 }
245 » while (x < -0.03125L) { 225 while (x < -0.03125L) {
246 » » z /= x; 226 z /= x;
247 » » x += 1.0; 227 x += 1.0;
248 » } 228 }
249 » if (x <= 0.03125L) 229 if (x <= 0.03125L)
250 » » goto small; 230 goto small;
251 » while (x < 2.0) { 231 while (x < 2.0) {
252 » » z /= x; 232 z /= x;
253 » » x += 1.0; 233 x += 1.0;
254 » } 234 }
255 » if (x == 2.0) 235 if (x == 2.0)
256 » » return z; 236 return z;
257 237
258 » x -= 2.0; 238 x -= 2.0;
259 » p = __polevll(x, P, 7); 239 p = __polevll(x, P, 7);
260 » q = __polevll(x, Q, 8); 240 q = __polevll(x, Q, 8);
261 » z = z * p / q; 241 z = z * p / q;
262 » return z; 242 return z;
263 243
264 small: 244 small:
265 » /* z==1 if x was originally +-0 */ 245 /* z==1 if x was originally +-0 */
266 » if (x == 0 && z != 1) 246 if (x == 0 && z != 1)
267 » » return x / x; 247 return x / x;
268 » if (x < 0.0) { 248 if (x < 0.0) {
269 » » x = -x; 249 x = -x;
270 » » q = z / (x * __polevll(x, SN, 8)); 250 q = z / (x * __polevll(x, SN, 8));
271 » } else 251 } else
272 » » q = z / (x * __polevll(x, S, 8)); 252 q = z / (x * __polevll(x, S, 8));
273 » return q; 253 return q;
274 } 254 }
275 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384 255 #elif LDBL_MANT_DIG == 113 && LDBL_MAX_EXP == 16384
276 // TODO: broken implementation to make things compile 256 // TODO: broken implementation to make things compile
277 long double tgammal(long double x) 257 long double tgammal(long double x) {
278 { 258 return tgamma(x);
279 » return tgamma(x);
280 } 259 }
281 #endif 260 #endif
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698