preven nans in d0 function when m^2 < 0

FlexibleSUSY · Mar 26, 2018 · 573a784 · 573a784
1 parent d5ed0b5
commit 573a784
Showing 1 changed file with 16 additions and 20 deletions.
diff --git a/src/pv2.cpp b/src/pv2.cpp
@@ -34,6 +34,7 @@ constexpr double TOL = 1e-4;
 constexpr double sqr(double a) noexcept { return a*a; }
 constexpr double pow3(double a) noexcept { return a*a*a; }
 constexpr double pow6(double a) noexcept { return a*a*a*a*a*a; }
+constexpr double log_abs(double a) noexcept { return std::log(std::abs(a)); }
 
 bool is_close(double m1, double m2, double tol) noexcept
 {
@@ -252,9 +253,6 @@ double h0(double p2, double m12, double m22, double q2) noexcept
 
 double c0(double m12, double m22, double m32) noexcept
 {
-   using std::log;
-   using std::abs;
-
    double ans = 0.;
 
    if (is_close(m12,0.,EPSTOL) && is_close(m22,0.,EPSTOL) && is_close(m32,0.,EPSTOL)) {
@@ -273,75 +271,73 @@ double c0(double m12, double m22, double m32) noexcept
       if (is_close(m22,m32,EPSTOL)) {
          ans = -1./m22;
       } else {
-         ans = log(abs(m32/m22))/(m22 - m32);
+         ans = log_abs(m32/m22)/(m22 - m32);
       }
    } else if (is_close(m22,0.,EPSTOL)) {
       if (is_close(m12,m32,EPSTOL)) {
          ans = -1./m12;
       } else {
-         ans = log(abs(m32/m12))/(m12 - m32);
+         ans = log_abs(m32/m12)/(m12 - m32);
       }
    } else if (is_close(m32,0.,EPSTOL)) {
       if (is_close(m12,m22,EPSTOL)) {
          ans = -1./m12;
       } else {
-         ans = log(abs(m22/m12))/(m12 - m22);
+         ans = log_abs(m22/m12)/(m12 - m22);
       }
    } else if (is_close(m22, m32, EPSTOL)) {
       if (is_close(m12, m22, EPSTOL)) {
          ans = ( - 0.5 / m22 );
       } else {
-         ans = ( m12 / sqr(m12-m22) * log(abs(m22/m12))
+         ans = ( m12 / sqr(m12-m22) * log_abs(m22/m12)
                  + 1.0 / (m12 - m22) );
       }
    } else if (is_close(m12, m22, EPSTOL)) {
-      ans = ( - ( 1.0 + m32 / (m22-m32) * log(abs(m32/m22)) )
+      ans = ( - (1.0 + m32 / (m22-m32) * log_abs(m32/m22))
               / (m22-m32) );
    } else if (is_close(m12, m32, EPSTOL)) {
-      ans = ( - (1.0 + m22 / (m32-m22) * log(abs(m22/m32)))
+      ans = ( - (1.0 + m22 / (m32-m22) * log_abs(m22/m32))
               / (m32-m22) );
    } else {
       ans = (1.0 / (m22 - m32) *
              (m22 / (m12 - m22) *
-              log(abs(m22 / m12)) -
+              log_abs(m22 / m12) -
               m32 / (m12 - m32) *
-              log(abs(m32 / m12))) );
+              log_abs(m32 / m12)) );
    }
 
    return ans;
 }
 
 double d0(double m12, double m22, double m32, double m42) noexcept
 {
-   using std::log;
-
    if (is_close(m12, m22, EPSTOL)) {
       if (is_close(m22,0.,EPSTOL)) {
          // d0 is undefined for m1 == m2 == 0
          return 0.;
       } else if (is_close(m32,0.,EPSTOL)) {
-         return (-m22 + m42 - m22 * log(m42/m22))/
+         return (-m22 + m42 - m22 * log_abs(m42/m22))/
             (m22 * sqr(m22 - m42));
       } else if (is_close(m42,0.,EPSTOL)) {
-         return (-m22 + m32 - m22 * log(m32/m22))/
+         return (-m22 + m32 - m22 * log_abs(m32/m22))/
             (m22 * sqr(m22 - m32));
       } else if (is_close(m22, m32, EPSTOL) && is_close(m22, m42, EPSTOL)) {
          return 1.0 / (6.0 * sqr(m22));
       } else if (is_close(m22, m32, EPSTOL)) {
-         return (sqr(m22) - sqr(m42) + 2.0 * m42 * m22 * log(m42 / m22)) /
+         return (sqr(m22) - sqr(m42) + 2.0 * m42 * m22 * log_abs(m42 / m22)) /
             (2.0 * m22 * sqr(m22 - m42) * (m22 - m42));
       } else if (is_close(m22, m42, EPSTOL)) {
-         return (sqr(m22) - sqr(m32) + 2.0 * m32 * m22 * log(m32 / m22)) /
+         return (sqr(m22) - sqr(m32) + 2.0 * m32 * m22 * log_abs(m32 / m22)) /
             (2.0 * m22 * sqr(m22 - m32) * (m22 - m32));
       } else if (is_close(m32, m42, EPSTOL)) {
          return -1.0 / sqr(m22 - m32) *
-            ((m22 + m32) / (m22 - m32) * log(m32 / m22) + 2.0);
+            ((m22 + m32) / (m22 - m32) * log_abs(m32 / m22) + 2.0);
       }
 
       return
-         (m42 / sqr(m22 - m42) * log(m42 / m22) +
+         (m42 / sqr(m22 - m42) * log_abs(m42 / m22) +
           m42 / (m22 * (m22 - m42)) -
-          m32 / sqr(m22 - m32) * log(m32 / m22) -
+          m32 / sqr(m22 - m32) * log_abs(m32 / m22) -
           m32 / (m22 * (m22 - m32))) / (m32 - m42);
    }
    return (c0(m12, m32, m42) - c0(m22, m32, m42)) / (m12 - m22);