| e_lgamma_r.c (271147) | e_lgamma_r.c (271651) |
|---|---|
| 1 2/* @(#)e_lgamma_r.c 1.3 95/01/18 */ 3/* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 * 13 */ 14 15#include <sys/cdefs.h> | 1 2/* @(#)e_lgamma_r.c 1.3 95/01/18 */ 3/* 4 * ==================================================== 5 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 6 * 7 * Developed at SunSoft, a Sun Microsystems, Inc. business. 8 * Permission to use, copy, modify, and distribute this 9 * software is freely granted, provided that this notice 10 * is preserved. 11 * ==================================================== 12 * 13 */ 14 15#include <sys/cdefs.h> |
| 16__FBSDID("$FreeBSD: head/lib/msun/src/e_lgamma_r.c 271147 2014-09-04 23:50:05Z kargl $"); | 16__FBSDID("$FreeBSD: head/lib/msun/src/e_lgamma_r.c 271651 2014-09-15 23:21:57Z kargl $"); |
| 17 18/* __ieee754_lgamma_r(x, signgamp) 19 * Reentrant version of the logarithm of the Gamma function 20 * with user provide pointer for the sign of Gamma(x). 21 * 22 * Method: 23 * 1. Argument Reduction for 0 < x <= 8 24 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may --- 53 unchanged lines hidden (view full) --- 78 * lgamma(1) = lgamma(2) = 0 79 * lgamma(x) ~ -log(|x|) for tiny x 80 * lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero 81 * lgamma(inf) = inf 82 * lgamma(-inf) = inf (bug for bug compatible with C99!?) 83 * 84 */ 85 | 17 18/* __ieee754_lgamma_r(x, signgamp) 19 * Reentrant version of the logarithm of the Gamma function 20 * with user provide pointer for the sign of Gamma(x). 21 * 22 * Method: 23 * 1. Argument Reduction for 0 < x <= 8 24 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may --- 53 unchanged lines hidden (view full) --- 78 * lgamma(1) = lgamma(2) = 0 79 * lgamma(x) ~ -log(|x|) for tiny x 80 * lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero 81 * lgamma(inf) = inf 82 * lgamma(-inf) = inf (bug for bug compatible with C99!?) 83 * 84 */ 85 |
| 86#include <float.h> 87 |
|
| 86#include "math.h" 87#include "math_private.h" 88 89static const volatile double vzero = 0; 90 91static const double 92zero= 0.00000000000000000000e+00, 93half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ --- 151 unchanged lines hidden (view full) --- 245 else {y=x-one;i=2;} 246 } 247 switch(i) { 248 case 0: 249 z = y*y; 250 p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); 251 p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); 252 p = y*p1+p2; | 88#include "math.h" 89#include "math_private.h" 90 91static const volatile double vzero = 0; 92 93static const double 94zero= 0.00000000000000000000e+00, 95half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ --- 151 unchanged lines hidden (view full) --- 247 else {y=x-one;i=2;} 248 } 249 switch(i) { 250 case 0: 251 z = y*y; 252 p1 = a0+z*(a2+z*(a4+z*(a6+z*(a8+z*a10)))); 253 p2 = z*(a1+z*(a3+z*(a5+z*(a7+z*(a9+z*a11))))); 254 p = y*p1+p2; |
| 253 r += (p-0.5*y); break; | 255 r += (p-y/2); break; |
| 254 case 1: 255 z = y*y; 256 w = z*y; 257 p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ 258 p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); 259 p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); 260 p = z*p1-(tt-w*(p2+y*p3)); 261 r += (tf + p); break; --- 6 unchanged lines hidden (view full) --- 268 else if(ix<0x40200000) { /* x < 8.0 */ 269 i = (int)x; 270 y = x-(double)i; 271 p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); 272 q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); 273 r = half*y+p/q; 274 z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ 275 switch(i) { | 256 case 1: 257 z = y*y; 258 w = z*y; 259 p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ 260 p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); 261 p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); 262 p = z*p1-(tt-w*(p2+y*p3)); 263 r += (tf + p); break; --- 6 unchanged lines hidden (view full) --- 270 else if(ix<0x40200000) { /* x < 8.0 */ 271 i = (int)x; 272 y = x-(double)i; 273 p = y*(s0+y*(s1+y*(s2+y*(s3+y*(s4+y*(s5+y*s6)))))); 274 q = one+y*(r1+y*(r2+y*(r3+y*(r4+y*(r5+y*r6))))); 275 r = half*y+p/q; 276 z = one; /* lgamma(1+s) = log(s) + lgamma(s) */ 277 switch(i) { |
| 276 case 7: z *= (y+6.0); /* FALLTHRU */ 277 case 6: z *= (y+5.0); /* FALLTHRU */ 278 case 5: z *= (y+4.0); /* FALLTHRU */ 279 case 4: z *= (y+3.0); /* FALLTHRU */ 280 case 3: z *= (y+2.0); /* FALLTHRU */ | 278 case 7: z *= (y+6); /* FALLTHRU */ 279 case 6: z *= (y+5); /* FALLTHRU */ 280 case 5: z *= (y+4); /* FALLTHRU */ 281 case 4: z *= (y+3); /* FALLTHRU */ 282 case 3: z *= (y+2); /* FALLTHRU */ |
| 281 r += __ieee754_log(z); break; 282 } 283 /* 8.0 <= x < 2**58 */ 284 } else if (ix < 0x43900000) { 285 t = __ieee754_log(x); 286 z = one/x; 287 y = z*z; 288 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); 289 r = (x-half)*(t-one)+w; 290 } else 291 /* 2**58 <= x <= inf */ 292 r = x*(__ieee754_log(x)-one); 293 if(hx<0) r = nadj - r; 294 return r; 295} | 283 r += __ieee754_log(z); break; 284 } 285 /* 8.0 <= x < 2**58 */ 286 } else if (ix < 0x43900000) { 287 t = __ieee754_log(x); 288 z = one/x; 289 y = z*z; 290 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); 291 r = (x-half)*(t-one)+w; 292 } else 293 /* 2**58 <= x <= inf */ 294 r = x*(__ieee754_log(x)-one); 295 if(hx<0) r = nadj - r; 296 return r; 297} |
| 298 299#if (LDBL_MANT_DIG == 53) 300__weak_reference(lgamma_r, lgammal_r); 301#endif 302 |
|