From d583081da9684da0d4bc5599218a13d7538f7f7b Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Wed, 6 Aug 2025 11:19:40 +0200 Subject: [PATCH] adding test and benchmark of Li2 against Pythia 8.315 --- test/alt/CMakeLists.txt | 2 ++ test/alt/alt.h | 2 ++ test/alt/pythia/CMakeLists.txt | 5 +++ test/alt/pythia/pythia.c | 65 ++++++++++++++++++++++++++++++++++ test/bench_Li.cpp | 6 ++++ test/test_Li2.cpp | 6 ++++ 6 files changed, 86 insertions(+) create mode 100644 test/alt/pythia/CMakeLists.txt create mode 100644 test/alt/pythia/pythia.c diff --git a/test/alt/CMakeLists.txt b/test/alt/CMakeLists.txt index ab70a6a5..39391477 100644 --- a/test/alt/CMakeLists.txt +++ b/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 diff --git a/test/alt/alt.h b/test/alt/alt.h index 35468c84..e031265f 100644 --- a/test/alt/alt.h +++ b/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); diff --git a/test/alt/pythia/CMakeLists.txt b/test/alt/pythia/CMakeLists.txt new file mode 100644 index 00000000..3554dbb4 --- /dev/null +++ b/test/alt/pythia/CMakeLists.txt @@ -0,0 +1,5 @@ +add_library(pythia + pythia.c +) + +target_include_directories(pythia PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}) diff --git a/test/alt/pythia/pythia.c b/test/alt/pythia/pythia.c new file mode 100644 index 00000000..0dca1713 --- /dev/null +++ b/test/alt/pythia/pythia.c @@ -0,0 +1,65 @@ +#include + +#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; +} diff --git a/test/bench_Li.cpp b/test/bench_Li.cpp index 3c602085..03d024ec 100644 --- a/test/bench_Li.cpp +++ b/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::max(), std::numeric_limits::epsilon()); }, values_d, + "pythia(max)", "double"); + bench_fn([&](long double x) { return polylogarithm::Li2(x); }, values_l, "polylogarithm C++", "long double"); diff --git a/test/test_Li2.cpp b/test/test_Li2.cpp index 3b0347d1..5c899888 100644 --- a/test/test_Li2.cpp +++ b/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::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::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); } }