Remove whitespace and 2 blank lines.
This commit is contained in:
parent
486de25d46
commit
3fbe079a34
Notes:
svn2git
2020-12-20 02:59:44 +00:00
svn path=/head/; revision=272457
@ -14,12 +14,12 @@
|
||||
__FBSDID("$FreeBSD$");
|
||||
|
||||
/* __ieee754_lgamma_r(x, signgamp)
|
||||
* Reentrant version of the logarithm of the Gamma function
|
||||
* with user provide pointer for the sign of Gamma(x).
|
||||
* Reentrant version of the logarithm of the Gamma function
|
||||
* with user provide pointer for the sign of Gamma(x).
|
||||
*
|
||||
* Method:
|
||||
* 1. Argument Reduction for 0 < x <= 8
|
||||
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
|
||||
* Since gamma(1+s)=s*gamma(s), for x in [0,8], we may
|
||||
* reduce x to a number in [1.5,2.5] by
|
||||
* lgamma(1+s) = log(s) + lgamma(s)
|
||||
* for example,
|
||||
@ -57,20 +57,20 @@ __FBSDID("$FreeBSD$");
|
||||
* by
|
||||
* 3 5 11
|
||||
* w = w0 + w1*z + w2*z + w3*z + ... + w6*z
|
||||
* where
|
||||
* where
|
||||
* |w - f(z)| < 2**-58.74
|
||||
*
|
||||
*
|
||||
* 4. For negative x, since (G is gamma function)
|
||||
* -x*G(-x)*G(x) = pi/sin(pi*x),
|
||||
* we have
|
||||
* G(x) = pi/(sin(pi*x)*(-x)*G(-x))
|
||||
* since G(-x) is positive, sign(G(x)) = sign(sin(pi*x)) for x<0
|
||||
* Hence, for x<0, signgam = sign(sin(pi*x)) and
|
||||
* Hence, for x<0, signgam = sign(sin(pi*x)) and
|
||||
* lgamma(x) = log(|Gamma(x)|)
|
||||
* = log(pi/(|x*sin(pi*x)|)) - lgamma(-x);
|
||||
* Note: one should avoid compute pi*(-x) directly in the
|
||||
* Note: one should avoid compute pi*(-x) directly in the
|
||||
* computation of sin(pi*(-x)).
|
||||
*
|
||||
*
|
||||
* 5. Special Cases
|
||||
* lgamma(2+s) ~ s*(1-Euler) for tiny s
|
||||
* lgamma(1) = lgamma(2) = 0
|
||||
@ -78,7 +78,6 @@ __FBSDID("$FreeBSD$");
|
||||
* lgamma(0) = lgamma(neg.integer) = inf and raise divide-by-zero
|
||||
* lgamma(inf) = inf
|
||||
* lgamma(-inf) = inf (bug for bug compatible with C99!?)
|
||||
*
|
||||
*/
|
||||
|
||||
#include <float.h>
|
||||
@ -187,9 +186,9 @@ sin_pi(double x)
|
||||
|
||||
switch (n) {
|
||||
case 0: y = __kernel_sin(pi*y,zero,0); break;
|
||||
case 1:
|
||||
case 1:
|
||||
case 2: y = __kernel_cos(pi*(0.5-y),zero); break;
|
||||
case 3:
|
||||
case 3:
|
||||
case 4: y = __kernel_sin(pi*(one-y),zero,0); break;
|
||||
case 5:
|
||||
case 6: y = -__kernel_cos(pi*(y-1.5),zero); break;
|
||||
@ -263,7 +262,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
|
||||
p3 = t2+w*(t5+w*(t8+w*(t11+w*t14)));
|
||||
p = z*p1-(tt-w*(p2+y*p3));
|
||||
r += (tf + p); break;
|
||||
case 2:
|
||||
case 2:
|
||||
p1 = y*(u0+y*(u1+y*(u2+y*(u3+y*(u4+y*u5)))));
|
||||
p2 = one+y*(v1+y*(v2+y*(v3+y*(v4+y*v5))));
|
||||
r += (-0.5*y + p1/p2);
|
||||
@ -291,7 +290,7 @@ __ieee754_lgamma_r(double x, int *signgamp)
|
||||
y = z*z;
|
||||
w = w0+z*(w1+y*(w2+y*(w3+y*(w4+y*(w5+y*w6)))));
|
||||
r = (x-half)*(t-one)+w;
|
||||
} else
|
||||
} else
|
||||
/* 2**58 <= x <= inf */
|
||||
r = x*(__ieee754_log(x)-one);
|
||||
if(hx<0) r = nadj - r;
|
||||
@ -301,4 +300,3 @@ __ieee754_lgamma_r(double x, int *signgamp)
|
||||
#if (LDBL_MANT_DIG == 53)
|
||||
__weak_reference(lgamma_r, lgammal_r);
|
||||
#endif
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user