RandomGenerator.h

/**************************************************************************
 * basf2 (Belle II Analysis Software Framework)                           *
 * Author: The Belle II Collaboration                                     *
 *                                                                        *
 * See git log for contributors and copyright holders.                    *
 * This file is licensed under LGPL-3.0, see LICENSE.md.                  *
 **************************************************************************/
 
#pragma once
 
#include <stdint.h>
#include <TRandom.h>
#include <vector>
 
namespace Belle2 {
  class RandomGenerator: public TRandom {
  public:
    enum EGeneratorMode {
      c_independent,
      c_runDependent,
      c_eventDependent
    };

    explicit RandomGenerator(const std::string& name = "Belle2 Random Generator");

    virtual ~RandomGenerator() {}

    void setSeed(const unsigned char* seed, unsigned int n);

    void setMode(EGeneratorMode mode) { m_mode = mode; }

    EGeneratorMode getMode() const { return m_mode; }

    const std::vector<unsigned char>& getSeed() const { return m_seed; }

    void initialize() { setState(0); }

    void barrier() { setState(m_barrier + 1); }

    void setBarrier(int barrierIndex) { setState(barrierIndex); }

    int getBarrier() const { return m_barrier; }

    uint64_t random64();

    uint32_t random32() { return random64() >> 32; }

    double random01();

    Double_t Rndm() override { return random01(); }
    Double_t Rndm(Int_t) override { return Rndm(); }
    void RndmArray(Int_t n, Float_t* array) override;
    void RndmArray(Int_t n, Double_t* array) override;
    void RndmArray(Int_t n, ULong64_t* array);

    void RndmArray(Int_t n, UInt_t* array);

    void RndmArray(Int_t n, Int_t* array) { RndmArray(n, (UInt_t*) array); }

    void RndmArray(Int_t n, Long64_t* array) { RndmArray(n, (ULong64_t*) array); }

    void RndmArray(Int_t n, unsigned char* array);
  private:
    void SetSeed(UInt_t) {}
    void SetSeed(ULong_t) override {}

    void setState(int barrier);

    uint64_t m_state[16];
    unsigned int m_index;
    int m_barrier;
    std::vector<unsigned char> m_seed;
    EGeneratorMode m_mode;
    ClassDefOverride(RandomGenerator, 2);
  };
 
  inline uint64_t RandomGenerator::random64()
  {
    //Generate random number using magic, taken from [arXiv:1402.6246] and only
    //changed to conform to naming scheme
    uint64_t s0 = m_state[ m_index ];
    uint64_t s1 = m_state[ m_index = (m_index + 1) & 15 ];
    s1 ^= s1 << 31;
    s1 ^= s1 >> 11;
    s0 ^= s0 >> 30;
    return (m_state[ m_index ] = s0 ^ s1) * 1181783497276652981LL;
  }
 
  inline void RandomGenerator::RndmArray(Int_t n, Float_t* array)
  {
    //We could optimize this more thoroughly since one 64bit random int is
    //enough to generate two floats but this would probably be rather academic
    for (int i = 0; i < n; ++i) array[i] = random01();
  }
 
  inline void RandomGenerator::RndmArray(Int_t n, Double_t* array)
  {
    //Fill the array, no optimization whatsoever necessary
    for (int i = 0; i < n; ++i) array[i] = random01();
  }
 
  inline void RandomGenerator::RndmArray(Int_t n, ULong64_t* array)
  {
    //Fill the array, no optimization whatsoever necessary
    for (int i = 0; i < n; ++i) array[i] = random64();
  }
 
  inline void RandomGenerator::RndmArray(Int_t n, UInt_t* array)
  {
    //Fill the most part of the array using 64bit numbers
    RndmArray(n / 2, (ULong64_t*)array);
    //Only for uneven number of elements we need to fill the last one using a
    //32bit number
    if (n % 2) array[n - 1] = random32();
  }
 
  inline double RandomGenerator::random01()
  {
    // There are two possibilities to generate a uniform double between 0 and
    // 1: multiply the integer by a constant or exploit the double
    // representation and use some bit shift magic. We have both implementation
    // here and they seem to produce the exact same output so we stick with the
    // more readable one but leave the bitshift solution just in case
#ifdef RANDOM_IEEE754
    //Generate a double in (0,1) using magic bit hackery with doubles: The
    //memory layout of a IEEE754 double precision floating point variable
    //is [sign(1)|exponent(11)|fraction(52)] with values in parentheses
    //being the number of bits. The actual value is then
    //-1^{sign} * (1.fraction) * 2^{exponent-1023}. Setting sign to 0 the
    //exponent to 1023 will thus return a value between 1 (inclusive) and 2
    //(exclusive). So we shift the 64bit integer to the right by 12 bits
    //(which gives as zeros for sign and exponent) and logical or this with
    //the correct binary representation of the exponent
 
    //To do this we use a union to modify the binary representation using
    //an integer and then return the double value
    union { uint64_t i; double d; } x;
    x.i = random64() >> 12;
    //This is a bit academic but we want (0,1) so if we happen to get
    //exactly zero we try again. Chance is 1 in 2^52
    if (x.i == 0) return random01();
    x.i |= 0x3FF0000000000000ULL;
    return x.d - 1.0;
#else
    //Generate a double (0,1) the traditional way by multiplying it with a
    //constant. As doubles only have a precision of 52 bits we need to
    //remove the 12 leading bits from our random int value
    const uint64_t x = random64() >> 12;
    //This is a bit academic but we want (0,1) so if we happen to get
    //exactly zero we try again. Chance is 1 in 2^52
    if (!x) return random01();
    //return x / 2^{52};
    return x * 2.220446049250313080847263336181640625e-16;
#endif
  }

}