| /* |
| * Mathlib : A C Library of Special Functions |
| * Copyright (C) 1998 Ross Ihaka |
| * Copyright (C) 2000-2001 The R Core Team |
| * |
| * This program is free software; you can redistribute it and/or modify |
| * it under the terms of the GNU General Public License as published by |
| * the Free Software Foundation; either version 2 of the License, or |
| * (at your option) any later version. |
| * |
| * This program is distributed in the hope that it will be useful, |
| * but WITHOUT ANY WARRANTY; without even the implied warranty of |
| * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the |
| * GNU General Public License for more details. |
| * |
| * You should have received a copy of the GNU General Public License |
| * along with this program; if not, a copy is available at |
| * https://www.R-project.org/Licenses/ |
| * |
| * SYNOPSIS |
| * |
| * #include <Rmath.h> |
| * double lgammacor(double x); |
| * |
| * DESCRIPTION |
| * |
| * Compute the log gamma correction factor for x >= 10 so that |
| * |
| * log(gamma(x)) = .5*log(2*pi) + (x-.5)*log(x) -x + lgammacor(x) |
| * |
| * [ lgammacor(x) is called Del(x) in other contexts (e.g. dcdflib)] |
| * |
| * NOTES |
| * |
| * This routine is a translation into C of a Fortran subroutine |
| * written by W. Fullerton of Los Alamos Scientific Laboratory. |
| * |
| * SEE ALSO |
| * |
| * Loader(1999)'s stirlerr() {in ./stirlerr.c} is *very* similar in spirit, |
| * is faster and cleaner, but is only defined "fast" for half integers. |
| */ |
| |
| #include "nmath.h" |
| |
| double attribute_hidden lgammacor(double x) |
| { |
| const static double algmcs[15] = { |
| +.1666389480451863247205729650822e+0, |
| -.1384948176067563840732986059135e-4, |
| +.9810825646924729426157171547487e-8, |
| -.1809129475572494194263306266719e-10, |
| +.6221098041892605227126015543416e-13, |
| -.3399615005417721944303330599666e-15, |
| +.2683181998482698748957538846666e-17, |
| -.2868042435334643284144622399999e-19, |
| +.3962837061046434803679306666666e-21, |
| -.6831888753985766870111999999999e-23, |
| +.1429227355942498147573333333333e-24, |
| -.3547598158101070547199999999999e-26, |
| +.1025680058010470912000000000000e-27, |
| -.3401102254316748799999999999999e-29, |
| +.1276642195630062933333333333333e-30 |
| }; |
| |
| double tmp; |
| |
| /* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 : |
| * xbig = 2 ^ 26.5 |
| * xmax = DBL_MAX / 48 = 2^1020 / 3 */ |
| #define nalgm 5 |
| #define xbig 94906265.62425156 |
| #define xmax 3.745194030963158e306 |
| |
| if (x < 10) |
| ML_ERR_return_NAN |
| else if (x >= xmax) { |
| ML_ERROR(ME_UNDERFLOW, "lgammacor"); |
| /* allow to underflow below */ |
| } |
| else if (x < xbig) { |
| tmp = 10 / x; |
| return chebyshev_eval(tmp * tmp * 2 - 1, algmcs, nalgm) / x; |
| } |
| return 1 / (x * 12); |
| } |