libs/random/performance/random_speed.cpp

   1 /* boost random_speed.cpp performance measurements
   2  *
   3  * Copyright Jens Maurer 2000
   4  * Distributed under the Boost Software License, Version 1.0. (See
   5  * accompanying file LICENSE_1_0.txt or copy at
   6  * http://www.boost.org/LICENSE_1_0.txt)
   7  *
   8  * $Id$
   9  */
  10
  11 #include <iostream>
  12 #include <cstdlib>
  13 #include <string>
  14 #include <boost/config.hpp>
  15 #include <boost/random.hpp>
  16 #include <boost/progress.hpp>
  17 #include <boost/shared_ptr.hpp>
  18
  19 /*
  20  * Configuration Section
  21  */
  22
  23 // define if your C library supports the non-standard drand48 family
  24 //#define HAVE_DRAND48
  25
  26 // define if you have the original mt19937int.c (with commented out main())
  27 #undef HAVE_MT19937INT_C
  28
  29 // define if you have the original mt19937ar.c (with commented out main())
  30 // #define HAVE_MT19937AR_C
  31
  32 // set to your CPU frequency
  33 static const double cpu_frequency = 1.87 * 1e9;
  34
  35 /*
  36  * End of Configuration Section
  37  */
  38
  39 /*
  40  * General portability note:
  41  * MSVC mis-compiles explicit function template instantiations.
  42  * For example, f<A>() and f<B>() are both compiled to call f<A>().
  43  * BCC is unable to implicitly convert a "const char *" to a std::string
  44  * when using explicit function template instantiations.
  45  *
  46  * Therefore, avoid explicit function template instantiations.
  47  */
  48
  49 // provides a run-time configurable linear congruential generator, just
  50 // for comparison
  51 template<class IntType>
  52 class linear_congruential
  53 {
  54 public:
  55   typedef IntType result_type;
  56
  57   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  58
  59   linear_congruential(IntType x0, IntType a, IntType c, IntType m)
  60     : _x(x0), _a(a), _c(c), _m(m) { }
  61   // compiler-generated copy ctor and assignment operator are fine
  62   void seed(IntType x0, IntType a, IntType c, IntType m)
  63     { _x = x0; _a = a; _c = c; _m = m; }
  64   void seed(IntType x0) { _x = x0; }
  65   result_type operator()() { _x = (_a*_x+_c) % _m; return _x; }
  66   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _c == 0 ? 1 : 0; }
  67   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return _m -1; }
  68
  69 private:
  70   IntType _x, _a, _c, _m;
  71 };
  72
  73
  74 // simplest "random" number generator possible, to check on overhead
  75 class counting
  76 {
  77 public:
  78   typedef int result_type;
  79
  80   BOOST_STATIC_CONSTANT(bool, has_fixed_range = false);
  81
  82   counting() : _x(0) { }
  83   result_type operator()() { return ++_x; }
  84   result_type min BOOST_PREVENT_MACRO_SUBSTITUTION () const { return 1; }
  85   result_type max BOOST_PREVENT_MACRO_SUBSTITUTION () const { return (std::numeric_limits<result_type>::max)(); }
  86
  87 private:
  88   int _x;
  89 };
  90
  91
  92 // decoration of variate_generator to make it runtime-exchangeable
  93 // for speed comparison
  94 template<class Ret>
  95 class RandomGenBase
  96 {
  97 public:
  98   virtual Ret operator()() = 0;
  99   virtual ~RandomGenBase() { }
 100 };
 101
 102 template<class URNG, class Dist, class Ret = typename Dist::result_type>
 103 class DynamicRandomGenerator
 104   : public RandomGenBase<Ret>
 105 {
 106 public:
 107   DynamicRandomGenerator(URNG& urng, const Dist& d) : _rng(urng, d) { }
 108   Ret operator()() { return _rng(); }
 109 private:
 110   boost::variate_generator<URNG&, Dist> _rng;
 111 };
 112
 113 template<class Ret>
 114 class GenericRandomGenerator
 115 {
 116 public:
 117   typedef Ret result_type;
 118
 119   GenericRandomGenerator() { };
 120   void set(boost::shared_ptr<RandomGenBase<Ret> > p) { _p = p; }
 121   // takes over ownership
 122   void set(RandomGenBase<Ret> * p) { _p.reset(p); }
 123   Ret operator()() { return (*_p)(); }
 124 private:
 125   boost::shared_ptr<RandomGenBase<Ret> > _p;
 126 };
 127
 128
 129 // start implementation of measuring timing
 130
 131 void show_elapsed(double end, int iter, const std::string & name)
 132 {
 133   double usec = end/iter*1e6;
 134   double cycles = usec * cpu_frequency/1e6;
 135   std::cout << name << ": "
 136             << usec*1e3 << " nsec/loop = "
 137             << cycles << " CPU cycles"
 138             << std::endl;
 139 }
 140
 141 #if 0
 142 template<class RNG>
 143 void timing(RNG & rng, int iter, const std::string& name)
 144 {
 145   // make sure we're not optimizing too much
 146   volatile typename RNG::result_type tmp;
 147   boost::timer t;
 148   for(int i = 0; i < iter; i++)
 149     tmp = rng();
 150   show_elapsed(t.elapsed(), iter, name);
 151 }
 152 #endif
 153
 154 // overload for using a copy, allows more concise invocation
 155 template<class RNG>
 156 void timing(RNG rng, int iter, const std::string& name)
 157 {
 158   // make sure we're not optimizing too much
 159   volatile typename RNG::result_type tmp;
 160   boost::timer t;
 161   for(int i = 0; i < iter; i++)
 162     tmp = rng();
 163   show_elapsed(t.elapsed(), iter, name);
 164 }
 165
 166 template<class RNG>
 167 void timing_sphere(RNG rng, int iter, const std::string & name)
 168 {
 169   boost::timer t;
 170   for(int i = 0; i < iter; i++) {
 171     // the special return value convention of uniform_on_sphere saves 20% CPU
 172     const std::vector<double> & tmp = rng();
 173     (void) tmp[0];
 174   }
 175   show_elapsed(t.elapsed(), iter, name);
 176 }
 177
 178 template<class RNG>
 179 void run(int iter, const std::string & name, RNG rng)
 180 {
 181   std::cout << (RNG::has_fixed_range ? "fixed-range " : "");
 182   // BCC has trouble with string autoconversion for explicit specializations
 183
 184   // make sure we're not optimizing too much
 185   volatile typename RNG::result_type tmp;
 186   boost::timer t;
 187   for(int i = 0; i < iter; i++)
 188     tmp = rng();
 189   show_elapsed(t.elapsed(), iter, name);
 190 }
 191
 192 #ifdef HAVE_DRAND48
 193 // requires non-standard C library support for srand48/lrand48
 194 struct lrand48_ {
 195     static const bool has_fixed_range = false;
 196     typedef long result_type;
 197     lrand48_() {
 198         using namespace std;
 199         srand48(1);
 200     }
 201     result_type operator()() {
 202         using namespace std;
 203         return lrand48();
 204     }
 205 };
 206 #endif
 207
 208 #ifdef HAVE_MT19937INT_C  // requires the original mt19937int.c
 209 extern "C" void sgenrand(unsigned long);
 210 extern "C" unsigned long genrand();
 211
 212 void run(int iter, const std::string & name, float)
 213 {
 214   sgenrand(4357);
 215   timing(genrand, iter, name, 0u);
 216 }
 217 #endif
 218
 219 #ifdef HAVE_MT19937AR_C
 220 extern "C" {
 221 void init_genrand(unsigned long s);
 222 unsigned long genrand_int32(void);
 223 }
 224 struct mt19937_c {
 225     static const bool has_fixed_range = false;
 226     mt19937_c() {
 227         init_genrand(5489);
 228     }
 229     typedef unsigned long result_type;
 230     result_type operator()() {
 231         return genrand_int32();
 232     }
 233 };
 234 #endif
 235
 236 template<class PRNG, class Dist>
 237 inline boost::variate_generator<PRNG&, Dist> make_gen(PRNG & rng, Dist d)
 238 {
 239   return boost::variate_generator<PRNG&, Dist>(rng, d);
 240 }
 241
 242 template<class Gen>
 243 void distrib(int iter, const std::string & name, const Gen &)
 244 {
 245   Gen gen;
 246
 247   timing(make_gen(gen, boost::random::uniform_int_distribution<>(-2, 4)),
 248          iter, name + " uniform_int");
 249
 250   timing(make_gen(gen, boost::random::uniform_smallint<>(-2, 4)),
 251          iter, name + " uniform_smallint");
 252
 253   timing(make_gen(gen, boost::random::bernoulli_distribution<>(0.5)),
 254          iter, name + " bernoulli");
 255
 256   timing(make_gen(gen, boost::random::geometric_distribution<>(0.5)),
 257          iter, name + " geometric");
 258
 259   timing(make_gen(gen, boost::random::binomial_distribution<int>(4, 0.8)),
 260          iter, name + " binomial");
 261
 262   timing(make_gen(gen, boost::random::negative_binomial_distribution<int>(4, 0.8)),
 263          iter, name + " negative_binomial");
 264
 265   timing(make_gen(gen, boost::random::poisson_distribution<>(1)),
 266          iter, name + " poisson");
 267
 268
 269   timing(make_gen(gen, boost::random::uniform_real_distribution<>(-5.3, 4.8)),
 270          iter, name + " uniform_real");
 271
 272   timing(make_gen(gen, boost::random::uniform_01<>()),
 273          iter, name + " uniform_01");
 274
 275   timing(make_gen(gen, boost::random::triangle_distribution<>(1, 2, 7)),
 276          iter, name + " triangle");
 277
 278   timing(make_gen(gen, boost::random::exponential_distribution<>(3)),
 279          iter, name + " exponential");
 280
 281   timing(make_gen(gen, boost::random::normal_distribution<>()),
 282                   iter, name + " normal polar");
 283
 284   timing(make_gen(gen, boost::random::lognormal_distribution<>()),
 285          iter, name + " lognormal");
 286
 287   timing(make_gen(gen, boost::random::chi_squared_distribution<>(4)),
 288          iter, name + " chi squared");
 289
 290   timing(make_gen(gen, boost::random::cauchy_distribution<>()),
 291          iter, name + " cauchy");
 292
 293   timing(make_gen(gen, boost::random::fisher_f_distribution<>(4, 5)),
 294          iter, name + " fisher f");
 295
 296   timing(make_gen(gen, boost::random::student_t_distribution<>(7)),
 297          iter, name + " student t");
 298
 299   timing(make_gen(gen, boost::random::gamma_distribution<>(2.8)),
 300          iter, name + " gamma");
 301
 302   timing(make_gen(gen, boost::random::weibull_distribution<>(3)),
 303          iter, name + " weibull");
 304
 305   timing(make_gen(gen, boost::random::extreme_value_distribution<>()),
 306          iter, name + " extreme value");
 307
 308   timing_sphere(make_gen(gen, boost::random::uniform_on_sphere<>(3)),
 309                 iter/10, name + " uniform_on_sphere");
 310 }
 311
 312 int main(int argc, char*argv[])
 313 {
 314   if(argc != 2) {
 315     std::cerr << "usage: " << argv[0] << " iterations" << std::endl;
 316     return 1;
 317   }
 318
 319   // okay, it's ugly, but it's only used here
 320   int iter =
 321 #ifndef BOOST_NO_STDC_NAMESPACE
 322     std::
 323 #endif
 324     atoi(argv[1]);
 325
 326 #if !defined(BOOST_NO_INT64_T) && \
 327     !defined(BOOST_NO_INCLASS_MEMBER_INITIALIZATION)
 328   run(iter, "rand48", boost::rand48());
 329   linear_congruential<boost::uint64_t>
 330     lcg48(boost::uint64_t(1)<<16 | 0x330e,
 331           boost::uint64_t(0xDEECE66DUL) | (boost::uint64_t(0x5) << 32), 0xB,
 332           boost::uint64_t(1)<<48);
 333   timing(lcg48, iter, "lrand48 run-time");
 334 #endif
 335
 336 #ifdef HAVE_DRAND48
 337   // requires non-standard C library support for srand48/lrand48
 338   run(iter, "lrand48", lrand48_());   // coded for lrand48()
 339 #endif
 340
 341   run(iter, "minstd_rand0", boost::minstd_rand0());
 342   run(iter, "minstd_rand", boost::minstd_rand());
 343   run(iter, "ecuyer combined", boost::ecuyer1988());
 344   run(iter, "kreutzer1986", boost::kreutzer1986());
 345   run(iter, "taus88", boost::taus88());
 346   run(iter, "knuth_b", boost::random::knuth_b());
 347
 348   run(iter, "hellekalek1995 (inversive)", boost::hellekalek1995());
 349
 350   run(iter, "mt11213b", boost::mt11213b());
 351   run(iter, "mt19937", boost::mt19937());
 352 #if !defined(BOOST_NO_INT64_T)
 353   run(iter, "mt19937_64", boost::mt19937_64());
 354 #endif
 355
 356   run(iter, "lagged_fibonacci607", boost::lagged_fibonacci607());
 357   run(iter, "lagged_fibonacci1279", boost::lagged_fibonacci1279());
 358   run(iter, "lagged_fibonacci2281", boost::lagged_fibonacci2281());
 359   run(iter, "lagged_fibonacci3217", boost::lagged_fibonacci3217());
 360   run(iter, "lagged_fibonacci4423", boost::lagged_fibonacci4423());
 361   run(iter, "lagged_fibonacci9689", boost::lagged_fibonacci9689());
 362   run(iter, "lagged_fibonacci19937", boost::lagged_fibonacci19937());
 363   run(iter, "lagged_fibonacci23209", boost::lagged_fibonacci23209());
 364   run(iter, "lagged_fibonacci44497", boost::lagged_fibonacci44497());
 365
 366   run(iter, "subtract_with_carry", boost::random::ranlux_base());
 367   run(iter, "subtract_with_carry_01", boost::random::ranlux_base_01());
 368   run(iter, "ranlux3", boost::ranlux3());
 369   run(iter, "ranlux4", boost::ranlux4());
 370   run(iter, "ranlux3_01", boost::ranlux3_01());
 371   run(iter, "ranlux4_01", boost::ranlux4_01());
 372   run(iter, "ranlux64_3", boost::ranlux3());
 373   run(iter, "ranlux64_4", boost::ranlux4());
 374   run(iter, "ranlux64_3_01", boost::ranlux3_01());
 375   run(iter, "ranlux64_4_01", boost::ranlux4_01());
 376   run(iter, "ranlux24", boost::ranlux3());
 377   run(iter, "ranlux48", boost::ranlux4());
 378
 379   run(iter, "counting", counting());
 380
 381 #ifdef HAVE_MT19937INT_C
 382   // requires the original mt19937int.c
 383   run<float>(iter, "mt19937 original");   // coded for sgenrand()/genrand()
 384 #endif
 385
 386 #ifdef HAVE_MT19937AR_C
 387   run(iter, "mt19937ar.c", mt19937_c());
 388 #endif
 389
 390   distrib(iter, "counting", counting());
 391
 392   distrib(iter, "minstd_rand", boost::minstd_rand());
 393
 394   distrib(iter, "kreutzer1986", boost::kreutzer1986());
 395
 396   distrib(iter, "mt19937", boost::mt19937());
 397
 398   distrib(iter, "lagged_fibonacci607", boost::lagged_fibonacci607());
 399 }