gm2_dilog.cpp

Go to the documentation of this file.

// ====================================================================
// This file is part of GM2Calc.
//
// GM2Calc is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// GM2Calc is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with GM2Calc.  If not, see
// <http://www.gnu.org/licenses/>.
// ====================================================================
 
#include "gm2_dilog.hpp"
#include <cmath>
#include <limits>
 
namespace gm2calc {
 
namespace {
 
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;
   }
 
   return { rz*a + b, iz*a };
}
 
/// 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;
 
   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)));
}
 
} // anonymous namespace
 
/**
 * @brief Real dilogarithm \f$\operatorname{Li}_2(x)\f$
 * @param x real argument
 * @return \f$\operatorname{Li}_2(x)\f$
 * @author Alexander Voigt
 * @note taken from the polylogarithm library 7.0.0
 *
 * Implemented as an economized Pade approximation with a
 * maximum error of 4.16e-18.
 * [[arXiv:2201.01678](https://arxiv.org/abs/2201.01678)].
 */
double dilog(double x) noexcept
{
   const double PI = 3.1415926535897932;
   const double P[] = {
      0.9999999999999999502e+0,
     -2.6883926818565423430e+0,
      2.6477222699473109692e+0,
     -1.1538559607887416355e+0,
      2.0886077795020607837e-1,
     -1.0859777134152463084e-2
   };
   const double Q[] = {
      1.0000000000000000000e+0,
     -2.9383926818565635485e+0,
      3.2712093293018635389e+0,
     -1.7076702173954289421e+0,
      4.1596017228400603836e-1,
     -3.9801343754084482956e-2,
      8.2743668974466659035e-4
   };
 
   double y = 0, r = 0, s = 1;
 
   // transform to [0, 1/2]
   if (x < -1) {
      const double l = std::log(1 - x);
      y = 1/(1 - x);
      r = -PI*PI/6 + l*(0.5*l - std::log(-x));
      s = 1;
   } else if (x == -1) {
      return -PI*PI/12;
   } else if (x < 0) {
      const double l = std::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 - std::log(x)*std::log1p(-x);
      s = -1;
   } else if (x == 1) {
      return PI*PI/6;
   } else if (x < 2) {
      const double l = std::log(x);
      y = 1 - 1/x;
      r = PI*PI/6 - l*(std::log(y) + 0.5*l);
      s = 1;
   } else {
      const double l = std::log(x);
      y = 1/x;
      r = PI*PI/3 - 0.5*l*l;
      s = -1;
   }
 
   const double y2 = y*y;
   const double y4 = y2*y2;
   const double p = P[0] + y * P[1] + y2 * (P[2] + y * P[3]) +
                    y4 * (P[4] + y * P[5]);
   const double q = Q[0] + y * Q[1] + y2 * (Q[2] + y * Q[3]) +
                    y4 * (Q[4] + y * Q[5] + y2 * Q[6]);
 
   return r + s*y*p/q;
}
 
/**
 * @brief Complex dilogarithm \f$\mathrm{Li}_2(z)\f$
 * @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
 * @note taken from the polylogarithm library 7.0.0
 */
std::complex<double> dilog(const std::complex<double>& z) noexcept
{
   const double PI = 3.1415926535897932;
 
   // bf[1..N-1] are the even Bernoulli numbers / (2 n + 1)!
   // generated by: Table[BernoulliB[2 n]/(2 n + 1)!, {n, 1, 9}]
   const double bf[10] = {
      - 1.0/4.0,
      + 1.0/36.0,
      - 1.0/3600.0,
      + 1.0/211680.0,
      - 1.0/10886400.0,
      + 1.0/526901760.0,
      - 4.0647616451442255e-11,
      + 8.9216910204564526e-13,
      - 1.9939295860721076e-14,
      + 4.5189800296199182e-16
   };
 
   const double rz = std::real(z);
   const double iz = std::imag(z);
 
   // special cases
   if (iz == 0) {
      if (rz <= 1) {
         return { dilog(rz), iz };
      }
      // rz > 1
      return { dilog(rz), -PI*std::log(rz) };
   }
 
   const double nz = std::norm(z);
 
   if (nz < std::numeric_limits<double>::epsilon()) {
      return z*(1.0 + 0.25*z);
   }
 
   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 (rz <= 0.5) {
      if (nz > 1) {
         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 = -log1p(-z);
         rest = 0;
         sgn = 1;
      }
   } else { // rz > 0.5
      if (nz <= 2*rz) {
         u = -std::log(z);
         rest = u*log1p(-z) + PI*PI/6;
         sgn = -1;
      } 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 auto u2(u*u);
 
   return sgn*(u + u2*(bf[0] + u*horner<1>(u2, bf))) + rest;
}
 
/**
 * @brief Clausen function \f$\mathrm{Cl}_2(\theta) = \mathrm{Im}(\mathrm{Li}_2(e^{i\theta}))\f$
 * @param x real angle
 * @return \f$\mathrm{Cl}_2(\theta)\f$
 * @author Alexander Voigt
 * @note Implemented as economized Padé approximation.
 * @note taken from the polylogarithm library 7.0.0
 */
double clausen_2(double x) noexcept
{
   const double PI = 3.14159265358979324;
   const double PI2 = 2*PI, PIH = PI/2, PI28 = PI*PI/8;
   double sgn = 1;
 
   if (x < 0) {
      x = -x;
      sgn = -1;
   }
 
   if (x >= PI2) {
      x = std::fmod(x, PI2);
   }
 
   if (x > PI) {
      const double p0 = 6.28125;
      const double p1 = 0.0019353071795864769253;
      x = (p0 - x) + p1;
      sgn = -sgn;
   }
 
   if (x == 0) {
      return x;
   } else if (x == PI) {
      return 0;
   }
 
   double h = 0;
 
   if (x < PIH) {
      const double P[] = {
         1.3888888888888889e-02, -4.3286930203743071e-04,
         3.2779814789973427e-06, -3.6001540369575084e-09
      };
      const double Q[] = {
         1.0000000000000000e+00, -3.6166589746694121e-02,
         3.6015827281202639e-04, -8.3646182842184428e-07
      };
      const double y = x*x;
      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);
   } else {
      const double P[] = {
         6.4005702446195512e-01, -2.0641655351338783e-01,
         2.4175305223497718e-02, -1.2355955287855728e-03,
         2.5649833551291124e-05, -1.4783829128773320e-07
      };
      const double Q[] = {
         1.0000000000000000e+00, -2.5299102015666356e-01,
         2.2148751048467057e-02, -7.8183920462457496e-04,
         9.5432542196310670e-06, -1.8184302880448247e-08
      };
      const double y = PI - x;
      const double z = y*y - PI28;
      const double z2 = z*z;
      const double z4 = z2*z2;
      const double p = P[0] + z * P[1] + z2 * (P[2] + z * P[3]) +
         z4 * (P[4] + z * P[5]);
      const double q = Q[0] + z * Q[1] + z2 * (Q[2] + z * Q[3]) +
         z4 * (Q[4] + z * Q[5]);
 
      h = y*p/q;
   }
 
   return sgn*h;
}
 
} // namespace gm2calc