25 constexpr double inf = std::numeric_limits<double>::infinity();
26 constexpr double nan = std::numeric_limits<double>::quiet_NaN();
27 constexpr double PI = 3.1415926535897932;
29 constexpr bool is_even(int64_t n)
noexcept {
return n % 2 == 0; }
31 constexpr bool is_finite(
const std::complex<double>& z)
noexcept
33 return std::isfinite(std::real(z)) && std::isfinite(std::imag(z));
37 std::complex<double>
clog(
const std::complex<double>& z)
noexcept
39 const double n = std::hypot(std::real(z), std::imag(z));
42 if (std::imag(z) == 0.0 && a < 0.0) {
51 std::complex<double>
Li_series(int64_t n,
const std::complex<double>& z)
noexcept
53 std::complex<double>
sum(0.0, 0.0), sum_old(0.0, 0.0), p = z;
60 p *= z*std::pow(k/(1.0 + k), n);
62 }
while (
sum != sum_old &&
63 k < std::numeric_limits<int64_t>::max() - 2);
69 std::complex<double>
Li_unity_pos(int64_t n,
const std::complex<double>& z)
noexcept
71 const std::complex<double> lnz =
clog(z);
72 const std::complex<double> lnz2 = lnz*lnz;
73 std::complex<double>
sum(
zeta(n), 0.0), p(1.0, 0.0);
75 for (int64_t j = 1; j < n - 1; ++j) {
76 p *= lnz/
static_cast<double>(j);
80 p *= lnz/
static_cast<double>(n - 1);
83 p *= lnz/
static_cast<double>(n);
86 p *= lnz/
static_cast<double>(n + 1);
89 for (int64_t j = (n + 3); j < std::numeric_limits<int64_t>::max() - 2; j += 2) {
90 p *= lnz2/
static_cast<double>((j - 1)*j);
91 const auto old_sum =
sum;
93 if (
sum == old_sum) {
break; }
100 std::complex<double>
stable_pow(
const std::complex<double>& z, int64_t n)
noexcept
102 if (std::imag(z) == 0) {
103 return { std::pow(std::real(z), n), 0.0 };
104 }
else if (std::real(z) == 0) {
105 const double p = std::pow(std::imag(z), n);
108 }
else if (n % 2 == 0) {
110 }
else if ((n - 1) % 4 == 0) {
115 return std::pow(z, n);
119 std::complex<double>
Li_unity_neg(int64_t n,
const std::complex<double>& z)
noexcept
121 const std::complex<double> lnz =
clog(z);
122 const std::complex<double> lnz2 = lnz*lnz;
123 std::complex<double>
sum = std::tgamma(1 - n)*
stable_pow(-lnz, n - 1);
124 std::complex<double> lnzk, sum_old, term;
143 }
while (
sum != sum_old);
149 std::complex<double>
Li_rest(int64_t n,
const std::complex<double>& z)
noexcept
151 const std::complex<double> lnz =
clog(-z);
152 const std::complex<double> lnz2 = lnz*lnz;
153 const int64_t kmax =
is_even(n) ? n/2 : (n - 1)/2;
154 std::complex<double> p =
is_even(n) ? 1.0 : lnz;
155 std::complex<double>
sum(0.0, 0.0), old_sum;
157 for (int64_t k = kmax; k != 0; --k) {
158 const double ifac =
inv_fac(n - 2*k);
159 if (ifac == 0) {
return 2.0*
sum; }
162 if (
sum == old_sum) {
break; }
181std::complex<double>
Li(int64_t n,
const std::complex<double>& z)
noexcept
183 if (std::isnan(std::real(z)) || std::isnan(std::imag(z))) {
185 }
else if (std::isinf(std::real(z)) || std::isinf(std::imag(z))) {
187 }
else if (z == 0.0) {
189 }
else if (z == 1.0) {
193 return {
zeta(n), 0.0};
194 }
else if (z == -1.0) {
200 if (4*
PI*
PI*nz < nl) {
202 }
else if (nl < 0.512*0.512*4*
PI*
PI) {
205 const auto sqrtz = std::sqrt(z);
206 return std::exp2(n - 1)*(
Li(n, sqrtz) +
Li(n, -sqrtz));
207 }
else if (n == -1) {
208 return z/((1.0 - z)*(1.0 - z));
212 return -
clog(1.0 - z);
226 const double sgn =
is_even(n) ? -1.0 : 1.0;
std::complex< double > Li_series(int64_t n, const std::complex< double > &z) noexcept
std::complex< double > clog(const std::complex< double > &z) noexcept
complex logarithm, converts -0.0 to 0.0
std::complex< double > stable_pow(const std::complex< double > &z, int64_t n) noexcept
returns z^n, treating Re(z) == 0 and Im(z) == 0 in a stable way
constexpr bool is_even(int64_t n) noexcept
std::complex< double > Li_unity_neg(int64_t n, const std::complex< double > &z) noexcept
Series expansion of Li_n(z) around z ~ 1, n < 0.
std::complex< double > Li_rest(int64_t n, const std::complex< double > &z) noexcept
returns remainder from inversion formula
std::complex< double > Li_unity_pos(int64_t n, const std::complex< double > &z) noexcept
Series expansion of Li_n(z) around z ~ 1, n > 0.
constexpr T norm(const Complex< T > &z) noexcept
double neg_eta(int64_t n) noexcept
negative Dirichlet eta function
std::complex< double > Li(int64_t n, const std::complex< double > &z) noexcept
Complex polylogarithm .
constexpr T arg(const Complex< T > &z) noexcept
double zeta(int64_t n) noexcept
Riemann zeta function for integer arguments.
std::complex< double > Li5(const std::complex< double > &z_) noexcept
Complex polylogarithm .
double Li4(double x) noexcept
Real 4-th order polylogarithm .
double Li2(double x) noexcept
Real dilogarithm .
bool is_finite(const gsl_vector *x)
Returns true if GSL vector contains only finite elements, false otherwise.
double Li3(double x) noexcept
Real trilogarithm .
Complex< T > log(const Complex< T > &z) noexcept
double harmonic(int64_t n) noexcept
harmonic number n
double inv_fac(int64_t n) noexcept
negative Dirichlet eta function
std::complex< double > Li6(const std::complex< double > &z_) noexcept
Complex polylogarithm .
auto sum(Idx ini, Idx fin, Function f) -> decltype(EvalEigenXpr< Idx >(ini, f))