r2rs-nmath 0.1.1

Statistics programming for Rust based on R's nmath package
Documentation 
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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99

// Translation of nmath's stirlerr
//
// Copyright (C) 2020 neek-sss
//
// 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 3 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, see <https://www.gnu.org/licenses/>.

use strafe_type::FloatConstraint;

use crate::func::lgamma;

/// 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!)
/// 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.c  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
pub fn stirlerr(n: f64) -> f64 {
    const S0: f64 = 0.083333333333333333333; /* 1/12 */
    const S1: f64 = 0.00277777777777777777778; /* 1/360 */
    const S2: f64 = 0.00079365079365079365079365; /* 1/1260 */
    const S3: f64 = 0.000595238095238095238095238; /* 1/1680 */
    const S4: f64 = 0.0008417508417508417508417508; /* 1/1188 */

    /*
      exact values for 0, 0.5, 1.0, 1.5, ..., 14.5, 15.0.
    */
    static sferr_halves: [f64; 31] = [
        0.0,                           /* n=0 - wrong, place holder only */
        0.1534264097200273452913848,   /* 0.5 */
        0.0810614667953272582196702,   /* 1.0 */
        0.0548141210519176538961390,   /* 1.5 */
        0.0413406959554092940938221,   /* 2.0 */
        0.03316287351993628748511048,  /* 2.5 */
        0.02767792568499833914878929,  /* 3.0 */
        0.02374616365629749597132920,  /* 3.5 */
        0.02079067210376509311152277,  /* 4.0 */
        0.01848845053267318523077934,  /* 4.5 */
        0.01664469118982119216319487,  /* 5.0 */
        0.01513497322191737887351255,  /* 5.5 */
        0.01387612882307074799874573,  /* 6.0 */
        0.01281046524292022692424986,  /* 6.5 */
        0.01189670994589177009505572,  /* 7.0 */
        0.01110455975820691732662991,  /* 7.5 */
        0.010411265261972096497478567, /* 8.0 */
        0.009799416126158803298389475, /* 8.5 */
        0.009255462182712732917728637, /* 9.0 */
        0.008768700134139385462952823, /* 9.5 */
        0.008330563433362871256469318, /* 10.0 */
        0.007934114564314020547248100, /* 10.5 */
        0.007573675487951840794972024, /* 11.0 */
        0.007244554301320383179543912, /* 11.5 */
        0.006942840107209529865664152, /* 12.0 */
        0.006665247032707682442354394, /* 12.5 */
        0.006408994188004207068439631, /* 13.0 */
        0.006171712263039457647532867, /* 13.5 */
        0.005951370112758847735624416, /* 14.0 */
        0.005746216513010115682023589, /* 14.5 */
        0.005554733551962801371038690, /* 15.0 */
    ];
    let mut nn = 0.0;

    if n <= 15.0 {
        nn = n + n;
        if nn == nn.trunc() {
            return sferr_halves[nn as usize];
        }
        return lgamma(n + 1.0).unwrap() - (n + 0.5) * n.ln() + n - strafe_consts::LN_SQRT_2TPI;
    }

    nn = n * n;
    if n > 500.0 {
        return (S0 - S1 / nn) / n;
    }
    if n > 80.0 {
        return (S0 - (S1 - S2 / nn) / nn) / n;
    }
    if n > 35.0 {
        return (S0 - (S1 - (S2 - S3 / nn) / nn) / nn) / n;
    }
    /* 15 < n <= 35 : */
    return (S0 - (S1 - (S2 - (S3 - S4 / nn) / nn) / nn) / nn) / n;
}