Natural logarithm of gamma function Input parameters: X - argument Result: logarithm of the absolute value of the Gamma(X). Output parameters: SgnGam - sign(Gamma(X)) Domain: 0 < X < 2.55e305 -2.55e305 < X < 0, X is not an integer. ACCURACY: arithmetic domain # trials peak rms IEEE 0, 3 28000 5.4e-16 1.1e-16
| 179 | Translated to AlgoPascal by Bochkanov Sergey (2005, 2006, 2007). |
| 180 | *************************************************************************/ |
| 181 | double lngamma(double x, double& sgngam) |
| 182 | { |
| 183 | #ifndef ALGLIB_INTERCEPTS_SPECFUNCS |
| 184 | double result; |
| 185 | double a; |
| 186 | double b; |
| 187 | double c; |
| 188 | double p; |
| 189 | double q; |
| 190 | double u; |
| 191 | double w; |
| 192 | double z; |
| 193 | int i; |
| 194 | double logpi; |
| 195 | double ls2pi; |
| 196 | double tmp; |
| 197 | |
| 198 | sgngam = 1; |
| 199 | logpi = 1.14472988584940017414; |
| 200 | ls2pi = 0.91893853320467274178; |
| 201 | if( ap::fp_less(x,-34.0) ) |
| 202 | { |
| 203 | q = -x; |
| 204 | w = lngamma(q, tmp); |
| 205 | p = ap::ifloor(q); |
| 206 | i = ap::round_f(p); |
| 207 | if( i%2==0 ) |
| 208 | { |
| 209 | sgngam = -1; |
| 210 | } |
| 211 | else |
| 212 | { |
| 213 | sgngam = 1; |
| 214 | } |
| 215 | z = q-p; |
| 216 | if( ap::fp_greater(z,0.5) ) |
| 217 | { |
| 218 | p = p+1; |
| 219 | z = p-q; |
| 220 | } |
| 221 | z = q*sin(ap::pi()*z); |
| 222 | result = logpi-log(z)-w; |
| 223 | return result; |
| 224 | } |
| 225 | if( ap::fp_less(x,13) ) |
| 226 | { |
| 227 | z = 1; |
| 228 | p = 0; |
| 229 | u = x; |
| 230 | while(ap::fp_greater_eq(u,3)) |
| 231 | { |
| 232 | p = p-1; |
| 233 | u = x+p; |
| 234 | z = z*u; |
| 235 | } |
| 236 | while(ap::fp_less(u,2)) |
| 237 | { |
| 238 | z = z/u; |
no test coverage detected