r2rs-nmath 0.1.1

use std::fmt::Debug;

use chrono::{DateTime, Utc};

/// # Random Number Generation
///
/// ## Description:
///
/// Implementors of this trait allow for the random generation of numbers.  By default, seeds are set
/// by the current time to allow for random and independent numbers, but a seed can be set by the
/// user.
///
/// ## Usage:
///
/// ```rust
/// # use r2rs_nmath::{rng::MarsagliaMulticarry, traits::RNG};
/// let mut rng = MarsagliaMulticarry::new();
/// rng.set_seed(1);
/// let state = rng.get_state();
/// let rand = rng.unif_rand();
/// ```
///
/// ## Details:
///
/// The currently available RNG kinds are given below.
///
/// * *Wichmann-Hill* The seed, ‘.Random.seed\[-1\] == r\[1:3\]’ is an
/// integer vector of length 3, where each ‘r\[i\]’ is in ‘1:(p\[i\]
/// - 1)’, where ‘p’ is the length 3 vector of primes, ‘p =
/// (30269, 30307, 30323)’.  The Wichmann-Hill generator has a
/// cycle length of 6.9536e12 (= ‘prod(p-1)/4’, see _Applied
/// Statistics_ (1984) *33*, 123 which corrects the original
/// article).  It exhibits 12 clear failures in the TestU01 Crush
/// suite and 22 in the BigCrush suite (_L'Ecuyer_, 2007).
///
/// * *Marsaglia-Multicarry*: A _multiply-with-carry_ RNG is used, as
/// recommended by George Marsaglia in his post to the mailing
/// list ‘sci.stat.math’.  It has a period of more than 2^60.
///
/// It exhibits 40 clear failures in L'Ecuyer's TestU01 Crush
/// suite.  Combined with Ahrens-Dieter or Kinderman-Ramage it
/// exhibits deviations from normality even for univariate
/// distribution generation.  See PR#18168 for a discussion.
///
/// The seed is two integers (all values allowed).
///
/// * *Super-Duper*: Marsaglia's famous Super-Duper from the 70's.
/// This is the original version which does _not_ pass the MTUPLE
/// test of the Diehard battery.  It has a period of about
/// 4.6*10^18 for most initial seeds.  The seed is two integers
/// (all values allowed for the first seed: the second must be
/// odd).
///
/// We use the implementation by Reeds _et al_ (1982-84).
///
/// The two seeds are the Tausworthe and congruence long
/// integers, respectively.  A one-to-one mapping to S's
/// ‘.Random.seed\[1:12\]’ is possible but we will not publish one,
/// not least as this generator is *not* exactly the same as that
/// in recent versions of S-PLUS.
///
/// It exhibits 25 clear failures in the TestU01 Crush suite
/// (_L'Ecuyer_, 2007).
///
/// * *Mersenne-Twister*: From Matsumoto and Nishimura (1998); code
/// updated in 2002.  A twisted GFSR with period 2^19937 - 1 and
/// equidistribution in 623 consecutive dimensions (over the
/// whole period).  The ‘seed’ is a 624-dimensional set of 32-bit
/// integers plus a current position in that set.
///
/// R uses its own initialization method due to B. D. Ripley and
/// is not affected by the initialization issue in the 1998 code
/// of Matsumoto and Nishimura addressed in a 2002 update.
///
/// It exhibits 2 clear failures in each of the TestU01 Crush and
/// the BigCrush suite (_L'Ecuyer_, 2007).
///
/// * *L'Ecuyer-CMRG*: A ‘combined multiple-recursive generator’ from
/// L'Ecuyer (1999), each element of which is a feedback
/// multiplicative generator with three integer elements: thus
/// the seed is a (signed) integer vector of length 6. The period
/// is around 2^191.
///
/// The 6 elements of the seed are internally regarded as 32-bit
/// unsigned integers.  Neither the first three nor the last
/// three should be all zero, and they are limited to less than
/// ‘4294967087’ and ‘4294944443’ respectively.
///
/// This is not particularly interesting of itself, but provides
/// the basis for the multiple streams used in package
/// ‘parallel’.
///
/// It exhibits 6 clear failures in each of the TestU01 Crush and
/// the BigCrush suite (_L'Ecuyer_, 2007).
///
/// ## Note:
///
/// Initially, there is no seed; a new one is created from the current
/// time.  Hence, different sessions will give different simulation results, by default.
/// However, the seed might be restored from a previous session if a
/// previously saved workspace is restored.
///
/// Do not rely on randomness of low-order bits from RNGs.  Most of
/// the supplied uniform generators return 32-bit integer values that
/// are converted to doubles, so they take at most 2^32 distinct
/// values and long runs will return duplicated values (Wichmann-Hill
/// is the exception, and all give at least 30 varying bits.)
///
/// ## Author(s):
///
/// of RNGkind: Martin Maechler. Current implementation, B. D. Ripley
/// with modifications by Duncan Murdoch.
///
/// ## References:
///
/// Ahrens, J. H. and Dieter, U. (1973).  Extensions of Forsythe's
/// method for random sampling from the normal distribution.
/// _Mathematics of Computation_, *27*, 927-937.
///
/// Becker, R. A., Chambers, J. M. and Wilks, A. R. (1988).  _The New
/// S Language_.  Wadsworth & Brooks/Cole.  (‘set.seed’, storing in
/// ‘.Random.seed’.)
///
/// Box, G. E. P. and Muller, M. E. (1958).  A note on the generation
/// of normal random deviates.  _Annals of Mathematical Statistics_,
/// *29*, 610-611.  doi:10.1214/aoms/1177706645
/// <https://doi.org/10.1214/aoms/1177706645>.
///
/// De Matteis, A. and Pagnutti, S. (1993).  Long-range Correlation
/// Analysis of the Wichmann-Hill Random Number Generator.
/// _Statistics and Computing_, *3*, 67-70.  doi:10.1007/BF00153065
/// <https://doi.org/10.1007/BF00153065>.
///
/// Kinderman, A. J. and Ramage, J. G. (1976).  Computer generation of
/// normal random variables.  _Journal of the American Statistical
/// Association_, *71*, 893-896.  doi:10.2307/2286857
/// <https://doi.org/10.2307/2286857>.
///
/// Knuth, D. E. (1997).  _The Art of Computer Programming_.  Volume
/// 2, third edition.
/// Source code at
/// <https://www-cs-faculty.stanford.edu/~knuth/taocp.html>.
///
/// Knuth, D. E. (2002).  _The Art of Computer Programming_.  Volume
/// 2, third edition, ninth printing.
///
/// L'Ecuyer, P. (1999).  Good parameters and implementations for
/// combined multiple recursive random number generators.  _Operations
/// Research_, *47*, 159-164.  doi:10.1287/opre.47.1.159
/// <https://doi.org/10.1287/opre.47.1.159>.
///
/// L'Ecuyer, P. and Simard, R. (2007).  TestU01: A C Library for
/// Empirical Testing of Random Number Generators _ACM Transactions on
/// Mathematical Software_, *33*, Article 22.
/// doi:10.1145/1268776.1268777
/// <https://doi.org/10.1145/1268776.1268777>.
///
/// The TestU01 C library is available from
/// <http://simul.iro.umontreal.ca/testu01/tu01.html> or also
/// <https://github.com/umontreal-simul/TestU01-2009>.
///
/// Marsaglia, G. (1997).  _A random number generator for C_.
/// Discussion paper, posting on Usenet newsgroup ‘sci.stat.math’ on
/// September 29, 1997.
///
/// Marsaglia, G. and Zaman, A. (1994).  Some portable
/// very-long-period random number generators.  _Computers in
/// Physics_, *8*, 117-121.  doi:10.1063/1.168514
/// <https://doi.org/10.1063/1.168514>.
///
/// Matsumoto, M. and Nishimura, T. (1998).  Mersenne Twister: A
/// 623-dimensionally equidistributed uniform pseudo-random number
/// generator, _ACM Transactions on Modeling and Computer Simulation_,
/// *8*, 3-30.
/// Source code formerly at
/// ‘http://www.math.keio.ac.jp/~matumoto/emt.html’.
/// Now see
/// <http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/VERSIONS/C-LANG/c-lang.html>.
///
/// Reeds, J., Hubert, S. and Abrahams, M. (1982-4).  C implementation
/// of SuperDuper, University of California at Berkeley.  (Personal
/// communication from Jim Reeds to Ross Ihaka.)
///
/// Wichmann, B. A.  and Hill, I. D. (1982).  Algorithm AS 183: An
/// Efficient and Portable Pseudo-random Number Generator.  _Applied
/// Statistics_, *31*, 188-190; Remarks: *34*, 198 and *35*, 89.
/// doi:10.2307/2347988 <https://doi.org/10.2307/2347988>.
///
/// ## Examples:
///
/// ```rust
/// # use r2rs_nmath::{rng::MarsagliaMulticarry, traits::RNG};
/// // Seed the current RNG, i.e., set the RNG status
/// let mut rng = MarsagliaMulticarry::new();
/// rng.set_seed(42);
/// let u1 = (0..30).map(|_| rng.unif_rand()).collect::<Vec<_>>();
///
/// rng.set_seed(42);
/// let u2 = (0..30).map(|_| rng.unif_rand()).collect::<Vec<_>>();
///
/// // These are the same because of identical RNG status
/// assert_eq!(u1, u2);
/// ```
///
/// ```rust
/// # use r2rs_nmath::{
/// #    rng::{MarsagliaMulticarry, SuperDuper},
/// #    traits::RNG,
/// # };
/// // Seed the current RNG, i.e., set the RNG status
/// let mut rng1 = MarsagliaMulticarry::new();
/// rng1.set_seed(42);
/// let u1 = (0..30).map(|_| rng1.unif_rand()).collect::<Vec<_>>();
///
/// let mut rng2 = SuperDuper::new();
/// rng2.set_seed(42);
/// let u2 = (0..30).map(|_| rng2.unif_rand()).collect::<Vec<_>>();
///
/// // These are different because of different RNG generators
/// assert_ne!(u1, u2);
/// ```
pub trait RNG: Clone + Debug + Sized {
    /// Sets the random generation seed.  By default, seeds are set from the time of object
    /// creation.
    fn set_seed(&mut self, seed: u32);

    /// Get the current random generation state for the generator.  This is the value that is
    /// internally used to create random numbers (created from the input seed).
    fn get_state(&self) -> Vec<i32>;

    /// Generate a random number from the generator.  This will then generate the next internal
    /// state to allow a new random number nxet time.
    fn unif_rand(&mut self) -> f64;

    /// Set the internal seed to the current time, effectively scrambling the generator.
    fn randomize(&mut self) {
        let time: DateTime<Utc> = Utc::now();
        let time_seed = time.timestamp_millis();
        self.set_seed(time_seed as u32);
    }

    /// Generate a random non-negative integer < 2 ^ bits in 16 bit chunks
    fn rbits(&mut self, bits: u32) -> f64 {
        let mut v: i64 = 0;
        let mut n = 0;
        while n <= bits {
            let v1 = (self.unif_rand() * 65536.0).floor();
            v = 65536 * v + v1 as i64;
            n += 16
        }
        // mask out the bits in the result that are not needed
        return (v & ((1 as i64) << bits as i64) - 1 as i64) as f64;
    }

    /// Rejection sampling from integers below the next larger power of two
    fn unif_index(&mut self, dn: f64) -> f64 {
        if dn <= 0.0 {
            return 0.0;
        }
        let bits = dn.log2().ceil();
        let mut dv = None;
        while dv.is_none() {
            let dv_temp = self.rbits(bits as u32);
            if dn > dv_temp {
                dv = Some(dv_temp);
            }
        }
        dv.unwrap()
    }
}