use B1 limit implementation which avoids numeric integration

The limit is triggered for p -> 0 and one mass 0 and the other not.

src/numerics.cpp

@@ -29,13 +29,11 @@
 
 #include "numerics.h"
 #include "numerics2.hpp"
-#include "rk.hpp"
 #ifdef USE_LOOPTOOLS
 #include "clooptools.h"
 #endif
 #include <algorithm>
 #include <cmath>
-#include <Eigen/Dense>
 
 namespace softsusy {
 
@@ -59,16 +57,6 @@ bool is_close(double m1, double m2, double tol) noexcept
   return (mmax - mmin <= max_tol);
 }
 
-double refnfn(double x, double p, double m1, double m2, double q) noexcept
-{
-  using flexiblesusy::fast_log;
-  const static std::complex<double> iEpsilon(0.0, TOL * 1.0e-20);
-
-  return std::real(x *
-    fast_log(((1 - x) * sqr(m1) + x * sqr(m2)
-              - x * (1 - x) * sqr(p) - iEpsilon) / sqr(q)));
-}
-
 /// returns a/b if a/b is finite, otherwise returns numeric_limits::max()
 template <typename T>
 T divide_finite(T a, T b) noexcept {
@@ -78,41 +66,6 @@ T divide_finite(T a, T b) noexcept {
    return result;
 }
 
-double integrandThreshbnr(double x, double p, double m1, double m2, double q) noexcept
-{
-  return refnfn(x, p, m1, m2, q);
-}
-
-Eigen::Array<double,1,1> dd(double x, double p, double m1, double m2, double q) noexcept
-{
-  Eigen::Array<double,1,1> dydx;
-  dydx(0) = -integrandThreshbnr(x, p, m1, m2, q);
-  return dydx;
-}
-
-// Returns real part of integral
-double bIntegral(double p, double m1, double m2, double q) noexcept
-{
-  using namespace flexiblesusy;
-
-  const double from = 0.0, to = 1.0, guess = 0.1, hmin = TOL * 1.0e-5;
-  const double eps = TOL * 1.0e-3;
-  Eigen::Array<double,1,1> v;
-  v(0) = 1.0;
-
-  try {
-     runge_kutta::integrateOdes(v, from, to, eps, guess, hmin,
-                                [p, m1, m2, q] (double x, const Eigen::Array<double,1,1>&) {
-                                   return dd(x, p, m1, m2, q);
-                                });
-  } catch (...) {
-     ERROR("B1 integral did not converge.");
-     v(0) = 1.;
-  }
-
-  return v(0) - 1.0;
-}
-
 double fB(const std::complex<double>& a) noexcept
 {
   using flexiblesusy::fast_log;
@@ -277,7 +230,10 @@ double b1(double p, double m1, double m2, double q) noexcept
           (24.*pow6(m12 - m22)) - 0.5*log(m22/q2);
     }
   } else {
-    ans = bIntegral(p, m1, m2, q);
+     if (m12 > m22)
+        ans = -0.5*log(m12/q2) + 0.75;
+     else
+        ans = -0.5*log(m22/q2) + 0.25;
   }
 
 #ifdef USE_LOOPTOOLS

test/test_pv.cpp

@@ -28,11 +28,63 @@
 
 #include <boost/test/unit_test.hpp>
 #include <boost/test/floating_point_comparison.hpp>
+#include "numerics2.hpp"
+#include "rk.hpp"
+#include <Eigen/Dense>
 
 using namespace std;
 using namespace flexiblesusy;
 using namespace flexiblesusy::passarino_veltman;
 
+constexpr double EPSTOL = 1.0e-11; ///< underflow accuracy
+constexpr double TOL = 1e-4;
+constexpr double sqr(double a) noexcept { return a*a; }
+
+double refnfn(double x, double p, double m1, double m2, double q) noexcept
+{
+  using flexiblesusy::fast_log;
+  const static std::complex<double> iEpsilon(0.0, TOL * 1.0e-20);
+
+  return std::real(x *
+    fast_log(((1 - x) * sqr(m1) + x * sqr(m2)
+              - x * (1 - x) * sqr(p) - iEpsilon) / sqr(q)));
+}
+
+double integrandThreshbnr(double x, double p, double m1, double m2, double q) noexcept
+{
+  return refnfn(x, p, m1, m2, q);
+}
+
+Eigen::Array<double,1,1> dd(double x, double p, double m1, double m2, double q) noexcept
+{
+  Eigen::Array<double,1,1> dydx;
+  dydx(0) = -integrandThreshbnr(x, p, m1, m2, q);
+  return dydx;
+}
+
+// Returns real part of integral
+double bIntegral(double p, double m1, double m2, double q) noexcept
+{
+  using namespace flexiblesusy;
+
+  const double from = 0.0, to = 1.0, guess = 0.1, hmin = TOL * 1.0e-5;
+  const double eps = TOL * 1.0e-3;
+  Eigen::Array<double,1,1> v;
+  v(0) = 1.0;
+
+  try {
+     runge_kutta::integrateOdes(v, from, to, eps, guess, hmin,
+                                [p, m1, m2, q] (double x, const Eigen::Array<double,1,1>&) {
+                                   return dd(x, p, m1, m2, q);
+                                });
+  } catch (...) {
+     ERROR("B1 integral did not converge.");
+     v(0) = 1.;
+  }
+
+  return v(0) - 1.0;
+}
+
 BOOST_AUTO_TEST_CASE(test_real_parts)
 {
     BOOST_CHECK_CLOSE_FRACTION(ReA0 (2,     1),  0.61370563888010943, 1e-14);
@@ -70,6 +122,56 @@ BOOST_AUTO_TEST_CASE( test_ReB1 )
 
    BOOST_CHECK_CLOSE(ReB1(p2, 0,0, scale2), -0.5*ReB0(p2, 0,0, scale2), 1e-10);
    BOOST_CHECK_CLOSE(ReB1(0,p2,p2, scale2), -0.5*ReB0(0,p2,p2, scale2), 1e-10);
+
+   BOOST_CHECK(ReB1(0,1e300,0,1e-10) != 0.);
+}
+
+struct Values {
+   Values(double p_, double m1_, double m2_, double q_)
+      : p(p_), m1(m1_), m2(m2_), q(q_) {}
+   double p{}, m1{}, m2{}, q{};
+};
+
+BOOST_AUTO_TEST_CASE( test_ReB1_integral )
+{
+   const std::vector<Values> vals = {
+      Values(0.  , 1.  , 0., 1.),
+      Values(1e-1, 1.  , 0., 1.),
+      Values(1e-2, 1.  , 0., 1.),
+      Values(1e-3, 1.  , 0., 1.),
+      Values(1e-4, 1.  , 0., 1.),
+      Values(1e-5, 1.  , 0., 1.),
+      Values(0.  , 1e20, 0., 1.),
+      Values(1e-1, 1e20, 0., 1.),
+      Values(1e-2, 1e20, 0., 1.),
+      Values(1e-3, 1e20, 0., 1.),
+      Values(1e-4, 1e20, 0., 1.),
+      Values(1e-5, 1e20, 0., 1.),
+      Values(0.  , 0.  , 1., 1.),
+      Values(1e-1, 0.  , 1., 1.),
+      Values(1e-2, 0.  , 1., 1.),
+      Values(1e-3, 0.  , 1., 1.),
+      Values(1e-4, 0.  , 1., 1.),
+      Values(1e-5, 0.  , 1., 1.),
+      Values(0.  , 1e20, 1., 1.),
+      Values(1e-1, 1e20, 1., 1.),
+      Values(1e-2, 1e20, 1., 1.),
+      Values(1e-3, 1e20, 1., 1.),
+      Values(1e-4, 1e20, 1., 1.),
+      Values(1e-5, 1e20, 1., 1.),
+      Values(1.  , 1.  , 1., 1.)
+   };
+
+   for (const auto v: vals) {
+      const auto p = v.p;
+      const auto m1 = v.m1;
+      const auto m2 = v.m2;
+      const auto q = v.q;
+
+      const auto v1 = softsusy::b1(p,m1,m2,q);
+      const auto v2 = bIntegral(p,m1,m2,q);
+      BOOST_CHECK_CLOSE_FRACTION(v1, v2, 8e-5);
+   }
 }
 
 BOOST_AUTO_TEST_CASE( test_ReB00 )