diff --git a/CMakeLists.txt b/CMakeLists.txt
index f3cafcc5c..3ab04837e 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -303,6 +303,7 @@ list(APPEND libopenmc_SOURCES
   src/plot.cpp
   src/position.cpp
   src/progress_bar.cpp
+  src/random_dist.cpp
   src/random_lcg.cpp
   src/reaction.cpp
   src/reaction_product.cpp
diff --git a/include/openmc/math_functions.h b/include/openmc/math_functions.h
index 79cd06325..0fd322556 100644
--- a/include/openmc/math_functions.h
+++ b/include/openmc/math_functions.h
@@ -8,9 +8,7 @@
 #include <complex>
 #include <cstdlib>
 
-#include "openmc/constants.h"
 #include "openmc/position.h"
-#include "openmc/random_lcg.h"
 
 
 namespace openmc {
@@ -138,71 +136,6 @@ extern "C" void rotate_angle_c(double uvw[3], double mu, const double* phi,
 Direction rotate_angle(Direction u, double mu, const double* phi,
   uint64_t* seed);
 
-//==============================================================================
-//! Samples an energy from the Maxwell fission distribution based on a direct
-//! sampling scheme.
-//!
-//! The probability distribution function for a Maxwellian is given as
-//! p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can be sampled using
-//! rule C64 in the Monte Carlo Sampler LA-9721-MS.
-//!
-//! \param T The tabulated function of the incoming energy
-//! \param seed A pointer to the pseudorandom seed
-//! \return The sampled outgoing energy
-//==============================================================================
-
-extern "C" double maxwell_spectrum(double T, uint64_t* seed);
-
-//==============================================================================
-//! Samples an energy from a Watt energy-dependent fission distribution.
-//!
-//! Although fitted parameters exist for many nuclides, generally the
-//! continuous tabular distributions (LAW 4) should be used in lieu of the Watt
-//! spectrum. This direct sampling scheme is an unpublished scheme based on the
-//! original Watt spectrum derivation (See F. Brown's MC lectures).
-//!
-//! \param a Watt parameter a
-//! \param b Watt parameter b
-//! \param seed A pointer to the pseudorandom seed
-//! \return The sampled outgoing energy
-//==============================================================================
-
-extern "C" double watt_spectrum(double a, double b, uint64_t* seed);
-
-//==============================================================================
-//! Samples an energy from the Gaussian energy-dependent fission distribution.
-//!
-//! Samples from a Normal distribution with a given mean and standard deviation
-//! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2
-//! Its sampled according to
-//! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf
-//! section 33.4.4
-//!
-//! @param mean mean of the Gaussian distribution
-//! @param std_dev standard deviation of the Gaussian distribution
-//! @param seed A pointer to the pseudorandom seed
-//! @result The sampled outgoing energy
-//==============================================================================
-
-extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed);
-
-//==============================================================================
-//! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
-//!
-//! This is another form of the Gaussian distribution but with more easily
-//! modifiable parameters
-//! https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
-//!
-//! @param e0 peak neutron energy [eV]
-//! @param m_rat ratio of the fusion reactants to AMU
-//! @param kt the ion temperature of the reactants [eV]
-//! @param seed A pointer to the pseudorandom seed
-//! @result The sampled outgoing energy
-//==============================================================================
-
-extern "C" double muir_spectrum(double e0, double m_rat, double kt,
-  uint64_t* seed);
-
 //==============================================================================
 //! Constructs a natural cubic spline.
 //!
diff --git a/include/openmc/random_dist.h b/include/openmc/random_dist.h
new file mode 100644
index 000000000..9012cd747
--- /dev/null
+++ b/include/openmc/random_dist.h
@@ -0,0 +1,74 @@
+#ifndef OPENMC_RANDOM_DIST_H
+#define OPENMC_RANDOM_DIST_H
+
+#include <cstdint> // for uint64_t
+
+namespace openmc {
+
+//==============================================================================
+//! Samples an energy from the Maxwell fission distribution based on a direct
+//! sampling scheme.
+//!
+//! The probability distribution function for a Maxwellian is given as
+//! p(x) = 2/(T*sqrt(pi))*sqrt(x/T)*exp(-x/T). This PDF can be sampled using
+//! rule C64 in the Monte Carlo Sampler LA-9721-MS.
+//!
+//! \param T The tabulated function of the incoming energy
+//! \param seed A pointer to the pseudorandom seed
+//! \return The sampled outgoing energy
+//==============================================================================
+
+extern "C" double maxwell_spectrum(double T, uint64_t* seed);
+
+//==============================================================================
+//! Samples an energy from a Watt energy-dependent fission distribution.
+//!
+//! Although fitted parameters exist for many nuclides, generally the
+//! continuous tabular distributions (LAW 4) should be used in lieu of the Watt
+//! spectrum. This direct sampling scheme is an unpublished scheme based on the
+//! original Watt spectrum derivation (See F. Brown's MC lectures).
+//!
+//! \param a Watt parameter a
+//! \param b Watt parameter b
+//! \param seed A pointer to the pseudorandom seed
+//! \return The sampled outgoing energy
+//==============================================================================
+
+extern "C" double watt_spectrum(double a, double b, uint64_t* seed);
+
+//==============================================================================
+//! Samples an energy from the Gaussian energy-dependent fission distribution.
+//!
+//! Samples from a Normal distribution with a given mean and standard deviation
+//! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2
+//! Its sampled according to
+//! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf
+//! section 33.4.4
+//!
+//! \param mean mean of the Gaussian distribution
+//! \param std_dev standard deviation of the Gaussian distribution
+//! \param seed A pointer to the pseudorandom seed
+//! \result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed);
+
+//==============================================================================
+//! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
+//!
+//! This is another form of the Gaussian distribution but with more easily
+//! modifiable parameters
+//! https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
+//!
+//! \param e0 peak neutron energy [eV]
+//! \param m_rat ratio of the fusion reactants to AMU
+//! \param kt the ion temperature of the reactants [eV]
+//! \param seed A pointer to the pseudorandom seed
+//! \result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed);
+
+} // namespace openmc
+
+#endif // OPENMC_RANDOM_DIST_H
diff --git a/src/distribution.cpp b/src/distribution.cpp
index 235bed231..68318c992 100644
--- a/src/distribution.cpp
+++ b/src/distribution.cpp
@@ -9,6 +9,7 @@
 
 #include "openmc/error.h"
 #include "openmc/math_functions.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/xml_interface.h"
 
diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp
index 7b32e9e4f..338486804 100644
--- a/src/distribution_energy.cpp
+++ b/src/distribution_energy.cpp
@@ -9,6 +9,7 @@
 #include "openmc/endf.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 #include "openmc/search.h"
 
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index 449c50467..1dd2726c7 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -2,6 +2,9 @@
 
 #include "Faddeeva.hh"
 
+#include "openmc/constants.h"
+#include "openmc/random_lcg.h"
+
 namespace openmc {
 
 //==============================================================================
@@ -669,49 +672,6 @@ Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* seed
 }
 
 
-double maxwell_spectrum(double T, uint64_t* seed) {
-  // Set the random numbers
-  double r1 = prn(seed);
-  double r2 = prn(seed);
-  double r3 = prn(seed);
-
-  // determine cosine of pi/2*r
-  double c = std::cos(PI / 2. * r3);
-
-  // Determine outgoing energy
-  double E_out = -T * (std::log(r1) + std::log(r2) * c * c);
-
-  return E_out;
-}
-
-
-double normal_variate(double mean, double standard_deviation, uint64_t* seed) {
-  // Sample a normal variate using Marsaglia's polar method
-  double x, y, r2;
-  do {
-    x = 2 * prn(seed) - 1.;
-    y = 2 * prn(seed) - 1.;
-    r2 = x*x + y*y;
-  } while (r2 > 1 || r2 == 0);
-  double z = std::sqrt(-2.0 * std::log(r2)/r2);
-  return mean + standard_deviation*z*x;
-}
-
-double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed) {
-  // https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
-  double sigma = std::sqrt(4.*e0*kt/m_rat);
-  return normal_variate(e0, sigma, seed);
-}
-
-
-double watt_spectrum(double a, double b, uint64_t* seed) {
-  double w = maxwell_spectrum(a, seed);
-  double E_out = w + 0.25 * a * a * b + (2. * prn(seed) - 1.) * std::sqrt(a * a * b * w);
-
-  return E_out;
-}
-
-
 void spline(int n, const double x[], const double y[], double z[])
 {
   std::vector<double> c_new(n-1);
diff --git a/src/random_dist.cpp b/src/random_dist.cpp
new file mode 100644
index 000000000..fe15131b2
--- /dev/null
+++ b/src/random_dist.cpp
@@ -0,0 +1,48 @@
+#include "openmc/random_dist.h"
+
+#include <cmath>
+
+#include "openmc/constants.h"
+#include "openmc/random_lcg.h"
+
+namespace openmc {
+
+double maxwell_spectrum(double T, uint64_t* seed) {
+  // Set the random numbers
+  double r1 = prn(seed);
+  double r2 = prn(seed);
+  double r3 = prn(seed);
+
+  // determine cosine of pi/2*r
+  double c = std::cos(PI / 2. * r3);
+
+  // Determine outgoing energy
+  return -T * (std::log(r1) + std::log(r2) * c * c);
+}
+
+double watt_spectrum(double a, double b, uint64_t* seed) {
+  double w = maxwell_spectrum(a, seed);
+  return w + 0.25 * a * a * b + (2. * prn(seed) - 1.) * std::sqrt(a * a * b * w);
+}
+
+
+double normal_variate(double mean, double standard_deviation, uint64_t* seed) {
+  // Sample a normal variate using Marsaglia's polar method
+  double x, y, r2;
+  do {
+    x = 2 * prn(seed) - 1.;
+    y = 2 * prn(seed) - 1.;
+    r2 = x*x + y*y;
+  } while (r2 > 1 || r2 == 0);
+  double z = std::sqrt(-2.0 * std::log(r2)/r2);
+  return mean + standard_deviation*z*x;
+}
+
+
+double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed) {
+  // https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
+  double sigma = std::sqrt(4.*e0*kt/m_rat);
+  return normal_variate(e0, sigma, seed);
+}
+
+} // namespace openmc
diff --git a/src/random_lcg.cpp b/src/random_lcg.cpp
index 707117d97..c3bd0d5ad 100644
--- a/src/random_lcg.cpp
+++ b/src/random_lcg.cpp
@@ -39,9 +39,7 @@ double prn(uint64_t* seed)
   uint64_t result = (word >> 43u) ^ word;
 
   // Convert output from unsigned integer to double
-  double result_d = ldexp(result, -64);
-
-  return result_d;
+  return ldexp(result, -64);
 }
 
 //==============================================================================
@@ -51,7 +49,7 @@ double prn(uint64_t* seed)
 double future_prn(int64_t n, uint64_t seed)
 {
   uint64_t fseed = future_seed(static_cast<uint64_t>(n), seed);
-  return prn( &fseed );
+  return prn(&fseed);
 }
 
 //==============================================================================
diff --git a/src/secondary_nbody.cpp b/src/secondary_nbody.cpp
index 8f6ba80e2..2220f37fb 100644
--- a/src/secondary_nbody.cpp
+++ b/src/secondary_nbody.cpp
@@ -5,6 +5,7 @@
 #include "openmc/constants.h"
 #include "openmc/hdf5_interface.h"
 #include "openmc/math_functions.h"
+#include "openmc/random_dist.h"
 #include "openmc/random_lcg.h"
 
 namespace openmc {