Sometimes a value gets negative zero causing a NaN on the sqrt.

Closes #61

1 files changed, 12 insertions, 1 deletions

diff --git a/src/mathfunc.cpp b/src/mathfunc.cpp
index fdd6a58..542093e 100644
--- a/src/mathfunc.cpp
+++ b/src/mathfunc.cpp
@@ -104,11 +104,22 @@ bool is_float(const std::string & s){
 }
 
 double safe_log(const double d) {
-  if (!is_legacy_mode())
+  if (!is_legacy_mode() && !is_check_mode())
     enforce_msg(d > 0.0, (std::string("Trying to take the log of ") + std::to_string(d)).c_str());
   return log(d);
 }
 
+double safe_sqrt(const double d) {
+  double d1 = d;
+  if (fabs(d < 0.001))
+    d1 = fabs(d);
+  if (!is_legacy_mode() && !is_check_mode())
+    enforce_msg(d1 >= 0.0, (std::string("Trying to take the sqrt of ") + std::to_string(d)).c_str());
+  if (d1 < 0.0 )
+    return nan("");
+  return sqrt(d1);
+}
+
 // calculate variance of a vector
 double VectorVar(const gsl_vector *v) {
   double d, m = 0.0, m2 = 0.0;