| e_lgamma_r.c (2117) | e_lgamma_r.c (8870) |
|---|---|
| 1/* @(#)er_lgamma.c 5.1 93/09/24 */ 2/* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this | 1/* @(#)er_lgamma.c 5.1 93/09/24 */ 2/* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunPro, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this |
| 8 * software is freely granted, provided that this notice | 8 * software is freely granted, provided that this notice |
| 9 * is preserved. 10 * ==================================================== 11 */ 12 13#ifndef lint | 9 * is preserved. 10 * ==================================================== 11 */ 12 13#ifndef lint |
| 14static char rcsid[] = "$Id: e_lgamma_r.c,v 1.5 1994/08/10 20:31:07 jtc Exp $"; | 14static char rcsid[] = "$Id: e_lgamma_r.c,v 1.1.1.1 1994/08/19 09:39:44 jkh Exp $"; |
| 15#endif 16 17/* __ieee754_lgamma_r(x, signgamp) | 15#endif 16 17/* __ieee754_lgamma_r(x, signgamp) |
| 18 * Reentrant version of the logarithm of the Gamma function 19 * with user provide pointer for the sign of Gamma(x). | 18 * Reentrant version of the logarithm of the Gamma function 19 * with user provide pointer for the sign of Gamma(x). |
| 20 * 21 * Method: 22 * 1. Argument Reduction for 0 < x <= 8 | 20 * 21 * Method: 22 * 1. Argument Reduction for 0 < x <= 8 |
| 23 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may | 23 * Since gamma(1+s)=s*gamma(s), for x in [0,8], we may |
| 24 * reduce x to a number in [1.5,2.5] by 25 * lgamma(1+s) = log(s) + lgamma(s) 26 * for example, 27 * lgamma(7.3) = log(6.3) + lgamma(6.3) 28 * = log(6.3*5.3) + lgamma(5.3) 29 * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) 30 * 2. Polynomial approximation of lgamma around its 31 * minimun ymin=1.461632144968362245 to maintain monotonicity. --- 21 unchanged lines hidden (view full) --- 53 * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+.... 54 * (better formula: 55 * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...) 56 * Let z = 1/x, then we approximation 57 * f(z) = lgamma(x) - (x-0.5)(log(x)-1) 58 * by 59 * 3 5 11 60 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z | 24 * reduce x to a number in [1.5,2.5] by 25 * lgamma(1+s) = log(s) + lgamma(s) 26 * for example, 27 * lgamma(7.3) = log(6.3) + lgamma(6.3) 28 * = log(6.3*5.3) + lgamma(5.3) 29 * = log(6.3*5.3*4.3*3.3*2.3) + lgamma(2.3) 30 * 2. Polynomial approximation of lgamma around its 31 * minimun ymin=1.461632144968362245 to maintain monotonicity. --- 21 unchanged lines hidden (view full) --- 53 * lgamma(x)~(x-0.5)log(x)-x+0.5*log(2pi)+1/(12x)-1/(360x**3)+.... 54 * (better formula: 55 * lgamma(x)~(x-0.5)*(log(x)-1)-.5*(log(2pi)-1) + ...) 56 * Let z = 1/x, then we approximation 57 * f(z) = lgamma(x) - (x-0.5)(log(x)-1) 58 * by 59 * 3 5 11 60 * w = w0 + w1*z + w2*z + w3*z + ... + w6*z |
| 61 * where | 61 * where |
| 62 * |w - f(z)| < 2**-58.74 | 62 * |w - f(z)| < 2**-58.74 |
| 63 * | 63 * |
| 64 * 4. For negative x, since (G is gamma function) 65 * -x*G(-x)*G(x) = pi/sin(pi*x), 66 * we have 67 * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) 68 * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 | 64 * 4. For negative x, since (G is gamma function) 65 * -x*G(-x)*G(x) = pi/sin(pi*x), 66 * we have 67 * G(x) = pi/(sin(pi*x)*(-x)*G(-x)) 68 * since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0 |
| 69 * Hence, for x<0, signgam = sign(sin(pi*x)) and | 69 * Hence, for x<0, signgam = sign(sin(pi*x)) and |
| 70 * lgamma(x) = log(|Gamma(x)|) 71 * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); | 70 * lgamma(x) = log(|Gamma(x)|) 71 * = log(pi/(|x*sin(pi*x)|)) - lgamma(-x); |
| 72 * Note: one should avoid compute pi*(-x) directly in the | 72 * Note: one should avoid compute pi*(-x) directly in the |
| 73 * computation of sin(pi*(-x)). | 73 * computation of sin(pi*(-x)). |
| 74 * | 74 * |
| 75 * 5. Special Cases 76 * lgamma(2+s) ~ s*(1-Euler) for tiny s 77 * lgamma(1)=lgamma(2)=0 78 * lgamma(x) ~ -log(x) for tiny x 79 * lgamma(0) = lgamma(inf) = inf 80 * lgamma(-integer) = +-inf | 75 * 5. Special Cases 76 * lgamma(2+s) ~ s*(1-Euler) for tiny s 77 * lgamma(1)=lgamma(2)=0 78 * lgamma(x) ~ -log(x) for tiny x 79 * lgamma(0) = lgamma(inf) = inf 80 * lgamma(-integer) = +-inf |
| 81 * | 81 * |
| 82 */ 83 84#include "math.h" 85#include "math_private.h" 86 87#ifdef __STDC__ | 82 */ 83 84#include "math.h" 85#include "math_private.h" 86 87#ifdef __STDC__ |
| 88static const double | 88static const double |
| 89#else | 89#else |
| 90static double | 90static double |
| 91#endif 92two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ 93half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 94one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 95pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ 96a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */ 97a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */ 98a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */ --- 96 unchanged lines hidden (view full) --- 195 GET_LOW_WORD(n,z); 196 n &= 1; 197 y = n; 198 n<<= 2; 199 } 200 } 201 switch (n) { 202 case 0: y = __kernel_sin(pi*y,zero,0); break; | 91#endif 92two52= 4.50359962737049600000e+15, /* 0x43300000, 0x00000000 */ 93half= 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */ 94one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */ 95pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */ 96a0 = 7.72156649015328655494e-02, /* 0x3FB3C467, 0xE37DB0C8 */ 97a1 = 3.22467033424113591611e-01, /* 0x3FD4A34C, 0xC4A60FAD */ 98a2 = 6.73523010531292681824e-02, /* 0x3FB13E00, 0x1A5562A7 */ --- 96 unchanged lines hidden (view full) --- 195 GET_LOW_WORD(n,z); 196 n &= 1; 197 y = n; 198 n<<= 2; 199 } 200 } 201 switch (n) { 202 case 0: y = __kernel_sin(pi*y,zero,0); break; |
| 203 case 1: | 203 case 1: |
| 204 case 2: y = __kernel_cos(pi*(0.5-y),zero); break; | 204 case 2: y = __kernel_cos(pi*(0.5-y),zero); break; |
| 205 case 3: | 205 case 3: |
| 206 case 4: y = __kernel_sin(pi*(one-y),zero,0); break; 207 case 5: 208 case 6: y = -__kernel_cos(pi*(y-1.5),zero); break; 209 default: y = __kernel_sin(pi*(y-2.0),zero,0); break; 210 } 211 return -y; 212} 213 --- 56 unchanged lines hidden (view full) --- 270 case 1: 271 z = y*y; 272 w = z*y; 273 p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ 274 p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); 275 p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); 276 p = z*p1-(tt-w*(p2+y*p3)); 277 r += (tf + p); break; | 206 case 4: y = __kernel_sin(pi*(one-y),zero,0); break; 207 case 5: 208 case 6: y = -__kernel_cos(pi*(y-1.5),zero); break; 209 default: y = __kernel_sin(pi*(y-2.0),zero,0); break; 210 } 211 return -y; 212} 213 --- 56 unchanged lines hidden (view full) --- 270 case 1: 271 z = y*y; 272 w = z*y; 273 p1 = t0+w*(t3+w*(t6+w*(t9 +w*t12))); /* parallel comp */ 274 p2 = t1+w*(t4+w*(t7+w*(t10+w*t13))); 275 p3 = t2+w*(t5+w*(t8+w*(t11+w*t14))); 276 p = z*p1-(tt-w*(p2+y*p3)); 277 r += (tf + p); break; |
| 278 case 2: | 278 case 2: |
| 279 p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); 280 p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); 281 r += (-0.5*y + p1/p2); 282 } 283 } 284 else if(ix<0x40200000) { /* x < 8.0 */ 285 i = (int)x; 286 t = zero; --- 12 unchanged lines hidden (view full) --- 299 } 300 /* 8.0 <= x < 2**58 */ 301 } else if (ix < 0x43900000) { 302 t = __ieee754_log(x); 303 z = one/x; 304 y = z*z; 305 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); 306 r = (x-half)*(t-one)+w; | 279 p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5))))); 280 p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5)))); 281 r += (-0.5*y + p1/p2); 282 } 283 } 284 else if(ix<0x40200000) { /* x < 8.0 */ 285 i = (int)x; 286 t = zero; --- 12 unchanged lines hidden (view full) --- 299 } 300 /* 8.0 <= x < 2**58 */ 301 } else if (ix < 0x43900000) { 302 t = __ieee754_log(x); 303 z = one/x; 304 y = z*z; 305 w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6))))); 306 r = (x-half)*(t-one)+w; |
| 307 } else | 307 } else |
| 308 /* 2**58 <= x <= inf */ 309 r = x*(__ieee754_log(x)-one); 310 if(hx<0) r = nadj - r; 311 return r; 312} | 308 /* 2**58 <= x <= inf */ 309 r = x*(__ieee754_log(x)-one); 310 if(hx<0) r = nadj - r; 311 return r; 312} |