exrandom_test.cpp

#include <iomanip>
#include <limits>
#include <cmath>
#include <exrandom/unit_normal_distribution.hpp>
#include <exrandom/unit_exponential_distribution.hpp>
#include <exrandom/unit_uniform_distribution.hpp>
#include <exrandom/discrete_normal_distribution.hpp>

int main() {
  static_assert(std::numeric_limits<double>::radix == 2 &&
                std::numeric_limits<double>::digits == 53,
                "tests require that doubles have precision of 53 binary bits");
  int retval = 0;
  std::cerr << std::fixed << std::setprecision(11);
  double tol = 0.5e-11;
  std::mt19937 g(5489u);        // Initialize URNG with default seed
  {
    unsigned long x, z = 4123659995;
    for (unsigned i = 0; i < 10000; ++i)
      x = g();
    if (x != z) {
      ++retval;
      std::cerr << "Error in std::mt19937:\n"
                << "  expected " << z << ", got " << x << "\n";
    }
  }
  {
    g.seed(1u);
    volatile double x = 0, y, z = -173.53065882716;
    exrandom::unit_uniform_distribution<double> U;
    for (unsigned i = 0; i < 1000000; ++i) {
      y = U(g) - 0.5;
      x += y;
    }
    if (std::abs(x - z) > tol) {
      ++retval;
      std::cerr << "Error in exrandom::unit_uniform_distribution:\n"
                << "  expected " << z << ", got " << x << "\n";
    }
  }
  {
    g.seed(2u);
    volatile double x = 0, y, z = 708.92395157383;
    exrandom::unit_exponential_distribution<double> E;
    for (unsigned i = 0; i < 1000000; ++i) {
      y = E(g) - 1.0;
      x += y;
    }
    if (std::abs(x - z) > tol) {
      ++retval;
      std::cerr << "Error in exrandom::unit_exponential_distribution:\n"
                << "  expected " << z << ", got " << x << "\n";
    }
  }
  {
    g.seed(3u);
    volatile double x = 0, z = 1598.67909684123;
    exrandom::unit_normal_distribution<double> N;
    for (unsigned i = 0; i < 1000000; ++i)
      x += N(g);
    if (std::abs(x - z) > tol) {
      ++retval;
      std::cerr << "Error in exrandom::unit_normal_distribution:\n"
                << "  expected " << z << ", got " << x << "\n";
    }
  }
  {
    g.seed(4u);
    long x = 0, z = 281125;
    exrandom::discrete_normal_distribution D(1,3,129,2);
    for (unsigned i = 0; i < 1000000; ++i)
      x += D(g);
    if (x != z) {
      ++retval;
      std::cerr << "Error in exrandom::discrete_normal_distribution:\n"
                << "  expected " << z << ", got " << x << "\n";
    }
  }
  if (retval > 0)
    std::cerr << "Some tests failed" << "\n";
  else
    std::cerr << "Tests completed successfully" << "\n";
  return retval;
}