1 //  (C) Copyright John Maddock 2005.

     2 //  (C) Copyright Henry S. Warren 2005.

     3 //  Use, modification and distribution are subject to the

     4 //  Boost Software License, Version 1.0. (See accompanying file

     5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

     6 

     7 #ifndef BOOST_TR1_RANDOM_HPP_INCLUDED

     8 #  define BOOST_TR1_RANDOM_HPP_INCLUDED

     9 #  include <boost/tr1/detail/config.hpp>

    10 

    11 #ifdef BOOST_HAS_TR1_RANDOM

    12 #  include BOOST_TR1_HEADER(random)

    13 #else

    14 // Boost.Random:

    15 #include <boost/random.hpp>

    16 #ifndef __SUNPRO_CC

    17     // Sunpros linker complains if we so much as include this...

    18 #   include <boost/nondet_random.hpp>

    19 #endif

    20 #include <boost/tr1/detail/functor2iterator.hpp>

    21 #include <boost/type_traits/is_fundamental.hpp>

    22 #include <boost/type_traits/is_same.hpp>

    23 

    24 namespace std { namespace tr1{

    25 

    26 using ::boost::variate_generator;

    27 

    28 template<class UIntType, UIntType a, UIntType c, UIntType m>

    29 class linear_congruential

    30 {

    31 private:

    32    typedef ::boost::random::linear_congruential<UIntType, a, c, m, 0> impl_type;

    33 public:

    34    // types

    35    typedef UIntType result_type;

    36    // parameter values

    37    BOOST_STATIC_CONSTANT(UIntType, multiplier = a);

    38    BOOST_STATIC_CONSTANT(UIntType, increment = c);

    39    BOOST_STATIC_CONSTANT(UIntType, modulus = m);

    40    // constructors and member function

    41    explicit linear_congruential(unsigned long x0 = 1)

    42       : m_gen(x0){}

    43    linear_congruential(const linear_congruential& that)

    44       : m_gen(that.m_gen){}

    45    template<class Gen> linear_congruential(Gen& g)

    46    {

    47       init1(g, ::boost::is_same<Gen,linear_congruential>());

    48    }

    49    void seed(unsigned long x0 = 1)

    50    { m_gen.seed(x0); }

    51    template<class Gen> void seed(Gen& g)

    52    { 

    53       init2(g, ::boost::is_fundamental<Gen>());

    54    }

    55    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

    56    { return (m_gen.min)(); }

    57    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

    58    { return (m_gen.max)(); }

    59    result_type operator()()

    60    {

    61       return m_gen(); 

    62    }

    63    bool operator==(const linear_congruential& that)const

    64    { return m_gen == that.m_gen; }

    65    bool operator!=(const linear_congruential& that)const

    66    { return m_gen != that.m_gen; }

    67 

    68 #if !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) && !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551))

    69   template<class CharT, class Traits>

    70   friend std::basic_ostream<CharT,Traits>&

    71   operator<<(std::basic_ostream<CharT,Traits>& os,

    72              const linear_congruential& lcg)

    73   {

    74     return os << lcg.m_gen; 

    75   }

    76 

    77   template<class CharT, class Traits>

    78   friend std::basic_istream<CharT,Traits>&

    79   operator>>(std::basic_istream<CharT,Traits>& is,

    80              linear_congruential& lcg)

    81   {

    82     return is >> lcg.m_gen;

    83   }

    84 #endif

    85 

    86 private:

    87    template <class Gen>

    88    void init1(Gen& g, const ::boost::true_type&)

    89    {

    90       m_gen = g.m_gen;

    91    }

    92    template <class Gen>

    93    void init1(Gen& g, const ::boost::false_type&)

    94    {

    95       init2(g, ::boost::is_fundamental<Gen>());

    96    }

    97    template <class Gen>

    98    void init2(Gen& g, const ::boost::true_type&)

    99    {

   100       m_gen.seed(static_cast<unsigned long>(g));

   101    }

   102    template <class Gen>

   103    void init2(Gen& g, const ::boost::false_type&)

   104    {

   105       //typedef typename Gen::result_type gen_rt;

   106       boost::tr1_details::functor2iterator<Gen, unsigned long> f1(g), f2;

   107       m_gen.seed(f1, f2);

   108    }

   109    impl_type m_gen;

   110 };

   111 

   112 template<class UIntType, int w, int n, int m, int r,

   113 UIntType a, int u, int s, UIntType b, int t, UIntType c, int l>

   114 class mersenne_twister

   115 {

   116    typedef ::boost::random::mersenne_twister

   117       <UIntType, w, n, m, r, a, u, s, b, t, c, l, 0> imp_type;

   118 public:

   119    // types

   120    typedef UIntType result_type;

   121    // parameter values

   122    BOOST_STATIC_CONSTANT(int, word_size = w);

   123    BOOST_STATIC_CONSTANT(int, state_size = n);

   124    BOOST_STATIC_CONSTANT(int, shift_size = m);

   125    BOOST_STATIC_CONSTANT(int, mask_bits = r);

   126    BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);

   127    BOOST_STATIC_CONSTANT(int, output_u = u);

   128    BOOST_STATIC_CONSTANT(int, output_s = s);

   129    BOOST_STATIC_CONSTANT(UIntType, output_b = b);

   130    BOOST_STATIC_CONSTANT(int, output_t = t);

   131    BOOST_STATIC_CONSTANT(UIntType, output_c = c);

   132    BOOST_STATIC_CONSTANT(int, output_l = l);

   133    // constructors and member function

   134    mersenne_twister(){}

   135    explicit mersenne_twister(unsigned long value)

   136       : m_gen(value == 0 ? 4357UL : value){}

   137    template<class Gen> mersenne_twister(Gen& g)

   138    {

   139       init1(g, ::boost::is_same<mersenne_twister,Gen>());

   140    }

   141    void seed()

   142    { m_gen.seed(); }

   143    void seed(unsigned long value)

   144    { m_gen.seed(value == 0 ? 5489UL : value); }

   145    template<class Gen> void seed(Gen& g)

   146    { init2(g, ::boost::is_fundamental<Gen>()); }

   147    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

   148    { return (m_gen.min)(); }

   149    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

   150    { return (m_gen.max)(); }

   151    result_type operator()()

   152    { return m_gen(); }

   153    bool operator==(const mersenne_twister& that)const

   154    { return m_gen == that.m_gen; }

   155    bool operator!=(const mersenne_twister& that)const

   156    { return m_gen != that.m_gen; }

   157 

   158 #if !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) && !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551))

   159    template<class CharT, class Traits>

   160    friend std::basic_ostream<CharT,Traits>&

   161    operator<<(std::basic_ostream<CharT,Traits>& os,

   162             const mersenne_twister& lcg)

   163    {

   164       return os << lcg.m_gen;

   165    }

   166 

   167    template<class CharT, class Traits>

   168    friend std::basic_istream<CharT,Traits>&

   169    operator>>(std::basic_istream<CharT,Traits>& is,

   170             mersenne_twister& lcg)

   171    {

   172       return is >> lcg.m_gen;

   173    }

   174 #endif

   175 private:

   176    template <class Gen>

   177    void init1(Gen& g, const ::boost::true_type&)

   178    {

   179       m_gen = g.m_gen;

   180    }

   181    template <class Gen>

   182    void init1(Gen& g, const ::boost::false_type&)

   183    {

   184       init2(g, ::boost::is_fundamental<Gen>());

   185    }

   186    template <class Gen>

   187    void init2(Gen& g, const ::boost::true_type&)

   188    {

   189       m_gen.seed(static_cast<unsigned long>(g == 0 ? 4357UL : g));

   190    }

   191    template <class Gen>

   192    void init2(Gen& g, const ::boost::false_type&)

   193    {

   194       m_gen.seed(g);

   195    }

   196    imp_type m_gen;

   197 };

   198 

   199 template<class IntType, IntType m, int s, int r>

   200 class subtract_with_carry

   201 {

   202 public:

   203    // types

   204    typedef IntType result_type;

   205    // parameter values

   206    BOOST_STATIC_CONSTANT(IntType, modulus = m);

   207    BOOST_STATIC_CONSTANT(int, long_lag = r);

   208    BOOST_STATIC_CONSTANT(int, short_lag = s);

   209 

   210    // constructors and member function

   211    subtract_with_carry(){}

   212    explicit subtract_with_carry(unsigned long value)

   213       : m_gen(value == 0 ? 19780503UL : value){}

   214    template<class Gen> subtract_with_carry(Gen& g)

   215    { init1(g, ::boost::is_same<Gen, subtract_with_carry<IntType, m, s, r> >()); }

   216    void seed(unsigned long value = 19780503ul)

   217    { m_gen.seed(value == 0 ? 19780503UL : value); }

   218    template<class Gen> void seed(Gen& g)

   219    { init2(g, ::boost::is_fundamental<Gen>()); }

   220    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

   221    { return (m_gen.min)(); }

   222    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

   223    { return (m_gen.max)(); }

   224    result_type operator()()

   225    { return m_gen(); }

   226    bool operator==(const subtract_with_carry& that)const

   227    { return m_gen == that.m_gen; }

   228    bool operator!=(const subtract_with_carry& that)const

   229    { return m_gen != that.m_gen; }

   230 

   231 #if !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) && !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551))

   232    template<class CharT, class Traits>

   233    friend std::basic_ostream<CharT,Traits>&

   234    operator<<(std::basic_ostream<CharT,Traits>& os,

   235             const subtract_with_carry& lcg)

   236    {

   237       return os << lcg.m_gen;

   238    }

   239 

   240    template<class CharT, class Traits>

   241    friend std::basic_istream<CharT,Traits>&

   242    operator>>(std::basic_istream<CharT,Traits>& is,

   243             subtract_with_carry& lcg)

   244    {

   245       return is >> lcg.m_gen;

   246    }

   247 #endif

   248 private:

   249    template <class Gen>

   250    void init1(Gen& g, const ::boost::true_type&)

   251    {

   252       m_gen = g.m_gen;

   253    }

   254    template <class Gen>

   255    void init1(Gen& g, const ::boost::false_type&)

   256    {

   257       init2(g, ::boost::is_fundamental<Gen>());

   258    }

   259    template <class Gen>

   260    void init2(Gen& g, const ::boost::true_type&)

   261    {

   262       m_gen.seed(static_cast<unsigned long>(g == 0 ? 19780503UL : g));

   263    }

   264    template <class Gen>

   265    void init2(Gen& g, const ::boost::false_type&)

   266    {

   267       m_gen.seed(g);

   268    }

   269    ::boost::random::subtract_with_carry<IntType, m, s, r, 0> m_gen;

   270 };

   271 

   272 template<class RealType, int w, int s, int r>

   273 class subtract_with_carry_01

   274 {

   275 public:

   276    // types

   277    typedef RealType result_type;

   278    // parameter values

   279    BOOST_STATIC_CONSTANT(int, word_size = w);

   280    BOOST_STATIC_CONSTANT(int, long_lag = r);

   281    BOOST_STATIC_CONSTANT(int, short_lag = s);

   282 

   283    // constructors and member function

   284    subtract_with_carry_01(){}

   285    explicit subtract_with_carry_01(unsigned long value)

   286       : m_gen(value == 0 ? 19780503UL : value){}

   287    template<class Gen> subtract_with_carry_01(Gen& g)

   288    { init1(g, ::boost::is_same<Gen, subtract_with_carry_01<RealType, w, s, r> >()); }

   289    void seed(unsigned long value = 19780503UL)

   290    { m_gen.seed(value == 0 ? 19780503UL : value); }

   291    template<class Gen> void seed(Gen& g)

   292    { init2(g, ::boost::is_fundamental<Gen>()); }

   293    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

   294    { return (m_gen.min)(); }

   295    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

   296    { return (m_gen.max)(); }

   297    result_type operator()()

   298    { return m_gen(); }

   299    bool operator==(const subtract_with_carry_01& that)const

   300    { return m_gen == that.m_gen; }

   301    bool operator!=(const subtract_with_carry_01& that)const

   302    { return m_gen != that.m_gen; }

   303 

   304 #if !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) && !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551))

   305    template<class CharT, class Traits>

   306    friend std::basic_ostream<CharT,Traits>&

   307    operator<<(std::basic_ostream<CharT,Traits>& os,

   308             const subtract_with_carry_01& lcg)

   309    {

   310       return os << lcg.m_gen;

   311    }

   312 

   313    template<class CharT, class Traits>

   314    friend std::basic_istream<CharT,Traits>&

   315    operator>>(std::basic_istream<CharT,Traits>& is,

   316             subtract_with_carry_01& lcg)

   317    {

   318       return is >> lcg.m_gen;

   319    }

   320 #endif

   321 private:

   322    template <class Gen>

   323    void init1(Gen& g, const ::boost::true_type&)

   324    {

   325       m_gen = g.m_gen;

   326    }

   327    template <class Gen>

   328    void init1(Gen& g, const ::boost::false_type&)

   329    {

   330       init2(g, ::boost::is_fundamental<Gen>());

   331    }

   332    template <class Gen>

   333    void init2(Gen& g, const ::boost::true_type&)

   334    {

   335       m_gen.seed(static_cast<unsigned long>(g == 0 ? 19780503UL : g));

   336    }

   337    template <class Gen>

   338    void init2(Gen& g, const ::boost::false_type&)

   339    {

   340       //typedef typename Gen::result_type gen_rt;

   341       boost::tr1_details::functor2iterator<Gen, unsigned long> f1(g), f2;

   342       m_gen.seed(f1, f2);

   343    }

   344    ::boost::random::subtract_with_carry_01<RealType, w, s, r, 0> m_gen;

   345 };

   346 

   347 using ::boost::random::discard_block;

   348 

   349 template<class UniformRandomNumberGenerator1, int s1, class UniformRandomNumberGenerator2, int s2>

   350 class xor_combine

   351 {

   352 public:

   353    // types

   354    typedef UniformRandomNumberGenerator1 base1_type;

   355    typedef UniformRandomNumberGenerator2 base2_type;

   356    typedef unsigned long result_type;

   357    // parameter values

   358    BOOST_STATIC_CONSTANT(int, shift1 = s1);

   359    BOOST_STATIC_CONSTANT(int, shift2 = s2);

   360    // constructors and member function

   361    xor_combine(){ init_minmax(); }

   362    xor_combine(const base1_type & rng1, const base2_type & rng2)

   363       : m_b1(rng1), m_b2(rng2) { init_minmax(); }

   364    xor_combine(unsigned long s)

   365       : m_b1(s), m_b2(s+1) { init_minmax(); }

   366    template<class Gen> xor_combine(Gen& g)

   367    { 

   368       init_minmax(); 

   369       init1(g, ::boost::is_same<Gen, xor_combine<UniformRandomNumberGenerator1, s1, UniformRandomNumberGenerator2, s2> >());

   370    }

   371    void seed()

   372    {

   373       m_b1.seed();

   374       m_b2.seed();

   375    }

   376    void seed(unsigned long s)

   377    {

   378       m_b1.seed(s);

   379       m_b2.seed(s+1);

   380    }

   381    template<class Gen> void seed(Gen& g)

   382    {

   383       init2(g, ::boost::is_fundamental<Gen>());

   384    }

   385 

   386    const base1_type& base1() const

   387    { return m_b1; }

   388    const base2_type& base2() const

   389    { return m_b2; }

   390    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

   391    { return m_min; }

   392    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

   393    { return m_max; }

   394    result_type operator()()

   395    { return (m_b1() << s1) ^ (m_b2() << s2); }

   396 

   397    bool operator == (const xor_combine& that)const

   398    { return (m_b1 == that.m_b1) && (m_b2 == that.m_b2); }

   399    bool operator != (const xor_combine& that)const

   400    { return !(*this == that); }

   401 

   402 #if !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) && !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551))

   403    template<class CharT, class Traits>

   404    friend std::basic_ostream<CharT,Traits>&

   405    operator<<(std::basic_ostream<CharT,Traits>& os,

   406             const xor_combine& lcg)

   407    {

   408       return os << lcg.m_b1 << " " << lcg.m_b2;

   409    }

   410 

   411    template<class CharT, class Traits>

   412    friend std::basic_istream<CharT,Traits>&

   413    operator>>(std::basic_istream<CharT,Traits>& is,

   414             xor_combine& lcg)

   415    {

   416       return is >> lcg.m_b1 >> lcg.m_b2;

   417    }

   418 #endif

   419 

   420 private:

   421    void init_minmax();

   422    base1_type m_b1;

   423    base2_type m_b2;

   424    result_type m_min;

   425    result_type m_max;

   426 

   427    template <class Gen>

   428    void init1(Gen& g, const ::boost::true_type&)

   429    {

   430       m_b1 = g.m_b1;

   431       m_b2 = g.m_b2;

   432    }

   433    template <class Gen>

   434    void init1(Gen& g, const ::boost::false_type&)

   435    {

   436       init2(g, ::boost::is_fundamental<Gen>());

   437    }

   438    template <class Gen>

   439    void init2(Gen& g, const ::boost::true_type&)

   440    {

   441       m_b1.seed(static_cast<unsigned long>(g));

   442       m_b2.seed(static_cast<unsigned long>(g));

   443    }

   444    template <class Gen>

   445    void init2(Gen& g, const ::boost::false_type&)

   446    {

   447       m_b1.seed(g);

   448       m_b2.seed(g);

   449    }

   450 };

   451 

   452 template<class UniformRandomNumberGenerator1, int s1, class UniformRandomNumberGenerator2, int s2>

   453 void xor_combine<UniformRandomNumberGenerator1, s1, UniformRandomNumberGenerator2, s2>::init_minmax()

   454 {

   455    //

   456    // The following code is based on that given in "Hacker's Delight"

   457    // by Henry S. Warren, (Addison-Wesley, 2003), and at 

   458    // http://www.hackersdelight.org/index.htm.

   459    // Used here by permission.

   460    //

   461    // calculation of minimum value:

   462    //

   463    result_type a = (m_b1.min)() << s1;

   464    result_type b = (m_b1.max)() << s1;

   465    result_type c = (m_b2.min)() << s2;

   466    result_type d = (m_b2.max)() << s2;

   467    result_type m, temp;

   468 

   469    m = 0x1uL << ((sizeof(result_type) * CHAR_BIT) - 1);

   470    while (m != 0) {

   471       if (~a & c & m) {

   472          temp = (a | m) & (static_cast<result_type>(0u) - m);

   473          if (temp <= b) a = temp;

   474       }

   475       else if (a & ~c & m) {

   476          temp = (c | m) & (static_cast<result_type>(0u) - m);

   477          if (temp <= d) c = temp;

   478       }

   479       m >>= 1;

   480    }

   481    m_min = a ^ c;

   482 

   483    //

   484    // calculation of maximum value:

   485    //

   486    if((((std::numeric_limits<result_type>::max)() >> s1) < (m_b1.max)())

   487       || ((((std::numeric_limits<result_type>::max)()) >> s2) < (m_b2.max)()))

   488    {

   489       m_max = (std::numeric_limits<result_type>::max)();

   490       return;

   491    }

   492    a = (m_b1.min)() << s1;

   493    b = (m_b1.max)() << s1;

   494    c = (m_b2.min)() << s2;

   495    d = (m_b2.max)() << s2;

   496 

   497    m = 0x1uL << ((sizeof(result_type) * CHAR_BIT) - 1);

   498 

   499    while (m != 0) {

   500       if (b & d & m) {

   501          temp = (b - m) | (m - 1);

   502          if (temp >= a) b = temp;

   503          else {

   504             temp = (d - m) | (m - 1);

   505             if (temp >= c) d = temp;

   506          }

   507       }

   508       m = m >> 1;

   509    }

   510    m_max = b ^ d;

   511 }

   512 

   513 typedef linear_congruential< ::boost::int32_t, 16807, 0, 2147483647> minstd_rand0;

   514 typedef linear_congruential< ::boost::int32_t, 48271, 0, 2147483647> minstd_rand;

   515 typedef mersenne_twister< ::boost::uint32_t, 32,624,397,31,0x9908b0df,11,7,0x9d2c5680,15,0xefc60000,18> mt19937;

   516 typedef subtract_with_carry_01<float, 24, 10, 24> ranlux_base_01;

   517 typedef subtract_with_carry_01<double, 48, 10, 24> ranlux64_base_01;

   518 typedef discard_block<subtract_with_carry< ::boost::int32_t, (1<<24), 10, 24>, 223, 24> ranlux3;

   519 typedef discard_block<subtract_with_carry< ::boost::int32_t, (1<<24), 10, 24>, 389, 24> ranlux4;

   520 typedef discard_block<subtract_with_carry_01<float, 24, 10, 24>, 223, 24> ranlux3_01;

   521 typedef discard_block<subtract_with_carry_01<float, 24, 10, 24>, 389, 24> ranlux4_01;

   522 

   523 #ifndef __SUNPRO_CC

   524 using ::boost::random_device;

   525 #endif

   526 using ::boost::uniform_int;

   527 

   528 class bernoulli_distribution

   529 {

   530 public:

   531    // types

   532    typedef int input_type;

   533    typedef bool result_type;

   534    // constructors and member function

   535    explicit bernoulli_distribution(double p = 0.5)

   536       : m_dist(p){}

   537    double p() const

   538    { return m_dist.p(); }

   539    void reset()

   540    { m_dist.reset(); }

   541    template<class UniformRandomNumberGenerator>

   542    result_type operator()(UniformRandomNumberGenerator& urng)

   543    {

   544       return m_dist(urng);

   545    }

   546 #if !defined(BOOST_NO_MEMBER_TEMPLATE_FRIENDS) && !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x551))

   547    template<class CharT, class Traits>

   548    friend std::basic_ostream<CharT,Traits>&

   549    operator<<(std::basic_ostream<CharT,Traits>& os,

   550             const bernoulli_distribution& lcg)

   551    {

   552       return os << lcg.m_dist;

   553    }

   554 

   555    template<class CharT, class Traits>

   556    friend std::basic_istream<CharT,Traits>&

   557    operator>>(std::basic_istream<CharT,Traits>& is,

   558             bernoulli_distribution& lcg)

   559    {

   560       return is >> lcg.m_dist;

   561    }

   562 #endif

   563 

   564 private:

   565    ::boost::bernoulli_distribution<double> m_dist;

   566 };

   567 //using ::boost::bernoulli_distribution;

   568 using ::boost::geometric_distribution;

   569 using ::boost::poisson_distribution;

   570 using ::boost::binomial_distribution;

   571 using ::boost::uniform_real;

   572 using ::boost::exponential_distribution;

   573 using ::boost::normal_distribution;

   574 using ::boost::gamma_distribution;

   575 

   576 } }

   577 

   578 #endif

   579 

   580 #endif

   581