CLHEP/Random/MixMaxRng.h

0001 //
0002 // -*- C++ -*-
0003 //
0004 // -----------------------------------------------------------------------
0005 //                          HEP Random
0006 //                       --- MixMaxRng ---
0007 //                       class header file
0008 // -----------------------------------------------------------------------
0009 //
0010 // This file interfaces the MixMax PseudoRandom Number Generator
0011 // proposed by:
0012 //
0013 // G.K.Savvidy and N.G.Ter-Arutyunian,
0014 //   On the Monte Carlo simulation of physical systems,
0015 //   J.Comput.Phys. 97, 566 (1991);
0016 //   Preprint EPI-865-16-86, Yerevan, Jan. 1986
0017 //   http://dx.doi.org/10.1016/0021-9991(91)90015-D
0018 //
0019 // K.Savvidy
0020 //   "The MIXMAX random number generator"
0021 //   Comp. Phys. Commun. (2015)
0022 //   http://dx.doi.org/10.1016/j.cpc.2015.06.003
0023 //
0024 // K.Savvidy and G.Savvidy
0025 //   "Spectrum and Entropy of C-systems. MIXMAX random number generator"
0026 //   Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33-38
0027 //   http://dx.doi.org/10.1016/j.chaos.2016.05.003
0028 //
0029 // =======================================================================
0030 // Implementation by Konstantin Savvidy - Copyright 2004-2023
0031 // July 2023 - Updated class structure upon suggestions from Marco Barbone
0032 // September 2023 - fix (re-)initialization from Gabriele Cosmo
0033 // =======================================================================
0034
0035 #ifndef MixMaxRng_h
0036 #define MixMaxRng_h 1
0037
0038 #include <array>
0039 #include <cstdint>
0040 #include "CLHEP/Random/defs.h"
0041 #include "CLHEP/Random/RandomEngine.h"
0042
0043 namespace CLHEP {
0044
0045 /**
0046   * @author  K.Savvidy
0047   * @ingroup random
0048   */
0049
0050 using myID_t = std::uint32_t;
0051 using myuint_t = std::uint64_t;
0052
0053 class
0054 #if (__cplusplus >= 201703L && __cpp_aligned_new >= 201606L)
0055   alignas(128)
0056 #endif
0057 MixMaxRng : public HepRandomEngine
0058 {
0059
0060   static const int N = 17;
0061
0062 public:
0063
0064   MixMaxRng(std::istream& is);
0065   MixMaxRng();
0066   MixMaxRng(long seed);
0067   ~MixMaxRng();
0068   // Constructors and destructor.
0069
0070   MixMaxRng(const MixMaxRng& rng);
0071   MixMaxRng& operator=(const MixMaxRng& rng);
0072   // Copy constructor and assignment operator.
0073
0074   inline double flat()
0075   {
0076     if (counter >= N) iterate();
0077     return INV_M61*static_cast<double>(V[counter++]);
0078   }
0079   // Returns a pseudo random number between 0 and 1
0080   // excluding the zero: in (0,1]
0081   // smallest number which it will give is approximately 10^-19
0082
0083   void flatArray (const int size, double* vect);
0084   // Fills the array "vect" of specified size with flat random values.
0085
0086   inline void setSeed(long longSeed, int = 0 /* extraSeed */)
0087   {
0088     seed_spbox(theSeed = longSeed);
0089   }
0090   // Sets the state of the algorithm according to seed.
0091
0092   void setSeeds(const long * seeds, int seedNum=0);
0093   // Sets the initial state of the engine according to the array of between one and four 32-bit seeds.
0094   // If the size of long is greater on the platform, only the lower 32-bits are used.
0095   // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely
0096   // to be independent and non-colliding for at least the next 10^100 random numbers
0097
0098   void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
0099   // Saves the the current engine state in the file given, by default MixMaxRngState.conf
0100
0101   void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
0102   // Reads a valid engine state from a given file, by default MixMaxRngState.conf
0103   // and restores it.
0104
0105   void showStatus() const;
0106   // Dumps the engine status on the screen.
0107
0108   inline operator double() { return flat(); }
0109   // Returns same as flat()
0110
0111   inline operator float() { return float( flat() ); }
0112   // less precise flat, faster if possible
0113
0114   inline operator unsigned int() { return static_cast<unsigned int>(get_next()); }
0115   // 32-bit flat. clhep_get_next() returns a 64-bit integer, of which
0116   // the lower 61 bits are random and upper 3 bits are zero
0117
0118   virtual std::ostream & put (std::ostream & os) const;
0119   virtual std::istream & get (std::istream & is);
0120   static  std::string beginTag ( );
0121   virtual std::istream & getState ( std::istream & is );
0122
0123   std::string name() const { return "MixMaxRng"; }
0124   static std::string engineName();
0125
0126   std::vector<unsigned long> put () const;
0127   bool get (const std::vector<unsigned long> & vec);
0128   bool getState (const std::vector<unsigned long> & vec);
0129
0130 private:
0131
0132   static constexpr long long int SPECIAL   = 0;
0133   static constexpr long long int SPECIALMUL= 36;
0134   static constexpr int BITS=61;
0135   static constexpr myuint_t M61=2305843009213693951ULL;
0136   static constexpr double INV_M61=0.43368086899420177360298E-18;
0137   static constexpr unsigned int VECTOR_STATE_SIZE = 2*N+4; // 2N+4 for MIXMAX
0138
0139   inline myuint_t MIXMAX_MOD_MERSENNE(myuint_t k)
0140   {
0141     return ((((k)) & M61) + (((k)) >> BITS) );
0142   }
0143
0144   static constexpr int rng_get_N();
0145   void seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );
0146   void seed_spbox(myuint_t seed);
0147   void print_state() const;
0148   myuint_t precalc();
0149   myuint_t get_next();
0150
0151   MixMaxRng Branch();
0152   void BranchInplace(int id);
0153
0154   MixMaxRng(myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );      // Constructor with four 32-bit seeds
0155   inline void seed64(myuint_t seedval)   // seed with one 64-bit seed
0156   {
0157     seed_uniquestream( 0, 0, (myID_t)(seedval>>32), (myID_t)seedval );
0158   }
0159
0160   inline void iterate()
0161   {
0162     myuint_t  tempP, tempV;
0163     V[0] = ( tempV = sumtot );
0164     myuint_t insumtot = V[0], ovflow = 0; // will keep a running sum of all new elements
0165     tempP = 0;                            // will keep a partial sum of all old elements
0166     myuint_t tempPO;
0167     tempPO = MULWU(tempP); tempP = modadd(tempP, V[1] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[1]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0168     tempPO = MULWU(tempP); tempP = modadd(tempP, V[2] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[2]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0169     tempPO = MULWU(tempP); tempP = modadd(tempP, V[3] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[3]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0170     tempPO = MULWU(tempP); tempP = modadd(tempP, V[4] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[4]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0171     tempPO = MULWU(tempP); tempP = modadd(tempP, V[5] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[5]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0172     tempPO = MULWU(tempP); tempP = modadd(tempP, V[6] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[6]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0173     tempPO = MULWU(tempP); tempP = modadd(tempP, V[7] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[7]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0174     tempPO = MULWU(tempP); tempP = modadd(tempP, V[8] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[8]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0175     tempPO = MULWU(tempP); tempP = modadd(tempP, V[9] ); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[9]  = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0176     tempPO = MULWU(tempP); tempP = modadd(tempP, V[10]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[10] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0177     tempPO = MULWU(tempP); tempP = modadd(tempP, V[11]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[11] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0178     tempPO = MULWU(tempP); tempP = modadd(tempP, V[12]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[12] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0179     tempPO = MULWU(tempP); tempP = modadd(tempP, V[13]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[13] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0180     tempPO = MULWU(tempP); tempP = modadd(tempP, V[14]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[14] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0181     tempPO = MULWU(tempP); tempP = modadd(tempP, V[15]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[15] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0182     tempPO = MULWU(tempP); tempP = modadd(tempP, V[16]); tempV = MIXMAX_MOD_MERSENNE(tempV+tempP+tempPO); V[16] = tempV; insumtot += tempV; if (insumtot < tempV) {++ovflow;}
0183     sumtot = MIXMAX_MOD_MERSENNE(MIXMAX_MOD_MERSENNE(insumtot) + (ovflow <<3 ));
0184
0185     counter=1;
0186   }
0187
0188   void state_init();
0189   inline myuint_t MULWU (myuint_t k)
0190   {
0191     return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) >> (BITS-SPECIALMUL))  );
0192   }
0193   myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld);
0194   myuint_t apply_bigskip(myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );
0195   inline myuint_t modadd(myuint_t xfoo, myuint_t xbar)
0196   {
0197     return MIXMAX_MOD_MERSENNE(xfoo+xbar);
0198   }
0199
0200 #if defined(__x86_64__)
0201   myuint_t mod128(__uint128_t s);
0202   myuint_t fmodmulM61(myuint_t cum, myuint_t a, myuint_t b);
0203 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows
0204   myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a);
0205 #endif
0206
0207   // Engine state
0208
0209   myuint_t V[N] = {0};
0210   myuint_t sumtot = 0;
0211   int counter = N;
0212 };
0213
0214 }  // namespace CLHEP
0215
0216 #endif