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