1

2
#include <math.h>
3
#include <stdint.h>
4

5
/*
6
 *  logfact[k] holds log(k!) for k = 0, 1, 2, ..., 125.
7
 */
8

9
static const double logfact[] = {
10
    0,
11
    0,
12
    0.69314718055994529,
13
    1.791759469228055,
14
    3.1780538303479458,
15
    4.7874917427820458,
16
    6.5792512120101012,
17
    8.5251613610654147,
18
    10.604602902745251,
19
    12.801827480081469,
20
    15.104412573075516,
21
    17.502307845873887,
22
    19.987214495661885,
23
    22.552163853123425,
24
    25.19122118273868,
25
    27.89927138384089,
26
    30.671860106080672,
27
    33.505073450136891,
28
    36.395445208033053,
29
    39.339884187199495,
30
    42.335616460753485,
31
    45.380138898476908,
32
    48.471181351835227,
33
    51.606675567764377,
34
    54.784729398112319,
35
    58.003605222980518,
36
    61.261701761002001,
37
    64.557538627006338,
38
    67.88974313718154,
39
    71.257038967168015,
40
    74.658236348830158,
41
    78.092223553315307,
42
    81.557959456115043,
43
    85.054467017581516,
44
    88.580827542197682,
45
    92.136175603687093,
46
    95.719694542143202,
47
    99.330612454787428,
48
    102.96819861451381,
49
    106.63176026064346,
50
    110.32063971475739,
51
    114.03421178146171,
52
    117.77188139974507,
53
    121.53308151543864,
54
    125.3172711493569,
55
    129.12393363912722,
56
    132.95257503561632,
57
    136.80272263732635,
58
    140.67392364823425,
59
    144.5657439463449,
60
    148.47776695177302,
61
    152.40959258449735,
62
    156.3608363030788,
63
    160.3311282166309,
64
    164.32011226319517,
65
    168.32744544842765,
66
    172.35279713916279,
67
    176.39584840699735,
68
    180.45629141754378,
69
    184.53382886144948,
70
    188.6281734236716,
71
    192.7390472878449,
72
    196.86618167289001,
73
    201.00931639928152,
74
    205.1681994826412,
75
    209.34258675253685,
76
    213.53224149456327,
77
    217.73693411395422,
78
    221.95644181913033,
79
    226.1905483237276,
80
    230.43904356577696,
81
    234.70172344281826,
82
    238.97838956183432,
83
    243.26884900298271,
84
    247.57291409618688,
85
    251.89040220972319,
86
    256.22113555000954,
87
    260.56494097186322,
88
    264.92164979855278,
89
    269.29109765101981,
90
    273.67312428569369,
91
    278.06757344036612,
92
    282.4742926876304,
93
    286.89313329542699,
94
    291.32395009427029,
95
    295.76660135076065,
96
    300.22094864701415,
97
    304.68685676566872,
98
    309.1641935801469,
99
    313.65282994987905,
100
    318.1526396202093,
101
    322.66349912672615,
102
    327.1852877037752,
103
    331.71788719692847,
104
    336.26118197919845,
105
    340.81505887079902,
106
    345.37940706226686,
107
    349.95411804077025,
108
    354.53908551944079,
109
    359.1342053695754,
110
    363.73937555556347,
111
    368.35449607240474,
112
    372.97946888568902,
113
    377.61419787391867,
114
    382.25858877306001,
115
    386.91254912321756,
116
    391.57598821732961,
117
    396.24881705179155,
118
    400.93094827891576,
119
    405.6222961611449,
120
    410.32277652693733,
121
    415.03230672824964,
122
    419.75080559954472,
123
    424.47819341825709,
124
    429.21439186665157,
125
    433.95932399501481,
126
    438.71291418612117,
127
    443.47508812091894,
128
    448.24577274538461,
129
    453.02489623849613,
130
    457.81238798127816,
131
    462.60817852687489,
132
    467.4121995716082,
133
    472.22438392698058,
134
    477.04466549258564,
135
    481.87297922988796
136
};
137

138
/*
139
 *  Compute log(k!)
140
 */
141

142 1
double logfactorial(int64_t k)
143
{
144 1
    const double halfln2pi = 0.9189385332046728;
145

146 1
    if (k < (int64_t) (sizeof(logfact)/sizeof(logfact[0]))) {
147
        /* Use the lookup table. */
148 1
        return logfact[k];
149
    }
150

151
    /*
152
     *  Use the Stirling series, truncated at the 1/k**3 term.
153
     *  (In a Python implementation of this approximation, the result
154
     *  was within 2 ULP of the best 64 bit floating point value for
155
     *  k up to 10000000.)
156
     */
157 1
    return (k + 0.5)*log(k) - k + (halfln2pi + (1.0/k)*(1/12.0 - 1/(360.0*k*k)));
158
}

Read our documentation on viewing source code .

Loading