boost/random/mixmax.hpp

0001 /* boost random/mixmax.hpp header file
0002  *
0003  * Copyright Kostas Savvidis 2008-2019
0004  *
0005  * Distributed under the Boost Software License, Version 1.0. (See
0006  * accompanying file LICENSE_1_0.txt or copy at
0007  * http://www.boost.org/LICENSE_1_0.txt)
0008  *
0009  * See http://www.boost.org for most recent version including documentation.
0010  *
0011  * $Id$
0012  *
0013  * Revision history
0014  *  2019-04-23 created
0015  */
0016
0017 #ifndef BOOST_RANDOM_MIXMAX_HPP
0018 #define BOOST_RANDOM_MIXMAX_HPP
0019
0020 #include <sstream>
0021 #include <boost/cstdint.hpp>
0022 #include <boost/array.hpp>
0023
0024 #include <boost/random/detail/seed.hpp>
0025 #include <boost/random/detail/seed_impl.hpp>
0026
0027 namespace boost {
0028 namespace random {
0029
0030 /**
0031  * Instantiations of class template mixmax_engine model,
0032  * \pseudo_random_number_generator .
0033  *  It uses the  MIXMAX generator algorithms from:
0034  *
0035  *  @blockquote
0036  *  G.K.Savvidy and N.G.Ter-Arutyunian,
0037  *  On the Monte Carlo simulation of physical systems,
0038  *  J.Comput.Phys. 97, 566 (1991);
0039  *  Preprint EPI-865-16-86, Yerevan, Jan. 1986
0040  *  http://dx.doi.org/10.1016/0021-9991(91)90015-D
0041  *
0042  *  K.Savvidy
0043  *  The MIXMAX random number generator
0044  *  Comp. Phys. Commun. 196 (2015), pp 161–165
0045  *  http://dx.doi.org/10.1016/j.cpc.2015.06.003
0046  *
0047  *  K.Savvidy and G.Savvidy
0048  *  Spectrum and Entropy of C-systems. MIXMAX random number generator
0049  *  Chaos, Solitons & Fractals, Volume 91, (2016) pp. 33–38
0050  *  http://dx.doi.org/10.1016/j.chaos.2016.05.003
0051  *  @endblockquote
0052  *
0053  * The generator crucially depends on the choice of the
0054  * parameters. The valid sets of parameters are from the published papers above.
0055  *
0056  */
0057
0058 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> // MIXMAX TEMPLATE PARAMETERS
0059 class mixmax_engine{
0060 public:
0061     // Interfaces required by C++11 std::random and boost::random
0062     typedef boost::uint64_t result_type ;
0063     BOOST_STATIC_CONSTANT(boost::uint64_t,mixmax_min=0);
0064     BOOST_STATIC_CONSTANT(boost::uint64_t,mixmax_max=((1ULL<<61)-1));
0065     BOOST_STATIC_CONSTEXPR result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() {return mixmax_min;}
0066     BOOST_STATIC_CONSTEXPR result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() {return mixmax_max;}
0067     static const bool has_fixed_range = false;
0068     BOOST_STATIC_CONSTANT(int,N=Ndim);     ///< The main internal parameter, size of the defining MIXMAX matrix
0069     // CONSTRUCTORS:
0070     explicit mixmax_engine();                       ///< Constructor, unit vector as initial state, acted on by A^2^512
0071     explicit mixmax_engine(boost::uint64_t);          ///< Constructor, one 64-bit seed
0072     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
0073     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
0074
0075 private: // DATATYPES
0076     struct rng_state_st{
0077         boost::array<boost::uint64_t, Ndim> V;
0078         boost::uint64_t sumtot;
0079         int counter;
0080     };
0081
0082     typedef struct rng_state_st rng_state_t;     // struct alias
0083     rng_state_t S;
0084
0085 public: // SEEDING FUNCTIONS
0086     template<class It> mixmax_engine(It& first, It last) { seed(first,last); }
0087     BOOST_RANDOM_DETAIL_SEED_SEQ_CONSTRUCTOR(mixmax_engine,  SeedSeq, seq){ seed(seq); }
0088
0089     /** Sets the state of the generator using values from an iterator range. */
0090     template<class It>
0091     void seed(It& first, It last){
0092         uint32_t v[4];
0093         detail::fill_array_int<32>(first, last, v);
0094         seed_uniquestream( &S, v[0], v[1], v[2], v[3]);
0095     }
0096     /** Sets the state of the generator using values from a seed_seq. */
0097     BOOST_RANDOM_DETAIL_SEED_SEQ_SEED(mixmax_engine, SeeqSeq, seq){
0098         uint32_t v[4];
0099         detail::seed_array_int<32>(seq, v);
0100         seed_uniquestream( &S, v[0], v[1], v[2], v[3]);
0101     }
0102
0103     /** return one uint64 between min=0 and max=2^61-1 */
0104     boost::uint64_t operator()(){
0105         if (S.counter<=(Ndim-1) ){
0106             return S.V[S.counter++];
0107         }else{
0108             S.sumtot = iterate_raw_vec(S.V.data(), S.sumtot);
0109             S.counter=2;
0110             return S.V[1];
0111         }
0112     }
0113
0114     /** Fills a range with random values */
0115     template<class Iter>
0116     void generate(Iter first, Iter last) { detail::generate_from_int(*this, first, last); }
0117
0118     void discard(boost::uint64_t nsteps) { for(boost::uint64_t j = 0; j < nsteps; ++j)  (*this)(); } ///< discard n steps, required in boost::random
0119
0120     /** save the state of the RNG to a stream */
0121     template<class CharT, class Traits>
0122     friend std::basic_ostream<CharT,Traits>&
0123     operator<< (std::basic_ostream<CharT,Traits>& ost, const mixmax_engine& me){
0124         ost << Ndim << " " << me.S.counter << " " << me.S.sumtot << " ";
0125         for (int j=0; (j< (Ndim) ); j++) {
0126         ost <<  (boost::uint64_t)me.S.V[j] << " ";
0127         }
0128         ost << "\n";
0129         ost.flush();
0130         return ost;
0131         }
0132
0133     /** read the state of the RNG from a stream */
0134     template<class CharT, class Traits>
0135     friend std::basic_istream<CharT,Traits>&
0136     operator>> (std::basic_istream<CharT,Traits> &in, mixmax_engine& me){
0137         // will set std::ios::failbit if the input format is not right
0138         boost::array<boost::uint64_t, Ndim> vec;
0139         boost::uint64_t sum=0, savedsum=0, counter=0;
0140         in >> counter >> std::ws;
0141         BOOST_ASSERT(counter==Ndim);
0142         in >> counter >> std::ws;
0143         in >> savedsum >> std::ws;
0144         for(int j=0;j<Ndim;j++) {
0145         in >> std::ws >> vec[j] ;
0146         sum=me.MOD_MERSENNE(sum+vec[j]);
0147     }
0148     if (sum == savedsum && counter>0 && counter<Ndim){
0149         me.S.V=vec; me.S.counter = counter; me.S.sumtot=savedsum;
0150     }else{
0151         in.setstate(std::ios::failbit);
0152     }
0153     return in;
0154     }
0155
0156 friend bool operator==(const mixmax_engine & x,
0157                        const mixmax_engine & y){return x.S.counter==y.S.counter && x.S.sumtot==y.S.sumtot && x.S.V==y.S.V ;}
0158 friend bool operator!=(const mixmax_engine & x,
0159                        const mixmax_engine & y){return !operator==(x,y);}
0160
0161
0162 private:
0163 BOOST_STATIC_CONSTANT(int, BITS=61);
0164 BOOST_STATIC_CONSTANT(boost::uint64_t, M61=2305843009213693951ULL);
0165 BOOST_STATIC_CONSTANT(boost::uint64_t, default_seed=1);
0166 inline boost::uint64_t MOD_MERSENNE(boost::uint64_t k) {return ((((k)) & M61) + (((k)) >> BITS) );}
0167 inline boost::uint64_t MULWU(boost::uint64_t k);
0168 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
0169 inline void seed_uniquestream( rng_state_t* Xin, uint32_t clusterID, uint32_t machineID, uint32_t runID, uint32_t  streamID );
0170 inline boost::uint64_t iterate_raw_vec(boost::uint64_t* Y, boost::uint64_t sumtotOld);
0171 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 );
0172 inline boost::uint64_t modadd(boost::uint64_t foo, boost::uint64_t bar);
0173 inline boost::uint64_t fmodmulM61(boost::uint64_t cum, boost::uint64_t s, boost::uint64_t a);
0174 };
0175
0176 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::mixmax_engine()
0177 ///< constructor, with no params, seeds with seed=0,  random numbers are as good as from any other seed
0178 {
0179     seed_uniquestream( &S, 0,  0, 0, default_seed);
0180 }
0181
0182 template <int Ndim, unsigned int SPECIALMUL, boost::int64_t SPECIAL> mixmax_engine  <Ndim, SPECIALMUL, SPECIAL> ::mixmax_engine(boost::uint64_t seedval){
0183     ///< constructor, one uint64_t seed, random numbers are statistically independent from any two distinct seeds, e.g. consecutive seeds are ok
0184     seed_uniquestream( &S, 0,  0,  (uint32_t)(seedval>>32), (uint32_t)seedval );
0185 }
0186
0187 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){
0188     // constructor, four 32-bit seeds for 128-bit seeding flexibility
0189     seed_uniquestream( &S, clusterID,  machineID,  runID,  streamID );
0190 }
0191
0192 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))  )  ;}
0193
0194 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){
0195     // operates with a raw vector, uses known sum of elements of Y
0196     boost::uint64_t  tempP=0, tempV=sumtotOld;
0197     Y[0] = tempV;
0198     boost::uint64_t sumtot = Y[0], ovflow = 0; // will keep a running sum of all new elements
0199     for (int i=1; i<Ndim; i++){
0200         boost::uint64_t tempPO = MULWU(tempP);
0201         tempV = (tempV+tempPO);
0202         tempP = modadd(tempP, Y[i]);
0203         tempV = modadd(tempV, tempP); // new Y[i] = old Y[i] + old partial * m
0204         Y[i] = tempV;
0205         sumtot += tempV; if (sumtot < tempV) {ovflow++;}
0206     }
0207     return MOD_MERSENNE(MOD_MERSENNE(sumtot) + (ovflow <<3 ));
0208 }
0209
0210 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){
0211     for (int i=0; i < Ndim; i++){
0212         X->V[i] = 0;
0213     }
0214     if (index<Ndim) { X->V[index] = 1; }else{ X->V[0]=1; }
0215     X->counter = Ndim;  // set the counter to Ndim if iteration should happen right away
0216     X->sumtot = 1;
0217 }
0218
0219
0220 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 ){
0221     seed_vielbein(Xin,0);
0222     Xin->sumtot = apply_bigskip(Xin->V.data(), Xin->V.data(),  clusterID,  machineID,  runID,   streamID );
0223     Xin->counter = 1;
0224 }
0225
0226
0227 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 ){
0228     /*
0229      makes a derived state vector, Vout, from the mother state vector Vin
0230      by skipping a large number of steps, determined by the given seeding ID's
0231
0232      it is mathematically guaranteed that the substreams derived in this way from the SAME (!!!) Vin will not collide provided
0233      1) at least one bit of ID is different
0234      2) less than 10^100 numbers are drawn from the stream
0235      (this is good enough : a single CPU will not exceed this in the lifetime of the universe, 10^19 sec,
0236      even if it had a clock cycle of Planck time, 10^44 Hz )
0237
0238      Caution: never apply this to a derived vector, just choose some mother vector Vin, for example the unit vector by seed_vielbein(X,0),
0239      and use it in all your runs, just change runID to get completely nonoverlapping streams of random numbers on a different day.
0240
0241      clusterID and machineID are provided for the benefit of large organizations who wish to ensure that a simulation
0242      which is running in parallel on a large number of  clusters and machines will have non-colliding source of random numbers.
0243
0244      did i repeat it enough times? the non-collision guarantee is absolute, not probabilistic
0245
0246      */
0247
0248
0249     const    boost::uint64_t skipMat17[128][17] =
0250 #include "boost/random/detail/mixmax_skip_N17.ipp"
0251     ;
0252
0253     const boost::uint64_t* skipMat[128];
0254     BOOST_ASSERT(Ndim==17);
0255     for (int i=0; i<128; i++) { skipMat[i] = skipMat17[i];}
0256
0257     uint32_t IDvec[4] = {streamID, runID, machineID, clusterID};
0258     boost::uint64_t Y[Ndim], cum[Ndim];
0259     boost::uint64_t sumtot=0;
0260
0261     for (int i=0; i<Ndim; i++) { Y[i] = Vin[i]; sumtot = modadd( sumtot, Vin[i]); } ;
0262     for (int IDindex=0; IDindex<4; IDindex++) { // go from lower order to higher order ID
0263         uint32_t id=IDvec[IDindex];
0264         int r = 0;
0265         while (id){
0266             if (id & 1) {
0267                 boost::uint64_t* rowPtr = (boost::uint64_t*)skipMat[r + IDindex*8*sizeof(uint32_t)];
0268                 for (int i=0; i<Ndim; i++){ cum[i] = 0; }
0269                 for (int j=0; j<Ndim; j++){              // j is lag, enumerates terms of the poly
0270                     // for zero lag Y is already given
0271                     boost::uint64_t coeff = rowPtr[j]; // same coeff for all i
0272                     for (int i =0; i<Ndim; i++){
0273                         cum[i] =  fmodmulM61( cum[i], coeff ,  Y[i] ) ;
0274                     }
0275                     sumtot = iterate_raw_vec(Y, sumtot);
0276                 }
0277                 sumtot=0;
0278                 for (int i=0; i<Ndim; i++){ Y[i] = cum[i]; sumtot = modadd( sumtot, cum[i]); } ;
0279             }
0280             id = (id >> 1); r++; // bring up the r-th bit in the ID
0281         }
0282     }
0283     sumtot=0;
0284     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
0285     return (sumtot) ;
0286 }
0287
0288 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){
0289     // works on all platforms, including 32-bit linux, PPC and PPC64, ARM and Windows
0290     const boost::uint64_t MASK32=0xFFFFFFFFULL;
0291     boost::uint64_t o,ph,pl,ah,al;
0292     o=(s)*a;
0293     ph = ((s)>>32);
0294     pl = (s) & MASK32;
0295     ah = a>>32;
0296     al = a & MASK32;
0297     o = (o & M61) + ((ph*ah)<<3) + ((ah*pl+al*ph + ((al*pl)>>32))>>29) ;
0298     o += cum;
0299     o = (o & M61) + ((o>>61));
0300     return o;
0301 }
0302
0303 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){
0304     return MOD_MERSENNE(foo+bar);
0305 }
0306
0307 /* @copydoc boost::random::detail::mixmax_engine_doc */
0308 /** Instantiation with a valid parameter set. */
0309 typedef mixmax_engine<17,36,0>          mixmax;
0310 }// namespace random
0311 }// namespace boost
0312
0313 #endif // BOOST_RANDOM_MIXMAX_HPP