source/include/SiPMRandom.h

0001 /** @class sipm::SiPMRandom SimSiPM/src/components/SiPMRandom.h SiPMRandom.h
0002  *
0003  * @brief Class for random number generation.
0004  *
0005  * Class used for random number generation. The simulation needs very fast
0006  * pseudo-random number generation, Xorhift256+ algorithm is used as it is one
0007  * of the fastest considering modern x86-64 architectures.
0008  *
0009  *  @author Edoardo Proserpio
0010  *  @date 2020
0011  */
0012
0013 #ifndef SIPM_RANDOM_H
0014 #define SIPM_RANDOM_H
0015
0016 #include <cmath>
0017 #include <cstring>
0018 #include <random>
0019 #include <stdint.h>
0020 #include <vector>
0021
0022 #ifdef __AVX2__
0023 #include <immintrin.h>
0024 #endif
0025
0026 #include "SiPMMath.h"
0027
0028 namespace sipm {
0029
0030 namespace SiPMRng {
0031
0032 /// @brief Implementation of xoshiro256++ 1.0 PRNG algorithm
0033 /** Written in 2019 by David Blackman and Sebastiano Vigna (vigna@acm.org)
0034  *
0035  * To the extent possible under law, the author has dedicated all copyright
0036  * and related and neighboring rights to this software to the public domain
0037  * worldwide. This software is distributed without any warranty.
0038  *
0039  * See <http://creativecommons.org/publicdomain/zero/1.0/>.
0040  * This is xoshiro256++ 1.0, one of our all-purpose, rock-solid generators.
0041  * It has excellent (sub-ns) speed, a state (256 bits) that is large
0042  * enough for any parallel application, and it passes all tests we are
0043  * aware of.
0044  *
0045  * For generating just floating-point numbers, xoshiro256+ is even faster.
0046  * The state must be seeded so that it is not everywhere zero. If you have
0047  * a 64-bit seed, we suggest to seed a splitmix64 generator and use its
0048  * output to fill s. */
0049 class Xorshift256plus {
0050 public:
0051   /// @brief Default contructor for Xorshift256plus
0052   Xorshift256plus() { seed(); }
0053   /// @brief Contructor for Xorshift256plus given a seed value
0054   Xorshift256plus(uint64_t aseed) { seed(aseed); }
0055   /// @brief Returns a pseud-random 64-bits intger
0056   inline uint64_t operator()() noexcept {
0057     const uint64_t result = s[0] + s[3];
0058
0059     const uint64_t t = s[1] << 17;
0060
0061     s[2] ^= s[0];
0062     s[3] ^= s[1];
0063     s[1] ^= s[2];
0064     s[0] ^= s[3];
0065
0066     s[2] ^= t;
0067
0068     s[3] = (s[3] << 45U) | (s[3] >> (64U - 45U));
0069     return result;
0070   }
0071   __attribute__((hot))
0072
0073   /// @brief Jump function for the alghoritm.
0074   /** Usefull in case the same generator is used in multiple instancies. The
0075    * jump function will make sure that pseudo-random values generated from the
0076    * different instancies are uncorrelated.
0077    */
0078   void
0079   jump();
0080   /// @brief Sets a random seed generated with rand()
0081   void seed();
0082   /// @brief Sets a new seed
0083   void seed(uint64_t);
0084
0085 private:
0086   alignas(64) uint64_t s[4];
0087 };
0088
0089 } // namespace SiPMRng
0090
0091 class SiPMRandom {
0092 public:
0093   SiPMRandom() = default;
0094   SiPMRandom(uint64_t aseed) noexcept { m_rng.seed(aseed); }
0095
0096   /** @brief Sets a seed for the rng.
0097    * @param aSeed Seed used to initialize the rng algorithm
0098    */
0099   void seed(const uint64_t aSeed) { m_rng.seed(aSeed); }
0100   /** @brief Sets a seed for the rng obtained from rand().*/
0101   void seed() { m_rng.seed(); }
0102   /** @brief This is the jump function for the generator. It is equivalent
0103    * to 2^128 calls to next(); it can be used to generate 2^128
0104    * non-overlapping subsequences for parallel computations.*/
0105   void jump() { m_rng.jump(); }
0106
0107   inline uint64_t operator()() { return m_rng(); }
0108
0109   // Uniform random in [0-1]
0110   inline double Rand() __attribute__((hot));
0111
0112   // Uniform integer in range [0-max)
0113   uint32_t randInteger(const uint32_t) __attribute__((hot));
0114   // Random gaussian given mean and sigma
0115   double randGaussian(const double, const double) __attribute__((hot));
0116   // Random exponential given mean
0117   double randExponential(const double);
0118   // Random poisson given mean
0119   uint32_t randPoisson(const double mu);
0120
0121   /** @brief Vector of random uniforms in [0-1] */
0122   std::vector<double> Rand(const uint32_t) __attribute__((hot));
0123   /** @brief Vector of random gaussian given mean an sigma */
0124   std::vector<double> randGaussian(const double, const double, const uint32_t) __attribute__((hot));
0125   /** @brief Vector of random integers in range [0-max) */
0126   std::vector<uint32_t> randInteger(const uint32_t max, const uint32_t n) __attribute__((hot));
0127   // Vector of random exponential given mean
0128   std::vector<double> randExponential(const double, const uint32_t);
0129
0130 private:
0131   SiPMRng::Xorshift256plus m_rng;
0132 };
0133
0134 /** Returns a uniform random in range (0,1)
0135  * Using getting highest 53 bits from a unit64 for the mantissa,
0136  * setting the exponent to get values in range (1-2), subtract 1 and type punning
0137  * to double.
0138  */
0139 inline double SiPMRandom::Rand() {
0140   double x;
0141   static constexpr uint64_t mask = 0x3ff0000000000000;
0142   const uint64_t u = (m_rng() >> 12) | mask;
0143   std::memcpy(&x, &u, 8);
0144   return x - 1;
0145 }
0146 } // namespace sipm
0147 #endif /* SIPM_RANDOM_H */