From 25d894a2f67b9fb988b2ada18d35b4de6015c0b8 Mon Sep 17 00:00:00 2001
From: Vinzent Steinberg <Vinzent.Steinberg@gmail.com>
Date: Wed, 22 May 2019 15:40:48 +0200
Subject: [PATCH] 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.
---
 Cargo.toml         |  1 -
 src/lib.rs         |  1 -
 src/moments/mod.rs | 48 +++++++++++++++++++++++++++++++++++++++++-----
 3 files changed, 43 insertions(+), 7 deletions(-)

diff --git a/Cargo.toml b/Cargo.toml
index b67c82e..c844ff1 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -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 }
diff --git a/src/lib.rs b/src/lib.rs
index 420ac05..83e60c7 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -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;
diff --git a/src/moments/mod.rs b/src/moments/mod.rs
index 95d65c8..0ef4bf3 100644
--- a/src/moments/mod.rs
+++ b/src/moments/mod.rs
@@ -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;