Adding refined implementation of Algorithm 490 (#44)

for testing purposes

Author: Expander

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;
+}

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);
 

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");
 

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);