diff --git a/test/alt/algorithm_490/algorithm_490.c b/test/alt/algorithm_490/algorithm_490.c index 8819246f..00562001 100644 --- a/test/alt/algorithm_490/algorithm_490.c +++ b/test/alt/algorithm_490/algorithm_490.c @@ -119,3 +119,82 @@ double algorithm_490(double X) return DILOG; } + + +/* + * Refined version of algorithm 490 without goto statements and with + * one division less. + * + * @author Refined by Alexander Voigt + */ +double algorithm_490_2(double x) +{ + const double PI = 3.1415926535897932; + // Table[1/(n(n+1)(n+2))^2, {n,2,25}] // N[#,17]& // CForm + const double C[24] = { + 0.0017361111111111111, 0.00027777777777777778, + 0.000069444444444444444, 0.000022675736961451247, + 8.8577097505668934e-6, 3.9367598891408415e-6, + 1.9290123456790123e-6, 1.0203040506070809e-6, + 5.7392102846648301e-7, 3.3959824169614379e-7, + 2.0964993492466020e-7, 1.3417595835178253e-7, + 8.8577097505668934e-8, 6.0073048827374087e-8, + 4.1717395019009783e-8, 2.9583526661680067e-8, + 2.1374098013063849e-8, 1.5703418948373440e-8, + 1.1712674050336388e-8, 8.8564643102732612e-9, + 6.7807304875529656e-9, 5.2509976895610166e-9, + 4.1091387245233399e-9, 3.2467268934505402e-9 + }; + + double y = 0, r = 0, s = 1; + + /* transform to [0, 1/2] */ + if (x < -1) { + const double l = log(1 - x); + y = 1/(1 - x); + r = -PI*PI/6 + l*(0.5*l - log(-x)); + s = 1; + } else if (x == -1) { + return -PI*PI/12; + } else if (x < 0) { + const double l = log1p(-x); + y = x/(x - 1); + r = -0.5*l*l; + s = -1; + } else if (x == 0) { + return x; + } else if (x < 0.5) { + y = x; + r = 0; + s = 1; + } else if (x < 1) { + y = 1 - x; + r = PI*PI/6 - log(x)*log1p(-x); + s = -1; + } else if (x == 1) { + return PI*PI/6; + } else if (x < 2) { + const double l = log(x); + y = 1 - 1/x; + r = PI*PI/6 - l*(log(y) + 0.5*l); + s = 1; + } else { + const double l = log(x); + y = 1/x; + r = PI*PI/3 - 0.5*l*l; + s = -1; + } + + double yn = 4*y*y*y; + double sum = yn/36; + + for (int n = 0; n < 24; ++n) { + const int d = (n + 2)*(n + 3)*(n + 4); + yn *= y; + sum += yn/(d*d); + } + + const double d = 1 + y*(4 + y); + + return r + s*(sum + y*(4 + 5.75*y) + 3*(1 + y)*(1 - y)*log(1 - y))/d; +} diff --git a/test/alt/alt.h b/test/alt/alt.h index dd4e5aed..35468c84 100644 --- a/test/alt/alt.h +++ b/test/alt/alt.h @@ -7,6 +7,7 @@ extern "C" { double algorithm_327(double x); double algorithm_490(double x); +double algorithm_490_2(double x); double babar_dilog(double x); diff --git a/test/bench_Li.cpp b/test/bench_Li.cpp index 8bc3dd8e..3c602085 100644 --- a/test/bench_Li.cpp +++ b/test/bench_Li.cpp @@ -311,6 +311,9 @@ int main() { bench_fn([&](double x) { return algorithm_490(x); }, values_d, "algorithm 490", "double"); + bench_fn([&](double x) { return algorithm_490_2(x); }, values_d, + "algorithm 490 2", "double"); + bench_fn([&](double x) { return hassani_dilog(x); }, values_d, "Hassani", "double"); diff --git a/test/test_Li2.cpp b/test/test_Li2.cpp index ff4b3590..3b0347d1 100644 --- a/test/test_Li2.cpp +++ b/test/test_Li2.cpp @@ -291,6 +291,7 @@ TEST_CASE("test_special_values") CHECK_CLOSE(Li2(std::real(z0)) , std::real(li0), eps); CHECK_CLOSE(algorithm_327(std::real(z0)), std::real(li0), eps); CHECK_CLOSE(algorithm_490(std::real(z0)), std::real(li0), eps); + CHECK_CLOSE(algorithm_490_2(std::real(z0)), std::real(li0), eps); #ifdef ENABLE_GSL CHECK_CLOSE(gsl_Li2(std::real(z0)) , std::real(li0), eps); #endif @@ -487,6 +488,7 @@ TEST_CASE("test_real_fixed_values") if (std::imag(z128) == 0.0L) { const auto li64_327 = algorithm_327(x64); const auto li64_490 = algorithm_490(x64); + const auto li64_490_2 = algorithm_490_2(x64); const auto li64_babar = babar_dilog(x64); const auto li64_cephes = cephes_dilog(x64); const auto li64_cephes_2 = cephes_dilog_2(x64); @@ -516,6 +518,7 @@ TEST_CASE("test_real_fixed_values") INFO("Li2(64) real = " << li64_expected << " (expected)"); INFO("Li2(64) real = " << li64_327 << " (algorithm 327)"); INFO("Li2(64) real = " << li64_490 << " (algorithm 490)"); + INFO("Li2(64) real = " << li64_490_2 << " (algorithm 490 2)"); INFO("Li2(64) real = " << li64_babar << " (BaBar)"); INFO("Li2(64) real = " << li64_cephes << " (cephes)"); INFO("Li2(64) real = " << li64_cephes_2 << " (cephes 2)"); @@ -541,6 +544,7 @@ TEST_CASE("test_real_fixed_values") CHECK_CLOSE(li32_poly_c , std::real(li32_expected) , 2*eps32); CHECK_CLOSE(li64_327 , std::real(li64_expected) , 10*eps64); CHECK_CLOSE(li64_490 , std::real(li64_expected) , 2*eps64); + CHECK_CLOSE(li64_490_2 , std::real(li64_expected) , 2*eps64); CHECK_CLOSE(li64_babar , std::real(li64_expected) , 100*eps64); CHECK_CLOSE(li64_cephes , std::real(li64_expected) , 2*eps64); CHECK_CLOSE(li64_cephes_2, std::real(li64_expected) , 2*eps64); @@ -669,6 +673,7 @@ TEST_CASE("test_real_random_values") #endif const double li2_327 = algorithm_327(v); const double li2_490 = algorithm_490(v); + const double li2_490_2 = algorithm_490_2(v); const double li2_babar = babar_dilog(v); const double li2_cephes = cephes_dilog(v); const double li2_cephes_2 = cephes_dilog_2(v); @@ -687,6 +692,7 @@ TEST_CASE("test_real_random_values") #endif INFO("Li2(64) real = " << li2_327 << " (algorithm 327)"); INFO("Li2(64) real = " << li2_490 << " (algorithm 490)"); + INFO("Li2(64) real = " << li2_490_2 << " (algorithm 490 2)"); INFO("Li2(64) real = " << li2_babar << " (BaBar)"); INFO("Li2(64) real = " << li2_cephes << " (cephes)"); INFO("Li2(64) real = " << li2_cephes_2 << " (cephes 2)"); @@ -703,6 +709,7 @@ TEST_CASE("test_real_random_values") #endif CHECK_CLOSE(li2, li2_327 , 10*eps64); CHECK_CLOSE(li2, li2_490 , eps64); + CHECK_CLOSE(li2, li2_490_2 , eps64); CHECK_CLOSE(li2, li2_babar , 100*eps64); CHECK_CLOSE(li2, li2_cephes , 2*eps64); CHECK_CLOSE(li2, li2_cephes_2, 2*eps64);