ceph/src/boost/boost/random/mixmax.hpp

   1 /* boost random/mixmax.hpp header file
   2  *
   3  * Copyright Kostas Savvidis 2008-2019
   4  *
   5  * Distributed under the Boost Software License, Version 1.0. (See
   6  * accompanying file LICENSE_1_0.txt or copy at
   7  * http://www.boost.org/LICENSE_1_0.txt)
   8  *
   9  * See http://www.boost.org for most recent version including documentation.
  10  *
  11  * $Id$
  12  *
  13  * Revision history
  14  *  2019-04-23 created
  15  */
  16
  17 #ifndef BOOST_RANDOM_MIXMAX_HPP
  18 #define BOOST_RANDOM_MIXMAX_HPP
  19
  20 #include <sstream>
  21 #include <boost/cstdint.hpp>
  22 #include <boost/array.hpp>
  23
  24 #include <boost/random/detail/seed.hpp>
  25 #include <boost/random/detail/seed_impl.hpp>
  26
  27 namespace boost {
  28 namespace random {
  29
  30 /**
  31  * Instantiations of class template mixmax_engine model,
  32  * \pseudo_random_number_generator .
  33  *  It uses the  MIXMAX generator algorithms from:
  34  *
  35  *  @blockquote
  36  *  G.K.Savvidy and N.G.Ter-Arutyunian,
  37  *  On the Monte Carlo simulation of physical systems,
  38  *  J.Comput.Phys. 97, 566 (1991);
  39  *  Preprint EPI-865-16-86, Yerevan, Jan. 1986
  40  *  http://dx.doi.org/10.1016/0021-9991(91)90015-D
  41  *
  42  *  K.Savvidy
  43  *  The MIXMAX random number generator
  44  *  Comp. Phys. Commun. 196 (2015), pp 161–165
  45  *  http://dx.doi.org/10.1016/j.cpc.2015.06.003
  46  *
  47  *  K.Savvidy and G.Savvidy
  48  *  Spectrum and Entropy of C-systems. MIXMAX random number generator
  49  *  Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33–38
  50  *  http://dx.doi.org/10.1016/j.chaos.2016.05.003
  51  *  @endblockquote
  52  *
  53  * The generator crucially depends on the choice of the
  54  * parameters. The valid sets of parameters are from the published papers above.
  55  *
  56  */
  57
  58 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> // MIXMAX TEMPLATE PARAMETERS
  59 class mixmax_engine{
  60 public:
  61     // Interfaces required by C++11 std::random and boost::random
  62     typedef boost::uint64_t result_type ;
  63     BOOST_STATIC_CONSTANT(boost::uint64_t,mixmax_min=0);
  64     BOOST_STATIC_CONSTANT(boost::uint64_t,mixmax_max=((1ULL<<61)-1));
  65     BOOST_STATIC_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() {return mixmax_min;}
  66     BOOST_STATIC_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() {return mixmax_max;}
  67     static const bool has_fixed_range = false;
  68     BOOST_STATIC_CONSTANT(int,N=Ndim);     ///< The main internal parameter, size of the defining MIXMAX matrix
  69     // CONSTRUCTORS:
  70     explicit mixmax_engine();                       ///< Constructor, unit vector as initial state, acted on by A^2^512
  71     explicit mixmax_engine(boost::uint64_t);          ///< Constructor, one 64-bit seed
  72     explicit mixmax_engine(uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID );  ///< Constructor, four 32-bit seeds for 128-bit seeding flexibility
  73     void seed(boost::uint64_t seedval=default_seed){seed_uniquestream( &S, 0, 0, (uint32_t)(seedval>>32), (uint32_t)seedval );} ///< seed with one 64-bit seed
  74
  75 private: // DATATYPES
  76     struct rng_state_st{
  77         boost::array<boost::uint64_t, Ndim> V;
  78         boost::uint64_t sumtot;
  79         int counter;
  80     };
  81
  82     typedef struct rng_state_st rng_state_t;     // struct alias
  83     rng_state_t S;
  84
  85 public: // SEEDING FUNCTIONS
  86     template<class It> mixmax_engine(It& first, It last) { seed(first,last); }
  87     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mixmax_engine,  SeedSeq, seq){ seed(seq); }
  88
  89     /** Sets the state of the generator using values from an iterator range. */
  90     template<class It>
  91     void seed(It& first, It last){
  92         uint32_t v[4];
  93         detail::fill_array_int<32>(first, last, v);
  94         seed_uniquestream( &S, v[0], v[1], v[2], v[3]);
  95     }
  96     /** Sets the state of the generator using values from a seed_seq. */
  97     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mixmax_engine, SeeqSeq, seq){
  98         uint32_t v[4];
  99         detail::seed_array_int<32>(seq, v);
 100         seed_uniquestream( &S, v[0], v[1], v[2], v[3]);
 101     }
 102
 103     /** return one uint64 between min=0 and max=2^61-1 */
 104     boost::uint64_t operator()(){
 105         if (S.counter<=(Ndim-1) ){
 106             return S.V[S.counter++];
 107         }else{
 108             S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot);
 109             S.counter=2;
 110             return S.V[1];
 111         }
 112     }
 113
 114     /** Fills a range with random values */
 115     template<class Iter>
 116     void generate(Iter first, Iter last) { detail::generate_from_int(*this, first, last); }
 117
 118     void discard(boost::uint64_t nsteps) { for(boost::uint64_t j = 0; j < nsteps; ++j)  (*this)(); } ///< discard n steps, required in boost::random
 119
 120     /** save the state of the RNG to a stream */
 121     template<class CharT, class Traits>
 122     friend std::basic_ostream<CharT,Traits>&
 123     operator<< (std::basic_ostream<CharT,Traits>& ost, const mixmax_engine& me){
 124         ost << Ndim << " " << me.S.counter << " " << me.S.sumtot << " ";
 125         for (int j=0; (j< (Ndim) ); j++) {
 126         ost <<  (boost::uint64_t)me.S.V[j] << " ";
 127         }
 128         ost << "\n";
 129         ost.flush();
 130         return ost;
 131         }
 132
 133     /** read the state of the RNG from a stream */
 134     template<class CharT, class Traits>
 135     friend std::basic_istream<CharT,Traits>&
 136     operator>> (std::basic_istream<CharT,Traits> &in, mixmax_engine& me){
 137         // will set std::ios::failbit if the input format is not right
 138         boost::array<boost::uint64_t, Ndim> vec;
 139         boost::uint64_t sum=0, savedsum=0, counter=0;
 140         in >> counter >> std::ws;
 141         BOOST_ASSERT(counter==Ndim);
 142         in >> counter >> std::ws;
 143         in >> savedsum >> std::ws;
 144         for(int j=0;j<Ndim;j++) {
 145         in >> std::ws >> vec[j] ;
 146         sum=me.MOD_MERSENNE(sum+vec[j]);
 147     }
 148     if (sum == savedsum && counter>0 && counter<Ndim){
 149         me.S.V=vec; me.S.counter = counter; me.S.sumtot=savedsum;
 150     }else{
 151         in.setstate(std::ios::failbit);
 152     }
 153     return in;
 154     }
 155
 156 friend bool operator==(const mixmax_engine & x,
 157                        const mixmax_engine & y){return x.S.counter==y.S.counter && x.S.sumtot==y.S.sumtot && x.S.V==y.S.V ;}
 158 friend bool operator!=(const mixmax_engine & x,
 159                        const mixmax_engine & y){return !operator==(x,y);}
 160
 161
 162 private:
 163 BOOST_STATIC_CONSTANT(int, BITS=61);
 164 BOOST_STATIC_CONSTANT(boost::uint64_t, M61=2305843009213693951ULL);
 165 BOOST_STATIC_CONSTANT(boost::uint64_t, default_seed=1);
 166 inline boost::uint64_t MOD_MERSENNE(boost::uint64_t k) {return ((((k)) & M61) + (((k)) >> BITS) );}
 167 inline boost::uint64_t MULWU(boost::uint64_t k);
 168 inline void seed_vielbein(rng_state_t* X, unsigned int i); // seeds with the i-th unit vector, i = 0..Ndim-1,  for testing only
 169 inline void seed_uniquestream( rng_state_t* Xin, uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID );
 170 inline boost::uint64_t iterate_raw_vec(boost::uint64_t* Y, boost::uint64_t sumtotOld);
 171 inline boost::uint64_t apply_bigskip(boost::uint64_t* Vout, boost::uint64_t* Vin, uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID );
 172 inline boost::uint64_t modadd(boost::uint64_t foo, boost::uint64_t bar);
 173 inline boost::uint64_t fmodmulM61(boost::uint64_t cum, boost::uint64_t s, boost::uint64_t a);
 174 };
 175
 176 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::mixmax_engine()
 177 ///< constructor, with no params, seeds with seed=0,  random numbers are as good as from any other seed
 178 {
 179     seed_uniquestream( &S, 0,  0, 0, default_seed);
 180 }
 181
 182 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::mixmax_engine(boost::uint64_t seedval){
 183     ///< constructor, one uint64_t seed, random numbers are statistically independent from any two distinct seeds, e.g. consecutive seeds are ok
 184     seed_uniquestream( &S, 0,  0,  (uint32_t)(seedval>>32), (uint32_t)seedval );
 185 }
 186
 187 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::mixmax_engine(uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID){
 188     // constructor, four 32-bit seeds for 128-bit seeding flexibility
 189     seed_uniquestream( &S, clusterID,  machineID,  runID,  streamID );
 190 }
 191
 192 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> uint64_t mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::MULWU (uint64_t k){ return (( (k)<<(SPECIALMUL) & M61) ^ ( (k) >> (BITS-SPECIALMUL))  )  ;}
 193
 194 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> boost::uint64_t mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::iterate_raw_vec(boost::uint64_t* Y, boost::uint64_t sumtotOld){
 195     // operates with a raw vector, uses known sum of elements of Y
 196     boost::uint64_t  tempP=0, tempV=sumtotOld;
 197     Y[0] = tempV;
 198     boost::uint64_t sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
 199     for (int i=1; i<Ndim; i++){
 200         boost::uint64_t tempPO = MULWU(tempP);
 201         tempV = (tempV+tempPO);
 202         tempP = modadd(tempP, Y[i]);
 203         tempV = modadd(tempV, tempP); // new Y[i] = old Y[i] + old partial * m
 204         Y[i] = tempV;
 205         sumtot += tempV; if (sumtot < tempV) {ovflow++;}
 206     }
 207     return MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
 208 }
 209
 210 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> void mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::seed_vielbein(rng_state_t* X, unsigned int index){
 211     for (int i=0; i < Ndim; i++){
 212         X->V[i] = 0;
 213     }
 214     if (index<Ndim) { X->V[index] = 1; }else{ X->V[0]=1; }
 215     X->counter = Ndim;  // set the counter to Ndim if iteration should happen right away
 216     X->sumtot = 1;
 217 }
 218
 219
 220 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> void mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::seed_uniquestream( rng_state_t* Xin, uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID ){
 221     seed_vielbein(Xin,0);
 222     Xin->sumtot = apply_bigskip(Xin->V.data(), Xin->V.data(),  clusterID,  machineID,  runID,   streamID );
 223     Xin->counter = 1;
 224 }
 225
 226
 227 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> boost::uint64_t mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::apply_bigskip( boost::uint64_t* Vout, boost::uint64_t* Vin, uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID ){
 228     /*
 229      makes a derived state vector, Vout, from the mother state vector Vin
 230      by skipping a large number of steps, determined by the given seeding ID's
 231
 232      it is mathematically guaranteed that the substreams derived in this way from the SAME (!!!) Vin will not collide provided
 233      1) at least one bit of ID is different
 234      2) less than 10^100 numbers are drawn from the stream
 235      (this is good enough : a single CPU will not exceed this in the lifetime of the universe, 10^19 sec,
 236      even if it had a clock cycle of Planck time, 10^44 Hz )
 237
 238      Caution: never apply this to a derived vector, just choose some mother vector Vin, for example the unit vector by seed_vielbein(X,0),
 239      and use it in all your runs, just change runID to get completely nonoverlapping streams of random numbers on a different day.
 240
 241      clusterID and machineID are provided for the benefit of large organizations who wish to ensure that a simulation
 242      which is running in parallel on a large number of  clusters and machines will have non-colliding source of random numbers.
 243
 244      did i repeat it enough times? the non-collision guarantee is absolute, not probabilistic
 245
 246      */
 247
 248
 249     const    boost::uint64_t skipMat17[128][17] =
 250 #include "boost/random/detail/mixmax_skip_N17.ipp"
 251     ;
 252
 253     const boost::uint64_t* skipMat[128];
 254     BOOST_ASSERT(Ndim==17);
 255     for (int i=0; i<128; i++) { skipMat[i] = skipMat17[i];}
 256
 257     uint32_t IDvec[4] = {streamID, runID, machineID, clusterID};
 258     boost::uint64_t Y[Ndim], cum[Ndim];
 259     boost::uint64_t sumtot=0;
 260
 261     for (int i=0; i<Ndim; i++) { Y[i] = Vin[i]; sumtot = modadd( sumtot, Vin[i]); } ;
 262     for (int IDindex=0; IDindex<4; IDindex++) { // go from lower order to higher order ID
 263         uint32_t id=IDvec[IDindex];
 264         int r = 0;
 265         while (id){
 266             if (id & 1) {
 267                 boost::uint64_t* rowPtr = (boost::uint64_t*)skipMat[r + IDindex*8*sizeof(uint32_t)];
 268                 for (int i=0; i<Ndim; i++){ cum[i] = 0; }
 269                 for (int j=0; j<Ndim; j++){              // j is lag, enumerates terms of the poly
 270                     // for zero lag Y is already given
 271                     boost::uint64_t coeff = rowPtr[j]; // same coeff for all i
 272                     for (int i =0; i<Ndim; i++){
 273                         cum[i] =  fmodmulM61( cum[i], coeff ,  Y[i] ) ;
 274                     }
 275                     sumtot = iterate_raw_vec(Y, sumtot);
 276                 }
 277                 sumtot=0;
 278                 for (int i=0; i<Ndim; i++){ Y[i] = cum[i]; sumtot = modadd( sumtot, cum[i]); } ;
 279             }
 280             id = (id >> 1); r++; // bring up the r-th bit in the ID
 281         }
 282     }
 283     sumtot=0;
 284     for (int i=0; i<Ndim; i++){ Vout[i] = Y[i]; sumtot = modadd( sumtot, Y[i]); } ;  // returns sumtot, and copy the vector over to Vout
 285     return (sumtot) ;
 286 }
 287
 288 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> inline boost::uint64_t mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::fmodmulM61(boost::uint64_t cum, boost::uint64_t s, boost::uint64_t a){
 289     // works on all platforms, including 32-bit linux, PPC and PPC64, ARM and Windows
 290     const boost::uint64_t MASK32=0xFFFFFFFFULL;
 291     boost::uint64_t o,ph,pl,ah,al;
 292     o=(s)*a;
 293     ph = ((s)>>32);
 294     pl = (s) & MASK32;
 295     ah = a>>32;
 296     al = a & MASK32;
 297     o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
 298     o += cum;
 299     o = (o & M61) + ((o>>61));
 300     return o;
 301 }
 302
 303 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> boost::uint64_t mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::modadd(boost::uint64_t foo, boost::uint64_t bar){
 304     return MOD_MERSENNE(foo+bar);
 305 }
 306
 307 /* @copydoc boost::random::detail::mixmax_engine_doc */
 308 /** Instantiation with a valid parameter set. */
 309 typedef mixmax_engine<17,36,0>          mixmax;
 310 }// namespace random
 311 }// namespace boost
 312
 313 #endif // BOOST_RANDOM_MIXMAX_HPP