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
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
#![allow(dead_code)]
#[cfg(test)]
#[macro_use]
extern crate approx;
extern crate num_bigint;
use num_bigint::BigUint;
use num_traits::cast::ToPrimitive;
use num_traits::{One, Zero};
use std::convert::TryFrom;
use std::f64;
pub fn bell(n: usize) -> BigUint {
let mut r1: Vec<BigUint> = vec![Zero::zero(); n];
let mut r2: Vec<BigUint> = vec![Zero::zero(); n];
r1[0] = One::one();
for k in 1..n {
r2[0] = r1[k - 1].clone();
for i in 1..(k + 1) {
r2[i] = r1[i - 1].clone() + &r2[i - 1];
}
let tmp = r1;
r1 = r2;
r2 = tmp;
}
r1[n - 1].clone()
}
pub fn lbell(n: usize) -> f64 {
let value = bell(n);
let n_bits = value.bits();
let threshold = 1022usize;
let log2 = if n_bits > threshold {
let n_shifted_bits = value.bits() - threshold;
let shifted_value = value >> n_shifted_bits;
if shifted_value.bits() > threshold {
return f64::INFINITY;
}
let y: f64 = shifted_value.to_f64().unwrap();
(n_shifted_bits as f64) + y.log2()
} else {
value.to_f64().unwrap().log2()
};
log2 / f64::consts::LOG2_E
}
#[doc(hidden)]
#[no_mangle]
pub extern "C" fn dahl_bellnumber__bell(n: i32) -> f64 {
if n < 0 {
return 0.0;
}
match usize::try_from(n) {
Ok(n) => match bell(n).to_f64() {
Some(x) => x,
None => f64::INFINITY,
},
Err(_) => f64::INFINITY,
}
}
#[doc(hidden)]
#[no_mangle]
pub extern "C" fn dahl_bellnumber__lbell(n: i32) -> f64 {
if n < 0 {
return 0.0;
}
match usize::try_from(n) {
Ok(n) => lbell(n),
Err(_) => f64::INFINITY,
}
}
#[cfg(test)]
mod tests {
use super::*;
#[test]
fn test_lbell() {
relative_eq!(lbell(220), 714.4033);
relative_eq!(bell(5).to_f64().unwrap(), 52.0);
}
}