[ExRandom lib-examples]

simple_normal

#include <iostream>
#include <random>
#include <exrandom/unit_normal_distribution.hpp>
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  exrandom::unit_normal_distribution<double> N;
  std::cout
    << "Sampling doubles exactly from the unit normal distribution:\n";
  for (int i = 0; i < 48; ++i)
    std::cout << N(g) << ((i + 1) % 8 ? ' ' : '\n');
  std::cout << "\n";
}

simple_uniform

从[0,1]范围 按照均匀分布采样

#include <iostream>
#include <random>
#include <exrandom/unit_uniform_distribution.hpp>
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  exrandom::unit_uniform_distribution<long double> N;
  std::cout
    << "Sampling long doubles exactly from the unit uniform distribution:\n";
  for (int i = 0; i < 48; ++i)
    std::cout << N(g) << ((i + 1) % 8 ? ' ' : '\n');
  std::cout << "\n";
}

simple_exponential

#include <iostream>
#include <random>
#include <exrandom/unit_exponential_distribution.hpp>
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  exrandom::unit_exponential_distribution<float> N;
  std::cout
    << "Sampling floats exactly from the unit exponential distribution:\n";
  for (int i = 0; i < 48; ++i)
    std::cout << N(g) << ((i + 1) % 8 ? ' ' : '\n');
  std::cout << "\n";
}

simple_discrete_normal

#include <iostream>
#include <random>
#include <exrandom/discrete_normal_distribution.hpp>
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  exrandom::discrete_normal_distribution N(1,3,129,2);
  std::cout
    << "Sampling exactly from the discrete normal distribution,\n"
    << "  mu = " << N.mu_num() << "/" << N.mu_den()
    << ", sigma = " << N.sigma_num() << "/" << N.sigma_den() << ":\n";
  for (int i = 0; i < 100; ++i)
    std::cout << N(g) << ((i + 1) % 20 ? ' ' : '\n');
  std::cout << "\n";
}

sample_uniform

#include <iostream>
#include <iomanip>
#include <random>
#include <exrandom/rand_digit.hpp>
#include <exrandom/unit_uniform_dist.hpp>
int main() {
  std::cout << "Sampling from the uniform distribution\n";
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  {
    const unsigned b = 16;
    std::cout << "Deviates rounded to floating point with base = " << b << "\n"
              << "columns are: initial u-rand, rounded float, final u-rand\n"
              << std::setprecision(std::numeric_limits<float>::digits10);
    exrandom::rand_digit<b> D;
    exrandom::unit_uniform_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
    for (int i = 0; i < 10; ++i) {
      N.generate(g, x);
      std::cout << x << " = ";
      std::cout << x.value<float>(g) << " = ";
      std::cout << x << "\n";
    }
    std::cout << "Total number of base " << b << " digits used = "
              << D.count() << "\n\n";
  }
  {
    const unsigned b = 10;
    int p = 6;
    std::cout << "Deviates rounded to fixed point with base = " << b
              << " and precision = " << p << "\n"
              << "columns are: initial u-rand, rounded fixed, final u-rand\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_uniform_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
    for (int i = 0; i < 10; ++i) {
      N.generate(g, x);
      std::cout << x << " = ";
      std::cout << x.print_fixed(g,p) << " = ";
      std::cout << x << "\n";
    }
    std::cout << "Total number of base " << b << " digits used = "
              << D.count() << "\n\n";
  }
}

sample_normal

#include <iostream>
#include <iomanip>
#include <random>
#include <exrandom/rand_digit.hpp>
#include <exrandom/unit_normal_dist.hpp>
int main() {
  std::cout << "Sampling from the normal distribution\n";
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  {
    const unsigned b = 16;
    std::cout << "Deviates rounded to floating point with base = " << b << "\n"
              << "columns are: initial u-rand, rounded float, final u-rand\n"
              << std::setprecision(std::numeric_limits<float>::digits10);
    exrandom::rand_digit<b> D;
    exrandom::unit_normal_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
    for (int i = 0; i < 10; ++i) {
      N.generate(g, x);
      std::cout << x << " = ";
      std::cout << x.value<float>(g) << " = ";
      std::cout << x << "\n";
    }
    std::cout << "Total number of base " << b << " digits used = "
              << D.count() << "\n\n";
  }
  {
    const unsigned b = 10;
    int p = 6;
    std::cout << "Deviates rounded to fixed point with base = " << b
              << " and precision = " << p << "\n"
              << "columns are: initial u-rand, rounded fixed, final u-rand\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_normal_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
    for (int i = 0; i < 10; ++i) {
      N.generate(g, x);
      std::cout << x << " = ";
      std::cout << x.print_fixed(g,p) << " = ";
      std::cout << x << "\n";
    }
    std::cout << "Total number of base " << b << " digits used = "
              << D.count() << "\n\n";
  }
}

sample_exponential

#include <iostream>
#include <iomanip>
#include <random>
#include <exrandom/rand_digit.hpp>
#include <exrandom/unit_exponential_dist.hpp>
int main() {
  std::cout << "Sampling from the exponential distribution\n";
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  {
    const unsigned b = 16;
    std::cout << "Deviates rounded to floating point with base = " << b << "\n"
              << "columns are: initial u-rand, rounded float, final u-rand\n"
              << std::setprecision(std::numeric_limits<float>::digits10);
    exrandom::rand_digit<b> D;
    exrandom::unit_exponential_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
    for (int i = 0; i < 10; ++i) {
      N.generate(g, x);
      std::cout << x << " = ";
      std::cout << x.value<float>(g) << " = ";
      std::cout << x << "\n";
    }
    std::cout << "Total number of base " << b << " digits used = "
              << D.count() << "\n\n";
  }
  {
    const unsigned b = 10;
    int p = 6;
    std::cout << "Deviates rounded to fixed point with base = " << b
              << " and precision = " << p << "\n"
              << "columns are: initial u-rand, rounded fixed, final u-rand\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_exponential_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
    for (int i = 0; i < 10; ++i) {
      N.generate(g, x);
      std::cout << x << " = ";
      std::cout << x.print_fixed(g,p) << " = ";
      std::cout << x << "\n";
    }
    std::cout << "Total number of base " << b << " digits used = "
              << D.count() << "\n\n";
  }
}

sample_distributions

#include <iostream>
#include <random>
#include <map>
#include <cmath>
#include <algorithm>
#include <string>
#include <exrandom/unit_normal_distribution.hpp>
#include <exrandom/unit_exponential_distribution.hpp>
#include <exrandom/unit_uniform_distribution.hpp>
#include <exrandom/discrete_normal_distribution.hpp>
template<typename Dist, typename Generator>
void histogram(long long num, Dist& d, Generator& g,
               double step, double offset) {
  std::map<int,unsigned> hist;
  const unsigned maxlength = 40; // Widest bar in histogram
  for (long long i = 0; i < num; ++i)
    ++hist[int(std::floor((double(d(g)) - offset) / step + 0.5))];
  unsigned m = 0;
  for (auto p : hist) m = (std::max)(m, p.second);
  m = (std::max)(m / maxlength, 1U);
  for (auto p : hist) {
    unsigned k = (p.second + m/2) / m;
    if (k)
      std::cout << typename Dist::result_type(p.first * step + offset) << '\t'
                << std::string(k, '*') << '\n';
  }
  std::cout << std::endl;
}
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  long long num = 1000000LL;
  {
    exrandom::unit_normal_distribution<double> d;
    std::cout << "Sampling the unit normal distribution\n";
    histogram(num, d, g, 0.2, 0.0);
  }
  {
    exrandom::unit_exponential_distribution<double> d;
    std::cout << "Sampling the unit exponential distribution\n";
    histogram(num, d, g, 0.2, 0.1);
  }
  {
    exrandom::unit_uniform_distribution<double> d;
    std::cout << "Sampling the unit uniform distribution\n";
    histogram(num, d, g, 0.1, 0.05);
  }
  {
    exrandom::discrete_normal_distribution d(6, 55, 10);
    std::cout << "Sampling the discrete normal distribution"
              << " mu = " << d.mu_num() << "/" << d.mu_den()
              << " sigma = " << d.sigma_num() << "/" << d.sigma_den() << "\n";
    histogram(num, d, g, 1, 0);
  }
}

sample_discrete_normal

#include <iostream>
#include <random>
#include <map>
#include <exrandom/rand_digit.hpp>
#include <exrandom/discrete_normal_dist.hpp>
int main() {
  std::cout << "Sampling from the discrete normal distribution\n";
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  const unsigned b = 2;
  exrandom::rand_digit<b> D;
  exrandom::discrete_normal_dist<exrandom::rand_digit<b>>
    N(D, 87331, 3, 16000051, 100);
  exrandom::i_rand<exrandom::rand_digit<b>> h(N.digit_generator());
  std::cout << "Parameters: "
            << " mu = " << N.mu_num() << "/" << N.mu_den()
            << " sigma = " << N.sigma_num() << "/" << N.sigma_den() << "\n"
            << "columns are: uniform range, final int\n";
  int num = 20;
  for (int i = 0; i < num; ++i) {
    N.generate(g, h);
    std::cout << h << " = ";
    std::cout << h(g) << "\n";
  }
    std::cout << "Total number of base " << b << " digits used = "
              << D.count() << "\n\n";
}

tabulate_normals

#include <iostream>
#include <string>
#include <vector>
#include <fstream>
#include <exrandom/table_gen.hpp>
#include <exrandom/rand_table.hpp>
#include <exrandom/unit_normal_dist.hpp>
int main() {
  std::vector<std::string> seed_table;
  // digits.txt is the the RAND table of random digits available from
  // http://www.rand.org/content/dam/rand/pubs/monograph_reports/MR1418/MR1418.digits.txt.zip
  std::ifstream seed_file("digits.txt");
  if (seed_file.good()) {
    std::string s;
    while (std::getline(seed_file, s)) {
      std::istringstream is(s);
      is >> s;
      std::string seed;
      for (int i = 0; i < 10; ++i) {
        is >> s;
        seed += s;
      }
      seed_table.push_back(seed);
    }
  } else {
    // lines 9077 - 9086, in case the full table is not available
    seed_table.push_back("91486866851718646640939841418208696361599740927130");
    seed_table.push_back("27085545979733960919843574531542387910433856492503");
    seed_table.push_back("50144629743871836759770203434855225597834197811609");
    seed_table.push_back("06513031977786096289011670052044879745609312522650");
    seed_table.push_back("27360659086699793904300006309496840459265099570954");
    seed_table.push_back("02069431357777418334914891893312167990996169079700");
    seed_table.push_back("32331134586674046873033847880723099262580390063311");
    seed_table.push_back("94820103089730893927855765220910774260125296553336");
    seed_table.push_back("13206851466217116627804451376819961678521211905382");
    seed_table.push_back("73582504716450282034788914850174375938706490445325");
  }
  exrandom::rand_table D;
  exrandom::u_rand<exrandom::rand_table> x(D);
  exrandom::unit_normal_dist<exrandom::rand_table> U(D);
  exrandom::table_gen g;
  for (size_t k = 0; k < seed_table.size(); ++k) {
    try {
      g.seed(seed_table[k]);
      long long c0 = D.count();
      U.generate(g,x);
      long long c1 = D.count();
      std::string x0 = x.print();
      std::string x1 = x.print_fixed(g, 6);
      std::string x2 = x.print();
      long long c2 = D.count();
      std::cout << seed_table[k].substr(0,int(c1-c0)) << "|"
                << seed_table[k].substr(int(c1-c0), int(c2-c1)) << "... -> ";
      std::cout << x0 << " = " << x2 << " ~ " << x1 << "\n";
    }
    catch (const std::exception&) {
      std::cerr << seed_table[k] << "... overflow\n";
    }
  }
}

multiprec_test

#include <iostream>
#include <random>
#if !defined(EXRANDOM_HAVE_MPREAL)
#define EXRANDOM_HAVE_MPREAL 0
#endif
#if !defined(EXRANDOM_HAVE_BOOST_DEC)
#define EXRANDOM_HAVE_BOOST_DEC 0
#endif
#if !defined(EXRANDOM_HAVE_BOOST_GMP)
#define EXRANDOM_HAVE_BOOST_GMP 0
#endif
#if !defined(EXRANDOM_HAVE_BOOST_MPFR)
#define EXRANDOM_HAVE_BOOST_MPFR 0
#endif
#if !defined(EXRANDOM_HAVE_BOOST_FLOAT128)
#define EXRANDOM_HAVE_BOOST_FLOAT128 0
#endif
#if EXRANDOM_HAVE_MPREAL
#include <mpreal.h>
#endif
#if EXRANDOM_HAVE_BOOST_MPFR
#include <boost/multiprecision/mpfr.hpp>
#endif
#if EXRANDOM_HAVE_BOOST_GMP
#include <boost/multiprecision/gmp.hpp>
#endif
#if EXRANDOM_HAVE_BOOST_DEC
#include <boost/multiprecision/cpp_dec_float.hpp>
#endif
#if EXRANDOM_HAVE_BOOST_FLOAT128
#include <boost/multiprecision/float128.hpp>
#endif
// N.B. This #include should come *after* the headers which define the
// floating point types.
#include <exrandom/unit_normal_distribution.hpp>
int main() {
#if EXRANDOM_HAVE_MPREAL || EXRANDOM_HAVE_BOOST_MPFR || \
  EXRANDOM_HAVE_BOOST_GMP || EXRANDOM_HAVE_BOOST_DEC || \
  EXRANDOM_HAVE_BOOST_FLOAT128
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
#endif
#if EXRANDOM_HAVE_MPREAL
  {
    typedef mpfr::mpreal real;
    real::set_default_prec(240);
    exrandom::unit_normal_distribution<real> N;
    std::cout
      << "Sampling mpreal from the unit normal distribution:\n"
      << std::setprecision(mpfr::bits2digits(real::get_default_prec()));
    for (int i = 0; i < 20; ++i)
      std::cout << N(g) << "\n";
    std::cout << "\n";
  }
#else
  std::cerr << "This example requires mpreal.h\n\n";
#endif
#if EXRANDOM_HAVE_BOOST_MPFR
  {
    typedef boost::multiprecision::mpfr_float_50 real;
    exrandom::unit_normal_distribution<real> N;
    std::cout
      << "Sampling mpfr_float_50 from the unit normal distribution:\n"
      << std::setprecision(50);
    for (int i = 0; i < 20; ++i)
      std::cout << N(g) << "\n";
    std::cout << "\n";
  }
#else
  std::cerr << "This example requires boost/multiprecision/mpfr.hpp\n\n";
#endif
#if EXRANDOM_HAVE_BOOST_GMP
  {
    typedef boost::multiprecision::mpf_float_50 real;
    exrandom::unit_normal_distribution<real> N;
    std::cout
      << "Sampling mpf_float_50 from the unit normal distribution:\n"
      << std::setprecision(50);
    for (int i = 0; i < 20; ++i)
      std::cout << N(g) << "\n";
    std::cout << "\n";
  }
#else
  std::cerr << "This example requires boost/multiprecision/gmp.hpp\n\n";
#endif
#if EXRANDOM_HAVE_BOOST_DEC
  {
    typedef
      boost::multiprecision::number<boost::multiprecision::cpp_dec_float<20U>>
      real;
    exrandom::unit_normal_distribution<real> N;
    std::cout <<
      "Sampling number<cpp_dec_float<20U>> from the unit normal distribution:\n"
      << std::setprecision(20);
    for (int i = 0; i < 20; ++i)
      std::cout << N(g) << "\n";
    std::cout << "\n";
  }
#else
  std::cerr
    << "This example requires boost/multiprecision/cpp_dec_float.hpp\n\n";
#endif
#if EXRANDOM_HAVE_BOOST_FLOAT128
  {
    typedef boost::multiprecision::float128 real;
    exrandom::unit_normal_distribution<real> N;
    std::cout
      << "Sampling float128 exactly from the unit normal distribution:\n"
      << std::setprecision(33);
    for (int i = 0; i < 20; ++i)
      std::cout << N(g) << "\n";
    std::cout << "\n";
  }
#else
  std::cerr << "This example requires boost/multiprecision/float128.hpp\n\n";
#endif
}

hist_data

#include <iostream>
#include <iomanip>
#include <random>
#include <map>
#include <exrandom/rand_digit.hpp>
#include <exrandom/u_rand.hpp>
#include <exrandom/unit_normal_dist.hpp>
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cerr << "Seed set to " << s << "\n"
            << "columns are l n prob(n/2^l)\n\n";
  const unsigned b = 2;
  exrandom::rand_digit<b> D;
  exrandom::unit_normal_dist<exrandom::rand_digit<b>> N(D);
  exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
  long long num = 5000000LL;
  unsigned maxl = 8, maxn = 3;
  std::map<std::pair<long long, long long>, long long> hist;
  for (long long i = 0; i < num; ++i) {
    N.generate(g, x);
    if (!(x.ndigits() <= maxl && x.integer() < maxn))
      continue;
    if (x.sign() < 1) x.negate();
    ++hist[x.rational()];
  }
  std::cout << std::setprecision(10);
  for (long l = 0, d = 1; l <= long(maxl); ++l, d *= 2)
    for (long n = 0; n < long(maxn)*d; ++n)
      std::cout << l << " " << n << " "
                << hist[std::pair<long long,long long>(n,d)]/double(num)
                << "\n";
}

Seed set to 2659261496
columns are l n prob(n/2^l)

0 0 0.398886
0 1 0
0 2 0
1 0 0.0995716
1 1 0.0028946
1 2 0.0907916
1 3 0.013058
1 4 0.0101286
1 5 0.0003244
2 0 0.0250126
2 1 0.0125894
2 2 0.0374668
2 3 0.0054856
2 4 0.0339686
2 5 0.0066612
2 6 0.023425
2 7 0.0078724
2 8 0.0082116
2 9 0.0010498
2 10 0.0022626
2 11 0.0005622
3 0 0.0062726
3 1 0.0046628
3 2 0.007737
3 3 0.0032878
3 4 0.0093202
3 5 0.002567
3 6 0.011083
3 7 0.0028742
3 8 0.009916
3 9 0.0025964
3 10 0.0092066
3 11 0.002766
3 12 0.0080532
3 13 0.0028834
3 14 0.0063594
3 15 0.0027322
3 16 0.003031
3 17 0.0004968
3 18 0.002043
3 19 0.000556
3 20 0.001284
3 21 0.0004462
3 22 0.0007146
3 23 0.0002926
4 0 0.0015626
4 1 0.001343
4 2 0.0017742
4 3 0.0011566
4 4 0.0019266
4 5 0.001005
4 6 0.0021374
4 7 0.000878
4 8 0.002337
4 9 0.0008076
4 10 0.0025204
4 11 0.0007606
4 12 0.002736
4 13 0.0008342
4 14 0.0029704
4 15 0.0009782
4 16 0.0027164
4 17 0.0007852
4 18 0.0026088
4 19 0.000789
4 20 0.0025176
4 21 0.0008134
4 22 0.002436
4 23 0.000829
4 24 0.0022638
4 25 0.0008578
4 26 0.0021186
4 27 0.0008554
4 28 0.0018986
4 29 0.0008566
4 30 0.001706
4 31 0.0008322
4 32 0.0009082
4 33 0.0001574
4 34 0.0007646
4 35 0.0001898
4 36 0.000643
4 37 0.0001792
4 38 0.0005344
4 39 0.0001746
4 40 0.000441
4 41 0.0001678
4 42 0.0003568
4 43 0.0001572
4 44 0.0002692
4 45 0.000119
4 46 0.0002048
4 47 9.88e-05
5 0 0.0003854
5 1 0.0003664
5 2 0.0004178
5 3 0.0003428
5 4 0.0004458
5 5 0.000311
5 6 0.0004526
5 7 0.0002892
5 8 0.0005004
5 9 0.0002944
5 10 0.0005306
5 11 0.0002478
5 12 0.0005344
5 13 0.000251
5 14 0.0005522
5 15 0.0002252
5 16 0.0005806
5 17 0.0002168
5 18 0.0006126
5 19 0.0002154
5 20 0.0006336
5 21 0.0002118
5 22 0.0006498
5 23 0.0002134
5 24 0.000698
5 25 0.0002324
5 26 0.0007056
5 27 0.000246
5 28 0.0007448
5 29 0.0002612
5 30 0.0007674
5 31 0.000278
5 32 0.000667
5 33 0.0002182
5 34 0.000684
5 35 0.00022
5 36 0.0006806
5 37 0.0002102
5 38 0.0006718
5 39 0.0002118
5 40 0.000639
5 41 0.0002174
5 42 0.000641
5 43 0.0002288
5 44 0.000625
5 45 0.0002278
5 46 0.0006266
5 47 0.0002372
5 48 0.0006008
5 49 0.0002328
5 50 0.0005752
5 51 0.0002412
5 52 0.000552
5 53 0.0002436
5 54 0.0005406
5 55 0.000233
5 56 0.0005074
5 57 0.00024
5 58 0.0004752
5 59 0.0002396
5 60 0.0004838
5 61 0.000216
5 62 0.0004584
5 63 0.0002322
5 64 0.0002296
5 65 4.94e-05
5 66 0.0002218
5 67 5.38e-05
5 68 0.0002056
5 69 5.06e-05
5 70 0.000209
5 71 5.66e-05
5 72 0.0001856
5 73 5.04e-05
5 74 0.0001632
5 75 4.62e-05
5 76 0.0001462
5 77 4.78e-05
5 78 0.0001416
5 79 4.82e-05
5 80 0.000125
5 81 4.62e-05
5 82 0.000115
5 83 4.42e-05
5 84 9.42e-05
5 85 4.74e-05
5 86 9.2e-05
5 87 3.72e-05
5 88 8.3e-05
5 89 4.3e-05
5 90 6.88e-05
5 91 3.98e-05
5 92 6.62e-05
5 93 4.02e-05
5 94 5.14e-05
5 95 2.96e-05
6 0 9.58e-05
6 1 9.14e-05
6 2 8.9e-05
6 3 8.68e-05
6 4 0.000101
6 5 9.1e-05
6 6 0.0001034
6 7 8.58e-05
6 8 0.000107
6 9 8.52e-05
6 10 0.0001272
6 11 8.06e-05
6 12 0.0001208
6 13 8.38e-05
6 14 0.0001262
6 15 7.34e-05
6 16 0.0001142
6 17 7.1e-05
6 18 0.0001272
6 19 6.32e-05
6 20 0.000131
6 21 6.86e-05
6 22 0.0001356
6 23 6.86e-05
6 24 0.0001392
6 25 6.52e-05
6 26 0.0001418
6 27 7.16e-05
6 28 0.0001424
6 29 5.78e-05
6 30 0.0001368
6 31 6.34e-05
6 32 0.0001466
6 33 5.98e-05
6 34 0.0001484
6 35 5.92e-05
6 36 0.000153
6 37 5.5e-05
6 38 0.0001504
6 39 5.76e-05
6 40 0.0001644
6 41 5.52e-05
6 42 0.0001682
6 43 5.22e-05
6 44 0.0001722
6 45 6.06e-05
6 46 0.0001656
6 47 5.7e-05
6 48 0.0001786
6 49 5.62e-05
6 50 0.0001736
6 51 5.92e-05
6 52 0.0001768
6 53 5.78e-05
6 54 0.0001746
6 55 6.38e-05
6 56 0.0001902
6 57 6.12e-05
6 58 0.000191
6 59 6.32e-05
6 60 0.000206
6 61 7.14e-05
6 62 0.0001982
6 63 7.56e-05
6 64 0.0001836
6 65 5.5e-05
6 66 0.0001696
6 67 5.4e-05
6 68 0.0001736
6 69 5.96e-05
6 70 0.000182
6 71 6.14e-05
6 72 0.0001786
6 73 5.1e-05
6 74 0.000157
6 75 5.12e-05
6 76 0.0001632
6 77 5.46e-05
6 78 0.0001668
6 79 5.48e-05
6 80 0.0001714
6 81 5.86e-05
6 82 0.00017
6 83 5.84e-05
6 84 0.0001594
6 85 6.38e-05
6 86 0.0001668
6 87 6.56e-05
6 88 0.0001586
6 89 6.36e-05
6 90 0.0001606
6 91 5.76e-05
6 92 0.0001566
6 93 5.6e-05
6 94 0.0001484
6 95 6.32e-05
6 96 0.0001644
6 97 5.66e-05
6 98 0.0001496
6 99 5.9e-05
6 100 0.000148
6 101 6.02e-05
6 102 0.0001458
6 103 6e-05
6 104 0.0001438
6 105 6.12e-05
6 106 0.0001456
6 107 6.5e-05
6 108 0.0001382
6 109 6.4e-05
6 110 0.0001322
6 111 6.36e-05
6 112 0.0001334
6 113 6.38e-05
6 114 0.0001366
6 115 5.92e-05
6 116 0.0001236
6 117 5.98e-05
6 118 0.0001248
6 119 5.98e-05
6 120 0.000128
6 121 5.72e-05
6 122 0.0001198
6 123 5.42e-05
6 124 0.0001158
6 125 6.12e-05
6 126 0.0001102
6 127 6.88e-05
6 128 5.98e-05
6 129 9.8e-06
6 130 6.44e-05
6 131 1.36e-05
6 132 5.3e-05
6 133 1.38e-05
6 134 5.88e-05
6 135 1.18e-05
6 136 5.28e-05
6 137 1.08e-05
6 138 5.52e-05
6 139 1.06e-05
6 140 4.96e-05
6 141 1.22e-05
6 142 4.52e-05
6 143 1.36e-05
6 144 4.76e-05
6 145 1.12e-05
6 146 4.26e-05
6 147 1.34e-05
6 148 4.66e-05
6 149 1.28e-05
6 150 4.02e-05
6 151 1.62e-05
6 152 4.2e-05
6 153 1.24e-05
6 154 3.68e-05
6 155 1.34e-05
6 156 3.56e-05
6 157 1.2e-05
6 158 3.82e-05
6 159 1.46e-05
6 160 3.62e-05
6 161 1.24e-05
6 162 3e-05
6 163 1.38e-05
6 164 2.98e-05
6 165 1.1e-05
6 166 3e-05
6 167 1.2e-05
6 168 3.2e-05
6 169 1.32e-05
6 170 2.8e-05
6 171 1.12e-05
6 172 2.54e-05
6 173 1.16e-05
6 174 2.04e-05
6 175 1.06e-05
6 176 2.22e-05
6 177 8.6e-06
6 178 1.62e-05
6 179 8e-06
6 180 1.9e-05
6 181 8.4e-06
6 182 2.14e-05
6 183 8.6e-06
6 184 1.64e-05
6 185 9.8e-06
6 186 1.46e-05
6 187 8.4e-06
6 188 1.94e-05
6 189 7.4e-06
6 190 1.58e-05
6 191 6.6e-06
7 0 2.16e-05
7 1 2.42e-05
7 2 2.12e-05
7 3 2.5e-05
7 4 2.4e-05
7 5 2.52e-05
7 6 2.96e-05
7 7 2.06e-05
7 8 2.72e-05
7 9 2.32e-05
7 10 2.66e-05
7 11 1.94e-05
7 12 3.26e-05
7 13 2.2e-05
7 14 2.54e-05
7 15 2.64e-05
7 16 2.86e-05
7 17 1.96e-05
7 18 3.26e-05
7 19 1.92e-05
7 20 3.04e-05
7 21 2.28e-05
7 22 2.8e-05
7 23 2.28e-05
7 24 2.98e-05
7 25 1.84e-05
7 26 3.58e-05
7 27 1.54e-05
7 28 2.86e-05
7 29 1.76e-05
7 30 2.96e-05
7 31 1.92e-05
7 32 2.56e-05
7 33 2.2e-05
7 34 2.94e-05
7 35 1.76e-05
7 36 2.74e-05
7 37 1.8e-05
7 38 3.6e-05
7 39 1.64e-05
7 40 2.94e-05
7 41 1.54e-05
7 42 3.38e-05
7 43 1.32e-05
7 44 3.72e-05
7 45 1.5e-05
7 46 3.12e-05
7 47 1.82e-05
7 48 3.02e-05
7 49 1.82e-05
7 50 3.5e-05
7 51 1.26e-05
7 52 3.18e-05
7 53 1.26e-05
7 54 3.64e-05
7 55 1.38e-05
7 56 3.58e-05
7 57 1.48e-05
7 58 2.96e-05
7 59 1.62e-05
7 60 3.68e-05
7 61 1.5e-05
7 62 3.5e-05
7 63 1.48e-05
7 64 3.34e-05
7 65 1.68e-05
7 66 3.58e-05
7 67 1.32e-05
7 68 4.12e-05
7 69 1.38e-05
7 70 3.56e-05
7 71 1.44e-05
7 72 3.98e-05
7 73 1.22e-05
7 74 3.94e-05
7 75 1.26e-05
7 76 3.62e-05
7 77 1.68e-05
7 78 4.02e-05
7 79 1.38e-05
7 80 3.98e-05
7 81 1.78e-05
7 82 3.84e-05
7 83 1.18e-05
7 84 4.22e-05
7 85 1.42e-05
7 86 3.82e-05
7 87 1.06e-05
7 88 4.62e-05
7 89 1.74e-05
7 90 3.86e-05
7 91 1.24e-05
7 92 4e-05
7 93 1.68e-05
7 94 4.6e-05
7 95 1.52e-05
7 96 4.56e-05
7 97 1.52e-05
7 98 4.48e-05
7 99 1.74e-05
7 100 4.34e-05
7 101 1.62e-05
7 102 4.36e-05
7 103 1.88e-05
7 104 4.4e-05
7 105 1.42e-05
7 106 4.58e-05
7 107 1.96e-05
7 108 4.56e-05
7 109 1.46e-05
7 110 4.94e-05
7 111 1.74e-05
7 112 4.58e-05
7 113 1.6e-05
7 114 4.92e-05
7 115 1.64e-05
7 116 4.48e-05
7 117 1.84e-05
7 118 5.18e-05
7 119 1.82e-05
7 120 4.6e-05
7 121 1.8e-05
7 122 5.3e-05
7 123 2.22e-05
7 124 5.22e-05
7 125 1.98e-05
7 126 4.78e-05
7 127 1.92e-05
7 128 4.62e-05
7 129 1.5e-05
7 130 3.84e-05
7 131 1.78e-05
7 132 4.3e-05
7 133 1.34e-05
7 134 4.62e-05
7 135 1e-05
7 136 4.84e-05
7 137 1.04e-05
7 138 4.5e-05
7 139 1.1e-05
7 140 4.18e-05
7 141 1.58e-05
7 142 4.56e-05
7 143 1.4e-05
7 144 3.94e-05
7 145 1.46e-05
7 146 4.28e-05
7 147 1.76e-05
7 148 4.64e-05
7 149 1.58e-05
7 150 4.24e-05
7 151 1.78e-05
7 152 4.74e-05
7 153 1.46e-05
7 154 4.52e-05
7 155 1.34e-05
7 156 4.12e-05
7 157 1.6e-05
7 158 4.9e-05
7 159 1.76e-05
7 160 4.14e-05
7 161 1.44e-05
7 162 4.7e-05
7 163 1.48e-05
7 164 4.36e-05
7 165 1.56e-05
7 166 4.64e-05
7 167 1.62e-05
7 168 4.22e-05
7 169 1.44e-05
7 170 4.16e-05
7 171 1.86e-05
7 172 3.66e-05
7 173 1.28e-05
7 174 4.4e-05
7 175 1.46e-05
7 176 4.38e-05
7 177 1.3e-05
7 178 4.12e-05
7 179 1.5e-05
7 180 4.38e-05
7 181 1.3e-05
7 182 4.04e-05
7 183 1.5e-05
7 184 4e-05
7 185 1.6e-05
7 186 4.16e-05
7 187 1.38e-05
7 188 3.76e-05
7 189 1.5e-05
7 190 4.26e-05
7 191 1.5e-05
7 192 4.02e-05
7 193 1.82e-05
7 194 3.3e-05
7 195 1.8e-05
7 196 3.42e-05
7 197 1.7e-05
7 198 3.98e-05
7 199 1.56e-05
7 200 3.52e-05
7 201 1.56e-05
7 202 3.66e-05
7 203 1.56e-05
7 204 3.48e-05
7 205 1.68e-05
7 206 3.88e-05
7 207 1.74e-05
7 208 3.76e-05
7 209 1.44e-05
7 210 3.7e-05
7 211 1.38e-05
7 212 3.5e-05
7 213 1.6e-05
7 214 3.86e-05
7 215 1.68e-05
7 216 3.22e-05
7 217 1.86e-05
7 218 3.14e-05
7 219 1.98e-05
7 220 3.88e-05
7 221 1.66e-05
7 222 3.38e-05
7 223 1.46e-05
7 224 3.32e-05
7 225 1.3e-05
7 226 3.26e-05
7 227 1.52e-05
7 228 3.54e-05
7 229 1.66e-05
7 230 2.98e-05
7 231 1.7e-05
7 232 3.46e-05
7 233 2.1e-05
7 234 3.32e-05
7 235 1.34e-05
7 236 3.02e-05
7 237 1.5e-05
7 238 3.58e-05
7 239 1.62e-05
7 240 3.12e-05
7 241 1.42e-05
7 242 2.86e-05
7 243 1.66e-05
7 244 2.98e-05
7 245 1.8e-05
7 246 2.96e-05
7 247 1.62e-05
7 248 3.02e-05
7 249 1.34e-05
7 250 2.96e-05
7 251 1.94e-05
7 252 2.94e-05
7 253 1.56e-05
7 254 2.52e-05
7 255 1.46e-05
7 256 1.64e-05
7 257 3.2e-06
7 258 1.66e-05
7 259 2.6e-06
7 260 1.28e-05
7 261 3e-06
7 262 1.44e-05
7 263 5.2e-06
7 264 1.34e-05
7 265 3e-06
7 266 1.2e-05
7 267 5e-06
7 268 1.32e-05
7 269 3.8e-06
7 270 1.44e-05
7 271 4.8e-06
7 272 1.24e-05
7 273 4.6e-06
7 274 1.22e-05
7 275 4.4e-06
7 276 1.38e-05
7 277 3.6e-06
7 278 1.46e-05
7 279 4.4e-06
7 280 1.18e-05
7 281 3.2e-06
7 282 1e-05
7 283 3.6e-06
7 284 1.32e-05
7 285 3.8e-06
7 286 1.3e-05
7 287 3.6e-06
7 288 1.36e-05
7 289 3.2e-06
7 290 1.3e-05
7 291 4.2e-06
7 292 1.16e-05
7 293 3.8e-06
7 294 1.04e-05
7 295 3.4e-06
7 296 1.18e-05
7 297 4e-06
7 298 9.6e-06
7 299 3.4e-06
7 300 1.26e-05
7 301 3.8e-06
7 302 1.08e-05
7 303 4.2e-06
7 304 1.02e-05
7 305 5.6e-06
7 306 1.3e-05
7 307 4.2e-06
7 308 1.1e-05
7 309 3.4e-06
7 310 1.08e-05
7 311 3.6e-06
7 312 8e-06
7 313 2.8e-06
7 314 1.04e-05
7 315 2.4e-06
7 316 8.2e-06
7 317 4e-06
7 318 7.6e-06
7 319 4.4e-06
7 320 9.2e-06
7 321 3.4e-06
7 322 6.4e-06
7 323 3.2e-06
7 324 5e-06
7 325 2.6e-06
7 326 7.4e-06
7 327 2.6e-06
7 328 7.8e-06
7 329 3e-06
7 330 8.8e-06
7 331 4e-06
7 332 8.2e-06
7 333 3.8e-06
7 334 7e-06
7 335 3.8e-06
7 336 7e-06
7 337 4.8e-06
7 338 5.6e-06
7 339 2.6e-06
7 340 5.4e-06
7 341 2.2e-06
7 342 5.6e-06
7 343 2.4e-06
7 344 5.8e-06
7 345 4.2e-06
7 346 5.6e-06
7 347 3.8e-06
7 348 6.4e-06
7 349 1.4e-06
7 350 5.6e-06
7 351 1.6e-06
7 352 4.8e-06
7 353 2.2e-06
7 354 4.4e-06
7 355 2.2e-06
7 356 6e-06
7 357 2.6e-06
7 358 5.4e-06
7 359 1.2e-06
7 360 4.8e-06
7 361 4.4e-06
7 362 7.4e-06
7 363 3e-06
7 364 4e-06
7 365 2.6e-06
7 366 3e-06
7 367 2.8e-06
7 368 3.8e-06
7 369 2.4e-06
7 370 2.2e-06
7 371 2.4e-06
7 372 3.8e-06
7 373 2.6e-06
7 374 3.8e-06
7 375 8e-07
7 376 2.8e-06
7 377 2e-06
7 378 4.2e-06
7 379 3e-06
7 380 3.8e-06
7 381 1e-06
7 382 3.4e-06
7 383 2.6e-06
8 0 5.4e-06
8 1 9.2e-06
8 2 5.6e-06
8 3 5.4e-06
8 4 7.2e-06
8 5 6e-06
8 6 7.2e-06
8 7 8.4e-06
8 8 6.6e-06
8 9 4e-06
8 10 6.6e-06
8 11 7.2e-06
8 12 5.6e-06
8 13 6.8e-06
8 14 5.6e-06
8 15 6.4e-06
8 16 7e-06
8 17 8.4e-06
8 18 5.8e-06
8 19 5.2e-06
8 20 7.4e-06
8 21 8.6e-06
8 22 7e-06
8 23 4.8e-06
8 24 5.8e-06
8 25 7.4e-06
8 26 8e-06
8 27 5.4e-06
8 28 6.4e-06
8 29 4.2e-06
8 30 8e-06
8 31 4e-06
8 32 6.4e-06
8 33 5.2e-06
8 34 7.4e-06
8 35 3.4e-06
8 36 5.4e-06
8 37 6.6e-06
8 38 5e-06
8 39 6e-06
8 40 7e-06
8 41 4.6e-06
8 42 5e-06
8 43 3.8e-06
8 44 5.6e-06
8 45 4e-06
8 46 7e-06
8 47 4e-06
8 48 6.4e-06
8 49 5e-06
8 50 8.2e-06
8 51 5.2e-06
8 52 8.4e-06
8 53 5.2e-06
8 54 7e-06
8 55 4.2e-06
8 56 6.6e-06
8 57 4.2e-06
8 58 4.4e-06
8 59 4.4e-06
8 60 6.8e-06
8 61 3.4e-06
8 62 6e-06
8 63 5.6e-06
8 64 8.4e-06
8 65 4.4e-06
8 66 8.2e-06
8 67 4.2e-06
8 68 5.4e-06
8 69 3.8e-06
8 70 7.4e-06
8 71 5.6e-06
8 72 8.2e-06
8 73 4.2e-06
8 74 8e-06
8 75 3.6e-06
8 76 9.6e-06
8 77 5e-06
8 78 5.6e-06
8 79 5.4e-06
8 80 9.2e-06
8 81 3.6e-06
8 82 8.2e-06
8 83 5.6e-06
8 84 8e-06
8 85 4.2e-06
8 86 9e-06
8 87 5.6e-06
8 88 7.8e-06
8 89 3.2e-06
8 90 7.6e-06
8 91 3.6e-06
8 92 1.04e-05
8 93 5.2e-06
8 94 8e-06
8 95 4.2e-06
8 96 9.8e-06
8 97 2.8e-06
8 98 7.8e-06
8 99 3.2e-06
8 100 8.2e-06
8 101 5e-06
8 102 7.4e-06
8 103 4e-06
8 104 9.2e-06
8 105 4.2e-06
8 106 1.02e-05
8 107 4.6e-06
8 108 8.4e-06
8 109 3.8e-06
8 110 7e-06
8 111 3.8e-06
8 112 9e-06
8 113 3.4e-06
8 114 9.2e-06
8 115 3.8e-06
8 116 9.6e-06
8 117 2.4e-06
8 118 1.1e-05
8 119 2.2e-06
8 120 9e-06
8 121 4.8e-06
8 122 7.8e-06
8 123 4.4e-06
8 124 6.6e-06
8 125 3.2e-06
8 126 1.08e-05
8 127 3.6e-06
8 128 8.8e-06
8 129 4e-06
8 130 8.2e-06
8 131 3.2e-06
8 132 1.02e-05
8 133 3.2e-06
8 134 8.8e-06
8 135 3.6e-06
8 136 8.4e-06
8 137 3e-06
8 138 7.2e-06
8 139 3.8e-06
8 140 8.4e-06
8 141 2.6e-06
8 142 1.1e-05
8 143 3.8e-06
8 144 9.2e-06
8 145 3.4e-06
8 146 8.4e-06
8 147 2.2e-06
8 148 1.06e-05
8 149 3e-06
8 150 8.2e-06
8 151 3.8e-06
8 152 1e-05
8 153 4e-06
8 154 1.16e-05
8 155 4e-06
8 156 1.18e-05
8 157 3.4e-06
8 158 1.22e-05
8 159 3.8e-06
8 160 9.4e-06
8 161 4.4e-06
8 162 1.16e-05
8 163 4e-06
8 164 8e-06
8 165 4.2e-06
8 166 9.6e-06
8 167 3.8e-06
8 168 1.12e-05
8 169 4.2e-06
8 170 1.04e-05
8 171 2.8e-06
8 172 8.6e-06
8 173 3.2e-06
8 174 1.18e-05
8 175 1.8e-06
8 176 8.8e-06
8 177 5.4e-06
8 178 9.8e-06
8 179 1.6e-06
8 180 1.04e-05
8 181 4.2e-06
8 182 9.2e-06
8 183 3.4e-06
8 184 1.16e-05
8 185 4.6e-06
8 186 9.6e-06
8 187 2.4e-06
8 188 1.06e-05
8 189 3.6e-06
8 190 1.08e-05
8 191 2.6e-06
8 192 9.4e-06
8 193 2.6e-06
8 194 1e-05
8 195 4.6e-06
8 196 1.12e-05
8 197 4.8e-06
8 198 1.14e-05
8 199 3.4e-06
8 200 1.24e-05
8 201 4.2e-06
8 202 1.12e-05
8 203 3.4e-06
8 204 1.18e-05
8 205 3.6e-06
8 206 1.22e-05
8 207 3.6e-06
8 208 1.02e-05
8 209 2e-06
8 210 1e-05
8 211 5.6e-06
8 212 9.2e-06
8 213 3.6e-06
8 214 1.18e-05
8 215 3.6e-06
8 216 9.4e-06
8 217 3.6e-06
8 218 8.8e-06
8 219 3.2e-06
8 220 1.1e-05
8 221 2.4e-06
8 222 1.4e-05
8 223 5.2e-06
8 224 1.28e-05
8 225 5.6e-06
8 226 1.06e-05
8 227 5.8e-06
8 228 1.14e-05
8 229 5.2e-06
8 230 1.28e-05
8 231 5e-06
8 232 1.44e-05
8 233 5e-06
8 234 1.4e-05
8 235 4.6e-06
8 236 1.24e-05
8 237 3.6e-06
8 238 1.16e-05
8 239 4.8e-06
8 240 1.02e-05
8 241 5.4e-06
8 242 1.7e-05
8 243 5.4e-06
8 244 1.14e-05
8 245 4.2e-06
8 246 8.6e-06
8 247 4.8e-06
8 248 1.06e-05
8 249 4.6e-06
8 250 1.32e-05
8 251 6.6e-06
8 252 1.24e-05
8 253 5.6e-06
8 254 1.2e-05
8 255 3e-06
8 256 9e-06
8 257 3.8e-06
8 258 1.06e-05
8 259 2.6e-06
8 260 1.08e-05
8 261 3e-06
8 262 9.8e-06
8 263 3.6e-06
8 264 1e-05
8 265 3.2e-06
8 266 1.08e-05
8 267 2.8e-06
8 268 1e-05
8 269 3.6e-06
8 270 9.8e-06
8 271 4.6e-06
8 272 1.12e-05
8 273 4.2e-06
8 274 1.16e-05
8 275 4e-06
8 276 1.22e-05
8 277 3e-06
8 278 9e-06
8 279 4e-06
8 280 1.26e-05
8 281 4.6e-06
8 282 1.16e-05
8 283 5.2e-06
8 284 1.1e-05
8 285 2.6e-06
8 286 1.06e-05
8 287 4e-06
8 288 1.12e-05
8 289 3e-06
8 290 9.8e-06
8 291 2.8e-06
8 292 1.16e-05
8 293 6.8e-06
8 294 1.08e-05
8 295 4e-06
8 296 1.14e-05
8 297 5.2e-06
8 298 1.24e-05
8 299 3.2e-06
8 300 1.12e-05
8 301 4.6e-06
8 302 9.8e-06
8 303 3.2e-06
8 304 1.08e-05
8 305 4e-06
8 306 9.6e-06
8 307 4e-06
8 308 1.24e-05
8 309 4.6e-06
8 310 1.08e-05
8 311 6e-06
8 312 1.06e-05
8 313 4e-06
8 314 8e-06
8 315 3.8e-06
8 316 1.08e-05
8 317 3.8e-06
8 318 8.8e-06
8 319 3.6e-06
8 320 1.54e-05
8 321 2.2e-06
8 322 9.4e-06
8 323 1.8e-06
8 324 1.12e-05
8 325 3.4e-06
8 326 9.6e-06
8 327 3e-06
8 328 9.4e-06
8 329 3.4e-06
8 330 1.1e-05
8 331 3.6e-06
8 332 9.8e-06
8 333 1.8e-06
8 334 1.14e-05
8 335 3e-06
8 336 8.4e-06
8 337 3.6e-06
8 338 1.04e-05
8 339 4e-06
8 340 9.2e-06
8 341 4e-06
8 342 9.2e-06
8 343 2.2e-06
8 344 9.8e-06
8 345 3.6e-06
8 346 1.06e-05
8 347 2.6e-06
8 348 1.18e-05
8 349 4e-06
8 350 1e-05
8 351 4.2e-06
8 352 1.36e-05
8 353 3.2e-06
8 354 1.2e-05
8 355 4.2e-06
8 356 1e-05
8 357 5e-06
8 358 1.16e-05
8 359 3.8e-06
8 360 8.6e-06
8 361 2.8e-06
8 362 1.08e-05
8 363 3.6e-06
8 364 9e-06
8 365 3.8e-06
8 366 1.02e-05
8 367 4.6e-06
8 368 9.6e-06
8 369 4e-06
8 370 8.8e-06
8 371 4.4e-06
8 372 9.2e-06
8 373 3.8e-06
8 374 9.8e-06
8 375 4.4e-06
8 376 1e-05
8 377 4e-06
8 378 1.1e-05
8 379 4.6e-06
8 380 9.2e-06
8 381 4.4e-06
8 382 8.6e-06
8 383 2.8e-06
8 384 1.08e-05
8 385 4.6e-06
8 386 1.06e-05
8 387 3.8e-06
8 388 1e-05
8 389 3.4e-06
8 390 1.22e-05
8 391 3.6e-06
8 392 1.14e-05
8 393 3.4e-06
8 394 9.4e-06
8 395 3.8e-06
8 396 1.22e-05
8 397 4.4e-06
8 398 9.8e-06
8 399 2.2e-06
8 400 8.2e-06
8 401 4.2e-06
8 402 9.4e-06
8 403 5.4e-06
8 404 1.1e-05
8 405 4e-06
8 406 1e-05
8 407 5e-06
8 408 1.08e-05
8 409 4.6e-06
8 410 8.6e-06
8 411 5e-06
8 412 7.6e-06
8 413 5.2e-06
8 414 1.18e-05
8 415 2.8e-06
8 416 8.8e-06
8 417 6.2e-06
8 418 8.6e-06
8 419 3.8e-06
8 420 7.6e-06
8 421 3.6e-06
8 422 9.2e-06
8 423 3.6e-06
8 424 9e-06
8 425 3.6e-06
8 426 1.04e-05
8 427 4.8e-06
8 428 1.02e-05
8 429 3.8e-06
8 430 8e-06
8 431 3.2e-06
8 432 9.8e-06
8 433 2.8e-06
8 434 8.8e-06
8 435 1.6e-06
8 436 8.8e-06
8 437 4.4e-06
8 438 9.6e-06
8 439 3.6e-06
8 440 7e-06
8 441 3e-06
8 442 8.2e-06
8 443 5.4e-06
8 444 7e-06
8 445 2.8e-06
8 446 7e-06
8 447 3.2e-06
8 448 9.8e-06
8 449 3.8e-06
8 450 8.4e-06
8 451 5.4e-06
8 452 9.2e-06
8 453 4.2e-06
8 454 8.6e-06
8 455 3.6e-06
8 456 7.4e-06
8 457 4.8e-06
8 458 7.4e-06
8 459 3.4e-06
8 460 7.8e-06
8 461 4.6e-06
8 462 6.8e-06
8 463 4.8e-06
8 464 9.2e-06
8 465 2.6e-06
8 466 7.4e-06
8 467 3.8e-06
8 468 7.8e-06
8 469 4.8e-06
8 470 6.8e-06
8 471 2.8e-06
8 472 1.08e-05
8 473 2.2e-06
8 474 8.8e-06
8 475 4e-06
8 476 7.4e-06
8 477 2.8e-06
8 478 8.8e-06
8 479 5.8e-06
8 480 7.4e-06
8 481 3.4e-06
8 482 7.2e-06
8 483 5.4e-06
8 484 1.12e-05
8 485 3.4e-06
8 486 5.6e-06
8 487 3.4e-06
8 488 1.04e-05
8 489 3.2e-06
8 490 9.8e-06
8 491 5e-06
8 492 6e-06
8 493 3.4e-06
8 494 9.4e-06
8 495 2.2e-06
8 496 1.02e-05
8 497 5e-06
8 498 8.2e-06
8 499 3.6e-06
8 500 7.4e-06
8 501 3e-06
8 502 8.8e-06
8 503 3.8e-06
8 504 7.2e-06
8 505 4.4e-06
8 506 6e-06
8 507 4e-06
8 508 6.2e-06
8 509 3e-06
8 510 7e-06
8 511 5.6e-06
8 512 3.6e-06
8 513 6e-07
8 514 3.4e-06
8 515 6e-07
8 516 4.2e-06
8 517 6e-07
8 518 4.6e-06
8 519 6e-07
8 520 2.6e-06
8 521 6e-07
8 522 4.6e-06
8 523 0
8 524 2.8e-06
8 525 2e-07
8 526 4.4e-06
8 527 1e-06
8 528 4e-06
8 529 1.4e-06
8 530 4.4e-06
8 531 8e-07
8 532 3.4e-06
8 533 1.2e-06
8 534 2.4e-06
8 535 2e-07
8 536 3.4e-06
8 537 1e-06
8 538 4e-06
8 539 1.4e-06
8 540 3.2e-06
8 541 1e-06
8 542 3.8e-06
8 543 1.4e-06
8 544 3.6e-06
8 545 4e-07
8 546 5.2e-06
8 547 6e-07
8 548 3.2e-06
8 549 1e-06
8 550 3.2e-06
8 551 1.8e-06
8 552 3e-06
8 553 6e-07
8 554 3e-06
8 555 8e-07
8 556 3.8e-06
8 557 6e-07
8 558 4.4e-06
8 559 1.6e-06
8 560 2.2e-06
8 561 8e-07
8 562 2.2e-06
8 563 1e-06
8 564 3.8e-06
8 565 2e-07
8 566 2.4e-06
8 567 6e-07
8 568 3.4e-06
8 569 8e-07
8 570 2.8e-06
8 571 1e-06
8 572 3.2e-06
8 573 1.2e-06
8 574 1.8e-06
8 575 1e-06
8 576 2.2e-06
8 577 1.2e-06
8 578 2.6e-06
8 579 6e-07
8 580 2.8e-06
8 581 6e-07
8 582 1.6e-06
8 583 1.6e-06
8 584 2.4e-06
8 585 1.2e-06
8 586 4.2e-06
8 587 2e-07
8 588 3.2e-06
8 589 6e-07
8 590 3.2e-06
8 591 1e-06
8 592 3.6e-06
8 593 1.2e-06
8 594 3.4e-06
8 595 4e-07
8 596 2.2e-06
8 597 0
8 598 2.8e-06
8 599 6e-07
8 600 3e-06
8 601 1.6e-06
8 602 3.2e-06
8 603 1e-06
8 604 3e-06
8 605 8e-07
8 606 1.6e-06
8 607 0
8 608 2.6e-06
8 609 1e-06
8 610 2.8e-06
8 611 1e-06
8 612 2.2e-06
8 613 4e-07
8 614 2.6e-06
8 615 6e-07
8 616 1.8e-06
8 617 1.4e-06
8 618 2.2e-06
8 619 1.4e-06
8 620 2e-06
8 621 1.6e-06
8 622 2e-06
8 623 6e-07
8 624 1.8e-06
8 625 0
8 626 3.2e-06
8 627 1.2e-06
8 628 2.4e-06
8 629 6e-07
8 630 2e-06
8 631 6e-07
8 632 1e-06
8 633 4e-07
8 634 3e-06
8 635 1.8e-06
8 636 1.8e-06
8 637 1e-06
8 638 1.4e-06
8 639 4e-07
8 640 1.4e-06
8 641 4e-07
8 642 3.2e-06
8 643 8e-07
8 644 2e-06
8 645 1.2e-06
8 646 2.2e-06
8 647 4e-07
8 648 2.2e-06
8 649 4e-07
8 650 1.8e-06
8 651 8e-07
8 652 1.4e-06
8 653 1.6e-06
8 654 1.8e-06
8 655 6e-07
8 656 2.4e-06
8 657 1.4e-06
8 658 2e-06
8 659 6e-07
8 660 1.4e-06
8 661 6e-07
8 662 1.8e-06
8 663 1.4e-06
8 664 1.6e-06
8 665 1.4e-06
8 666 1.2e-06
8 667 6e-07
8 668 3e-06
8 669 8e-07
8 670 1.4e-06
8 671 0
8 672 8e-07
8 673 6e-07
8 674 1e-06
8 675 4e-07
8 676 1.8e-06
8 677 8e-07
8 678 1.6e-06
8 679 1.6e-06
8 680 2.6e-06
8 681 1e-06
8 682 2e-06
8 683 2e-07
8 684 2.4e-06
8 685 6e-07
8 686 2.4e-06
8 687 8e-07
8 688 1.2e-06
8 689 2e-07
8 690 1.6e-06
8 691 8e-07
8 692 1.6e-06
8 693 8e-07
8 694 1.2e-06
8 695 4e-07
8 696 6e-07
8 697 4e-07
8 698 1.8e-06
8 699 1.4e-06
8 700 1.8e-06
8 701 1.6e-06
8 702 1.2e-06
8 703 2e-07
8 704 1.4e-06
8 705 8e-07
8 706 6e-07
8 707 6e-07
8 708 1.8e-06
8 709 1.4e-06
8 710 1.4e-06
8 711 1e-06
8 712 2.2e-06
8 713 1.2e-06
8 714 8e-07
8 715 4e-07
8 716 1.2e-06
8 717 6e-07
8 718 1.2e-06
8 719 2e-07
8 720 8e-07
8 721 1e-06
8 722 1.2e-06
8 723 6e-07
8 724 1.2e-06
8 725 4e-07
8 726 8e-07
8 727 1.2e-06
8 728 1e-06
8 729 6e-07
8 730 1.4e-06
8 731 1.2e-06
8 732 1.8e-06
8 733 4e-07
8 734 8e-07
8 735 8e-07
8 736 1.6e-06
8 737 1.2e-06
8 738 2e-06
8 739 8e-07
8 740 6e-07
8 741 8e-07
8 742 1e-06
8 743 6e-07
8 744 1e-06
8 745 2e-07
8 746 1.2e-06
8 747 4e-07
8 748 4e-07
8 749 1.8e-06
8 750 1.8e-06
8 751 0
8 752 1.8e-06
8 753 4e-07
8 754 6e-07
8 755 2e-07
8 756 6e-07
8 757 2e-07
8 758 8e-07
8 759 4e-07
8 760 2e-07
8 761 8e-07
8 762 6e-07
8 763 6e-07
8 764 1.6e-06
8 765 0
8 766 1.6e-06
8 767 1.4e-06

discrete_count_bits

#include <iostream>
#include <random>
#include <exrandom/rand_digit.hpp>
#include <exrandom/i_rand.hpp>
#include <exrandom/discrete_normal_dist.hpp>
#include <exrandom/aux_info.hpp>
template<typename Generator>
void count_bits(Generator& g, long long num,
                int mu_num, int mu_den, int sigma_num, int sigma_den) {
  const unsigned b = 2;
  exrandom::rand_digit<b> D;
  exrandom::discrete_normal_dist<exrandom::rand_digit<b>>
    N(D, mu_num, mu_den, sigma_num, sigma_den);
  exrandom::i_rand<exrandom::rand_digit<b>> h(N.digit_generator());
  long long counta = 0, countb = 0, fractotal = 0;
  for (long long i = 0; i < num; ++i) {
    long long c0 = D.count();
    N.generate(g, h);
    counta += D.count() - c0;
    fractotal += h.entropy();
    h(g);
    countb += D.count() - c0;
  }
  double H;
  exrandom::aux_info::discrete_normal_norm(N.param(), H);
  H /= std::log(2.0);
  double B = counta/double(num),
    L = fractotal/double(num),
    F = countb/double(num),
    T = F - H;
  std::cout << "mu = " << N.mu_num() << "/" << N.mu_den()
            << ", sigma = " << N.sigma_num() << "/" << N.sigma_den()
            << ":\n";
  std::cout << " <B> = " << B
            << ", <L> = " << L
            << ", <F> = <B> + <L> = " << F
            << ",\n H = " << H
            << ", T = <F> - H = " << T << "" << std::endl;
}
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  long long num = 2000000LL;
  std::cout
    << "Sampling from the discrete normal distribution (Algorithm D)\n";
  count_bits(g, num,  1, 7, 16, 100);
  count_bits(g, num,  2, 7, 16, 10);
  count_bits(g, num,  3, 7, 16, 1);
  count_bits(g, num,  0, 7, 160, 1);
  count_bits(g, num, -1, 7, 1600, 1);
  count_bits(g, num, -2, 7, 16000, 1);
  count_bits(g, num, -3, 7, 160000, 1);
  count_bits(g, num,  1, 7, 1600000, 1);
  count_bits(g, num,  0, 7, (1<<17)-1, 1);
  count_bits(g, num,  0, 7, (1<<17),   1);
  count_bits(g, num,  0, 1, (1<<17)+1, 1);
  count_bits(g, num,  0, 1, (1<<18)-1, 1);
  count_bits(g, num,  0, 1, (1<<18),   1);
  count_bits(g, num,  0, 1, (1<<18)+1, 1);
}

count_bits

#include <iostream>
#include <random>
#include <map>
#include <limits>
#include <exrandom/rand_digit.hpp>
#include <exrandom/u_rand.hpp>
#include <exrandom/unit_normal_dist.hpp>
#include <exrandom/unit_normal_kahn.hpp>
#include <exrandom/unit_exponential_dist.hpp>
#include <exrandom/unit_uniform_dist.hpp>
#include <exrandom/aux_info.hpp>
void print_hist(const std::map<int, long long>& hist) {
  const int maxlength = 70;     // Widest bar in histogram
  long long m = 0;
  for (auto p = hist.begin(); p != hist.end(); ++p)
    m = (std::max)(m, p->second);
  m = (std::max)(m / maxlength, 1LL);
  for (auto p = hist.begin(); p != hist.end(); ++p) {
    unsigned k = unsigned((p->second + m/2) / m);
    if (k) {
      std::cout << p->first << '\t';
      for (unsigned i = 0; i < k; ++i)
        std::cout << '*';
      std::cout << '\n';
    }
  }
  std::cout << std::endl;
}
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n" << std::endl;
  const unsigned b = 2;
  long long num = 1000000LL;
  std::cout << std::fixed << std::setprecision(5);
  {
    std::cout << "Sampling from unit normal distribution\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_normal_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>>
      x(N.digit_generator());
    std::map<int, long long> bitcount;
    std::map<int, long long> fraccount;
    long long counta = 0, countb = 0, fractotal = 0;
    for (long long i = 0; i < num; ++i) {
      long long c0 = D.count();
      N.generate(g, x);
      counta += D.count() - c0;
      fractotal += x.ndigits();
      ++fraccount[int(x.ndigits())];
      x.value<double>(g);
      countb += D.count() - c0;
      ++bitcount[int(D.count() - c0)];
    }
    double H = exrandom::aux_info::normal_entropy() / std::log(2.0),
      Q = exrandom::aux_info::normal_mean_exponent(),
      B = counta/double(num),
      L = fractotal/double(num),
      C = B - L, T = C - H,
      F = countb/double(num);
    int d = std::numeric_limits<double>::digits + 1;
    std::cout << "<B> = " << B
              << ", <L> = " << L
              << ", C = <B> - <L> = " << C
              << ",\nH = " << H
              << ", T = C - H = " << T
              << ",\nQ = " << Q
              << ", F = C - Q + " << d << " = " << C - Q + d
              << " (predicted) = " << F << " (measured)\n" << std::endl;
    std::cout << "Histogram of number of bits in fraction\n";
    print_hist(fraccount);
    std::cout << "Histogram of number of bits needed to generate doubles\n";
    print_hist(bitcount);
  }
  {
    std::cout << "Sampling from unit normal distribution (Kahn method)\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_normal_kahn<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>>
      x(N.digit_generator());
    std::map<int, long long> bitcount;
    std::map<int, long long> fraccount;
    long long counta = 0, countb = 0, fractotal = 0;
    for (long long i = 0; i < num; ++i) {
      long long c0 = D.count();
      N.generate(g, x);
      counta += D.count() - c0;
      fractotal += x.ndigits();
      ++fraccount[int(x.ndigits())];
      x.value<double>(g);
      countb += D.count() - c0;
      ++bitcount[int(D.count() - c0)];
    }
    double H = exrandom::aux_info::normal_entropy() / std::log(2.0),
      Q = exrandom::aux_info::normal_mean_exponent(),
      B = counta/double(num),
      L = fractotal/double(num),
      C = B - L, T = C - H,
      F = countb/double(num);
    int d = std::numeric_limits<double>::digits + 1;
    std::cout << "<B> = " << B
              << ", <L> = " << L
              << ", C = <B> - <L> = " << C
              << ",\nH = " << H
              << ", T = C - H = " << T
              << ",\nQ = " << Q
              << ", F = C - Q + " << d << " = " << C - Q + d
              << " (predicted) = " << F << " (measured)\n" << std::endl;
    std::cout << "Histogram of number of bits in fraction\n";
    print_hist(fraccount);
    std::cout << "Histogram of number of bits needed to generate doubles\n";
    print_hist(bitcount);
  }
  {
    std::cout
      << "Sampling from unit exponential distribution (Algorithm V)\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_exponential_dist<exrandom::rand_digit<b>,false> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>>
      x(N.digit_generator());
    std::map<int, long long> bitcount;
    std::map<int, long long> fraccount;
    long long counta = 0, countb = 0, fractotal = 0;
    for (long long i = 0; i < num; ++i) {
      long long c0 = D.count();
      N.generate(g, x);
      counta += D.count() - c0;
      fractotal += x.ndigits();
      ++fraccount[int(x.ndigits())];
      x.value<double>(g);
      countb += D.count() - c0;
      ++bitcount[int(D.count() - c0)];
    }
    double H = exrandom::aux_info::exponential_entropy() / std::log(2.0),
      Q = exrandom::aux_info::exponential_mean_exponent(),
      B = counta/double(num),
      L = fractotal/double(num),
      C = B - L, T = C - H,
      F = countb/double(num);
    int d = std::numeric_limits<double>::digits + 1;
    std::cout << "<B> = " << B
              << ", <L> = " << L
              << ", C = <B> - <L> = " << C
              << ",\nH = " << H
              << ", T = C - H = " << T
              << ",\nQ = " << Q
              << ", F = C - Q + " << d << " = " << C - Q + d
              << " (predicted) = " << F << " (measured)\n" << std::endl;
  }
  {
    std::cout
      << "Sampling from unit exponential distribution (Algorithm E)\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_exponential_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>>
      x(N.digit_generator());
    std::map<int, long long> bitcount;
    std::map<int, long long> fraccount;
    long long counta = 0, countb = 0, fractotal = 0;
    for (long long i = 0; i < num; ++i) {
      long long c0 = D.count();
      N.generate(g, x);
      counta += D.count() - c0;
      fractotal += x.ndigits();
      ++fraccount[int(x.ndigits())];
      x.value<double>(g);
      countb += D.count() - c0;
      ++bitcount[int(D.count() - c0)];
    }
    double H = exrandom::aux_info::exponential_entropy() / std::log(2.0),
      Q = exrandom::aux_info::exponential_mean_exponent(),
      B = counta/double(num),
      L = fractotal/double(num),
      C = B - L, T = C - H,
      F = countb/double(num);
    int d = std::numeric_limits<double>::digits + 1;
    std::cout << "<B> = " << B
              << ", <L> = " << L
              << ", C = <B> - <L> = " << C
              << ",\nH = " << H
              << ", T = C - H = " << T
              << ",\nQ = " << Q
              << ", F = C - Q + " << d << " = " << C - Q + d
              << " (predicted) = " << F << " (measured)\n" << std::endl;
    std::cout << "Histogram of number of bits in fraction\n";
    print_hist(fraccount);
    std::cout << "Histogram of number of bits needed to generate doubles\n";
    print_hist(bitcount);
  }
  {
    std::cout
      << "Sampling from unit uniform distribution\n";
    exrandom::rand_digit<b> D;
    exrandom::unit_uniform_dist<exrandom::rand_digit<b>> N(D);
    exrandom::u_rand<exrandom::rand_digit<b>> x(N.digit_generator());
    std::map<int, long long> bitcount;
    std::map<int, long long> fraccount;
    long long counta = 0, countb = 0, fractotal = 0;
    for (long long i = 0; i < num; ++i) {
      long long c0 = D.count();
      N.generate(g, x);
      counta += D.count() - c0;
      fractotal += x.ndigits();
      ++fraccount[int(x.ndigits())];
      x.value<double>(g);
      countb += D.count() - c0;
      ++bitcount[int(D.count() - c0)];
    }
    double H = exrandom::aux_info::uniform_entropy() / std::log(2.0),
      Q = exrandom::aux_info::uniform_mean_exponent(),
      B = counta/double(num),
      L = fractotal/double(num),
      C = B - L, T = C - H,
      F = countb/double(num);
    int d = std::numeric_limits<double>::digits + 1;
    std::cout << "<B> = " << B
              << ", <L> = " << L
              << ", C = <B> - <L> = " << C
              << ",\nH = " << H
              << ", T = C - H = " << T
              << ",\nQ = " << Q
              << ", F = C - Q + " << d << " = " << C - Q + d
              << " (predicted) = " << F << " (measured)\n" << std::endl;
    std::cout << "Histogram of number of bits in fraction\n";
    print_hist(fraccount);
    std::cout << "Histogram of number of bits needed to generate doubles\n";
    print_hist(bitcount);
  }
}

Seed set to 1699814630

Sampling from unit normal distribution
= 24.05239, = 1.44300, C = - = 22.60939,
H = 2.04710, T = C - H = 20.56229,
Q = -0.41664, F = C - Q + 54 = 77.02603 (predicted) = 77.02741 (measured)

Histogram of number of bits in fraction
0 **********************************************************************
1 **************************************
2 *****************************
3 ******************
4 **********
5 *****
6 ***
7 *
8 *

Histogram of number of bits needed to generate doubles
58 **************************************************************
59 ***********************************************
60 **********************************************************************
61 *******************************************
62 ********************************************************
63 *************************************
64 *****************************************
65 ******************************
66 *********************************
67 ***************************
68 ***************************
69 ************************
70 ************************
71 *********************
72 *********************
73 ********************
74 *******************
75 ******************
76 *****************
77 ****************
78 ***************
79 **************
80 **************
81 *************
82 *************
83 ************
84 ************
85 ***********
86 ***********
87 **********
88 **********
89 *********
90 *********
91 *********
92 ********
93 ********
94 *******
95 *******
96 *******
97 ******
98 ******
99 ******
100 ******
101 *****
102 *****
103 *****
104 *****
105 ****
106 ****
107 ****
108 ****
109 ****
110 ****
111 ***
112 ***
113 ***
114 ***
115 ***
116 ***
117 ***
118 ***
119 **
120 **
121 **
122 **
123 **
124 **
125 **
126 **
127 **
128 **
129 **
130 *
131 *
132 *
133 *
134 *
135 *
136 *
137 *
138 *
139 *
140 *
141 *
142 *
143 *
144 *
145 *
146 *
147 *
148 *
149 *
150 *
151 *
152 *

Sampling from unit normal distribution (Kahn method)
= 21.45269, = 2.01914, C = - = 19.43355,
H = 2.04710, T = C - H = 17.38646,
Q = -0.41664, F = C - Q + 54 = 73.85019 (predicted) = 73.84885 (measured)

Histogram of number of bits in fraction
1 **********************************************************************
2 *************************
3 ***************
4 *********
5 *****
6 **
7 *
8 *

Histogram of number of bits needed to generate doubles
60 ******************************************************************
61 **********************************************************************
62 ******************************************************************
63 **********************************************************
64 ******************************************************
65 ************************************************
66 **********************************************
67 ********************************************
68 ******************************************
69 ***************************************
70 ************************************
71 **********************************
72 *******************************
73 *****************************
74 ***************************
75 *************************
76 ***********************
77 *********************
78 ********************
79 ******************
80 *****************
81 ****************
82 ***************
83 **************
84 *************
85 ************
86 ***********
87 **********
88 *********
89 *********
90 ********
91 ********
92 *******
93 *******
94 ******
95 ******
96 ******
97 *****
98 *****
99 *****
100 ****
101 ****
102 ****
103 ****
104 ***
105 ***
106 ***
107 ***
108 ***
109 **
110 **
111 **
112 **
113 **
114 **
115 **
116 **
117 *
118 *
119 *
120 *
121 *
122 *
123 *
124 *
125 *
126 *
127 *
128 *
129 *
130 *
131 *
132 *
133 *

Sampling from unit exponential distribution (Algorithm V)
= 9.31982, = 2.05436, C = - = 7.26546,
H = 1.44270, T = C - H = 5.82277,
Q = -0.33275, F = C - Q + 54 = 61.59821 (predicted) = 61.59797 (measured)

Sampling from unit exponential distribution (Algorithm E)
= 7.23335, = 1.74352, C = - = 5.48983,
H = 1.44270, T = C - H = 4.04714,
Q = -0.33275, F = C - Q + 54 = 59.82258 (predicted) = 59.82374 (measured)

Histogram of number of bits in fraction
1 **********************************************************************
2 *******************
3 **********
4 *****
5 ***
6 *
7 *

Histogram of number of bits needed to generate doubles
56 **********************************************************************
57 ***********************************
58 ***************************
59 **************
60 **************
61 **********
62 **********
63 *******
64 ******
65 *****
66 ****
67 ***
68 ***
69 **
70 **
71 **
72 *
73 *
74 *
75 *
76 *
77 *
78 *

Sampling from unit uniform distribution
= 0.00000, = 0.00000, C = - = 0.00000,
H = 0.00000, T = C - H = 0.00000,
Q = -1.00000, F = C - Q + 54 = 55.00000 (predicted) = 54.99764 (measured)

Histogram of number of bits in fraction
0 **********************************************************************

Histogram of number of bits needed to generate doubles
54 **********************************************************************
55 ***********************************
56 *****************
57 *********
58 ****
59 **
60 *
61 *

compare_times

#include <iostream>
#include <iomanip>
#include <random>
#include <chrono>
#include <exrandom/unit_normal_distribution.hpp>
#include <exrandom/unit_exponential_distribution.hpp>
#include <exrandom/discrete_normal_distribution.hpp>
template<typename Dist, typename Generator>
double timer(long long num, Dist& d, Generator& g) {
  typename Dist::result_type sum = 0;
  auto t0 = std::chrono::high_resolution_clock::now();
  for (long long i = 0; i < num; ++i)
    sum += d(g);
  auto t1 = std::chrono::high_resolution_clock::now();
  double dt = double(std::chrono::duration_cast<std::chrono::nanoseconds>
                       (t1 - t0).count());
  return dt/num;
}
template<typename Generator>
void discrete_timer(Generator& g, long long num,
                    int mu_num, int mu_den, int sigma_num, int sigma_den) {
  exrandom::discrete_normal_distribution d(mu_num, mu_den,
                                           sigma_num, sigma_den);
  double t = timer(num, d, g);
  std::cout << "  time with mu = " << d.mu_num() << "/" << d.mu_den()
            << ", sigma = " << d.sigma_num() << "/" << d.sigma_den()
            << ": " << t << " ns" << std::endl;
}
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "Seed set to " << s << "\n\n";
  long long num = 5000000LL;
  std::cout << std::fixed << std::setprecision(1);
  std::cout << "Compare times to sample from the normal distribution\n";
  {
    typedef float real;
    std::normal_distribution<real> d1;
    double t1 = timer(num, d1, g);
    exrandom::unit_normal_distribution<real> d2;
    double t2 = timer(num, d2, g);
    std::cout << "  time with prec "
              << std::numeric_limits<real>::digits
              << ": C++11/random = " << t1
              << " ns; exrandom = " << t2 << " ns" << std::endl;
  }
  {
    typedef double real;
    std::normal_distribution<real> d1;
    double t1 = timer(num, d1, g);
    exrandom::unit_normal_distribution<real> d2;
    double t2 = timer(num, d2, g);
    std::cout << "  time with prec "
              << std::numeric_limits<real>::digits
              << ": C++11/random = " << t1
              << " ns; exrandom = " << t2 << " ns" << std::endl;
  }
  {
    // For Visual Studio long double is the same as double
    typedef long double real;
    std::normal_distribution<real> d1;
    double t1 = timer(num, d1, g);
    exrandom::unit_normal_distribution<real> d2;
    double t2 = timer(num, d2, g);
    std::cout << "  time with prec "
              << std::numeric_limits<real>::digits
              << ": C++11/random = " << t1
              << " ns; exrandom = " << t2 << " ns" << std::endl;
  }
  std::cout << "Compare times to sample from the exponential distribution\n";
  {
    typedef float real;
    std::exponential_distribution<real> d1;
    double t1 = timer(num, d1, g);
    exrandom::unit_exponential_distribution<real> d2;
    double t2 = timer(num, d2, g);
    std::cout << "  time with prec "
              << std::numeric_limits<real>::digits
              << ": C++11/random = " << t1
              << " ns; exrandom = " << t2 << " ns" << std::endl;
  }
  {
    typedef double real;
    std::exponential_distribution<real> d1;
    double t1 = timer(num, d1, g);
    exrandom::unit_exponential_distribution<real> d2;
    double t2 = timer(num, d2, g);
    std::cout << "  time with prec "
              << std::numeric_limits<real>::digits
              << ": C++11/random = " << t1
              << " ns; exrandom = " << t2 << " ns" << std::endl;
  }
  {
    // For Visual Studio long double is the same as double
    typedef long double real;
    std::exponential_distribution<real> d1;
    double t1 = timer(num, d1, g);
    exrandom::unit_exponential_distribution<real> d2;
    double t2 = timer(num, d2, g);
    std::cout << "  time with prec "
              << std::numeric_limits<real>::digits
              << ": C++11/random = " << t1
              << " ns; exrandom = " << t2 << " ns" << std::endl;
  }
  std::cout << "Times to sample from the discrete normal distribution\n";
  discrete_timer(g, num,  1, 7, 16, 100);
  discrete_timer(g, num,  2, 7, 16, 10);
  discrete_timer(g, num,  3, 7, 16, 1);
  discrete_timer(g, num,  0, 7, 160, 1);
  discrete_timer(g, num, -1, 7, 1600, 1);
  discrete_timer(g, num, -2, 7, 16000, 1);
  discrete_timer(g, num, -3, 7, 160000, 1);
  discrete_timer(g, num,  1, 7, 1600000, 1);
  discrete_timer(g, num,  0, 7, (1<<17)-1, 1);
  discrete_timer(g, num,  0, 7, (1<<17),   1);
  discrete_timer(g, num,  0, 1, (1<<17)+1, 1);
  discrete_timer(g, num,  0, 1, (1<<18)-1, 1);
  discrete_timer(g, num,  0, 1, (1<<18),   1);
  discrete_timer(g, num,  0, 1, (1<<18)+1, 1);
}

chisq_test

#include <iostream>
#include <iomanip>
#include <cmath>
#include <vector>
#include <map>
#include <exrandom/unit_normal_distribution.hpp>
#include <exrandom/unit_exponential_distribution.hpp>
#include <exrandom/unit_uniform_distribution.hpp>
#include <exrandom/discrete_normal_distribution.hpp>
#include <exrandom/aux_info.hpp>
// Compute probs for num bins [x0+n*dx, x0+(n+1)*dx] for n=[0,num) plus an
// (num+1)th bin for everything else
std::vector<double> probs(double (*f)(double),
                          double x0, double dx, int nbins) {
  std::vector<double> r(nbins+1, 0);
  double s = 0;
  for (int n = 0; n < nbins; ++n) {
    r[n] = f(x0 + (n + 1) * dx) - f(x0 + n * dx);
    s += r[n];
  }
  r[nbins] = 1 - s;
  return r;
}
// Compute probs for discrete normal using num bins [x0+n*dx, x0+(n+1)*dx) for
// i=[0,num) plus an (num+1)th bin for everything else
template<typename param_type>
std::vector<double> discrete_normal_probs(const param_type& p,
                                          int x0, int dx, int nbins) {
  double H;
  double norm = exrandom::aux_info::discrete_normal_norm(p, H);
  std::vector<double> r(nbins+1, 0);
  double s = 0;
  for (int n = 0; n < nbins; ++n) {
    for (int j = 0; j < dx; ++j)
      r[n] +=
        exrandom::aux_info::discrete_normal_prob(p, x0 + dx * n + j, norm);
    s += r[n];
  }
  r[nbins] = 1 - s;
  return r;
}
double chisqf(const std::vector<double>& probs,
              std::map<int, long long>& counts) {
  int nbins = int(probs.size()) - 1;
  if (nbins < 1) return 0;
  long long num = 0;
  for (auto p = counts.begin(); p != counts.end(); ++p)
    num += p->second;
  double v = 0, x;
  long long s = 0;
  for (int n = 0; n < nbins; ++n) {
    s += counts[n];
    x = counts[n] - num * probs[n];
    v += x * x / (num * probs[n]);
  }
  x = (num - s) - num * probs[nbins];
  v += x * x / (num * probs[nbins]);
  return v;
}
template<typename Generator>
void discrete_chisq(Generator& g, long long num,
                    int mu_num, int mu_den, int sigma_num, int sigma_den,
                    int x0, int dx, int DOF) {
    exrandom::discrete_normal_distribution
      d(mu_num, mu_den, sigma_num, sigma_den);
    double mu = d.mu_num() / double(d.mu_den()),
      sigma = d.sigma_num() / double(d.sigma_den());
    std::map<int, long long> hist;
    for (long long i = 0; i < num; ++i)
      ++hist[int(std::floor((d(g) - x0 + 0.5)/dx))];
    double chisq = chisqf(discrete_normal_probs(d.param(), x0, dx, DOF),
                          hist);
    std::cout << "discrete_normal_distribution (mu = " << mu
              << ", sigma = " << sigma << "):\n            samples = "
              << num << ", DOF = " << DOF
              << ", chi-squared = " << chisq << std::endl;
  }
int main() {
  unsigned s = std::random_device()(); // Set seed from random_device
  std::mt19937 g(s);                   // Initialize URNG
  std::cout << "With 50 degrees of freedom, chisq should lie\n"
            << "   between 29.71 and 76.15, 98% of the time\n"
            << "   between 34.76 and 67.50, 90% of the time\n"
            << "   between 42.94 and 56.33, 50% of the time\n"
            << "See Knuth TAOCP, Vol 2, Sec. 3.3.1\n"
            << "Seed set to " << s << "\n\n";
  std::cout << std::fixed << std::setprecision(2);
  long long num =  5000000LL;
  {
    exrandom::unit_normal_distribution<double> d;
    int DOF = 50;
    double x0 = -4, dx = 0.16;  // 50 bins in [-4,4]
    std::map<int, long long> hist;
    for (long long i = 0; i < num; ++i)
      ++hist[int(std::floor( (d(g) - x0) / dx ))];
    double chisq = chisqf(probs(exrandom::aux_info::cumulative_normal,
                                x0, dx, DOF),
                          hist);
    std::cout << "unit_normal_distribution: samples = "
              << num << ", DOF = " << DOF
              << ", chi-squared = " << chisq << std::endl;
  }
  {
    exrandom::unit_exponential_distribution<double> d;
    int DOF = 50;
    double x0 = 0, dx = 0.16;   // 50 bins in [0,8]
    std::map<int, long long> hist;
    for (long long i = 0; i < num; ++i)
      ++hist[int(std::floor( (d(g) - x0) / dx ))];
    double chisq = chisqf(probs(exrandom::aux_info::cumulative_exponential,
                                x0, dx, DOF),
                          hist);
    std::cout << "unit_exponential_distribution: samples = "
              << num << ", DOF = " << DOF
              << ", chi-squared = " << chisq << std::endl;
  }
  {
    exrandom::unit_uniform_distribution<double> d;
    int DOF = 50;
    double x0 = 0, dx = 1/51.0; // 50 bins in [0,50/51]
    std::map<int, long long> hist;
    for (long long i = 0; i < num; ++i)
      ++hist[int(std::floor( (d(g) - x0) / dx ))];
    double chisq = chisqf(probs(exrandom::aux_info::cumulative_uniform,
                                x0, dx, DOF),
                          hist);
    std::cout << "unit_uniform_distribution: samples = "
              << num << ", DOF = " << DOF
              << ", chi-squared = " << chisq << std::endl;
  }
  discrete_chisq(g, num, 0, 1, 6, 1, -24, 1, 50);
  discrete_chisq(g, num, 1, 3, 6, 1, -24, 1, 50);
  discrete_chisq(g, num, 0, 1, 61, 10, -24, 1, 50);
  discrete_chisq(g, num, -5, 3, 69, 10, -24, 1, 50);
  discrete_chisq(g, num, 201, 7, 1301, 2, -2500, 100, 50);
}