File randomgaussian.h

File List > code_source > templat > include > TempLat > util > random > randomgaussian.h

Go to the documentation of this file

#ifndef TEMPLAT_UTIL_RANDOM_RANDOMGAUSSIAN_H
#define TEMPLAT_UTIL_RANDOM_RANDOMGAUSSIAN_H

/* This file is part of TempLat, available at https://cosmolattice.github.io/templat .
   Copyright 2021-2026 The TempLat authors, see AUTHORS.md.
   Released under the MIT license, see LICENSE.md. */

// File info: Main contributor(s): Wessel Valkenburg, Franz R. Sattler, Year: 2025

#include <iomanip>
#include <sstream>

#include "TempLat/util/random/randomuniform.h"

namespace TempLat
{
  MakeException(RandomGaussianWrongCallOrderException);
  template <typename T> class RandomGaussian
  {
    using INT = typename RandomUniform<T>::IntegerType;

  public:
    // Put public methods here. These should change very little over time.

    using IntegerType = INT;

    RandomGaussian(const std::string &seed) : mRandomUniform(seed) {}

    auto getSeed() const { return mRandomUniform.getSeed(); }
    auto getSeedString() const { return mRandomUniform.getSeedString(); }

    DEVICE_INLINE_FUNCTION
    auto getPair(INT r, INT c, INT g, bool real = false, bool unitary = false) const
    { // Even if this is not completely consistent with the name, it is convenient to be able to use this class to
      // generate numbers with a real gaussian distribution or uniformly on the unit disk.
      return getNextGaussianPair(r, c, g, real, unitary);
    }

    DEVICE_INLINE_FUNCTION
    T get(INT r, INT c, INT g) const { return getPair(r, c, g, false, false)[0]; }

    friend std::ostream &operator<<(std::ostream &ostream, const RandomGaussian &pr)
    {
      ostream << "RandomGaussian - seed string: \"" << pr.mRandomUniform.getSeedString()
              << "\" - seed value: " << pr.mRandomUniform.getSeed();
      return ostream;
    }

    std::string saveState() const { return mRandomUniform.saveState(); }

    void loadState(const std::string &state) { mRandomUniform.loadState(state); }

  private:
    /* Put all member variables and private methods here. These may change arbitrarily. */
    RandomUniform<T> mRandomUniform;

    static constexpr T cTwoPi =
        6.2831853071795864769252867665590057683943387987502116419498891846156328125724179972560696506842341359642961730265646132941876892;

    DEVICE_INLINE_FUNCTION
    device::array<T, 2u> getNextGaussianPair(INT r, INT c, INT g, bool real = false, bool unitary = false) const
    { // Even if this is not completely consistent with the name, it is convenient to be able to use this class to
      // generate numbers with a real gaussian distribution or uniformly on the unit disk.

      const auto result = mRandomUniform.getPair(r, c, g);

      const T &r0 = result[0];
      const T &r1 = result[1];

      const T boxMullerR = unitary ? 1 : (r0 == 0 ? std::numeric_limits<T>::max() : device::sqrt(-2 * device::log(r0)));
      const T boxMullerTheta = real ? 0 : cTwoPi * r1;

      return device::array<T, 2u>{boxMullerR * device::cos(boxMullerTheta), boxMullerR * device::sin(boxMullerTheta)};
    }
  };
} // namespace TempLat

#endif