use same dilogarithm function in FORTRAN and C++ module for 2-loop co…

…rrections

This modification

1) fixes the following compilation error on MacOS when creating
   libDSZ.dynlib:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Undefined symbols for architecture x86_64:

  "himalaya::dilog(double)", referenced from:

      himalaya::mssm_twoloophiggs::delta_ma2_2loop_at_as_mst1_eq_mst2(double, double, double, double, double, double, double, double, double, double, double) in DSZHiggs.cpp.o

      himalaya::mssm_twoloophiggs::delta_mh2_2loop_at_as(double, double, double, double, double, double, double, double, double, double, double, int) in DSZHiggs.cpp.o

ld: symbol(s) not found for architecture x86_64
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

2) ensures that the same precision is used to evaluate the 2-loop
   corrections both in the FORTRAN and in the C++ module

CMakeLists.txt

@@ -93,7 +93,8 @@ set(HIM_LIBSOURCES
 set(DSZ_LIBSOURCES
   ${SOURCE_PATH}/mh2l/DSZHiggs.cpp
   ${SOURCE_PATH}/mh2l/DSZHiggs.f
-  ${SOURCE_PATH}/mh2l/functs.f)
+  ${SOURCE_PATH}/mh2l/functs.f
+  ${SOURCE_PATH}/mh2l/li2.cpp)
 
 # build the himalaya and DSZ library
 add_library(DSZ_static STATIC ${DSZ_LIBSOURCES})

source/include/mh2l/li2.hpp

@@ -0,0 +1,18 @@
+// ====================================================================
+// This file is part of Himalaya.
+//
+// Himalaya is licenced under the GNU General Public License (GNU GPL)
+// version 3.
+// ====================================================================
+
+#include <complex>
+
+namespace himalaya {
+
+/// real dilogarithm from FORTRAN module
+double li2(double);
+
+/// complex dilogarithm from FORTRAN module
+std::complex<double> li2(const std::complex<double>&);
+
+} // namespace himalaya

source/mh2l/DSZHiggs.cpp

@@ -7,8 +7,8 @@
 
 #include "DSZHiggs.hpp"
 #include "DSZHiggs.h"
-#include "dilog.hpp"
 #include "EFTFlags.hpp"
+#include "li2.hpp"
 #include <cmath>
 #include <limits>
 #include <mutex>
@@ -44,7 +44,6 @@ Eigen::Matrix<double, 2, 2> delta_mh2_2loop_at_as_st_0_mst1_eq_mst2(
    using std::log;
    using std::sin;
    using std::pow;
-   using himalaya::dilog;
 
    constexpr double Pi2 = M_PI * M_PI;
    constexpr double Pi4 = M_PI * M_PI * M_PI * M_PI;
@@ -76,8 +75,8 @@ Eigen::Matrix<double, 2, 2> delta_mh2_2loop_at_as_st_0_mst1_eq_mst2(
       g*(-1 + sqrtabs(del/g2))))))/ g2 + (2*(g + t - T)*(Pi2 -
       3*logabs(t/g)*logabs(T/g) + 6*logabs((g + t - T -
       g*sqrtabs(del/g2))/(2.*g))*logabs((g - t + T -
-      g*sqrtabs(del/g2))/(2.*g)) - 6*dilog((g + t - T -
-      g*sqrtabs(del/g2))/(2.*g)) - 6*dilog((g - t + T -
+      g*sqrtabs(del/g2))/(2.*g)) - 6*li2((g + t - T -
+      g*sqrtabs(del/g2))/(2.*g)) - 6*li2((g - t + T -
       g*sqrtabs(del/g2))/(2.*g))))/ (3.*g2)))/(2.*pow(fabs(del)/g2,1.5))))/
       (8.*Pi4*T);
 
@@ -89,17 +88,17 @@ Eigen::Matrix<double, 2, 2> delta_mh2_2loop_at_as_st_0_mst1_eq_mst2(
       + T) + 2*g*t*(t + 5*T))*logabs(t/g) +
       4*g2*T*logabs(T/g)))/(del*T) + 2*logabs(T/q) -
       8*sqr(logabs(T/q)) + 5*logabs(Tsq/tsq) + sqr(logabs(Tsq/tsq)) +
-      (4*g2*(g + t - T)*(Pi2 - 6*dilog((g + t - T -
-      g*sqrtabs(del/g2))/ (2.*g)) - 6*dilog((g - t + T -
+      (4*g2*(g + t - T)*(Pi2 - 6*li2((g + t - T -
+      g*sqrtabs(del/g2))/ (2.*g)) - 6*li2((g - t + T -
       g*sqrtabs(del/g2))/(2.*g)) - 3*logabs(t/g)*logabs(T/g) +
       6*logabs((-rdel + g + t - T)/(2.*g))* logabs((-rdel + g - t +
       T)/(2.*g))))/(3.*pow(fabs(del),1.5)) + (4*g*(g + t - T)*(Pi2 -
-      6*dilog((g + t - T - g*sqrtabs(del/g2))/(2.*g)) - 6*dilog((g - t
+      6*li2((g + t - T - g*sqrtabs(del/g2))/(2.*g)) - 6*li2((g - t
       + T - g*sqrtabs(del/g2))/ (2.*g)) - 3*logabs(t/g)*logabs(T/g) +
       6*logabs((g + t - T - g*sqrtabs(del/g2))/(2.*g))* logabs((g - t
       + T - g*sqrtabs(del/g2))/(2.*g))))/ (3.*del*sqrtabs(del/g2)) +
       (8*mg*mu*(1/T - logabs(t/q)/T + (g*(g + t - T)* (Pi2 -
-      6*dilog((g + t - T - g*sqrtabs(del/g2))/(2.*g)) - 6*dilog((g - t
+      6*li2((g + t - T - g*sqrtabs(del/g2))/(2.*g)) - 6*li2((g - t
       + T - g*sqrtabs(del/g2))/ (2.*g)) - 3*logabs(t/g)*logabs(T/g) +
       6*logabs((-rdel + g + t - T)/(2.*g))*logabs((-rdel + g - t +
       T)/(2.*g))))/ (3.*pow(fabs(del),1.5)) + (g*(-(logabs(t/g)/T) -
@@ -167,7 +166,7 @@ double phi(double x, double y, double z)
    const double xm = 0.5 * (1 - (u - v) - lambda);
 
    return 1./lambda * (2*logabs(xp)*logabs(xm) - logabs(u)*logabs(v) -
-                       2*(dilog(xp) + dilog(xm)) + M_PI*M_PI/3.);
+                       2*(li2(xp) + li2(xm)) + M_PI*M_PI/3.);
 }
 
 /// First derivative of phi[t,T,g] w.r.t. T
@@ -183,8 +182,8 @@ double dphi_010(double t, double T, double g)
    const double abbr = (-4*t*T)/g2 + sqr(1 - t/g - T/g);
    const double rabbr = sqrtabs(abbr);
 
-   return ((g + t - T)*(Pi2 - 6*dilog((g - rabbr*g + t - T)/(2.*g)) -
-      6*dilog((g - rabbr*g - t + T)/(2.*g)) -
+   return ((g + t - T)*(Pi2 - 6*li2((g - rabbr*g + t - T)/(2.*g)) -
+      6*li2((g - rabbr*g - t + T)/(2.*g)) -
       3*logabs(t/g)*logabs(T/g) + 6*logabs((g - rabbr*g + t -
       T)/(2.*g))*logabs((g - rabbr*g - t + T)/(2.*g))) + (3*rabbr*g* (
       rabbr*g*((-1 + rabbr)*g + t - T)*logabs(t/g) +

source/mh2l/functs.f

@@ -105,3 +105,26 @@ FUNCTION pCSPEN(ZZ)
      *               +.5D0*CDLOG(Z-1D0)**2-CDLOG(Z)*CDLOG(1D0-Z)
 c        write(*,*) 'cspen:', cspen
       END
+
+*
+***********************************************************************
+*
+
+!>    @brief C wrapper for complex dilogarithm
+      subroutine li2c(re_in, im_in, re_out, im_out)
+     $ bind(C, name="li2c_")
+
+      implicit none
+
+      double precision, intent(in) :: re_in, im_in
+      double precision, intent(out) :: re_out, im_out
+      complex*16 z, l, pCSPEN
+
+      z = dcmplx(re_in, im_in)
+      l = pCSPEN(z)
+
+      re_out = dble(l)
+      im_out = dimag(l)
+
+      return
+      end

source/mh2l/li2.cpp

@@ -0,0 +1,36 @@
+// ====================================================================
+// This file is part of Himalaya.
+//
+// Himalaya is licenced under the GNU General Public License (GNU GPL)
+// version 3.
+// ====================================================================
+
+#include "li2.hpp"
+#include <complex>
+
+/// complex dilogarithm from FORTRAN module
+extern "C" int li2c_(double*, double*, double*, double*);
+
+namespace himalaya {
+
+/// real dilogarithm from FORTRAN module
+double li2(double x)
+{
+   const std::complex<double> z(x, 0.);
+   return std::real(li2(z));
+}
+
+/// complex dilogarithm from FORTRAN module
+std::complex<double> li2(const std::complex<double>& z)
+{
+   double re_in = z.real();
+   double im_in = z.imag();
+   double re_out = 0.;
+   double im_out = 0.;
+
+   li2c_(&re_in, &im_in, &re_out, &im_out);
+
+   return std::complex<double>(re_out, im_out);
+}
+
+} // namespace himalaya

test/CMakeLists.txt

@@ -19,6 +19,7 @@ add_himalaya_test(test_EFT_expressions)
 add_himalaya_test(test_FO_expressions)
 add_himalaya_test(test_FO_vs_EFT_expressions)
 add_himalaya_test(test_Himalaya_interface)
+add_himalaya_test(test_li2)
 add_himalaya_test(test_delta_lambda)
 
 if(Mathematica_FOUND)