/src/r-source/src/nmath/lgamma.c

Source
/*
 *  Mathlib : A C Library of Special Functions
 *  Copyright (C) 2000-2020 The R Core Team
 *  Copyright (C) 1998 Ross Ihaka
 *
 *  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 lgammafn_sign(double x, int *sgn);
 *    double lgammafn(double x);
 *
 *  DESCRIPTION
 *
 *    The function lgammafn computes log|gamma(x)|.  The function
 *    lgammafn_sign in addition assigns the sign of the gamma function
 *    to the address in the second argument if this is not NULL.
 *
 *  NOTES
 *
 *    This routine is a translation into C of a Fortran subroutine
 *    by W. Fullerton of Los Alamos Scientific Laboratory.
 *
 *    The accuracy of this routine compares (very) favourably
 *    with those of the Sun Microsystems portable mathematical
 *    library.
 *
 *  ./toms708.c  has  gamln()
 */

#include "nmath.h"

double lgammafn_sign(double x, int *sgn)
{
    double ans, y, sinpiy;

#ifdef NOMORE_FOR_THREADS
    static double xmax = 0.;
    static double dxrel = 0.;

    if (xmax == 0) {/* initialize machine dependent constants _ONCE_ */
  xmax = d1mach(2)/log(d1mach(2));/* = 2.533 e305  for IEEE double */
  dxrel = sqrt (d1mach(4));/* sqrt(Eps) ~ 1.49 e-8  for IEEE double */
    }
#else
/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
   xmax  = DBL_MAX / log(DBL_MAX) = 2^1024 / (1024 * log(2)) = 2^1014 / log(2)
   dxrel = sqrt(DBL_EPSILON) = 2^-26 = 5^26 * 1e-26 (is *exact* below !)
 */
#define xmax  2.5327372760800758e+305
#define dxrel 1.490116119384765625e-8
#endif

    if (sgn != NULL) *sgn = 1;

#ifdef IEEE_754
    if(ISNAN(x)) return x;
#endif

    if (sgn != NULL && x < 0 && fmod(floor(-x), 2.) == 0)
  *sgn = -1;

    if (x <= 0 && x == trunc(x)) { /* Negative integer argument */
  // No warning: this is the best answer; was  ML_WARNING(ME_RANGE, "lgamma");
  return ML_POSINF;/* +Inf, since lgamma(x) = log|gamma(x)| */
    }

    y = fabs(x);

    if (y < 1e-306) return -log(y); // denormalized range, R change
    if (y <= 10) return log(fabs(gammafn(x)));
    /*
      ELSE  y = |x| > 10 ---------------------- */

    if (y > xmax) {
  // No warning: +Inf is the best answer
  return ML_POSINF;
    }

    if (x > 0) { /* i.e. y = x > 10 */
#ifdef IEEE_754
  if(x > 1e17)
      return(x*(log(x) - 1.));
  else if(x > 4934720.)
      return(M_LN_SQRT_2PI + (x - 0.5) * log(x) - x);
  else
#endif
      return M_LN_SQRT_2PI + (x - 0.5) * log(x) - x + lgammacor(x);
    }
    /* else: x < -10; y = -x > 10 */
    sinpiy = fabs(sinpi(y));

    if (sinpiy == 0) { /* Negative integer argument ===
        Now UNNECESSARY: caught above */
  MATHLIB_WARNING(" ** should NEVER happen! *** [lgamma.c: Neg.int, y=%g]\n",y);
  ML_WARN_return_NAN;
    }

    ans = M_LN_SQRT_PId2 + (x - 0.5) * log(y) - x - log(sinpiy) - lgammacor(y);

    if(fabs((x - trunc(x - 0.5)) * ans / x) < dxrel) {

  /* The answer is less than half precision because
   * the argument is too near a negative integer; e.g. for  lgamma(1e-7 - 11) */

  ML_WARNING(ME_PRECISION, "lgamma");
    }

    return ans;
}

double lgammafn(double x)
{
    return lgammafn_sign(x, NULL);
}

Coverage Report

Created: 2026-10-28 06:40

Line	Count	Source
1		/*
2		* Mathlib : A C Library of Special Functions
3		* Copyright (C) 2000-2020 The R Core Team
4		* Copyright (C) 1998 Ross Ihaka
5		*
6		* This program is free software; you can redistribute it and/or modify
7		* it under the terms of the GNU General Public License as published by
8		* the Free Software Foundation; either version 2 of the License, or
9		* (at your option) any later version.
10		*
11		* This program is distributed in the hope that it will be useful,
12		* but WITHOUT ANY WARRANTY; without even the implied warranty of
13		* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14		* GNU General Public License for more details.
15		*
16		* You should have received a copy of the GNU General Public License
17		* along with this program; if not, a copy is available at
18		* https://www.R-project.org/Licenses/
19		*
20		* SYNOPSIS
21		*
22		* #include <Rmath.h>
23		* double lgammafn_sign(double x, int *sgn);
24		* double lgammafn(double x);
25		*
26		* DESCRIPTION
27		*
28		* The function lgammafn computes log\|gamma(x)\|. The function
29		* lgammafn_sign in addition assigns the sign of the gamma function
30		* to the address in the second argument if this is not NULL.
31		*
32		* NOTES
33		*
34		* This routine is a translation into C of a Fortran subroutine
35		* by W. Fullerton of Los Alamos Scientific Laboratory.
36		*
37		* The accuracy of this routine compares (very) favourably
38		* with those of the Sun Microsystems portable mathematical
39		* library.
40		*
41		* ./toms708.c has gamln()
42		*/
43
44		#include "nmath.h"
45
46		double lgammafn_sign(double x, int *sgn)
47	0	{
48	0	double ans, y, sinpiy;
49
50		#ifdef NOMORE_FOR_THREADS
51		static double xmax = 0.;
52		static double dxrel = 0.;
53
54		if (xmax == 0) {/* initialize machine dependent constants _ONCE_ */
55		xmax = d1mach(2)/log(d1mach(2));/* = 2.533 e305 for IEEE double */
56		dxrel = sqrt (d1mach(4));/* sqrt(Eps) ~ 1.49 e-8 for IEEE double */
57		}
58		#else
59		/* For IEEE double precision DBL_EPSILON = 2^-52 = 2.220446049250313e-16 :
60		xmax = DBL_MAX / log(DBL_MAX) = 2^1024 / (1024 * log(2)) = 2^1014 / log(2)
61		dxrel = sqrt(DBL_EPSILON) = 2^-26 = 5^26 * 1e-26 (is exact below !)
62		*/
63	0	#define xmax 2.5327372760800758e+305
64	0	#define dxrel 1.490116119384765625e-8
65	0	#endif
66
67	0	if (sgn != NULL) *sgn = 1;
68
69	0	#ifdef IEEE_754
70	0	if(ISNAN(x)) return x;
71	0	#endif
72
73	0	if (sgn != NULL && x < 0 && fmod(floor(-x), 2.) == 0)
74	0	*sgn = -1;
75
76	0	if (x <= 0 && x == trunc(x)) { /* Negative integer argument */
77		// No warning: this is the best answer; was ML_WARNING(ME_RANGE, "lgamma");
78	0	return ML_POSINF;/* +Inf, since lgamma(x) = log\|gamma(x)\| */
79	0	}
80
81	0	y = fabs(x);
82
83	0	if (y < 1e-306) return -log(y); // denormalized range, R change
84	0	if (y <= 10) return log(fabs(gammafn(x)));
85		/*
86		ELSE y = \|x\| > 10 ---------------------- */
87
88	0	if (y > xmax) {
89		// No warning: +Inf is the best answer
90	0	return ML_POSINF;
91	0	}
92
93	0	if (x > 0) { /* i.e. y = x > 10 */
94	0	#ifdef IEEE_754
95	0	if(x > 1e17)
96	0	return(x*(log(x) - 1.));
97	0	else if(x > 4934720.)
98	0	return(M_LN_SQRT_2PI + (x - 0.5) * log(x) - x);
99	0	else
100	0	#endif
101	0	return M_LN_SQRT_2PI + (x - 0.5) * log(x) - x + lgammacor(x);
102	0	}
103		/* else: x < -10; y = -x > 10 */
104	0	sinpiy = fabs(sinpi(y));
105
106	0	if (sinpiy == 0) { /* Negative integer argument ===
107		Now UNNECESSARY: caught above */
108	0	MATHLIB_WARNING(" should NEVER happen! * [lgamma.c: Neg.int, y=%g]\n",y);
109	0	ML_WARN_return_NAN;
110	0	}
111
112	0	ans = M_LN_SQRT_PId2 + (x - 0.5) * log(y) - x - log(sinpiy) - lgammacor(y);
113
114	0	if(fabs((x - trunc(x - 0.5)) * ans / x) < dxrel) {
115
116		/* The answer is less than half precision because
117		* the argument is too near a negative integer; e.g. for lgamma(1e-7 - 11) */
118
119	0	ML_WARNING(ME_PRECISION, "lgamma");
120	0	}
121
122	0	return ans;
123	0	}
124
125		double lgammafn(double x)
126	0	{
127	0	return lgammafn_sign(x, NULL);
128	0	}