From 4d5bd5b9edb8b9ddcb39b0df6cba83f08a138bbe Mon Sep 17 00:00:00 2001 From: Alexander Voigt Date: Mon, 18 Mar 2024 12:52:33 +0100 Subject: [PATCH] update dilog and clausen_2 (#17) Co-authored-by: Alexander Voigt --- source/himalaya/misc/Li2.cpp | 115 +++++++++++++++------------ source/himalaya/misc/complex.hpp | 129 ------------------------------- 2 files changed, 67 insertions(+), 177 deletions(-) delete mode 100644 source/himalaya/misc/complex.hpp diff --git a/source/himalaya/misc/Li2.cpp b/source/himalaya/misc/Li2.cpp index 7b596086..71094131 100644 --- a/source/himalaya/misc/Li2.cpp +++ b/source/himalaya/misc/Li2.cpp @@ -6,7 +6,6 @@ // ==================================================================== #include "himalaya/misc/Li2.hpp" -#include "himalaya/misc/complex.hpp" #include #include @@ -20,24 +19,41 @@ namespace himalaya { namespace { - template - Complex horner(const Complex& z, const T (&coeffs)[N]) noexcept - { - static_assert(0 <= Nstart && Nstart < N && N >= 2, "invalid array bounds"); +template +std::complex horner(const std::complex& 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 +std::complex log1p(const std::complex& z) noexcept +{ + const std::complex u = T(1) + z; - return Complex(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 /** @@ -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; @@ -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 dilog(const std::complex& z_) noexcept +std::complex dilog(const std::complex& z) noexcept { const double PI = 3.1415926535897932; - const Complex 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}] @@ -148,50 +163,53 @@ std::complex dilog(const std::complex& 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::epsilon()) { return z*(1.0 + 0.25*z); } - Complex u(0.0, 0.0), rest(0.0, 0.0); + std::complex 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 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 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 u2(u*u); + const auto u2(u*u); return sgn*(u + u2*(bf[0] + u*horner<1>(u2, bf))) + rest; } @@ -225,7 +243,9 @@ double clausen_2(double x) noexcept sgn = -sgn; } - if (x == 0 || x == PI) { + if (x == 0) { + return x; + } else if (x == PI) { return 0; } @@ -233,20 +253,19 @@ double clausen_2(double x) noexcept 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, diff --git a/source/himalaya/misc/complex.hpp b/source/himalaya/misc/complex.hpp deleted file mode 100644 index 07f7a6bc..00000000 --- a/source/himalaya/misc/complex.hpp +++ /dev/null @@ -1,129 +0,0 @@ -// ==================================================================== -// This file is part of Himalaya. -// -// Himalaya is licenced under the GNU General Public License (GNU GPL) -// version 3. -// ==================================================================== - -#pragma once - -#include -#include - -namespace himalaya { - -template -struct Complex { - constexpr Complex(T re_ = T{}, T im_ = T{}) : re(re_), im(im_) {} - operator std::complex() const noexcept { return std::complex(re, im); } - T re{}; - T im{}; -}; - -/// converts -0.0 to 0.0 -template -Complex pos(const Complex& z) noexcept -{ - const T rz = z.re == T(0) ? std::abs(z.re) : z.re; - const T iz = z.im == T(0) ? std::abs(z.im) : z.im; - return { rz, iz }; -} - -template -constexpr T arg(const Complex& z) noexcept -{ - return std::atan2(z.im, z.re); -} - -template -constexpr Complex conj(const Complex& z) noexcept -{ - return { z.re, -z.im }; -} - -template -Complex log(const Complex& z_) noexcept -{ - const Complex z = pos(z_); - return { 0.5*std::log(norm_sqr(z)), arg(z) }; -} - -template -constexpr T norm_sqr(const Complex& z) noexcept -{ - return z.re*z.re + z.im*z.im; -} - -template -constexpr Complex operator+(const Complex& a, const Complex& b) noexcept -{ - return { a.re + b.re, a.im + b.im }; -} - -template -constexpr Complex operator+(const Complex& z, T x) noexcept -{ - return { z.re + x, z.im }; -} - -template -constexpr Complex operator+(T x, const Complex& z) noexcept -{ - return { x + z.re, z.im }; -} - -template -constexpr Complex operator-(const Complex& a, const Complex& b) noexcept -{ - return { a.re - b.re, a.im - b.im }; -} - -template -constexpr Complex operator-(T x, const Complex& z) noexcept -{ - return { x - z.re, -z.im }; -} - -template -constexpr Complex operator-(const Complex& z, T x) noexcept -{ - return { z.re - x, z.im }; -} - -template -constexpr Complex operator-(const Complex& z) noexcept -{ - return { -z.re, -z.im }; -} - -template -constexpr Complex operator*(const Complex& a, const Complex& b) noexcept -{ - return { a.re*b.re - a.im*b.im, a.re*b.im + a.im*b.re }; -} - -template -constexpr Complex operator*(T x, const Complex& z) noexcept -{ - return { x*z.re, x*z.im }; -} - -template -constexpr Complex operator*(const Complex& z, T x) noexcept -{ - return x*z; -} - -template -constexpr Complex operator/(T x, const Complex& z) noexcept -{ - return x*conj(z)/norm_sqr(z); -} - -template -constexpr Complex operator/(const Complex& z, T x) noexcept -{ - return { z.re/x, z.im/x }; -} - -} // namespace himalaya