1 | module math |
2 | |
3 | // factorial calculates the factorial of the provided value. |
4 | pub fn factorial(n f64) f64 { |
5 | // For a large positive argument (n >= factorials_table.len) return max_f64 |
6 | if n >= factorials_table.len { |
7 | return max_f64 |
8 | } |
9 | // Otherwise return n!. |
10 | if n == f64(i64(n)) && n >= 0.0 { |
11 | return factorials_table[i64(n)] |
12 | } |
13 | return gamma(n + 1.0) |
14 | } |
15 | |
16 | // log_factorial calculates the log-factorial of the provided value. |
17 | pub fn log_factorial(n f64) f64 { |
18 | // For a large positive argument (n < 0) return max_f64 |
19 | if n < 0 { |
20 | return -max_f64 |
21 | } |
22 | // If n < N then return ln(n!). |
23 | if n != f64(i64(n)) { |
24 | return log_gamma(n + 1) |
25 | } else if n < log_factorials_table.len { |
26 | return log_factorials_table[i64(n)] |
27 | } |
28 | // Otherwise return asymptotic expansion of ln(n!). |
29 | return log_factorial_asymptotic_expansion(int(n)) |
30 | } |
31 | |
32 | fn log_factorial_asymptotic_expansion(n int) f64 { |
33 | m := 6 |
34 | mut term := []f64{} |
35 | xx := f64((n + 1) * (n + 1)) |
36 | mut xj := f64(n + 1) |
37 | log_factorial := log_sqrt_2pi - xj + (xj - 0.5) * log(xj) |
38 | mut i := 0 |
39 | for i = 0; i < m; i++ { |
40 | term << bernoulli[i] / xj |
41 | xj *= xx |
42 | } |
43 | mut sum := term[m - 1] |
44 | for i = m - 2; i >= 0; i-- { |
45 | if abs(sum) <= abs(term[i]) { |
46 | break |
47 | } |
48 | sum = term[i] |
49 | } |
50 | for i >= 0 { |
51 | sum += term[i] |
52 | i-- |
53 | } |
54 | return log_factorial + sum |
55 | } |
56 | |
57 | // factoriali returns 1 for n <= 0 and -1 if the result is too large for a 64 bit integer |
58 | pub fn factoriali(n int) i64 { |
59 | if n <= 0 { |
60 | return i64(1) |
61 | } |
62 | |
63 | if n < 21 { |
64 | return i64(factorials_table[n]) |
65 | } |
66 | |
67 | return i64(-1) |
68 | } |