rng.cc

Go to the documentation of this file.

#define TF_RNG_CC_VERSION \
  "TF_RNG_CC   V1.1"

#include <iostream>
#include <tfxx/rng.h>
#include <time.h>
#include <gsl/gsl_rng.h>
#include <string>

namespace tfxx {

  namespace numeric {

  /* ---------------------------------------------------------------------- */
  // class RNGgaussian
      
    // constructor
    RNGgaussian::RNGgaussian(const double& std,
                             const double& mean,
                             const char* type):
      Mstd(std), Mmean(mean)
    {
      if (std::string(type) == "default")
      {
        gsl_rng_env_setup();
        const gsl_rng_type* T=gsl_rng_default;
        MR=gsl_rng_alloc(T);
        TFXX_assert(MR!=NULL, "gsl_rng_alloc(T) failed to initialize RNG"); 
      }
      else
      {
        const gsl_rng_type* T=NULL;
        const gsl_rng_type **t, **t0;
        t0 = gsl_rng_types_setup ();

        for (t = t0; *t != 0; t++)
        {
          if (std::string((*t)->name)==std::string(type)) { T=(*t); }
        }
        TFXX_Xassert(T!=NULL, type, UTException);
        MR=gsl_rng_alloc(T);
        TFXX_assert(MR!=NULL, "gsl_rng_alloc(T) failed to initialize RNG"); 
      }
      this->set(gsl_rng_default_seed);
    }

  /* ---------------------------------------------------------------------- */

    // destructor
    RNGgaussian::~RNGgaussian()
    {
      gsl_rng_free(MR);
    }

  /* ---------------------------------------------------------------------- */

    // return value
    double RNGgaussian::value() const
    { return(Mmean+Mstd*gsl_ran_ugaussian(MR)); }

  /* ---------------------------------------------------------------------- */

    // return value
    double RNGgaussian::operator()() const
    { return(this->value()); }

  /* ---------------------------------------------------------------------- */

    // set seed
    void RNGgaussian::set(const unsigned long int& seed) 
    { 
      gsl_rng_set(MR, seed); 
      Mseed=seed;
    }

  /* ---------------------------------------------------------------------- */

    // set seed using time
    void RNGgaussian::set() 
    { this->set(time(0)); }

  /* ---------------------------------------------------------------------- */

    void RNGgaussian::rng_list_types(std::ostream& os)
    {
      const gsl_rng_type **t, **t0;

      t0 = gsl_rng_types_setup ();

      for (t = t0; *t != 0; t++)
      {
        os << std::string((*t)->name) << std::endl;
      }

    } // RNGgaussian::rng_list_types(std::ostream* os)

  /* ====================================================================== */

    RNGgaussian::UTException::UTException(const char* type, 
                                          const char* file,
                                          const int& line,
                                          const char* condition):
      TBase("unkown RNG type requested", file, line, condition),
      Mtype(type) 
    {
      if (this->report_on_construct_is_true()) 
      { 
        std::cerr << "  Requested RNG type: " << Mtype << std::endl;
      }
    } // RNGgaussian::UTException::UTException

  /* ====================================================================== */

    const char* RNGgaussian::comment_gsl_rng_ranlux=
    {
"The recommended random number generator type is \"ranlxd2\"\n"
"\n"
"Quotation from \n"
"https://www.gnu.org/software/gsl/doc/html/index.html\n"
"https://www.gnu.org/software/gsl/doc/html/rng.html#random-number-generator-algorithms\n"
"\n"
"The following generators are recommended for use in simulation. They have\n"
"extremely long periods, low correlation and pass most statistical tests.\n"
"For the most reliable source of uncorrelated numbers, the second-generation\n"
"RANLUX generators have the strongest proof of randomness.\n"
"\n"
"gsl_rng_ranlxd1\n"
"gsl_rng_ranlxd2\n"
"\n"
"    These generators produce double precision output (48 bits) from the\n"
"    RANLXS generator. The library provides two luxury levels ranlxd1 and\n"
"    ranlxd2, in increasing order of strength.\n"
"\n"
"gsl_rng_ranlux\n"
"gsl_rng_ranlux389\n"
"\n"
"    The ranlux generator is an implementation of the original algorithm\n"
"    developed by Luscher. It uses a lagged-fibonacci-with-skipping\n"
"    algorithm to produce “luxury random numbers”. It is a 24-bit\n"
"    generator, originally designed for single-precision IEEE floating\n"
"    point numbers. This implementation is based on integer arithmetic,\n"
"    while the second-generation versions RANLXS and RANLXD described\n"
"    above provide floating-point implementations which will be faster\n"
"    on many platforms. The period of the generator is about 10^{171}.\n"
"    The algorithm has mathematically proven properties and it can\n"
"    provide truly decorrelated numbers at a known level of randomness.\n"
"    The default level of decorrelation recommended by Luscher is\n"
"    provided by gsl_rng_ranlux, while gsl_rng_ranlux389 gives the\n"
"    highest level of randomness, with all 24 bits decorrelated. Both\n"
"    types of generator use 24 words of state per generator.\n"
"\n"
"    For more information see,\n"
"\n"
"    *  M. Luscher, “A portable high-quality random number generator for\n"
"       lattice field theory calculations”, Computer Physics Communications,\n"
"       79 (1994) 100–110.\n"
"       DOI: 10.1016/0010-4655(94)90232-1.\n"
"    *  F. James, “RANLUX: A Fortran implementation of the high-quality\n"
"       pseudo-random number generator of Luscher”, Computer Physics\n"
"       Communications, 79 (1994) 111–114\n"
"       DOI: 10.1016/0010-4655(94)90233-X\n"
    }; // RNGgaussian::comment_gsl_rng_ranlux

  } // namespace numeric

} // namespace tfxx

/* ----- END OF rng.cc ----- */