1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
use crate::lgamma::*;
use crate::nmath::*;
const S0: f64 = 0.083_333_333_333_333_33; /* 1/12 */
const S1: f64 = 0.002_777_777_777_777_778; /* 1/360 */
const S2: f64 = 0.000_793_650_793_650_793_7; /* 1/1260 */
const S3: f64 = 0.000_595_238_095_238_095_3; /* 1/1680 */
const S4: f64 = 0.000_841_750_841_750_841_7; /* 1/1188 */
/// exact values for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0.
const SFERR_HALVES: [f64; 31] = [
0.0, /* n=0 - wrong, place holder only */
0.153_426_409_720_027_36, /* 0.5 */
0.081_061_466_795_327_26, /* 1.0 */
0.054_814_121_051_917_65, /* 1.5 */
0.041_340_695_955_409_3, /* 2.0 */
0.033_162_873_519_936_29, /* 2.5 */
0.027_677_925_684_998_34, /* 3.0 */
0.023_746_163_656_297_496, /* 3.5 */
0.020_790_672_103_765_093, /* 4.0 */
0.018_488_450_532_673_187, /* 4.5 */
0.016_644_691_189_821_193, /* 5.0 */
0.015_134_973_221_917_378, /* 5.5 */
0.013_876_128_823_070_748, /* 6.0 */
0.012_810_465_242_920_227, /* 6.5 */
0.011_896_709_945_891_77, /* 7.0 */
0.011_104_559_758_206_917, /* 7.5 */
0.010_411_265_261_972_096, /* 8.0 */
0.009_799_416_126_158_804, /* 8.5 */
0.009_255_462_182_712_733, /* 9.0 */
0.008_768_700_134_139_386, /* 9.5 */
0.008_330_563_433_362_87, /* 10.0 */
0.007_934_114_564_314_02, /* 10.5 */
0.007_573_675_487_951_841, /* 11.0 */
0.007_244_554_301_320_383, /* 11.5 */
0.006_942_840_107_209_53, /* 12.0 */
0.006_665_247_032_707_682, /* 12.5 */
0.006_408_994_188_004_207, /* 13.0 */
0.006_171_712_263_039_458, /* 13.5 */
0.005_951_370_112_758_847_5, /* 14.0 */
0.005_746_216_513_010_115_5, /* 14.5 */
0.005_554_733_551_962_801, /* 15.0 */
];
/// Computes the log of the error term in Stirling's formula.
///
/// For n > 15, uses the series 1/12n - 1/360n^3 + ...
/// For n <=15, integers or half-integers, uses stored values.
/// For other n < 15, uses lgamma directly (don't use this to
/// write lgamma!)
pub fn stirlerr(n: f64) -> f64 {
// stirlerr(n) = log(n!) - log( sqrt(2*pi*n)*(n/e)^n )
// = log Gamma(n+1) - 1/2 * [log(2*pi) + log(n)] - n*[log(n) - 1]
// = log Gamma(n+1) - (n + 1/2) * log(n) + n - log(2*pi)/2
//
// see also lgammacor() in ./lgammacor.rs which computes almost the same!
//
// NB: stirlerr(n/2) is called from dt() *and* gamma(n/2) when n is integer and n/2 <= 50
if n <= 15.0 {
let nn = n + n;
if nn == nn.floor() {
return SFERR_HALVES[nn as usize];
}
lgammafn(n + 1.0) - (n + 0.5) * n.ln() + n - M_LN_SQRT_2PI
} else {
let nn = n * n;
if n > 500.0 {
(S0 - S1 / nn) / n
} else if n > 80.0 {
(S0 - (S1 - S2 / nn) / nn) / n
} else if n > 35.0 {
(S0 - (S1 - (S2 - S3 / nn) / nn) / nn) / n
} else {
// 15 < n <= 35
(S0 - (S1 - (S2 - (S3 - S4 / nn) / nn) / nn) / nn) / n
}
}
}