Commit da600570 authored by Jan-Oliver Kaiser's avatar Jan-Oliver Kaiser

Add memoization for log_fac_stirling.

parent 849ad30a
...@@ -10,9 +10,37 @@ pub mod stirling_approximation{ ...@@ -10,9 +10,37 @@ pub mod stirling_approximation{
} }
pub fn log_fac_stirling(n: f64) -> f64{ use std;
let result: f64 = n*n.ln() - n; use atomic_f64::AtomicF64;
result / f64::consts::LN_2
const MAX_LOG_FAC_MEMOIZATION: usize = 1024*1024;
lazy_static!{
static ref LOG_FAC : Box<[AtomicF64; MAX_LOG_FAC_MEMOIZATION]> = {
let mut result : Box<[AtomicF64; MAX_LOG_FAC_MEMOIZATION]> = unsafe { box std::mem::uninitialized() };
for i in 0..MAX_LOG_FAC_MEMOIZATION {
let uninit = std::mem::replace(&mut result[i], AtomicF64::new(0.0));
std::mem::forget(uninit)
};
result
};
}
#[inline]
pub fn log_fac_stirling(n: usize) -> f64{
let mut result : f64;
if n < MAX_LOG_FAC_MEMOIZATION {
let result = LOG_FAC[n-1].get();
if result != 0.0 {
return result
}
}
result = ((n*n) as f64).ln() - (n as f64);
result /= f64::consts::LN_2;
if n < MAX_LOG_FAC_MEMOIZATION {
LOG_FAC[n-1].set(result)
}
result
} }
...@@ -25,8 +53,8 @@ pub mod stirling_approximation{ ...@@ -25,8 +53,8 @@ pub mod stirling_approximation{
if (n == 0) || (k == 0) { if (n == 0) || (k == 0) {
res = f64::INFINITY; res = f64::INFINITY;
} else { } else {
let n = n as f64; // let n = n as f64;
let k = k as f64; // let k = k as f64;
res = log_fac_stirling(n) - log_fac_stirling(k) - log_fac_stirling(n - k); res = log_fac_stirling(n) - log_fac_stirling(k) - log_fac_stirling(n - k);
} }
} }
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment