From 31aad1431af3c38752eb7a473aa50315f2cf44f8 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 14:27:25 +0200 Subject: [PATCH 1/9] rand_distr: Check bounds before returning Poisson sample Before, we might trigger undefined behaviour if the sample gets too large. --- rand_distr/src/poisson.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index 96bc8b9b6b4..e32641eaebf 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -101,6 +101,8 @@ impl Distribution for Poisson { // now the result is a random variable greater than 0 with Cauchy distribution // the result should be an integer value result = result.floor(); + assert!(result >= 0.); + assert!(result <= ::core::u64::MAX as f64); int_result = result as u64; // this is the ratio of the Poisson distribution to the comparison distribution From c03e2c861339c005372d0fd2c4222fea012c1b94 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 14:36:41 +0200 Subject: [PATCH 2/9] rand_distr: Add support for `f64` to `Poisson` Internally, we generate `f64`, so it makes sense to let the user access that result directly, instead of forcing a conversion from `u64`. --- rand_distr/src/poisson.rs | 36 +++++++++++++++++++++--------------- 1 file changed, 21 insertions(+), 15 deletions(-) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index e32641eaebf..3b4799ce5ea 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -24,7 +24,7 @@ use crate::utils::log_gamma; /// use rand_distr::{Poisson, Distribution}; /// /// let poi = Poisson::new(2.0).unwrap(); -/// let v = poi.sample(&mut rand::thread_rng()); +/// let v: u64 = poi.sample(&mut rand::thread_rng()); /// println!("{} is from a Poisson(2) distribution", v); /// ``` #[derive(Clone, Copy, Debug)] @@ -62,30 +62,28 @@ impl Poisson { } } -impl Distribution for Poisson { - fn sample(&self, rng: &mut R) -> u64 { +impl Distribution for Poisson { + fn sample(&self, rng: &mut R) -> f64 { // using the algorithm from Numerical Recipes in C // for low expected values use the Knuth method if self.lambda < 12.0 { - let mut result = 0; + let mut result = 0.; let mut p = 1.0; while p > self.exp_lambda { p *= rng.gen::(); - result += 1; + result += 1.; } - result - 1 + result - 1. } // high expected values - rejection method else { - let mut int_result: u64; - // we use the Cauchy distribution as the comparison distribution // f(x) ~ 1/(1+x^2) let cauchy = Cauchy::new(0.0, 1.0).unwrap(); + let mut result; loop { - let mut result; let mut comp_dev; loop { @@ -101,9 +99,6 @@ impl Distribution for Poisson { // now the result is a random variable greater than 0 with Cauchy distribution // the result should be an integer value result = result.floor(); - assert!(result >= 0.); - assert!(result <= ::core::u64::MAX as f64); - int_result = result as u64; // this is the ratio of the Poisson distribution to the comparison distribution // the magic value scales the distribution function to a range of approximately 0-1 @@ -117,11 +112,20 @@ impl Distribution for Poisson { break; } } - int_result + result } } } +impl Distribution for Poisson { + fn sample(&self, rng: &mut R) -> u64 { + let result: f64 = self.sample(rng); + assert!(result >= 0.); + assert!(result <= ::core::u64::MAX as f64); + result as u64 + } +} + #[cfg(test)] mod test { use crate::Distribution; @@ -133,7 +137,8 @@ mod test { let mut rng = crate::test::rng(123); let mut sum = 0; for _ in 0..1000 { - sum += poisson.sample(&mut rng); + let s: u64 = poisson.sample(&mut rng); + sum += s; } let avg = (sum as f64) / 1000.0; println!("Poisson average: {}", avg); @@ -147,7 +152,8 @@ mod test { let mut rng = crate::test::rng(123); let mut sum = 0; for _ in 0..1000 { - sum += poisson.sample(&mut rng); + let s: u64 = poisson.sample(&mut rng); + sum += s; } let avg = (sum as f64) / 1000.0; println!("Poisson average: {}", avg); From 638b6be0a585154d5119db8003baf88e9ce944c0 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 15:57:59 +0200 Subject: [PATCH 3/9] rand_distr: Make Poisson generic It now supports `f32` as well. --- rand_distr/src/poisson.rs | 64 ++++++++++++++++++++------------------- rand_distr/src/utils.rs | 48 ++++++++++++++++++++--------- 2 files changed, 67 insertions(+), 45 deletions(-) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index 3b4799ce5ea..c6351db11e1 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -10,8 +10,8 @@ //! The Poisson distribution. use rand::Rng; -use crate::{Distribution, Cauchy}; -use crate::utils::log_gamma; +use crate::{Distribution, Cauchy, Standard}; +use crate::utils::{log_gamma, Float}; /// The Poisson distribution `Poisson(lambda)`. /// @@ -28,13 +28,13 @@ use crate::utils::log_gamma; /// println!("{} is from a Poisson(2) distribution", v); /// ``` #[derive(Clone, Copy, Debug)] -pub struct Poisson { - lambda: f64, +pub struct Poisson { + lambda: N, // precalculated values - exp_lambda: f64, - log_lambda: f64, - sqrt_2lambda: f64, - magic_val: f64, + exp_lambda: N, + log_lambda: N, + sqrt_2lambda: N, + magic_val: N, } /// Error type returned from `Poisson::new`. @@ -44,11 +44,11 @@ pub enum Error { ShapeTooSmall, } -impl Poisson { +impl Poisson { /// Construct a new `Poisson` with the given shape parameter /// `lambda`. - pub fn new(lambda: f64) -> Result { - if !(lambda > 0.0) { + pub fn new(lambda: N) -> Result, Error> { + if !(lambda > N::from(0.0)) { return Err(Error::ShapeTooSmall); } let log_lambda = lambda.ln(); @@ -56,31 +56,33 @@ impl Poisson { lambda, exp_lambda: (-lambda).exp(), log_lambda, - sqrt_2lambda: (2.0 * lambda).sqrt(), - magic_val: lambda * log_lambda - log_gamma(1.0 + lambda), + sqrt_2lambda: (N::from(2.0) * lambda).sqrt(), + magic_val: lambda * log_lambda - log_gamma(N::from(1.0) + lambda), }) } } -impl Distribution for Poisson { - fn sample(&self, rng: &mut R) -> f64 { +impl Distribution for Poisson +where Standard: Distribution +{ + fn sample(&self, rng: &mut R) -> N { // using the algorithm from Numerical Recipes in C // for low expected values use the Knuth method - if self.lambda < 12.0 { - let mut result = 0.; - let mut p = 1.0; + if self.lambda < N::from(12.0) { + let mut result = N::from(0.); + let mut p = N::from(1.0); while p > self.exp_lambda { - p *= rng.gen::(); - result += 1.; + p *= rng.gen::(); + result += N::from(1.); } - result - 1. + result - N::from(1.) } // high expected values - rejection method else { // we use the Cauchy distribution as the comparison distribution // f(x) ~ 1/(1+x^2) - let cauchy = Cauchy::new(0.0, 1.0).unwrap(); + let cauchy = Cauchy::new(N::from(0.0), N::from(1.0)).unwrap(); let mut result; loop { @@ -92,7 +94,7 @@ impl Distribution for Poisson { // shift the peak of the comparison ditribution result = self.sqrt_2lambda * comp_dev + self.lambda; // repeat the drawing until we are in the range of possible values - if result >= 0.0 { + if result >= N::from(0.0) { break; } } @@ -104,11 +106,11 @@ impl Distribution for Poisson { // the magic value scales the distribution function to a range of approximately 0-1 // since it is not exact, we multiply the ratio by 0.9 to avoid ratios greater than 1 // this doesn't change the resulting distribution, only increases the rate of failed drawings - let check = 0.9 * (1.0 + comp_dev * comp_dev) - * (result * self.log_lambda - log_gamma(1.0 + result) - self.magic_val).exp(); + let check = N::from(0.9) * (N::from(1.0) + comp_dev * comp_dev) + * (result * self.log_lambda - log_gamma(N::from(1.0) + result) - self.magic_val).exp(); // check with uniform random value - if below the threshold, we are within the target distribution - if rng.gen::() <= check { + if rng.gen::() <= check { break; } } @@ -117,12 +119,12 @@ impl Distribution for Poisson { } } -impl Distribution for Poisson { +impl Distribution for Poisson +where Standard: Distribution +{ fn sample(&self, rng: &mut R) -> u64 { - let result: f64 = self.sample(rng); - assert!(result >= 0.); - assert!(result <= ::core::u64::MAX as f64); - result as u64 + let result: N = self.sample(rng); + result.into_ui().unwrap() } } diff --git a/rand_distr/src/utils.rs b/rand_distr/src/utils.rs index e5f107f23e5..70ed5cfa1c4 100644 --- a/rand_distr/src/utils.rs +++ b/rand_distr/src/utils.rs @@ -33,9 +33,13 @@ pub trait Float: Copy + Sized + cmp::PartialOrd fn pi() -> Self; /// Support approximate representation of a f64 value fn from(x: f64) -> Self; + /// Support converting to an unsigned integer. + fn into_ui(self) -> Option; /// Take the absolute value of self fn abs(self) -> Self; + /// Take the largest integer less than or equal to self + fn floor(self) -> Self; /// Take the exponential of self fn exp(self) -> Self; @@ -53,8 +57,16 @@ pub trait Float: Copy + Sized + cmp::PartialOrd impl Float for f32 { fn pi() -> Self { core::f32::consts::PI } fn from(x: f64) -> Self { x as f32 } + fn into_ui(self) -> Option { + if self >= 0. && self <= ::core::u64::MAX as f32 { + Some(self as u64) + } else { + None + } + } fn abs(self) -> Self { self.abs() } + fn floor(self) -> Self { self.floor() } fn exp(self) -> Self { self.exp() } fn ln(self) -> Self { self.ln() } @@ -67,8 +79,16 @@ impl Float for f32 { impl Float for f64 { fn pi() -> Self { core::f64::consts::PI } fn from(x: f64) -> Self { x } + fn into_ui(self) -> Option { + if self >= 0. && self <= ::core::u64::MAX as f64 { + Some(self as u64) + } else { + None + } + } fn abs(self) -> Self { self.abs() } + fn floor(self) -> Self { self.floor() } fn exp(self) -> Self { self.exp() } fn ln(self) -> Self { self.ln() } @@ -91,33 +111,33 @@ impl Float for f64 { /// `Ag(z)` is an infinite series with coefficients that can be calculated /// ahead of time - we use just the first 6 terms, which is good enough /// for most purposes. -pub(crate) fn log_gamma(x: f64) -> f64 { +pub(crate) fn log_gamma(x: N) -> N { // precalculated 6 coefficients for the first 6 terms of the series - let coefficients: [f64; 6] = [ - 76.18009172947146, - -86.50532032941677, - 24.01409824083091, - -1.231739572450155, - 0.1208650973866179e-2, - -0.5395239384953e-5, + let coefficients: [N; 6] = [ + N::from(76.18009172947146), + N::from(-86.50532032941677), + N::from(24.01409824083091), + N::from(-1.231739572450155), + N::from(0.1208650973866179e-2), + N::from(-0.5395239384953e-5), ]; // (x+0.5)*ln(x+g+0.5)-(x+g+0.5) - let tmp = x + 5.5; - let log = (x + 0.5) * tmp.ln() - tmp; + let tmp = x + N::from(5.5); + let log = (x + N::from(0.5)) * tmp.ln() - tmp; // the first few terms of the series for Ag(x) - let mut a = 1.000000000190015; + let mut a = N::from(1.000000000190015); let mut denom = x; - for coeff in &coefficients { - denom += 1.0; + for &coeff in &coefficients { + denom += N::from(1.0); a += coeff / denom; } // get everything together // a is Ag(x) // 2.5066... is sqrt(2pi) - log + (2.5066282746310005 * a / x).ln() + log + (N::from(2.5066282746310005) * a / x).ln() } /// Sample a random number using the Ziggurat method (specifically the From eec9bba1e737e2f05362cabcbbd756b6d0b24a5c Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 16:08:34 +0200 Subject: [PATCH 4/9] rand_distr: Add tests for f64 sampling from Poisson --- rand_distr/src/poisson.rs | 36 ++++++++++++++++++++++++------------ 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index c6351db11e1..7acd08d0cd4 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -137,14 +137,20 @@ mod test { fn test_poisson_10() { let poisson = Poisson::new(10.0).unwrap(); let mut rng = crate::test::rng(123); - let mut sum = 0; + let mut sum_u64 = 0; + let mut sum_f64 = 0.; for _ in 0..1000 { - let s: u64 = poisson.sample(&mut rng); - sum += s; + let s_u64: u64 = poisson.sample(&mut rng); + let s_f64: f64 = poisson.sample(&mut rng); + sum_u64 += s_u64; + sum_f64 += s_f64; + } + let avg_u64 = (sum_u64 as f64) / 1000.0; + let avg_f64 = sum_f64 / 1000.0; + println!("Poisson averages: {} (u64) {} (f64)", avg_u64, avg_f64); + for &avg in &[avg_u64, avg_f64] { + assert!((avg - 10.0).abs() < 0.5); // not 100% certain, but probable enough } - let avg = (sum as f64) / 1000.0; - println!("Poisson average: {}", avg); - assert!((avg - 10.0).abs() < 0.5); // not 100% certain, but probable enough } #[test] @@ -152,14 +158,20 @@ mod test { // Take the 'high expected values' path let poisson = Poisson::new(15.0).unwrap(); let mut rng = crate::test::rng(123); - let mut sum = 0; + let mut sum_u64 = 0; + let mut sum_f64 = 0.; for _ in 0..1000 { - let s: u64 = poisson.sample(&mut rng); - sum += s; + let s_u64: u64 = poisson.sample(&mut rng); + let s_f64: f64 = poisson.sample(&mut rng); + sum_u64 += s_u64; + sum_f64 += s_f64; + } + let avg_u64 = (sum_u64 as f64) / 1000.0; + let avg_f64 = sum_f64 / 1000.0; + println!("Poisson average: {} (u64) {} (f64)", avg_u64, avg_f64); + for &avg in &[avg_u64, avg_f64] { + assert!((avg - 15.0).abs() < 0.5); // not 100% certain, but probable enough } - let avg = (sum as f64) / 1000.0; - println!("Poisson average: {}", avg); - assert!((avg - 15.0).abs() < 0.5); // not 100% certain, but probable enough } #[test] From 90833cae3123cbcab8820815b7838ece07738723 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 16:13:57 +0200 Subject: [PATCH 5/9] rand_distr: Add tests for f32 sampling from Poisson --- rand_distr/src/poisson.rs | 41 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index 7acd08d0cd4..069767c57bd 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -174,6 +174,47 @@ mod test { } } + #[test] + fn test_poisson_10_f32() { + let poisson = Poisson::new(10.0f32).unwrap(); + let mut rng = crate::test::rng(123); + let mut sum_u64 = 0; + let mut sum_f32 = 0.; + for _ in 0..1000 { + let s_u64: u64 = poisson.sample(&mut rng); + let s_f32: f32 = poisson.sample(&mut rng); + sum_u64 += s_u64; + sum_f32 += s_f32; + } + let avg_u64 = (sum_u64 as f32) / 1000.0; + let avg_f32 = sum_f32 / 1000.0; + println!("Poisson averages: {} (u64) {} (f32)", avg_u64, avg_f32); + for &avg in &[avg_u64, avg_f32] { + assert!((avg - 10.0).abs() < 0.5); // not 100% certain, but probable enough + } + } + + #[test] + fn test_poisson_15_f32() { + // Take the 'high expected values' path + let poisson = Poisson::new(15.0f32).unwrap(); + let mut rng = crate::test::rng(123); + let mut sum_u64 = 0; + let mut sum_f32 = 0.; + for _ in 0..1000 { + let s_u64: u64 = poisson.sample(&mut rng); + let s_f32: f32 = poisson.sample(&mut rng); + sum_u64 += s_u64; + sum_f32 += s_f32; + } + let avg_u64 = (sum_u64 as f32) / 1000.0; + let avg_f32 = sum_f32 / 1000.0; + println!("Poisson average: {} (u64) {} (f32)", avg_u64, avg_f32); + for &avg in &[avg_u64, avg_f32] { + assert!((avg - 15.0).abs() < 0.5); // not 100% certain, but probable enough + } + } + #[test] #[should_panic] fn test_poisson_invalid_lambda_zero() { From 2553cb5ddaf4f9663888e8815db808240ea71eb6 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 16:18:33 +0200 Subject: [PATCH 6/9] rand_distr: Encourage inlining of Float methods --- rand_distr/src/utils.rs | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/rand_distr/src/utils.rs b/rand_distr/src/utils.rs index 70ed5cfa1c4..3d540a9ad1a 100644 --- a/rand_distr/src/utils.rs +++ b/rand_distr/src/utils.rs @@ -55,8 +55,11 @@ pub trait Float: Copy + Sized + cmp::PartialOrd } impl Float for f32 { + #[inline] fn pi() -> Self { core::f32::consts::PI } + #[inline] fn from(x: f64) -> Self { x as f32 } + #[inline] fn into_ui(self) -> Option { if self >= 0. && self <= ::core::u64::MAX as f32 { Some(self as u64) @@ -65,20 +68,30 @@ impl Float for f32 { } } + #[inline] fn abs(self) -> Self { self.abs() } + #[inline] fn floor(self) -> Self { self.floor() } + #[inline] fn exp(self) -> Self { self.exp() } + #[inline] fn ln(self) -> Self { self.ln() } + #[inline] fn sqrt(self) -> Self { self.sqrt() } + #[inline] fn powf(self, power: Self) -> Self { self.powf(power) } + #[inline] fn tan(self) -> Self { self.tan() } } impl Float for f64 { + #[inline] fn pi() -> Self { core::f64::consts::PI } + #[inline] fn from(x: f64) -> Self { x } + #[inline] fn into_ui(self) -> Option { if self >= 0. && self <= ::core::u64::MAX as f64 { Some(self as u64) @@ -87,14 +100,21 @@ impl Float for f64 { } } + #[inline] fn abs(self) -> Self { self.abs() } + #[inline] fn floor(self) -> Self { self.floor() } + #[inline] fn exp(self) -> Self { self.exp() } + #[inline] fn ln(self) -> Self { self.ln() } + #[inline] fn sqrt(self) -> Self { self.sqrt() } + #[inline] fn powf(self, power: Self) -> Self { self.powf(power) } + #[inline] fn tan(self) -> Self { self.tan() } } From 15b9a39ee2968521a2ee4eb66fe681adda1f3bd4 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 17:02:31 +0200 Subject: [PATCH 7/9] Address review feedback --- rand_distr/src/poisson.rs | 6 ++++-- rand_distr/src/utils.rs | 6 +++--- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index 069767c57bd..500cc82c088 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -44,7 +44,9 @@ pub enum Error { ShapeTooSmall, } -impl Poisson { +impl Poisson +where Standard: Distribution +{ /// Construct a new `Poisson` with the given shape parameter /// `lambda`. pub fn new(lambda: N) -> Result, Error> { @@ -124,7 +126,7 @@ where Standard: Distribution { fn sample(&self, rng: &mut R) -> u64 { let result: N = self.sample(rng); - result.into_ui().unwrap() + result.to_u64().unwrap() } } diff --git a/rand_distr/src/utils.rs b/rand_distr/src/utils.rs index 3d540a9ad1a..b0f3105ec95 100644 --- a/rand_distr/src/utils.rs +++ b/rand_distr/src/utils.rs @@ -34,7 +34,7 @@ pub trait Float: Copy + Sized + cmp::PartialOrd /// Support approximate representation of a f64 value fn from(x: f64) -> Self; /// Support converting to an unsigned integer. - fn into_ui(self) -> Option; + fn to_u64(self) -> Option; /// Take the absolute value of self fn abs(self) -> Self; @@ -60,7 +60,7 @@ impl Float for f32 { #[inline] fn from(x: f64) -> Self { x as f32 } #[inline] - fn into_ui(self) -> Option { + fn to_u64(self) -> Option { if self >= 0. && self <= ::core::u64::MAX as f32 { Some(self as u64) } else { @@ -92,7 +92,7 @@ impl Float for f64 { #[inline] fn from(x: f64) -> Self { x } #[inline] - fn into_ui(self) -> Option { + fn to_u64(self) -> Option { if self >= 0. && self <= ::core::u64::MAX as f64 { Some(self as u64) } else { From 84d89a7c77206d92eccff1b70c42f7074d6eea3b Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Wed, 15 May 2019 17:20:42 +0200 Subject: [PATCH 8/9] Always calculate log_gamma using f64 --- rand_distr/src/poisson.rs | 6 +++--- rand_distr/src/utils.rs | 37 ++++++++++++++++++++++++------------- 2 files changed, 27 insertions(+), 16 deletions(-) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index 500cc82c088..8f5d287374c 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -11,7 +11,7 @@ use rand::Rng; use crate::{Distribution, Cauchy, Standard}; -use crate::utils::{log_gamma, Float}; +use crate::utils::Float; /// The Poisson distribution `Poisson(lambda)`. /// @@ -59,7 +59,7 @@ where Standard: Distribution exp_lambda: (-lambda).exp(), log_lambda, sqrt_2lambda: (N::from(2.0) * lambda).sqrt(), - magic_val: lambda * log_lambda - log_gamma(N::from(1.0) + lambda), + magic_val: lambda * log_lambda - (N::from(1.0) + lambda).log_gamma(), }) } } @@ -109,7 +109,7 @@ where Standard: Distribution // since it is not exact, we multiply the ratio by 0.9 to avoid ratios greater than 1 // this doesn't change the resulting distribution, only increases the rate of failed drawings let check = N::from(0.9) * (N::from(1.0) + comp_dev * comp_dev) - * (result * self.log_lambda - log_gamma(N::from(1.0) + result) - self.magic_val).exp(); + * (result * self.log_lambda - (N::from(1.0) + result).log_gamma() - self.magic_val).exp(); // check with uniform random value - if below the threshold, we are within the target distribution if rng.gen::() <= check { diff --git a/rand_distr/src/utils.rs b/rand_distr/src/utils.rs index b0f3105ec95..c275d0e5cba 100644 --- a/rand_distr/src/utils.rs +++ b/rand_distr/src/utils.rs @@ -52,6 +52,8 @@ pub trait Float: Copy + Sized + cmp::PartialOrd /// Take the tangent of self fn tan(self) -> Self; + /// Take the logarithm of the gamma function of self + fn log_gamma(self) -> Self; } impl Float for f32 { @@ -84,6 +86,13 @@ impl Float for f32 { #[inline] fn tan(self) -> Self { self.tan() } + #[inline] + fn log_gamma(self) -> Self { + let result = log_gamma(self as f64); + assert!(result <= ::core::f32::MAX as f64); + assert!(result >= ::core::f32::MIN as f64); + result as f32 + } } impl Float for f64 { @@ -116,6 +125,8 @@ impl Float for f64 { #[inline] fn tan(self) -> Self { self.tan() } + #[inline] + fn log_gamma(self) -> Self { log_gamma(self) } } /// Calculates ln(gamma(x)) (natural logarithm of the gamma @@ -131,33 +142,33 @@ impl Float for f64 { /// `Ag(z)` is an infinite series with coefficients that can be calculated /// ahead of time - we use just the first 6 terms, which is good enough /// for most purposes. -pub(crate) fn log_gamma(x: N) -> N { +pub(crate) fn log_gamma(x: f64) -> f64 { // precalculated 6 coefficients for the first 6 terms of the series - let coefficients: [N; 6] = [ - N::from(76.18009172947146), - N::from(-86.50532032941677), - N::from(24.01409824083091), - N::from(-1.231739572450155), - N::from(0.1208650973866179e-2), - N::from(-0.5395239384953e-5), + let coefficients: [f64; 6] = [ + 76.18009172947146, + -86.50532032941677, + 24.01409824083091, + -1.231739572450155, + 0.1208650973866179e-2, + -0.5395239384953e-5, ]; // (x+0.5)*ln(x+g+0.5)-(x+g+0.5) - let tmp = x + N::from(5.5); - let log = (x + N::from(0.5)) * tmp.ln() - tmp; + let tmp = x + 5.5; + let log = (x + 0.5) * tmp.ln() - tmp; // the first few terms of the series for Ag(x) - let mut a = N::from(1.000000000190015); + let mut a = 1.000000000190015; let mut denom = x; for &coeff in &coefficients { - denom += N::from(1.0); + denom += 1.0; a += coeff / denom; } // get everything together // a is Ag(x) // 2.5066... is sqrt(2pi) - log + (N::from(2.5066282746310005) * a / x).ln() + log + (2.5066282746310005 * a / x).ln() } /// Sample a random number using the Ziggurat method (specifically the From ec99801390226f31c3e3bc53264147ef484a9702 Mon Sep 17 00:00:00 2001 From: Vinzent Steinberg Date: Thu, 16 May 2019 13:46:00 +0200 Subject: [PATCH 9/9] rand_distr: Inline Poisson sampling --- rand_distr/src/poisson.rs | 2 ++ 1 file changed, 2 insertions(+) diff --git a/rand_distr/src/poisson.rs b/rand_distr/src/poisson.rs index 8f5d287374c..4f4a0b7a3d9 100644 --- a/rand_distr/src/poisson.rs +++ b/rand_distr/src/poisson.rs @@ -67,6 +67,7 @@ where Standard: Distribution impl Distribution for Poisson where Standard: Distribution { + #[inline] fn sample(&self, rng: &mut R) -> N { // using the algorithm from Numerical Recipes in C @@ -124,6 +125,7 @@ where Standard: Distribution impl Distribution for Poisson where Standard: Distribution { + #[inline] fn sample(&self, rng: &mut R) -> u64 { let result: N = self.sample(rng); result.to_u64().unwrap()