From 1319031a0fe8a951e4487530d978b68f55773b9a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 29 Jun 2020 09:23:15 -0500 Subject: [PATCH] Prevent potential divide-by-zero in ContinuousTabular::sample --- src/distribution_energy.cpp | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/distribution_energy.cpp b/src/distribution_energy.cpp index 5d08cc4d2c..ba1e9f696c 100644 --- a/src/distribution_energy.cpp +++ b/src/distribution_energy.cpp @@ -238,12 +238,14 @@ double ContinuousTabular::sample(double E, uint64_t* seed) const double E_l_k1 = distribution_[l].e_out[k+1]; double p_l_k1 = distribution_[l].p[k+1]; - double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k); - if (frac == 0.0) { - E_out = E_l_k + (r1 - c_k)/p_l_k; - } else { - E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k + - 2.0*frac*(r1 - c_k))) - p_l_k)/frac; + if (E_l_k != E_l_k1) { + double frac = (p_l_k1 - p_l_k)/(E_l_k1 - E_l_k); + if (frac == 0.0) { + E_out = E_l_k + (r1 - c_k)/p_l_k; + } else { + E_out = E_l_k + (std::sqrt(std::max(0.0, p_l_k*p_l_k + + 2.0*frac*(r1 - c_k))) - p_l_k)/frac; + } } } else { throw std::runtime_error{"Unexpected interpolation for continuous energy "