From fd091b04316db9dc5fafadbd6bdbe60b127408a9 Mon Sep 17 00:00:00 2001 From: Daniel Mueller Date: Thu, 2 Jan 2020 08:32:06 -0800 Subject: 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 --- rand/rand_distr/src/utils.rs | 234 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 234 insertions(+) create mode 100644 rand/rand_distr/src/utils.rs (limited to 'rand/rand_distr/src/utils.rs') 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 or the MIT license +// , 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 + + ops::Add + + ops::Sub + + ops::Mul + + ops::Div + + 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; + + /// 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 { + 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 { + 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( + 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::() < pdf(x) { + return x; + } + } +} -- cgit v1.2.1