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 #  ifdef BOOST_HAS_INCLUDE_NEXT

    13 #     include_next BOOST_TR1_HEADER(random)

    14 #  else

    15 #     include <boost/tr1/detail/config_all.hpp>

    16 #     include BOOST_TR1_STD_HEADER(BOOST_TR1_PATH(random))

    17 #  endif

    18 #else

    19 // Boost.Random:

    20 #include <boost/random.hpp>

    21 #ifndef __SUNPRO_CC

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

    23 #   include <boost/nondet_random.hpp>

    24 #endif

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

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

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

    28 

    29 namespace std { namespace tr1{

    30 

    31 using ::boost::variate_generator;

    32 

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

    34 class linear_congruential

    35 {

    36 private:

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

    38 public:

    39    // types

    40    typedef UIntType result_type;

    41    // parameter values

    42    BOOST_STATIC_CONSTANT(UIntType, multiplier = a);

    43    BOOST_STATIC_CONSTANT(UIntType, increment = c);

    44    BOOST_STATIC_CONSTANT(UIntType, modulus = m);

    45    // constructors and member function

    46    explicit linear_congruential(unsigned long x0 = 1)

    47       : m_gen(x0){}

    48    linear_congruential(const linear_congruential& that)

    49       : m_gen(that.m_gen){}

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

    51    {

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

    53    }

    54    void seed(unsigned long x0 = 1)

    55    { m_gen.seed(x0); }

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

    57    { 

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

    59    }

    60    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

    62    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

    64    result_type operator()()

    65    {

    66       return m_gen(); 

    67    }

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

    69    { return m_gen == that.m_gen; }

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

    71    { return m_gen != that.m_gen; }

    72 

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

    74   template<class CharT, class Traits>

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

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

    77              const linear_congruential& lcg)

    78   {

    79     return os << lcg.m_gen; 

    80   }

    81 

    82   template<class CharT, class Traits>

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

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

    85              linear_congruential& lcg)

    86   {

    87     return is >> lcg.m_gen;

    88   }

    89 #endif

    90 

    91 private:

    92    template <class Gen>

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

    94    {

    95       m_gen = g.m_gen;

    96    }

    97    template <class Gen>

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

    99    {

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

   101    }

   102    template <class Gen>

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

   104    {

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

   106    }

   107    template <class Gen>

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

   109    {

   110       //typedef typename Gen::result_type gen_rt;

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

   112       m_gen.seed(f1, f2);

   113    }

   114    impl_type m_gen;

   115 };

   116 

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

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

   119 class mersenne_twister

   120 {

   121    typedef ::boost::random::mersenne_twister

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

   123 public:

   124    // types

   125    typedef UIntType result_type;

   126    // parameter values

   127    BOOST_STATIC_CONSTANT(int, word_size = w);

   128    BOOST_STATIC_CONSTANT(int, state_size = n);

   129    BOOST_STATIC_CONSTANT(int, shift_size = m);

   130    BOOST_STATIC_CONSTANT(int, mask_bits = r);

   131    BOOST_STATIC_CONSTANT(UIntType, parameter_a = a);

   132    BOOST_STATIC_CONSTANT(int, output_u = u);

   133    BOOST_STATIC_CONSTANT(int, output_s = s);

   134    BOOST_STATIC_CONSTANT(UIntType, output_b = b);

   135    BOOST_STATIC_CONSTANT(int, output_t = t);

   136    BOOST_STATIC_CONSTANT(UIntType, output_c = c);

   137    BOOST_STATIC_CONSTANT(int, output_l = l);

   138    // constructors and member function

   139    mersenne_twister(){}

   140    explicit mersenne_twister(unsigned long value)

   141       : m_gen(value == 0 ? 5489UL : value){}

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

   143    {

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

   145    }

   146    void seed()

   147    { m_gen.seed(); }

   148    void seed(unsigned long value)

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

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

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

   152    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

   154    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

   156    result_type operator()()

   157    { return m_gen(); }

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

   159    { return m_gen == that.m_gen; }

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

   161    { return m_gen != that.m_gen; }

   162 

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

   164    template<class CharT, class Traits>

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

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

   167             const mersenne_twister& lcg)

   168    {

   169       return os << lcg.m_gen;

   170    }

   171 

   172    template<class CharT, class Traits>

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

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

   175             mersenne_twister& lcg)

   176    {

   177       return is >> lcg.m_gen;

   178    }

   179 #endif

   180 private:

   181    template <class Gen>

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

   183    {

   184       m_gen = g.m_gen;

   185    }

   186    template <class Gen>

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

   188    {

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

   190    }

   191    template <class Gen>

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

   193    {

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

   195    }

   196    template <class Gen>

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

   198    {

   199       m_gen.seed(g);

   200    }

   201    imp_type m_gen;

   202 };

   203 

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

   205 class subtract_with_carry

   206 {

   207 public:

   208    // types

   209    typedef IntType result_type;

   210    // parameter values

   211    BOOST_STATIC_CONSTANT(IntType, modulus = m);

   212    BOOST_STATIC_CONSTANT(int, long_lag = r);

   213    BOOST_STATIC_CONSTANT(int, short_lag = s);

   214 

   215    // constructors and member function

   216    subtract_with_carry(){}

   217    explicit subtract_with_carry(unsigned long value)

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

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

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

   221    void seed(unsigned long value = 19780503ul)

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

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

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

   225    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

   227    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

   229    result_type operator()()

   230    { return m_gen(); }

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

   232    { return m_gen == that.m_gen; }

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

   234    { return m_gen != that.m_gen; }

   235 

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

   237    template<class CharT, class Traits>

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

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

   240             const subtract_with_carry& lcg)

   241    {

   242       return os << lcg.m_gen;

   243    }

   244 

   245    template<class CharT, class Traits>

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

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

   248             subtract_with_carry& lcg)

   249    {

   250       return is >> lcg.m_gen;

   251    }

   252 #endif

   253 private:

   254    template <class Gen>

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

   256    {

   257       m_gen = g.m_gen;

   258    }

   259    template <class Gen>

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

   261    {

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

   263    }

   264    template <class Gen>

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

   266    {

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

   268    }

   269    template <class Gen>

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

   271    {

   272       m_gen.seed(g);

   273    }

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

   275 };

   276 

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

   278 class subtract_with_carry_01

   279 {

   280 public:

   281    // types

   282    typedef RealType result_type;

   283    // parameter values

   284    BOOST_STATIC_CONSTANT(int, word_size = w);

   285    BOOST_STATIC_CONSTANT(int, long_lag = r);

   286    BOOST_STATIC_CONSTANT(int, short_lag = s);

   287 

   288    // constructors and member function

   289    subtract_with_carry_01(){}

   290    explicit subtract_with_carry_01(unsigned long value)

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

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

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

   294    void seed(unsigned long value = 19780503UL)

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

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

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

   298    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

   300    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

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

   302    result_type operator()()

   303    { return m_gen(); }

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

   305    { return m_gen == that.m_gen; }

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

   307    { return m_gen != that.m_gen; }

   308 

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

   310    template<class CharT, class Traits>

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

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

   313             const subtract_with_carry_01& lcg)

   314    {

   315       return os << lcg.m_gen;

   316    }

   317 

   318    template<class CharT, class Traits>

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

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

   321             subtract_with_carry_01& lcg)

   322    {

   323       return is >> lcg.m_gen;

   324    }

   325 #endif

   326 private:

   327    template <class Gen>

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

   329    {

   330       m_gen = g.m_gen;

   331    }

   332    template <class Gen>

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

   334    {

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

   336    }

   337    template <class Gen>

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

   339    {

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

   341    }

   342    template <class Gen>

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

   344    {

   345       //typedef typename Gen::result_type gen_rt;

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

   347       m_gen.seed(f1, f2);

   348    }

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

   350 };

   351 

   352 using ::boost::random::discard_block;

   353 

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

   355 class xor_combine

   356 {

   357 public:

   358    // types

   359    typedef UniformRandomNumberGenerator1 base1_type;

   360    typedef UniformRandomNumberGenerator2 base2_type;

   361    typedef unsigned long result_type;

   362    // parameter values

   363    BOOST_STATIC_CONSTANT(int, shift1 = s1);

   364    BOOST_STATIC_CONSTANT(int, shift2 = s2);

   365    // constructors and member function

   366    xor_combine(){ init_minmax(); }

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

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

   369    xor_combine(unsigned long s)

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

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

   372    { 

   373       init_minmax(); 

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

   375    }

   376    void seed()

   377    {

   378       m_b1.seed();

   379       m_b2.seed();

   380    }

   381    void seed(unsigned long s)

   382    {

   383       m_b1.seed(s);

   384       m_b2.seed(s+1);

   385    }

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

   387    {

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

   389    }

   390 

   391    const base1_type& base1() const

   392    { return m_b1; }

   393    const base2_type& base2() const

   394    { return m_b2; }

   395    result_type min BOOST_PREVENT_MACRO_SUBSTITUTION() const

   396    { return m_min; }

   397    result_type max BOOST_PREVENT_MACRO_SUBSTITUTION() const

   398    { return m_max; }

   399    result_type operator()()

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

   401 

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

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

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

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

   406 

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

   408    template<class CharT, class Traits>

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

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

   411             const xor_combine& lcg)

   412    {

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

   414    }

   415 

   416    template<class CharT, class Traits>

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

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

   419             xor_combine& lcg)

   420    {

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

   422    }

   423 #endif

   424 

   425 private:

   426    void init_minmax();

   427    base1_type m_b1;

   428    base2_type m_b2;

   429    result_type m_min;

   430    result_type m_max;

   431 

   432    template <class Gen>

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

   434    {

   435       m_b1 = g.m_b1;

   436       m_b2 = g.m_b2;

   437    }

   438    template <class Gen>

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

   440    {

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

   442    }

   443    template <class Gen>

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

   445    {

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

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

   448    }

   449    template <class Gen>

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

   451    {

   452       m_b1.seed(g);

   453       m_b2.seed(g);

   454    }

   455 };

   456 

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

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

   459 {

   460    //

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

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

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

   464    // Used here by permission.

   465    //

   466    // calculation of minimum value:

   467    //

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

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

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

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

   472    result_type m, temp;

   473 

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

   475    while (m != 0) {

   476       if (~a & c & m) {

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

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

   479       }

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

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

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

   483       }

   484       m >>= 1;

   485    }

   486    m_min = a ^ c;

   487 

   488    //

   489    // calculation of maximum value:

   490    //

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

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

   493    {

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

   495       return;

   496    }

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

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

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

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

   501 

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

   503 

   504    while (m != 0) {

   505       if (b & d & m) {

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

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

   508          else {

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

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

   511          }

   512       }

   513       m = m >> 1;

   514    }

   515    m_max = b ^ d;

   516 }

   517 

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

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

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

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

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

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

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

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

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

   527 

   528 #ifndef __SUNPRO_CC

   529 using ::boost::random_device;

   530 #endif

   531 using ::boost::uniform_int;

   532 

   533 class bernoulli_distribution

   534 {

   535 public:

   536    // types

   537    typedef int input_type;

   538    typedef bool result_type;

   539    // constructors and member function

   540    explicit bernoulli_distribution(double p = 0.5)

   541       : m_dist(p){}

   542    double p() const

   543    { return m_dist.p(); }

   544    void reset()

   545    { m_dist.reset(); }

   546    template<class UniformRandomNumberGenerator>

   547    result_type operator()(UniformRandomNumberGenerator& urng)

   548    {

   549       return m_dist(urng);

   550    }

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

   552    template<class CharT, class Traits>

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

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

   555             const bernoulli_distribution& lcg)

   556    {

   557       return os << lcg.m_dist;

   558    }

   559 

   560    template<class CharT, class Traits>

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

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

   563             bernoulli_distribution& lcg)

   564    {

   565       return is >> lcg.m_dist;

   566    }

   567 #endif

   568 

   569 private:

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

   571 };

   572 //using ::boost::bernoulli_distribution;

   573 using ::boost::geometric_distribution;

   574 using ::boost::poisson_distribution;

   575 using ::boost::binomial_distribution;

   576 using ::boost::uniform_real;

   577 using ::boost::exponential_distribution;

   578 using ::boost::normal_distribution;

   579 using ::boost::gamma_distribution;

   580 

   581 } }

   582 

   583 #endif

   584 

   585 #endif

   586