Coverage Report

Created: 2026-09-14 06:56

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/r-source/src/nmath/stirlerr.c
Line
Count
Source
1
/*
2
 *  AUTHOR
3
 *    Catherine Loader, catherine@research.bell-labs.com.
4
 *    October 23, 2000.
5
 *
6
 *  Merge in to R:
7
 *  Copyright (C) 2000-2025, The R Core Team
8
 *
9
 *  This program is free software; you can redistribute it and/or modify
10
 *  it under the terms of the GNU General Public License as published by
11
 *  the Free Software Foundation; either version 2 of the License, or
12
 *  (at your option) any later version.
13
 *
14
 *  This program is distributed in the hope that it will be useful,
15
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
16
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17
 *  GNU General Public License for more details.
18
 *
19
 *  You should have received a copy of the GNU General Public License
20
 *  along with this program; if not, a copy is available at
21
 *  https://www.R-project.org/Licenses/
22
 *
23
 *
24
 *  DESCRIPTION
25
 *
26
 *    Computes the log of the error term in Stirling's formula.
27
 *      For n > 15, uses the series 1/12n - 1/360n^3 + ...
28
 *      For n <=15, integers or half-integers, uses stored values.
29
 *      For other n < 15, uses lgamma directly (don't use this to
30
 *        write lgamma!)
31
 *
32
 * Merge in to R:
33
 * Copyright (C) 2000, The R Core Team
34
 * R has lgammafn, and lgamma is not part of ISO C
35
 */
36
37
#include "nmath.h"
38
39
/* stirlerr(n) = log(n!) - log( sqrt(2*pi*n)*(n/e)^n )
40
 *             = log Gamma(n+1) - 1/2 * [log(2*pi) + log(n)] - n*[log(n) - 1]
41
 *             = log Gamma(n+1) - (n + 1/2) * log(n) + n - log(2*pi)/2
42
 *
43
 * see also lgammacor() in ./lgammacor.c  which computes almost the same!
44
 *
45
 * NB: stirlerr(n/2) & stirlerr((n+1)/2) are called from dt(x,n) for all real n > 0 ;
46
 *     stirlerr(x) from gammafn(x) when |x| > 10,  2|x| is integer, but |x| is *not* in {11:50}
47
 *     stirlerr(x) from dpois_raw(x, lam) for any x > 0  which itself is called by many,
48
 *                                        including pgamma(), hence ppois(), ..
49
50
 *     stirlerr(n), stirlerr(x), stirlerr(n-x) from binom_raw(x, n, ..) for all possible 0 < x < n
51
 */
52
53
attribute_hidden double stirlerr(double n)
54
0
{
55
56
0
#define S0 0.083333333333333333333       /* 1/12 */
57
0
#define S1 0.00277777777777777777778     /* 1/360 */
58
0
#define S2 0.00079365079365079365079365  /* 1/1260 */
59
0
#define S3 0.000595238095238095238095238 /* 1/1680 */
60
0
#define S4 0.0008417508417508417508417508/* 1/1188 */
61
0
#define S5 0.0019175269175269175269175262 // 691/360360
62
0
#define S6 0.0064102564102564102564102561 // 1/156
63
0
#define S7 0.029550653594771241830065352  // 3617/122400
64
0
#define S8 0.17964437236883057316493850   // 43867/244188
65
0
#define S9 1.3924322169059011164274315    // 174611/125400
66
0
#define S10 13.402864044168391994478957 // 77683/5796
67
0
#define S11 156.84828462600201730636509 // 236364091/1506960
68
0
#define S12 2193.1033333333333333333333 // 657931/300
69
0
#define S13 36108.771253724989357173269 // 3392780147/93960
70
0
#define S14 691472.26885131306710839498 // 1723168255201/2492028
71
0
#define S15 15238221.539407416192283370 // 7709321041217/505920
72
0
#define S16 382900751.39141414141414141 // 151628697551/396
73
/* #define S17 10882266035.784391089015145 // 26315271553053477373/2418179400 */
74
75
/*
76
  exact values for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0.
77
*/
78
0
    const static double sferr_halves[31] = {
79
0
  0.0, /* n=0 - wrong, place holder only */
80
0
  0.1534264097200273452913848,  /* 0.5 */
81
0
  0.0810614667953272582196702,  /* 1.0 */
82
0
  0.0548141210519176538961390,  /* 1.5 */
83
0
  0.0413406959554092940938221,  /* 2.0 */
84
0
  0.03316287351993628748511048, /* 2.5 */
85
0
  0.02767792568499833914878929, /* 3.0 */
86
0
  0.02374616365629749597132920, /* 3.5 */
87
0
  0.02079067210376509311152277, /* 4.0 */
88
0
  0.01848845053267318523077934, /* 4.5 */
89
0
  0.01664469118982119216319487, /* 5.0 */
90
0
  0.01513497322191737887351255, /* 5.5 */
91
0
  0.01387612882307074799874573, /* 6.0 */
92
0
  0.01281046524292022692424986, /* 6.5 */
93
0
  0.01189670994589177009505572, /* 7.0 */
94
0
  0.01110455975820691732662991, /* 7.5 */
95
0
  0.010411265261972096497478567, /* 8.0 */
96
0
  0.009799416126158803298389475, /* 8.5 */
97
0
  0.009255462182712732917728637, /* 9.0 */
98
0
  0.008768700134139385462952823, /* 9.5 */
99
0
  0.008330563433362871256469318, /* 10.0 */
100
0
  0.007934114564314020547248100, /* 10.5 */
101
0
  0.007573675487951840794972024, /* 11.0 */
102
0
  0.007244554301320383179543912, /* 11.5 */
103
0
  0.006942840107209529865664152, /* 12.0 */
104
0
  0.006665247032707682442354394, /* 12.5 */
105
0
  0.006408994188004207068439631, /* 13.0 */
106
0
  0.006171712263039457647532867, /* 13.5 */
107
0
  0.005951370112758847735624416, /* 14.0 */
108
0
  0.005746216513010115682023589, /* 14.5 */
109
0
  0.005554733551962801371038690  /* 15.0 */
110
0
    };
111
0
    double nn;
112
113
0
    if (n <= 23.5) {
114
0
  nn = n + n;
115
0
  if (n <= 15. && (nn == (int)nn)) return sferr_halves[(int)nn];
116
  // else:
117
0
  if (n <= 5.25) {
118
0
      if(n >= 1.) { // "MM2"; slightly more accurate than direct form
119
0
    double l_n = log(n);        // ldexp(u, -1) == u/2
120
0
    return lgamma(n) + n*(1 - l_n) + ldexp(l_n - M_LN_2PI, -1);
121
0
      }
122
0
      else // n < 1
123
0
    return lgamma1p(n) - (n + 0.5)*log(n) + n - M_LN_SQRT_2PI;
124
0
  }
125
  // else  5.25 < n <= 23.5
126
0
  nn = n*n;
127
0
  if (n > 12.8) return (S0-(S1-(S2-(S3-(S4-(S5 -S6/nn)/nn)/nn)/nn)/nn)/nn)/n;   // k = 7
128
0
  if (n > 12.3) return (S0-(S1-(S2-(S3-(S4-(S5-(S6 -S7/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; // k = 8
129
0
  if (n >  8.9) return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7 -S8/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; // k = 9
130
  /* if (n >  7.9)  return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8 -S9/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; skip k = 10 */
131
0
  if (n >  7.3) return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-S10/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; // 11
132
  /* if (n >  6.5)  return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-S11/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; skip k=12*/
133
0
  if (n >  6.6) return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-S12/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;
134
  /* if (n >  5.7)  return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-S13/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; skip k= 14 */
135
0
  if (n >  6.1) return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-S14/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; // k = 15
136
  /* ....   return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-(S14-S15/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;
137
   * skip order k=16 : never "good" for double prec */
138
  /* 6.1 >= n > 5.25 */
139
0
  return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-(S14-(S15-S16/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n;
140
  /* return (S0-(S1-(S2-(S3-(S4-(S5-(S6-(S7-(S8-(S9-(S10-(S11-(S12-(S13-(S14-(S15-(S16-S17/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/nn)/n; */
141
142
0
     } else { // n > 23.5
143
0
  nn = n*n;
144
0
  if (n > 15.7e6) return S0/n;
145
0
  if (n > 6180) return (S0 -S1/nn)/n;
146
0
  if (n > 205)    return (S0-(S1 -S2/nn)/nn)/n;
147
0
  if (n > 86) return (S0-(S1-(S2 -S3/nn)/nn)/nn)/n;
148
0
  if (n > 27) return (S0-(S1-(S2-(S3 -S4/nn)/nn)/nn)/nn)/n;
149
  /* 23.5 < n <= 27 */
150
0
  return (S0-(S1-(S2-(S3-(S4 -S5/nn)/nn)/nn)/nn)/nn)/n;
151
0
    }
152
0
}