Skip to content

Commit

Permalink
update dilog and clausen_2
Browse files Browse the repository at this point in the history
  • Loading branch information
@Expander
Expander committed Mar 18, 2024
1 parent ce0336b commit fbddbbf
Show file tree
Hide file tree
Showing 2 changed files with 67 additions and 177 deletions.
115 changes: 67 additions & 48 deletions source/himalaya/misc/Li2.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
// ====================================================================

#include "himalaya/misc/Li2.hpp"
#include "himalaya/misc/complex.hpp"
#include <cmath>
#include <limits>

Expand All @@ -20,24 +19,41 @@ namespace himalaya {

namespace {

template <int Nstart, typename T, int N>
Complex<T> horner(const Complex<T>& z, const T (&coeffs)[N]) noexcept
{
static_assert(0 <= Nstart && Nstart < N && N >= 2, "invalid array bounds");
template <int Nstart, typename T, int N>
std::complex<T> horner(const std::complex<T>& z, const T (&coeffs)[N]) noexcept
{
static_assert(0 <= Nstart && Nstart < N && N >= 2, "invalid array bounds");

const T rz = std::real(z);
const T iz = std::imag(z);
const T r = rz + rz;
const T s = std::norm(z);
T a = coeffs[N - 1], b = coeffs[N - 2];

for (int i = N - 3; i >= Nstart; --i) {
const T t = a;
a = b + r*a;
b = coeffs[i] - s*t;
}

const T r = z.re + z.re;
const T s = z.re * z.re + z.im * z.im;
T a = coeffs[N - 1], b = coeffs[N - 2];
return { rz*a + b, iz*a };
}

for (int i = N - 3; i >= Nstart; --i) {
const T t = a;
a = b + r * a;
b = coeffs[i] - s * t;
}
/// returns log(1 + z) for complex z
template <typename T>
std::complex<T> log1p(const std::complex<T>& z) noexcept
{
const std::complex<T> u = T(1) + z;

return Complex<T>(z.re*a + b, z.im*a);
if (std::real(u) == T(1) && std::imag(u) == T(0)) {
return z;
} else if (std::real(u) <= T(0)) {
return std::log(u);
}

return std::log(u)*(z/(u - T(1)));
}

} // namespace

/**
Expand Down Expand Up @@ -87,14 +103,14 @@ double dilog(double x) noexcept
r = -0.5*l*l;
s = -1;
} else if (x == 0) {
return 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 - std::log(x)*std::log(y);
r = PI*PI/6 - std::log(x)*std::log1p(-x);
s = -1;
} else if (x == 1) {
return PI*PI/6;
Expand Down Expand Up @@ -122,16 +138,15 @@ double dilog(double x) noexcept

/**
* @brief Complex dilogarithm \f$\mathrm{Li}_2(z)\f$
* @param z_ complex argument
* @param z complex argument
* @return \f$\mathrm{Li}_2(z)\f$
* @note Implementation translated from SPheno to C++
* @author Werner Porod
* @note translated to C++ by Alexander Voigt
*/
std::complex<double> dilog(const std::complex<double>& z_) noexcept
std::complex<double> dilog(const std::complex<double>& z) noexcept
{
const double PI = 3.1415926535897932;
const Complex<double> z = { std::real(z_), std::imag(z_) };

// bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
// generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 9}]
Expand All @@ -148,50 +163,53 @@ std::complex<double> dilog(const std::complex<double>& z_) noexcept
+ 4.5189800296199182e-16
};

const double rz = std::real(z);
const double iz = std::imag(z);

// special cases
if (z.im == 0) {
if (z.re <= 1) {
return dilog(z.re);
if (iz == 0) {
if (rz <= 1) {
return { dilog(rz), iz };
}
// z.re > 1
return { dilog(z.re), -PI*std::log(z.re) };
// rz > 1
return { dilog(rz), -PI*std::log(rz) };
}

const double nz = norm_sqr(z);
const double nz = std::norm(z);

if (nz < std::numeric_limits<double>::epsilon()) {
return z*(1.0 + 0.25*z);
}

Complex<double> u(0.0, 0.0), rest(0.0, 0.0);
std::complex<double> u(0.0, 0.0), rest(0.0, 0.0);
double sgn = 1;

// transformation to |z|<1, Re(z)<=0.5
if (z.re <= 0.5) {
if (rz <= 0.5) {
if (nz > 1) {
const Complex<double> lz = log(-z);
u = -log(1.0 - 1.0 / z);
const auto lz = std::log(-z);
u = -log1p(-1.0/z);
rest = -0.5*lz*lz - PI*PI/6;
sgn = -1;
} else { // nz <= 1
u = -log(1.0 - z);
u = -log1p(-z);
rest = 0;
sgn = 1;
}
} else { // z.re > 0.5
if (nz <= 2*z.re) {
u = -log(z);
rest = u*log(1.0 - z) + PI*PI/6;
} else { // rz > 0.5
if (nz <= 2*rz) {
u = -std::log(z);
rest = u*log1p(-z) + PI*PI/6;
sgn = -1;
} else { // nz > 2*z.re
const Complex<double> lz = log(-z);
u = -log(1.0 - 1.0 / z);
} else { // nz > 2*rz
const auto lz = std::log(-z);
u = -log1p(-1.0/z);
rest = -0.5*lz*lz - PI*PI/6;
sgn = -1;
}
}

const Complex<double> u2(u*u);
const auto u2(u*u);

return sgn*(u + u2*(bf[0] + u*horner<1>(u2, bf))) + rest;
}
Expand Down Expand Up @@ -225,28 +243,29 @@ double clausen_2(double x) noexcept
sgn = -sgn;
}

if (x == 0 || x == PI) {
if (x == 0) {
return x;
} else if (x == PI) {
return 0;
}

double h = 0;

if (x < PIH) {
const double P[] = {
2.7951565822419270e-02, -8.8865360514541522e-04,
6.8282348222485902e-06, -7.5276232403566808e-09
1.3888888888888889e-02, -4.3286930203743071e-04,
3.2779814789973427e-06, -3.6001540369575084e-09
};
const double Q[] = {
1.0000000000000000e+00, -3.6904397961160525e-02,
3.7342870576106476e-04, -8.7460760866531179e-07
1.0000000000000000e+00, -3.6166589746694121e-02,
3.6015827281202639e-04, -8.3646182842184428e-07
};
const double y = x*x;
const double z = y - PI28;
const double z2 = z*z;
const double p = P[0] + z * P[1] + z2 * (P[2] + z * P[3]);
const double q = Q[0] + z * Q[1] + z2 * (Q[2] + z * Q[3]);
const double y2 = y*y;
const double p = P[0] + y * P[1] + y2 * (P[2] + y * P[3]);
const double q = Q[0] + y * Q[1] + y2 * (Q[2] + y * Q[3]);

h = x*(1 - std::log(x) + y*p/q/2);
h = x*(1 - std::log(x) + y*p/q);
} else {
const double P[] = {
6.4005702446195512e-01, -2.0641655351338783e-01,
Expand Down
129 changes: 0 additions & 129 deletions source/himalaya/misc/complex.hpp
View file

This file was deleted.

0 comments on commit fbddbbf

Please sign in to comment.