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-2017

0031 // =======================================================================

0032
0033 #ifndef MixMaxRng_h
0034 #define MixMaxRng_h 1
0035
0036 #include <array>
0037 #include <cstdint>
0038 #include "CLHEP/Random/defs.h"
0039 #include "CLHEP/Random/RandomEngine.h"
0040
0041 namespace CLHEP {
0042
0043 /**

0044   * @author  K.Savvidy

0045   * @ingroup random

0046   */
0047
0048 using myID_t = std::uint32_t;
0049 using myuint_t = unsigned long long int;
0050
0051 class MixMaxRng: public HepRandomEngine {
0052
0053   static const int N = 17;
0054
0055 public:
0056
0057   MixMaxRng(std::istream& is);
0058   MixMaxRng();
0059   MixMaxRng(long seed);
0060   ~MixMaxRng();
0061   // Constructors and destructor.

0062
0063   MixMaxRng(const MixMaxRng& rng);
0064   MixMaxRng& operator=(const MixMaxRng& rng);
0065   // Copy constructor and assignment operator.

0066
0067   double flat() { return (S.counter<=(N-1)) ? generate(S.counter):iterate(); }
0068   // Returns a pseudo random number between 0 and 1

0069   // (excluding the zero: in (0,1] )

0070   // smallest number which it will give is approximately 10^-19

0071
0072   void flatArray (const int size, double* vect);
0073   // Fills the array "vect" of specified size with flat random values.

0074
0075   void setSeed(long seed, int dum=0);
0076   // Sets the state of the algorithm according to seed.

0077
0078   void setSeeds(const long * seeds, int seedNum=0);
0079   // Sets the initial state of the engine according to the array of between one and four 32-bit seeds.

0080   // If the size of long is greater on the platform, only the lower 32-bits are used.

0081   // Streams created from seeds differing by at least one bit somewhere are guaranteed absolutely

0082   // to be independent and non-colliding for at least the next 10^100 random numbers

0083
0084   void saveStatus( const char filename[] = "MixMaxRngState.conf" ) const;
0085   // Saves the the current engine state in the file given, by default MixMaxRngState.conf

0086
0087   void restoreStatus( const char filename[] = "MixMaxRngState.conf" );
0088   // Reads a valid engine state from a given file, by default MixMaxRngState.conf

0089   // and restores it.

0090
0091   void showStatus() const;
0092   // Dumps the engine status on the screen.

0093
0094   operator double();
0095   // Returns same as flat()

0096   operator float();
0097   // less precise flat, faster if possible

0098   operator unsigned int();
0099   // 32-bit flat

0100
0101   virtual std::ostream & put (std::ostream & os) const;
0102   virtual std::istream & get (std::istream & is);
0103   static  std::string beginTag ( );
0104   virtual std::istream & getState ( std::istream & is );
0105
0106   std::string name() const { return "MixMaxRng"; }
0107   static std::string engineName();
0108
0109   std::vector<unsigned long> put () const;
0110   bool get (const std::vector<unsigned long> & v);
0111   bool getState (const std::vector<unsigned long> & v);
0112
0113 private:
0114
0115   static constexpr long long int SPECIAL   = ((N==17)? 0 : ((N==240)? 487013230256099140ULL:0) ); // etc...

0116   static constexpr long long int SPECIALMUL= ((N==17)? 36: ((N==240)? 51                   :53) ); // etc...

0117   // Note the potential for confusion...

0118   static constexpr int BITS=61;
0119   static constexpr myuint_t M61=2305843009213693951ULL;
0120   static constexpr double INV_M61=0.43368086899420177360298E-18;
0121   static constexpr unsigned int VECTOR_STATE_SIZE = 2*N+4; // 2N+4 for MIXMAX

0122
0123   #define MIXMAX_MOD_MERSENNE(k) ((((k)) & M61) + (((k)) >> BITS) )
0124
0125   static constexpr int rng_get_N();
0126   static constexpr long long int rng_get_SPECIAL();
0127   static constexpr int rng_get_SPECIALMUL();
0128   void seed_uniquestream( myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );
0129   void seed_spbox(myuint_t);
0130   void print_state() const;
0131   myuint_t precalc();
0132   myuint_t get_next() ;
0133   inline double get_next_float() { return get_next_float_packbits(); }
0134   // Returns a random double with all 52 bits random, in the range (0,1]

0135
0136   MixMaxRng Branch();
0137   void BranchInplace(int id);
0138
0139   MixMaxRng(myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );      // Constructor with four 32-bit seeds

0140   inline void seed64(myuint_t seedval) { seed_uniquestream( 0, 0, (myID_t)(seedval>>32), (myID_t)seedval ); } // seed with one 64-bit seed

0141
0142   double generate(int i);
0143   double iterate();
0144
0145   double get_next_float_packbits();
0146 #if defined __GNUC__
0147 #pragma GCC diagnostic push
0148 #pragma GCC diagnostic ignored "-Wstrict-aliasing"
0149 #pragma GCC diagnostic ignored "-Wuninitialized"
0150 #endif
0151   inline double convert1double(myuint_t u)
0152   {
0153     const double one = 1;
0154     const myuint_t onemask = *(myuint_t*)&one;
0155     myuint_t tmp = (u>>9) | onemask; // bits between 52 and 62 dont affect the result!

0156     double d = *(double*)&tmp;
0157     return d-1.0;
0158   }
0159 #if defined __GNUC__
0160 #pragma GCC diagnostic pop
0161 #endif
0162   myuint_t MOD_MULSPEC(myuint_t k);
0163   myuint_t MULWU(myuint_t k);
0164   void seed_vielbein( unsigned int i); // seeds with the i-th unit vector, i = 0..N-1,  for testing only

0165   myuint_t iterate_raw_vec(myuint_t* Y, myuint_t sumtotOld);
0166   myuint_t apply_bigskip(myuint_t* Vout, myuint_t* Vin, myID_t clusterID, myID_t machineID, myID_t runID, myID_t  streamID );
0167   myuint_t modadd(myuint_t foo, myuint_t bar);
0168 #if defined(__x86_64__)
0169   myuint_t mod128(__uint128_t s);
0170   myuint_t fmodmulM61(myuint_t cum, myuint_t a, myuint_t b);
0171 #else // on all other platforms, including 32-bit linux, PPC and PPC64, ARM and all Windows

0172   myuint_t fmodmulM61(myuint_t cum, myuint_t s, myuint_t a);
0173 #endif
0174
0175 private:
0176
0177   struct rng_state_st
0178   {
0179       std::array<myuint_t, N> V;
0180       myuint_t sumtot;
0181       int counter;
0182   };
0183
0184   typedef struct rng_state_st rng_state_t;     // struct alias

0185   rng_state_t S;
0186 };
0187
0188 }  // namespace CLHEP

0189
0190 #endif