Skip to content

Adding refined implementation of Algorithm 490 #44

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 1 commit into from
Mar 20, 2025
Merged

Adding refined implementation of Algorithm 490 #44

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
79 changes: 79 additions & 0 deletions test/alt/algorithm_490/algorithm_490.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
}
1 change: 1 addition & 0 deletions test/alt/alt.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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);

Expand Down
3 changes: 3 additions & 0 deletions test/bench_Li.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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");

Expand Down
7 changes: 7 additions & 0 deletions test/test_Li2.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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)");
Expand All @@ -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);
Expand Down Expand Up @@ -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);
Expand All @@ -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)");
Expand All @@ -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);
Expand Down