dahl_bellnumber/
lib.rs

1//! # Bell Number
2//!
3//! `dahl_bellnumber` is a collection of functions related to the [Bell number](https://en.wikipedia.org/wiki/Bell_number),
4//! which gives the number of partitions of a set.
5//!
6//!
7
8#![allow(dead_code)]
9
10#[cfg(test)]
11#[macro_use]
12extern crate approx;
13extern crate num_bigint;
14
15use num_bigint::BigUint;
16use num_traits::cast::ToPrimitive;
17use num_traits::{One, Zero};
18use std::convert::TryFrom;
19use std::f64;
20
21/// Compute the [Bell number](https://en.wikipedia.org/wiki/Bell_number).
22///
23/// # Examples
24///
25/// ```
26/// let answer = dahl_bellnumber::bell(5);
27///
28/// use std::convert::TryFrom;
29/// use num_traits::cast::ToPrimitive;
30///
31/// assert_eq!(answer, num_bigint::BigUint::try_from(52_u32).unwrap());
32/// assert_eq!(answer.to_f64().unwrap(), 52.0);
33/// ```
34pub fn bell(n: usize) -> BigUint {
35    let mut r1: Vec<BigUint> = vec![Zero::zero(); n];
36    let mut r2: Vec<BigUint> = vec![Zero::zero(); n];
37    r1[0] = One::one();
38    for k in 1..n {
39        r2[0] = r1[k - 1].clone();
40        for i in 1..(k + 1) {
41            r2[i] = r1[i - 1].clone() + &r2[i - 1];
42        }
43        let tmp = r1;
44        r1 = r2;
45        r2 = tmp;
46    }
47    r1[n - 1].clone()
48}
49
50/// Compute the natural logarithm of the [Bell number](https://en.wikipedia.org/wiki/Bell_number).
51///
52/// # Examples
53///
54/// ```
55/// let answer = dahl_bellnumber::lbell(5);
56///
57/// assert!( (answer - 52.0_f64.ln()).abs() < 0.00000001 );
58/// ```
59pub fn lbell(n: usize) -> f64 {
60    let value = bell(n);
61    let n_bits = value.bits();
62    let threshold = 1022_u64;
63    let log2 = if n_bits > threshold {
64        let n_shifted_bits = value.bits() - threshold;
65        let shifted_value = value >> n_shifted_bits;
66        if shifted_value.bits() > threshold {
67            return f64::INFINITY;
68        }
69        let y: f64 = shifted_value.to_f64().unwrap();
70        (n_shifted_bits as f64) + y.log2()
71    } else {
72        value.to_f64().unwrap().log2()
73    };
74    log2 / f64::consts::LOG2_E
75}
76
77/// C-friendly wrapper over the `bell` function.
78///
79#[doc(hidden)]
80#[no_mangle]
81pub extern "C" fn dahl_bellnumber__bell(n: i32) -> f64 {
82    if n < 0 {
83        return 0.0;
84    }
85    match usize::try_from(n) {
86        Ok(n) => match bell(n).to_f64() {
87            Some(x) => x,
88            None => f64::INFINITY,
89        },
90        Err(_) => f64::INFINITY,
91    }
92}
93
94/// C-friendly wrapper over the `lbell` function.
95///
96#[doc(hidden)]
97#[no_mangle]
98pub extern "C" fn dahl_bellnumber__lbell(n: i32) -> f64 {
99    if n < 0 {
100        return 0.0;
101    }
102    match usize::try_from(n) {
103        Ok(n) => lbell(n),
104        Err(_) => f64::INFINITY,
105    }
106}
107
108#[cfg(test)]
109mod tests {
110    use super::*;
111
112    #[test]
113    fn test_lbell() {
114        assert_relative_eq!(lbell(220), 714.4032630589774);
115        assert_relative_eq!(bell(5).to_f64().unwrap(), 52.0);
116    }
117
118}