native_digamma_math.hpp

Go to the documentation of this file.

#pragma once
 
#include "native_pi_math.hpp"
 
namespace uni20::detail::native_math
{
template <NativeReal R> R positive_digamma(R x)
{
  constexpr int bits = numeric_limits<R>::digits();
  constexpr int cutoff = bits <= 24 ? 8 : bits <= 64 ? 16 : 32;
  constexpr int terms = bits <= 24 ? 4 : bits <= 53 ? 8 : bits <= 64 ? 10 : 16;
  // B_(2k)/(2k), generated by the exact rational recurrence in
  // scripts/generate-digamma-coefficients.py. All integer literals fit int64.
  constexpr R coefficients[] = {R(1LL) / R(12LL),
                                R(-1LL) / R(120LL),
                                R(1LL) / R(252LL),
                                R(-1LL) / R(240LL),
                                R(1LL) / R(132LL),
                                R(-691LL) / R(32760LL),
                                R(1LL) / R(12LL),
                                R(-3617LL) / R(8160LL),
                                R(43867LL) / R(14364LL),
                                R(-174611LL) / R(6600LL),
                                R(77683LL) / R(276LL),
                                R(-236364091LL) / R(65520LL),
                                R(657931LL) / R(12LL),
                                R(-3392780147LL) / R(3480LL),
                                R(1723168255201LL) / R(85932LL),
                                R(-7709321041217LL) / R(16320LL)};
  int shifts = x < R(cutoff) ? static_cast<int>(std::ceil(R(cutoff) - x)) : 0;
  R y = x + R(shifts), inverse = R(1) / y, square = inverse * inverse;
  R polynomial = coefficients[terms - 1];
  for (int k = terms - 2; k >= 0; --k)
    polynomial = std::fma(polynomial, square, coefficients[k]);
 
  R sum = std::log(y), low = 0;
  auto add = [&](R term) {
    R next = sum + term;
    low += std::abs(sum) >= std::abs(term) ? (sum - next) + term : (term - next) + sum;
    sum = next;
  };
  add(-inverse / 2);
  add(-square * polynomial);
  // Form each shifted argument from the original x, without accumulating the
  // rounding error of repeated x += 1. Retain cancellation residuals in low.
  for (int k = 0; k < shifts; ++k)
    add(-R(1) / (x + R(k)));
  return sum + low;
}
 
template <NativeReal R> R digamma(R x)
{
  using L = numeric_limits<R>;
  if (std::isnan(x) || x == -L::infinity()) return L::quiet_NaN();
  if (x == L::infinity()) return x;
  // The omitted term is O(x), negligible relative to the reciprocal here.
  // This also preserves signed-zero poles and the final overflow boundary,
  // without passing tiny x through a subnormal trigonometric denominator.
  if (std::abs(x) <= L::epsilon()) return -R(1) / x - euler_gamma<R>;
  if (x > 0) return native_math::positive_digamma(x);
  if (x == std::trunc(x)) return L::quiet_NaN();
  return native_math::positive_digamma(R(1) - x) - pi<R> / native_math::tanpi(x);
}
} // namespace uni20::detail::native_math