Skip to content

Commit

Permalink
Replace Marsaglia polar method with Box-muller to generate a normally…
Browse files Browse the repository at this point in the history
… distributed random number (kokkos#6556)

* Kokkos Random: Replace Marsaglia polar method with Box-muller to generate a normally distributed random number

* Apply clang-formatting

* Add const qualifier to some internal variables

* Update Kokkos_Random.hpp
  • Loading branch information
Shihab-Shahriar committed Nov 9, 2023
1 parent 91ee4e1 commit 0a83695
Showing 1 changed file with 16 additions and 18 deletions.
34 changes: 16 additions & 18 deletions algorithms/src/Kokkos_Random.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -849,18 +849,17 @@ class Random_XorShift64 {
return drand(end - start) + start;
}

// Marsaglia polar method for drawing a standard normal distributed random
// Box-muller method for drawing a standard normal distributed random
// number
KOKKOS_INLINE_FUNCTION
double normal() {
double S = 2.0;
double U;
while (S >= 1.0) {
U = 2.0 * drand() - 1.0;
const double V = 2.0 * drand() - 1.0;
S = U * U + V * V;
}
return U * std::sqrt(-2.0 * std::log(S) / S);
constexpr auto two_pi = 2 * Kokkos::numbers::pi_v<double>;

const double u = drand();
const double v = drand();
const double r = Kokkos::sqrt(-2.0 * Kokkos::log(u));
const double theta = v * two_pi;
return r * Kokkos::cos(theta);
}

KOKKOS_INLINE_FUNCTION
Expand Down Expand Up @@ -1094,18 +1093,17 @@ class Random_XorShift1024 {
return drand(end - start) + start;
}

// Marsaglia polar method for drawing a standard normal distributed random
// Box-muller method for drawing a standard normal distributed random
// number
KOKKOS_INLINE_FUNCTION
double normal() {
double S = 2.0;
double U;
while (S >= 1.0) {
U = 2.0 * drand() - 1.0;
const double V = 2.0 * drand() - 1.0;
S = U * U + V * V;
}
return U * std::sqrt(-2.0 * std::log(S) / S);
constexpr auto two_pi = 2 * Kokkos::numbers::pi_v<double>;

const double u = drand();
const double v = drand();
const double r = Kokkos::sqrt(-2.0 * Kokkos::log(u));
const double theta = v * two_pi;
return r * Kokkos::cos(theta);
}

KOKKOS_INLINE_FUNCTION
Expand Down

0 comments on commit 0a83695

Please sign in to comment.