| Index: fusl/src/math/tgamma.c
|
| diff --git a/fusl/src/math/tgamma.c b/fusl/src/math/tgamma.c
|
| index feb0616ec674db8028211c684652d0b5ae9b99f6..466acf0d8a400631f9ee357bb736fb0b9fada947 100644
|
| --- a/fusl/src/math/tgamma.c
|
| +++ b/fusl/src/math/tgamma.c
|
| @@ -190,51 +190,50 @@ double tgamma(double x) {
|
| }
|
|
|
| #if 0
|
| -double __lgamma_r(double x, int *sign)
|
| -{
|
| - double r, absx;
|
| -
|
| - *sign = 1;
|
| -
|
| - /* special cases */
|
| - if (!isfinite(x))
|
| - /* lgamma(nan)=nan, lgamma(+-inf)=inf */
|
| - return x*x;
|
| -
|
| - /* integer arguments */
|
| - if (x == floor(x) && x <= 2) {
|
| - /* n <= 0: lgamma(n)=inf with divbyzero */
|
| - /* n == 1,2: lgamma(n)=0 */
|
| - if (x <= 0)
|
| - return 1/0.0;
|
| - return 0;
|
| - }
|
| -
|
| - absx = fabs(x);
|
| -
|
| - /* lgamma(x) ~ -log(|x|) for tiny |x| */
|
| - if (absx < 0x1p-54) {
|
| - *sign = 1 - 2*!!signbit(x);
|
| - return -log(absx);
|
| - }
|
| -
|
| - /* use tgamma for smaller |x| */
|
| - if (absx < 128) {
|
| - x = tgamma(x);
|
| - *sign = 1 - 2*!!signbit(x);
|
| - return log(fabs(x));
|
| - }
|
| -
|
| - /* second term (log(S)-g) could be more precise here.. */
|
| - /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
|
| - r = (absx-0.5)*(log(absx+gmhalf)-1) + (log(S(absx)) - (gmhalf+0.5));
|
| - if (x < 0) {
|
| - /* reflection formula for negative x */
|
| - x = sinpi(absx);
|
| - *sign = 2*!!signbit(x) - 1;
|
| - r = log(pi/(fabs(x)*absx)) - r;
|
| - }
|
| - return r;
|
| +double __lgamma_r(double x, int* sign) {
|
| + double r, absx;
|
| +
|
| + *sign = 1;
|
| +
|
| + /* special cases */
|
| + if (!isfinite(x))
|
| + /* lgamma(nan)=nan, lgamma(+-inf)=inf */
|
| + return x * x;
|
| +
|
| + /* integer arguments */
|
| + if (x == floor(x) && x <= 2) {
|
| + /* n <= 0: lgamma(n)=inf with divbyzero */
|
| + /* n == 1,2: lgamma(n)=0 */
|
| + if (x <= 0)
|
| + return 1 / 0.0;
|
| + return 0;
|
| + }
|
| +
|
| + absx = fabs(x);
|
| +
|
| + /* lgamma(x) ~ -log(|x|) for tiny |x| */
|
| + if (absx < 0x1p-54) {
|
| + *sign = 1 - 2 * !!signbit(x);
|
| + return -log(absx);
|
| + }
|
| +
|
| + /* use tgamma for smaller |x| */
|
| + if (absx < 128) {
|
| + x = tgamma(x);
|
| + *sign = 1 - 2 * !!signbit(x);
|
| + return log(fabs(x));
|
| + }
|
| +
|
| + /* second term (log(S)-g) could be more precise here.. */
|
| + /* or with stirling: (|x|-0.5)*(log(|x|)-1) + poly(1/|x|) */
|
| + r = (absx - 0.5) * (log(absx + gmhalf) - 1) + (log(S(absx)) - (gmhalf + 0.5));
|
| + if (x < 0) {
|
| + /* reflection formula for negative x */
|
| + x = sinpi(absx);
|
| + *sign = 2 * !!signbit(x) - 1;
|
| + r = log(pi / (fabs(x) * absx)) - r;
|
| + }
|
| + return r;
|
| }
|
|
|
| weak_alias(__lgamma_r, lgamma_r);
|
|
|