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

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

Issue 1706293003: [fusl] Remove some more tabs (Closed) Base URL: git@github.com:domokit/mojo.git@master
Patch Set: 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 /* 1 /*
2 "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964) 2 "A Precision Approximation of the Gamma Function" - Cornelius Lanczos (1964)
3 "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001) 3 "Lanczos Implementation of the Gamma Function" - Paul Godfrey (2001)
4 "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004) 4 "An Analysis of the Lanczos Gamma Approximation" - Glendon Ralph Pugh (2004)
5 5
6 approximation method: 6 approximation method:
7 7
8 (x - 0.5) S(x) 8 (x - 0.5) S(x)
9 Gamma(x) = (x + g - 0.5) * ---------------- 9 Gamma(x) = (x + g - 0.5) * ----------------
10 exp(x + g - 0.5) 10 exp(x + g - 0.5)
(...skipping 172 matching lines...) Expand 10 before | Expand all | Expand 10 after
183 dy = -dy; 183 dy = -dy;
184 z = -z; 184 z = -z;
185 } 185 }
186 r += dy * (gmhalf + 0.5) * r / y; 186 r += dy * (gmhalf + 0.5) * r / y;
187 z = pow(y, 0.5 * z); 187 z = pow(y, 0.5 * z);
188 y = r * z * z; 188 y = r * z * z;
189 return y; 189 return y;
190 } 190 }
191 191
192 #if 0 192 #if 0
193 double __lgamma_r(double x, int *sign) 193 double __lgamma_r(double x, int* sign) {
194 { 194 double r, absx;
195 » double r, absx;
196 195
197 » *sign = 1; 196 *sign = 1;
198 197
199 » /* special cases */ 198 /* special cases */
200 » if (!isfinite(x)) 199 if (!isfinite(x))
201 » » /* lgamma(nan)=nan, lgamma(+-inf)=inf */ 200 /* lgamma(nan)=nan, lgamma(+-inf)=inf */
202 » » return x*x; 201 return x * x;
203 202
204 » /* integer arguments */ 203 /* integer arguments */
205 » if (x == floor(x) && x <= 2) { 204 if (x == floor(x) && x <= 2) {
206 » » /* n <= 0: lgamma(n)=inf with divbyzero */ 205 /* n <= 0: lgamma(n)=inf with divbyzero */
207 » » /* n == 1,2: lgamma(n)=0 */ 206 /* n == 1,2: lgamma(n)=0 */
208 » » if (x <= 0) 207 if (x <= 0)
209 » » » return 1/0.0; 208 return 1 / 0.0;
210 » » return 0; 209 return 0;
211 » } 210 }
212 211
213 » absx = fabs(x); 212 absx = fabs(x);
214 213
215 » /* lgamma(x) ~ -log(|x|) for tiny |x| */ 214 /* lgamma(x) ~ -log(|x|) for tiny |x| */
216 » if (absx < 0x1p-54) { 215 if (absx < 0x1p-54) {
217 » » *sign = 1 - 2*!!signbit(x); 216 *sign = 1 - 2 * !!signbit(x);
218 » » return -log(absx); 217 return -log(absx);
219 » } 218 }
220 219
221 » /* use tgamma for smaller |x| */ 220 /* use tgamma for smaller |x| */
222 » if (absx < 128) { 221 if (absx < 128) {
223 » » x = tgamma(x); 222 x = tgamma(x);
224 » » *sign = 1 - 2*!!signbit(x); 223 *sign = 1 - 2 * !!signbit(x);
225 » » return log(fabs(x)); 224 return log(fabs(x));
226 » } 225 }
227 226
228 » /* second term (log(S)-g) could be more precise here.. */ 227 /* second term (log(S)-g) could be more precise here.. */
229 » /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */ 228 /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
230 » r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5)); 229 r = (absx - 0.5) * (log(absx + gmhalf) - 1) + (log(S(absx)) - (gmhalf + 0.5));
231 » if (x < 0) { 230 if (x < 0) {
232 » » /* reflection formula for negative x */ 231 /* reflection formula for negative x */
233 » » x = sinpi(absx); 232 x = sinpi(absx);
234 » » *sign = 2*!!signbit(x) - 1; 233 *sign = 2 * !!signbit(x) - 1;
235 » » r = log(pi/(fabs(x)*absx)) - r; 234 r = log(pi / (fabs(x) * absx)) - r;
236 » } 235 }
237 » return r; 236 return r;
238 } 237 }
239 238
240 weak_alias(__lgamma_r, lgamma_r); 239 weak_alias(__lgamma_r, lgamma_r);
241 #endif 240 #endif
OLDNEW

Powered by Google App Engine
This is Rietveld 408576698