unit_normal_dist.hpp

Go to the documentation of this file.

/**
 * @file unit_normal_dist.hpp
 * @author Charles Karney <charles.karney@sri.com>
 * @brief Definition of unit_normal_dist
 *
 * Copyright (c) Charles Karney (2014-2020) and licensed under the MIT/X11
 * License.  For more information, see http://exrandom.sourceforge.net/
 */

#if !defined(EXRANDOM_UNIT_NORMAL_DIST_HPP)
#define EXRANDOM_UNIT_NORMAL_DIST_HPP 1

#include <algorithm>            // for std::min, std::max

#include <exrandom/u_rand.hpp>

namespace exrandom {

  /**
   * @brief Sample u-rands exactly from the unit normal distribution.
   *
   * This samples from the unit normal distribution
   * P(x) = exp(&minus; x<sup>2</sup>/2) / &radic;(2&pi;).
   * This implements Algorithm N with improvements due to Du et al. (2020).
   *
   * @tparam digit_gen the type of digit generator.
   *
   * This class allows a u-rand to be returned via the
   * unit_normal_dist::generate member function or a floating point result via
   * the unit_normal_dist::value<RealType> member function.
   *
   * See unit_normal_distribution for a simpler interface to sampling normal
   * deviates.  (But, you can't obtain u-rands with this class and the
   * constructor locks you into a specific floating point type.)
   */
  template<typename digit_gen> class unit_normal_dist {
  public:
    /**
     * The constructor.
     *
     * @param D a reference to the digit generator to be used.
     */
    unit_normal_dist(digit_gen& D) : _D(D), _y(_D), _z(_D), _x(_D) {}

    /**
     * Generate the next deviate as a u_rand.
     *
     * @tparam Generator the type of g.
     * @param g the random generator engine.
     * @param[out] x the u_rand to set.
     */
    template<typename Generator>
    void generate(Generator& g, u_rand<digit_gen>& x) {
      for (;;) {
        int k = G(g);                                // step 1
        k = S(k); if (k < 0) continue;               // step 2
        x.init();                                    // step 3
        int j = k; while (j-- && E(g, x)) {};        // step 4
        if (!(j < 0)) continue;
        if (!B(g, x)) continue;                      // step 5
        x.set_integer(k);                            // step 6
        if (_y.init().less_than_half(g)) x.negate(); // step 7
        return;                                      // step 8
      }
    }

    /**
     * Generate the next deviate and round it to a floating point number.
     *
     * @tparam RealType the floating point type of the result.
     * @tparam Generator the type of g.
     * @param g the random generator engine.
     * @return the normal deviate.
     */
    template<typename RealType, typename Generator>
    RealType value(Generator& g) {
      generate(g, _x);
      return _x.template value<RealType>(g);
    }

    /**
     * @return a reference to the digit generator used in the constructor.
     */
    digit_gen& digit_generator() const { return _D; }

    /**
     * The base of the digit generator.
     */
    static const uint_t base = digit_gen::base;
  private:
    static const int bits = digit_gen::bits;
    static const bool power_of_two = digit_gen::power_of_two;
    // Disable copy assignment
    unit_normal_dist& operator=(const unit_normal_dist&);
    digit_gen& _D;
    u_rand<digit_gen> _y, _z, _x; // temporary storage

    // Algorithm H: true with probability exp(-1/2).
    template<typename Generator>
    bool H(Generator& g) {
      if (!_y.init().less_than_half(g)) return true;
      for (;;) {
        if (!_z.init().less_than(g, _y)) return false;
        if (!_y.init().less_than(g, _z)) return true;
      }
    }

    // Step N1: return n >= 0 with prob. exp(-n/2) * (1 - exp(-1/2)).
    template<typename Generator>
    int G(Generator& g)
    { int n = 0; while (H(g)) ++n; return n; }

    // Step N2: return the square root of n >= 0 if it's a perfect square else -1
    int S(int n) const {
      for (int k = 0, k2 = 0; k2 <= n; ++k, k2 += 2*k - 1) {
        // Here k2 = k * k; note that k^2 - (k - 1)^2 = 2*k - 1
        if (n == k2) return k;
      }
      return -1;
    }

    // Algorithm E: true with probability exp(-x) for x in (0, 1).
    template<typename Generator>
    bool E(Generator& g, u_rand<digit_gen>& x) {
      if (!_y.init().less_than(g, x)) return true;
      for (;;) {
        if (!_z.init().less_than(g, _y)) return false;
        if (!_y.init().less_than(g, _z)) return true;
      }
    }

    // Algorithm B: true with prob exp(-x^2 / 2).
    template<typename Generator>
    bool B(Generator& g, u_rand<digit_gen>& x) {
      int n = 0;
      for (;; ++n) {
        if (_z.init().less_than_half(g)) break;
        if (!_z.init().less_than(g, n ? _y : x)) break;
        if (!_y.init().less_than(g, x)) break;
        _y.swap(_z);            // an efficient way of doing y = z
      }
      return (n % 2) == 0;
    }

  };

}

#endif  // EXRANDOM_UNIT_NORMAL_DIST_HPP