Update nitrokey crate to 0.4.0

This change finally updates the version of the nitrokey crate that we
consume to 0.4.0. Along with that we update rand_core, one of its
dependencies, to 0.5.1. Further more we add cfg-if in version 0.1.10 and
getrandom in version 0.1.13, both of which are now new (non-development)
dependencies.

Import subrepo nitrokey/:nitrokey at e81057037e9b4f370b64c0a030a725bc6bdfb870
Import subrepo cfg-if/:cfg-if at 4484a6faf816ff8058088ad857b0c6bb2f4b02b2
Import subrepo getrandom/:getrandom at d661aa7e1b8cc80b47dabe3d2135b3b47d2858af
Import subrepo rand/:rand at d877ed528248b52d947e0484364a4e1ae59ca502

1 files changed, 234 insertions, 0 deletions

diff --git a/rand/rand_distr/src/utils.rs b/rand/rand_distr/src/utils.rs
new file mode 100644
index 0000000..75b3500
--- /dev/null
+++ b/rand/rand_distr/src/utils.rs
@@ -0,0 +1,234 @@
+// Copyright 2018 Developers of the Rand project.
+//
+// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
+// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
+// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
+// option. This file may not be copied, modified, or distributed
+// except according to those terms.
+
+//! Math helper functions
+
+use rand::Rng;
+use crate::ziggurat_tables;
+use rand::distributions::hidden_export::IntoFloat;
+use core::{cmp, ops};
+
+/// Trait for floating-point scalar types
+/// 
+/// This allows many distributions to work with `f32` or `f64` parameters and is
+/// potentially extensible. Note however that the `Exp1` and `StandardNormal`
+/// distributions are implemented exclusively for `f32` and `f64`.
+/// 
+/// The bounds and methods are based purely on internal
+/// requirements, and will change as needed.
+pub trait Float: Copy + Sized + cmp::PartialOrd
+    + ops::Neg<Output = Self>
+    + ops::Add<Output = Self>
+    + ops::Sub<Output = Self>
+    + ops::Mul<Output = Self>
+    + ops::Div<Output = Self>
+    + ops::AddAssign + ops::SubAssign + ops::MulAssign + ops::DivAssign
+{
+    /// The constant π
+    fn pi() -> Self;
+    /// Support approximate representation of a f64 value
+    fn from(x: f64) -> Self;
+    /// Support converting to an unsigned integer.
+    fn to_u64(self) -> Option<u64>;
+    
+    /// 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;
+    /// Take the natural logarithm of self
+    fn ln(self) -> Self;
+    /// Take square root of self
+    fn sqrt(self) -> Self;
+    /// Take self to a floating-point power
+    fn powf(self, power: Self) -> Self;
+    
+    /// 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 {
+    #[inline]
+    fn pi() -> Self { core::f32::consts::PI }
+    #[inline]
+    fn from(x: f64) -> Self { x as f32 }
+    #[inline]
+    fn to_u64(self) -> Option<u64> {
+        if self >= 0. && self <= ::core::u64::MAX as f32 {
+            Some(self as u64)
+        } else {
+            None
+        }
+    }
+    
+    #[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() }
+    #[inline]
+    fn log_gamma(self) -> Self {
+        let result = log_gamma(self.into());
+        assert!(result <= ::core::f32::MAX.into());
+        assert!(result >= ::core::f32::MIN.into());
+        result as f32
+    }
+}
+
+impl Float for f64 {
+    #[inline]
+    fn pi() -> Self { core::f64::consts::PI }
+    #[inline]
+    fn from(x: f64) -> Self { x }
+    #[inline]
+    fn to_u64(self) -> Option<u64> {
+        if self >= 0. && self <= ::core::u64::MAX as f64 {
+            Some(self as u64)
+        } else {
+            None
+        }
+    }
+    
+    #[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() }
+    #[inline]
+    fn log_gamma(self) -> Self { log_gamma(self) }
+}
+
+/// Calculates ln(gamma(x)) (natural logarithm of the gamma
+/// function) using the Lanczos approximation.
+///
+/// The approximation expresses the gamma function as:
+/// `gamma(z+1) = sqrt(2*pi)*(z+g+0.5)^(z+0.5)*exp(-z-g-0.5)*Ag(z)`
+/// `g` is an arbitrary constant; we use the approximation with `g=5`.
+///
+/// Noting that `gamma(z+1) = z*gamma(z)` and applying `ln` to both sides:
+/// `ln(gamma(z)) = (z+0.5)*ln(z+g+0.5)-(z+g+0.5) + ln(sqrt(2*pi)*Ag(z)/z)`
+///
+/// `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 {
+    // 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,
+    ];
+
+    // (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;
+
+    // the first few terms of the series for Ag(x)
+    let mut a = 1.000000000190015;
+    let mut denom = x;
+    for &coeff in &coefficients {
+        denom += 1.0;
+        a += coeff / denom;
+    }
+
+    // get everything together
+    // a is Ag(x)
+    // 2.5066... is sqrt(2pi)
+    log + (2.5066282746310005 * a / x).ln()
+}
+
+/// Sample a random number using the Ziggurat method (specifically the
+/// ZIGNOR variant from Doornik 2005). Most of the arguments are
+/// directly from the paper:
+///
+/// * `rng`: source of randomness
+/// * `symmetric`: whether this is a symmetric distribution, or one-sided with P(x < 0) = 0.
+/// * `X`: the $x_i$ abscissae.
+/// * `F`: precomputed values of the PDF at the $x_i$, (i.e. $f(x_i)$)
+/// * `F_DIFF`: precomputed values of $f(x_i) - f(x_{i+1})$
+/// * `pdf`: the probability density function
+/// * `zero_case`: manual sampling from the tail when we chose the
+///    bottom box (i.e. i == 0)
+
+// the perf improvement (25-50%) is definitely worth the extra code
+// size from force-inlining.
+#[inline(always)]
+pub(crate) fn ziggurat<R: Rng + ?Sized, P, Z>(
+            rng: &mut R,
+            symmetric: bool,
+            x_tab: ziggurat_tables::ZigTable,
+            f_tab: ziggurat_tables::ZigTable,
+            mut pdf: P,
+            mut zero_case: Z)
+            -> f64 where P: FnMut(f64) -> f64, Z: FnMut(&mut R, f64) -> f64 {
+    loop {
+        // As an optimisation we re-implement the conversion to a f64.
+        // From the remaining 12 most significant bits we use 8 to construct `i`.
+        // This saves us generating a whole extra random number, while the added
+        // precision of using 64 bits for f64 does not buy us much.
+        let bits = rng.next_u64();
+        let i = bits as usize & 0xff;
+
+        let u = if symmetric {
+            // Convert to a value in the range [2,4) and substract to get [-1,1)
+            // We can't convert to an open range directly, that would require
+            // substracting `3.0 - EPSILON`, which is not representable.
+            // It is possible with an extra step, but an open range does not
+            // seem neccesary for the ziggurat algorithm anyway.
+            (bits >> 12).into_float_with_exponent(1) - 3.0
+        } else {
+            // Convert to a value in the range [1,2) and substract to get (0,1)
+            (bits >> 12).into_float_with_exponent(0)
+            - (1.0 - std::f64::EPSILON / 2.0)
+        };
+        let x = u * x_tab[i];
+
+        let test_x = if symmetric { x.abs() } else {x};
+
+        // algebraically equivalent to |u| < x_tab[i+1]/x_tab[i] (or u < x_tab[i+1]/x_tab[i])
+        if test_x < x_tab[i + 1] {
+            return x;
+        }
+        if i == 0 {
+            return zero_case(rng, u);
+        }
+        // algebraically equivalent to f1 + DRanU()*(f0 - f1) < 1
+        if f_tab[i + 1] + (f_tab[i] - f_tab[i + 1]) * rng.gen::<f64>() < pdf(x) {
+            return x;
+        }
+    }
+}