From 5f6f89e5c7903da9fd6984dfc95300ccd0af630b Mon Sep 17 00:00:00 2001
From: Vinzent Steinberg <Vinzent.Steinberg@gmail.com>
Date: Mon, 25 Feb 2019 14:05:00 +0100
Subject: [PATCH] Binomial: Use integer exponentation

This increase performance significantly, see below. However, it requires
shipping our own integer exponentiation implementation, because `u64`
exponents are not supported by `std` or `num`.

Before:
```
test distr_binomial       ... bench:  88,365 ns/iter (+/- 2,215) =90 MB/s
test distr_binomial_small ... bench:  33,775 ns/iter (+/- 5,494) =236 MB/s
test misc_binomial_1      ... bench:      13 ns/iter (+/- 2)
test misc_binomial_10     ... bench:      72 ns/iter (+/- 1)
test misc_binomial_100    ... bench:      71 ns/iter (+/- 11)
test misc_binomial_1000   ... bench:     624 ns/iter (+/- 19)
test misc_binomial_1e12   ... bench:     681 ns/iter (+/- 18)
```

After:
```
test distr_binomial       ... bench:  44,354 ns/iter (+/- 3,518) = 180 MB/s
test distr_binomial_small ... bench:  22,736 ns/iter (+/- 514) = 351 MB/s
test misc_binomial_1      ... bench:      10 ns/iter (+/- 0)
test misc_binomial_10     ... bench:      23 ns/iter (+/- 0)
test misc_binomial_100    ... bench:      27 ns/iter (+/- 4)
test misc_binomial_1000   ... bench:     621 ns/iter (+/- 15)
test misc_binomial_1e12   ... bench:     686 ns/iter (+/- 20)
```
---
 src/distributions/binomial.rs | 29 ++++++++++++++++++++++++++---
 1 file changed, 26 insertions(+), 3 deletions(-)

diff --git a/src/distributions/binomial.rs b/src/distributions/binomial.rs
index 907155b24ed..2bad03b78bd 100644
--- a/src/distributions/binomial.rs
+++ b/src/distributions/binomial.rs
@@ -47,6 +47,31 @@ impl Binomial {
     }
 }
 
+/// Raise a `base` to the power of `exp`, using exponentiation by squaring.
+///
+/// This implementation is based on the one in the `num_traits` crate. It is
+/// slightly modified to accept `u64` exponents.
+fn pow(mut base: f64, mut exp: u64) -> f64 {
+    if exp == 0 {
+        return 1.;
+    }
+
+    while exp & 1 == 0 {
+        base *= base;
+        exp >>= 1;
+    }
+
+    let mut acc = base;
+    while exp > 1 {
+        exp >>= 1;
+        base *= base;
+        if exp & 1 == 1 {
+            acc *= base;
+        }
+    }
+    acc
+}
+
 impl Distribution<u64> for Binomial {
     fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> u64 {
         // Handle these values directly.
@@ -77,9 +102,7 @@ impl Distribution<u64> for Binomial {
             let q = 1. - p;
             let s = p / q;
             let a = ((self.n + 1) as f64) * s;
-            let mut r = q.powf(self.n as f64);
-            // FIXME: Using integer exponentiation should be faster, but it is not available for
-            // `u64` exponents, neither in `std` or `num`.
+            let mut r = pow(q, self.n);
             let mut u: f64 = rng.gen();
             let mut x = 0;
             while u > r as f64 {