Skip to content

Commit

Permalink
Bugfix: imprecision of complex li2 for small arguments (#39)
Browse files Browse the repository at this point in the history 
Note: This change leads to a performance decrease due to the use of `std::log(const std::complex<T>&)`, which is more precise than `std::log(|z|) + I*std::arg(z)` for values of z close to unity.
  • Loading branch information
@Expander
Expander authored Mar 14, 2024
1 parent d47291c commit a6ea863
Show file tree
Hide file tree
Showing 34 changed files with 656 additions and 665 deletions.
1 change: 1 addition & 0 deletions src/c/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@ add_library(polylogarithm_c
Li4.c
Li5.c
Li6.c
log.c
)

target_include_directories(polylogarithm_c PUBLIC ${CMAKE_CURRENT_SOURCE_DIR})
Expand Down
44 changes: 22 additions & 22 deletions src/c/Li2.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include <complex.h>
#include <float.h>
#include <math.h>
#include "fast_clog.h"
#include "log.h"


static long double hornerl(long double x, const long double* c, int len)
Expand Down Expand Up @@ -331,23 +331,23 @@ float _Complex cli2f(float _Complex z)
/* transformation to |z|<1, Re(z)<=0.5f */
if (rz <= 0.5f) {
if (nz > 1.0f) {
const float _Complex lz = fast_clogf(-z);
u = -fast_clogf(1.0f - 1.0f / z);
const float _Complex lz = clogf(-z);
u = -clogf(1.0f - 1.0f / z);
rest = -0.5f*lz*lz - PI*PI/6;
sgn = -1;
} else { /* nz <= 1 */
u = -fast_clogf(1.0f - z);
u = -clogf(1.0f - z);
rest = 0;
sgn = 1;
}
} else { /* rz > 0.5f */
if (nz <= 2*rz) {
u = -fast_clogf(z);
rest = u*fast_clogf(1.0f - z) + PI*PI/6;
u = -clogf(z);
rest = u*clogf(1.0f - z) + PI*PI/6;
sgn = -1;
} else { /* nz > 2*rz */
const float _Complex lz = fast_clogf(-z);
u = -fast_clogf(1.0f - 1.0f / z);
const float _Complex lz = clogf(-z);
u = -clogf(1.0f - 1.0f / z);
rest = -0.5f*lz*lz - PI*PI/6;
sgn = -1;
}
Expand Down Expand Up @@ -412,23 +412,23 @@ double _Complex cli2(double _Complex z)
/* transformation to |z|<1, Re(z)<=0.5 */
if (rz <= 0.5) {
if (nz > 1.0) {
const double _Complex lz = fast_clog(-z);
u = -fast_clog(1.0 - 1.0 / z);
const double _Complex lz = clog(-z);
u = -clog1p(-1.0/z);
rest = -0.5*lz*lz - PI*PI/6;
sgn = -1;
} else { /* nz <= 1 */
u = -fast_clog(1.0 - z);
u = -clog1p(-z);
rest = 0;
sgn = 1;
}
} else { /* rz > 0.5 */
if (nz <= 2*rz) {
u = -fast_clog(z);
rest = u*fast_clog(1.0 - z) + PI*PI/6;
u = -clog(z);
rest = u*clog1p(-z) + PI*PI/6;
sgn = -1;
} else { /* nz > 2*rz */
const double _Complex lz = fast_clog(-z);
u = -fast_clog(1.0 - 1.0 / z);
const double _Complex lz = clog(-z);
u = -clog1p(-1.0/z);
rest = -0.5*lz*lz - PI*PI/6;
sgn = -1;
}
Expand Down Expand Up @@ -517,23 +517,23 @@ long double _Complex cli2l(long double _Complex z)
/* transformation to |z|<1, Re(z)<=0.5 */
if (rz <= 0.5L) {
if (nz > 1.0L) {
const long double _Complex lz = fast_clogl(-z);
u = -fast_clogl(1.0L - 1.0L/z);
const long double _Complex lz = clogl(-z);
u = -clog1pl(-1.0L/z);
rest = -0.5L*lz*lz - PI*PI/6;
sgn = -1;
} else { /* nz <= 1 */
u = -fast_clogl(1.0L - z);
u = -clog1pl(-z);
rest = 0;
sgn = 1;
}
} else { /* rz > 0.5L */
if (nz <= 2*rz) {
u = -fast_clogl(z);
rest = u * fast_clogl(1.0L - z) + PI*PI/6;
u = -clogl(z);
rest = u * clog1pl(-z) + PI*PI/6;
sgn = -1;
} else { /* nz > 2*rz */
const long double _Complex lz = fast_clogl(-z);
u = -fast_clogl(1.0L - 1.0L/z);
const long double _Complex lz = clogl(-z);
u = -clog1pl(-1.0L/z);
rest = -0.5L*lz*lz - PI*PI/6;
sgn = -1;
}
Expand Down
14 changes: 7 additions & 7 deletions src/c/Li3.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include <complex.h>
#include <float.h>
#include <math.h>
#include "fast_clog.h"
#include "log.h"

/// Li_3(x) for x in [-1,0]
static double li3_neg(double x)
Expand Down Expand Up @@ -145,7 +145,7 @@ double _Complex cli3(double _Complex z)
const double _Complex u4 = u2*u2;
const double _Complex u8 = u4*u4;
const double _Complex c0 = zeta3 + u*(zeta2 - u2/12.0);
const double _Complex c1 = 0.25 * (3.0 - 2.0*fast_pos_clog(-u));
const double _Complex c1 = 0.25 * (3.0 - 2.0*pos_clog(-u));

const double cs[7] = {
-3.4722222222222222e-03, 1.1574074074074074e-05,
Expand All @@ -165,11 +165,11 @@ double _Complex cli3(double _Complex z)
double _Complex u = 0.0, rest = 0.0;

if (nz <= 1.0) {
u = -fast_pos_clog(1.0 - z);
u = -pos_clog(1.0 - z);
} else { // nz > 1
const double arg = pz > 0.0 ? pz - PI : pz + PI;
const double _Complex lmz = lnz + arg*I; // clog(-z)
u = -fast_pos_clog(1.0 - 1.0/z);
u = -pos_clog(1.0 - 1.0/z);
rest = -lmz*(lmz*lmz/6.0 + zeta2);
}

Expand Down Expand Up @@ -275,7 +275,7 @@ long double _Complex cli3l(long double _Complex z)
const long double _Complex u = lnz + pz*I; // clog(z)
const long double _Complex u2 = u*u;
const long double _Complex c0 = zeta3 + u*(zeta2 - u2/12.0L);
const long double _Complex c1 = 0.25L * (3.0L - 2.0L*fast_pos_clogl(-u));
const long double _Complex c1 = 0.25L * (3.0L - 2.0L*pos_clogl(-u));

const long double cs[] = {
-3.47222222222222222222222222222222222e-03L,
Expand Down Expand Up @@ -317,11 +317,11 @@ long double _Complex cli3l(long double _Complex z)
long double _Complex u = 0.0L, rest = 0.0L;

if (nz <= 1.0L) {
u = -fast_pos_clogl(1.0L - z);
u = -pos_clogl(1.0L - z);
} else { // nz > 1
const long double arg = pz > 0.0L ? pz - PI : pz + PI;
const long double _Complex lmz = lnz + arg*I; // clog(-z)
u = -fast_pos_clogl(1.0L - 1.0L/z);
u = -pos_clogl(1.0L - 1.0L/z);
rest = -lmz*(lmz*lmz/6.0L + zeta2);
}

Expand Down
14 changes: 7 additions & 7 deletions src/c/Li4.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include <complex.h>
#include <float.h>
#include <math.h>
#include "fast_clog.h"
#include "log.h"

/// Li_4(x) for x in [-1,0]
static double li4_neg(double x)
Expand Down Expand Up @@ -200,7 +200,7 @@ double _Complex cli4(double _Complex z)
const double _Complex u8 = u4*u4;
const double c1 = 1.2020569031595943; // zeta(3)
const double c2 = 0.82246703342411322;
const double _Complex c3 = (11.0/6.0 - fast_pos_clog(-u))/6.0;
const double _Complex c3 = (11.0/6.0 - pos_clog(-u))/6.0;
const double c4 = -1.0/48.0;

const double cs[7] = {
Expand All @@ -224,12 +224,12 @@ double _Complex cli4(double _Complex z)
double sgn = 1;

if (nz <= 1.0) {
u = -fast_pos_clog(1.0 - z);
u = -pos_clog(1.0 - z);
} else { // nz > 1
const double arg = pz > 0.0 ? pz - PI : pz + PI;
const double _Complex lmz = lnz + arg*I; // clog(-z)
const double _Complex lmz2 = lmz*lmz;
u = -fast_pos_clog(1.0 - 1.0/z);
u = -pos_clog(1.0 - 1.0/z);
r = 1.0/360.0*(-7*PI4 + lmz2*(-30.0*PI2 - 15.0*lmz2));
sgn = -1;
}
Expand Down Expand Up @@ -336,7 +336,7 @@ long double _Complex cli4l(long double _Complex z)
const long double _Complex u2 = u*u;
const long double c1 = 1.20205690315959428539973816151144999L; // zeta(3)
const long double c2 = 0.822467033424113218236207583323012595L;
const long double _Complex c3 = (11.0L/6.0L - fast_pos_clogl(-u))/6.0L;
const long double _Complex c3 = (11.0L/6.0L - pos_clogl(-u))/6.0L;
const long double c4 = -1.0L/48.0L;

const long double cs[] = {
Expand Down Expand Up @@ -377,12 +377,12 @@ long double _Complex cli4l(long double _Complex z)
long double sgn = 1;

if (nz <= 1.0L) {
u = -fast_pos_clogl(1.0L - z);
u = -pos_clogl(1.0L - z);
} else { // nz > 1
const long double arg = pz > 0.0 ? pz - PI : pz + PI;
const long double _Complex lmz = lnz + arg*I; // clog(-z)
const long double _Complex lmz2 = lmz*lmz;
u = -fast_pos_clogl(1.0L - 1.0L/z);
u = -pos_clogl(1.0L - 1.0L/z);
r = 1.0L/360.0L*(-7*PI4 + lmz2*(-30.0L*PI2 - 15.0L*lmz2));
sgn = -1;
}
Expand Down
14 changes: 7 additions & 7 deletions src/c/Li5.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include <complex.h>
#include <float.h>
#include <math.h>
#include "fast_clog.h"
#include "log.h"


/**
Expand Down Expand Up @@ -61,7 +61,7 @@ double _Complex cli5(double _Complex z)
const double c1 = 1.0823232337111382; // zeta(4)
const double c2 = 0.60102845157979714; // zeta(3)/2
const double c3 = 0.27415567780803774;
const double _Complex c4 = (25.0/12.0 - fast_pos_clog(-u))/24.0;
const double _Complex c4 = (25.0/12.0 - pos_clog(-u))/24.0;
const double c5 = -1.0/240.0;

const double cs[6] = {
Expand All @@ -84,12 +84,12 @@ double _Complex cli5(double _Complex z)
double _Complex u = 0.0, rest = 0.0;

if (nz <= 1.0) {
u = -fast_pos_clog(1.0 - z);
u = -pos_clog(1.0 - z);
} else { // nz > 1
const double arg = pz > 0.0 ? pz - PI : pz + PI;
const double _Complex lmz = lnz + arg*I; // clog(-z)
const double _Complex lmz2 = lmz*lmz;
u = -fast_pos_clog(1.0 - 1.0/z);
u = -pos_clog(1.0 - 1.0/z);
rest = -1.0/360.0*lmz*(7*PI4 + lmz2*(10.0*PI2 + 3.0*lmz2));
}

Expand Down Expand Up @@ -196,7 +196,7 @@ long double _Complex cli5l(long double _Complex z)
const long double c1 = 1.08232323371113819151600369654116790L; // zeta(4)
const long double c2 = 0.601028451579797142699869080755724995L; // zeta(3)/2
const long double c3 = 0.274155677808037739412069194441004198L;
const long double _Complex c4 = (25.0L/12.0L - fast_pos_clogl(-u))/24.0L;
const long double _Complex c4 = (25.0L/12.0L - pos_clogl(-u))/24.0L;
const long double c5 = -1.0L/240.0L;

const long double cs[] = {
Expand Down Expand Up @@ -236,12 +236,12 @@ long double _Complex cli5l(long double _Complex z)
long double _Complex u = 0.0L, rest = 0.0L;

if (nz <= 1.0L) {
u = -fast_pos_clogl(1.0L - z);
u = -pos_clogl(1.0L - z);
} else { // nz > 1
const long double arg = pz > 0.0 ? pz - PI : pz + PI;
const long double _Complex lmz = lnz + arg*I; // clog(-z)
const long double _Complex lmz2 = lmz*lmz;
u = -fast_pos_clogl(1.0L - 1.0L/z);
u = -pos_clogl(1.0L - 1.0L/z);
rest = -1.0L/360.0L*lmz*(7*PI4 + lmz2*(10.0L*PI2 + 3.0L*lmz2));
}

Expand Down
14 changes: 7 additions & 7 deletions src/c/Li6.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
#include <complex.h>
#include <float.h>
#include <math.h>
#include "fast_clog.h"
#include "log.h"


/**
Expand Down Expand Up @@ -62,7 +62,7 @@ double _Complex cli6(double _Complex z)
const double c2 = 0.54116161685556910;
const double c3 = 0.20034281719326571;
const double c4 = 0.068538919452009435;
const double _Complex c5 = (137.0/60.0 - fast_pos_clog(-u))/120.0;
const double _Complex c5 = (137.0/60.0 - pos_clog(-u))/120.0;
const double c6 = -1.0/1440.0;

const double cs[5] = {
Expand All @@ -86,12 +86,12 @@ double _Complex cli6(double _Complex z)
double sgn = 1;

if (nz <= 1.0) {
u = -fast_pos_clog(1.0 - z);
u = -pos_clog(1.0 - z);
} else { // nz > 1
const double arg = pz > 0.0 ? pz - PI : pz + PI;
const double _Complex lmz = lnz + arg*I; // clog(-z)
const double _Complex lmz2 = lmz*lmz;
u = -fast_pos_clog(1.0 - 1.0/z);
u = -pos_clog(1.0 - 1.0/z);
r = -31.0*PI6/15120.0 + lmz2*(-7.0/720.0*PI4 + lmz2*(-1.0/144.0*PI2 - 1.0/720.0*lmz2));
sgn = -1;
}
Expand Down Expand Up @@ -202,7 +202,7 @@ long double _Complex cli6l(long double _Complex z)
const long double c2 = 0.541161616855569095758001848270583951L;
const long double c3 = 0.200342817193265714233289693585241665L;
const long double c4 = 0.0685389194520094348530172986102510496L;
const long double _Complex c5 = (137.0L/60.0L - fast_pos_clogl(-u))/120.0L;
const long double _Complex c5 = (137.0L/60.0L - pos_clogl(-u))/120.0L;
const long double c6 = -1.0L/1440.0L;

const long double cs[] = {
Expand Down Expand Up @@ -247,12 +247,12 @@ long double _Complex cli6l(long double _Complex z)
long double sgn = 1;

if (nz <= 1.0L) {
u = -fast_pos_clogl(1.0L - z);
u = -pos_clogl(1.0L - z);
} else { // nz > 1
const long double arg = pz > 0.0 ? pz - PI : pz + PI;
const long double _Complex lmz = lnz + arg*I; // clog(-z)
const long double _Complex lmz2 = lmz*lmz;
u = -fast_pos_clogl(1.0L - 1.0L/z);
u = -pos_clogl(1.0L - 1.0L/z);
r = -31.0L*PI6/15120.0L
+ lmz2*(-7.0L/720.0L*PI4 +
lmz2*(-1.0L/144.0L*PI2 - 1.0L/720.0L*lmz2));
Expand Down
Loading

0 comments on commit a6ea863

Please sign in to comment.