Use custom implementation of binomial coefficient

This lets us get rid of the num-integer dependency and makes the
performance of the code generated by `define_moments` close to that of
`Kurtosis`. Before, it was several times slower.

However, the custom implementation is more vulnerable to integer
overflow. In practise, this should not matter, since it does not make
sense to calculate moments of very high order.
This commit is contained in:
Vinzent Steinberg 2019-05-22 15:40:48 +02:00
parent 5b8ace3fb4
commit 25d894a2f6
3 changed files with 43 additions and 7 deletions

View File

@ -27,7 +27,6 @@ name = "kurtosis"
[dependencies]
num-traits = "0.2"
num-integer = "0.1"
float-ord = "0.2"
serde = { version = "1", optional = true }
serde_derive = { version = "1", optional = true }

View File

@ -99,7 +99,6 @@ extern crate serde;
#[cfg(feature = "serde1")]
#[macro_use] extern crate serde_big_array;
extern crate num_traits;
extern crate num_integer;
#[macro_use] mod macros;
#[macro_use] mod moments;

View File

@ -14,9 +14,12 @@ pub type MeanWithError = Variance;
/// Define an estimator of all moments up to a number given at compile time.
///
/// This uses a [general algorithm][paper] and is less efficient than the
/// specialized implementations (such as [`Mean`], [`Variance`], [`Skewness`]
/// and [`Kurtosis`]), but it works for any number of moments >= 4.
/// This uses a [general algorithm][paper] and is slightly less efficient than
/// the specialized implementations (such as [`Mean`], [`Variance`],
/// [`Skewness`] and [`Kurtosis`]), but it works for any number of moments >= 4.
///
/// (In practise, there is an upper limit due to integer overflow and possibly
/// numerical issues.)
///
/// [paper]: https://doi.org/10.1007/s00180-015-0637-z.
/// [`Mean`]: ./struct.Mean.html
@ -30,7 +33,6 @@ pub type MeanWithError = Variance;
/// ```
/// # extern crate core;
/// # extern crate conv;
/// # extern crate num_integer;
/// # extern crate num_traits;
/// #[cfg(feature = "serde1")]
/// extern crate serde;
@ -61,7 +63,43 @@ macro_rules! define_moments {
($name:ident, $MAX_MOMENT:expr) => (
use ::conv::ApproxFrom;
use ::num_traits::pow;
use ::num_integer::IterBinomial;
/// An iterator over binomial coefficients.
struct IterBinomial {
a: u64,
n: u64,
k: u64,
}
impl IterBinomial {
/// For a given n, iterate over all binomial coefficients binomial(n, k), for k=0...n.
#[inline]
pub fn new(n: u64) -> IterBinomial {
IterBinomial {
k: 0,
a: 1,
n: n,
}
}
}
impl Iterator for IterBinomial {
type Item = u64;
#[inline]
fn next(&mut self) -> Option<u64> {
if self.k > self.n {
return None;
}
self.a = if !(self.k == 0) {
self.a * (self.n - self.k + 1) / self.k
} else {
1
};
self.k += 1;
Some(self.a)
}
}
/// The maximal order of the moment to be calculated.
const MAX_MOMENT: usize = $MAX_MOMENT;