boost/boost/tr1/random.hpp

   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