Adding test and benchmark of Li2 against Pythia 8.315

test/alt/CMakeLists.txt

@@ -10,6 +10,7 @@ add_subdirectory(koelbig)
 add_subdirectory(feynhiggs)
 add_subdirectory(morris)
 add_subdirectory(pade)
+add_subdirectory(pythia)
 add_subdirectory(sherpa)
 add_subdirectory(spheno)
 add_subdirectory(tsil)
@@ -35,6 +36,7 @@ target_link_libraries(alt
     feynhiggs
     morris
     pade
+    pythia
     sherpa
     spheno
     tsil

test/alt/alt.h

@@ -25,6 +25,8 @@ void feynhiggs_dilog(long double re, long double im, long double* res_re, long d
 
 double morris_dilog(double x);
 
+double pythia_dilog(const double x, const double kmax, const double xerr);
+
 void hdecay_dilog(double re, double im, double* res_re, double* res_im);
 
 void hollik_dilog(double re, double im, double* res_re, double* res_im);

test/alt/pythia/CMakeLists.txt

@@ -0,0 +1,5 @@
+add_library(pythia
+  pythia.c
+)
+
+target_include_directories(pythia PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})

test/alt/pythia/pythia.c

@@ -0,0 +1,65 @@
+#include <math.h>
+
+#ifndef M_PI
+#define M_PI 3.1415926535897932
+#endif
+
+/*
+  Dilogarithm
+
+  Implementation from Pythia 8.315.
+  Translated to C by Alexander Voigt.
+ */
+double pythia_dilog(const double x, const double kmax, const double xerr)
+{
+   if (x < 0.0) {
+      return 0.5*pythia_dilog(x*x, kmax, xerr) - pythia_dilog(-x, kmax, xerr);
+   }
+
+   if (x <= 0.5) {
+      double sum = x, term = x;
+
+      for (int k = 2; k < kmax; k++) {
+         double rk = (k-1.0)/k;
+         term *= x*rk*rk;
+         sum += term;
+         if (fabs(term/sum) < xerr) return sum;
+      }
+
+      /* cout << "Maximum number of iterations exceeded in pythia_dilog" << endl; */
+      return sum;
+   }
+
+   if (x < 1.0) {
+      return M_PI*M_PI/6.0 - pythia_dilog(1.0 - x, kmax, xerr) - log(x)*log(1.0 - x);
+   }
+
+   if (x == 1.0) {
+      return M_PI*M_PI/6.0;
+   }
+
+   if (x <= 1.01) {
+      const double eps = x - 1.0,
+         lne = log(eps),
+         c0 = M_PI*M_PI/6.0,
+         c1 =  1.0 - lne,
+         c2 = -(1.0 - 2.0*lne)/4.0,
+         c3 = (1.0 - 3.0*lne)/9.0,
+         c4 = -(1.0 - 4.0*lne)/16.0,
+         c5 = (1.0 - 5.0*lne)/25.0,
+         c6 = -(1.0 - 6.0*lne)/36.0,
+         c7 = (1.0 - 7.0*lne)/49.0,
+         c8 = -(1.0 - 8.0*lne)/64.0;
+
+      return c0 + eps*(c1 + eps*(c2 + eps*(c3 + eps*(c4 + eps*(c5 + eps*(c6 + eps*(c7 + eps*c8)))))));
+   }
+
+   double logx = log(x);
+
+   if (x <= 2.0) {
+      return M_PI*M_PI/6.0 + pythia_dilog(1.0 - 1.0/x, kmax, xerr) -
+         logx*(log(1.0 - 1.0/x) + 0.5*logx);
+   }
+
+   return M_PI*M_PI/3.0 - pythia_dilog(1.0/x, kmax, xerr) - 0.5*logx*logx;
+}

test/bench_Li.cpp

@@ -323,6 +323,12 @@ int main() {
    bench_fn([&](double x) { return morris_dilog(x); }, values_d,
             "Morris", "double");
 
+   bench_fn([&](double x) { return pythia_dilog(x, 100.0, 1e-9); }, values_d,
+            "pythia(default)", "double");
+
+   bench_fn([&](double x) { return pythia_dilog(x, std::numeric_limits<double>::max(), std::numeric_limits<double>::epsilon()); }, values_d,
+            "pythia(max)", "double");
+
    bench_fn([&](long double x) { return polylogarithm::Li2(x); }, values_l,
             "polylogarithm C++", "long double");
 

test/test_Li2.cpp

@@ -496,6 +496,7 @@ TEST_CASE("test_real_fixed_values")
          const auto li64_koelbig  = koelbig_dilog(x64);
          const auto li128_koelbig = koelbig_dilogl(x128);
          const auto li64_morris   = morris_dilog(x64);
+         const auto li64_pythia   = pythia_dilog(x64, 100.0, std::numeric_limits<double>::epsilon());
 #ifdef ENABLE_GSL
          const auto li64_gsl      = gsl_Li2(x64);
 #endif
@@ -525,6 +526,7 @@ TEST_CASE("test_real_fixed_values")
          INFO("Li2(64)  real = " << li64_hassani   << " (hassani)");
          INFO("Li2(64)  real = " << li64_koelbig   << " (koelbig)");
          INFO("Li2(64)  real = " << li64_morris    << " (morris)");
+         INFO("Li2(64)  real = " << li64_pythia    << " (pythia(default))");
 #ifdef ENABLE_GSL
          INFO("Li2(64)  real = " << li64_gsl       << " (GSL)");
 #endif
@@ -551,6 +553,7 @@ TEST_CASE("test_real_fixed_values")
          CHECK_CLOSE(li64_hassani , std::real(li64_expected) , 100*eps64);
          CHECK_CLOSE(li64_koelbig , std::real(li64_expected) , 2*eps64);
          CHECK_CLOSE(li64_morris  , std::real(li64_expected) , 2*eps64);
+         CHECK_CLOSE(li64_pythia  , std::real(li64_expected) , 2*eps64);
 #ifdef ENABLE_GSL
          CHECK_CLOSE(li64_gsl     , std::real(li64_expected) , 2*eps64);
 #endif
@@ -680,6 +683,7 @@ TEST_CASE("test_real_random_values")
       const double li2_hassani = hassani_dilog(v);
       const double li2_koelbig = koelbig_dilog(v);
       const double li2_morris = morris_dilog(v);
+      const double li2_pythia = pythia_dilog(v, 100.0, std::numeric_limits<double>::epsilon());
 
       INFO("x = " << v);
       INFO("Li2(64) real = " << li2          << " (polylogarithm C++)");
@@ -699,6 +703,7 @@ TEST_CASE("test_real_random_values")
       INFO("Li2(64) real = " << li2_hassani  << " (Hassani)");
       INFO("Li2(64) real = " << li2_koelbig  << " (Koelbig)");
       INFO("Li2(64) real = " << li2_morris   << " (Morris)");
+      INFO("Li2(64) real = " << li2_pythia   << " (pythia)");
 
       CHECK_CLOSE(li2, li2_c       , eps64);
 #ifdef ENABLE_FORTRAN
@@ -716,6 +721,7 @@ TEST_CASE("test_real_random_values")
       CHECK_CLOSE(li2, li2_hassani , 100*eps64);
       CHECK_CLOSE(li2, li2_koelbig , 2*eps64);
       CHECK_CLOSE(li2, li2_morris  , eps64);
+      CHECK_CLOSE(li2, li2_pythia  , 2*eps64);
    }
 }