native_upper_gamma_guard.hpp

Go to the documentation of this file.

#pragma once
 
#include <uni20/core/detail/native_math.hpp>
 
namespace uni20::detail::native_math::upper_gamma_guard
{
// Two native words retain guard bits through cancellation and final range
// decisions. No MPFR or wider scalar type is used. See native_upper_gamma.md.
template <NativeReal R> struct two_word
{
    R hi = 0, lo = 0;
    two_word() = default;
    two_word(R x) : hi(x) {}
    two_word(R h, R l) : hi(h), lo(l) {}
    static two_word sum(R a, R b)
    {
      R s = a + b, v = s - a;
      return {s, (a - (s - v)) + (b - v)};
    }
    friend two_word operator+(two_word a, two_word b)
    {
      auto s = two_word::sum(a.hi, b.hi);
      auto t = two_word::sum(a.lo, b.lo);
      auto u = two_word::sum(s.lo, t.hi);
      auto v = two_word::sum(s.hi, u.hi);
      return two_word::sum(v.hi, v.lo + u.lo + t.lo);
    }
    friend two_word operator-(two_word a) { return {-a.hi, -a.lo}; }
    friend two_word operator-(two_word a, two_word b) { return a + -b; }
    friend two_word operator*(two_word a, two_word b)
    {
      R p = a.hi * b.hi;
      R e = std::fma(a.hi, b.hi, -p);
      e = std::fma(a.hi, b.lo, e);
      e = std::fma(a.lo, b.hi, e);
      e = std::fma(a.lo, b.lo, e);
      return two_word::sum(p, e);
    }
    friend two_word operator/(two_word a, two_word b)
    {
      R q = a.hi / b.hi;
      two_word result(q), r = a - b * result;
      R t = r.hi / b.hi;
      result = result + two_word(t);
      r = a - b * result;
      return result + two_word(r.hi / b.hi);
    }
    two_word scaled(int n) const { return {std::ldexp(hi, n), std::ldexp(lo, n)}; }
    R rounded() const { return hi + lo; }
};
// Leave room for accumulated two-word roundoff: 40/85/96/145 working bits.
template <NativeReal R> R epsilon()
{
  constexpr int p = numeric_limits<R>::digits();
  return std::ldexp(R(1), -p - std::min(p - 8, 32));
}
template <NativeReal R> two_word<R> log_two()
{
  // Generated by scripts/generate-upper-gamma-bounds.py. Each 16-bit chunk
  // converts exactly even at native float32 precision.
  constexpr unsigned chunks[] = {0xb172, 0x17f7, 0xd1cf, 0x79ab, 0xc9e3, 0xb398, 0x03f2, 0xf6af,
                                 0x40f3, 0x4326, 0x7298, 0xb62d, 0x8a0d, 0x175b, 0x8baa, 0xfa2b};
  two_word<R> s;
  constexpr int n = (2 * numeric_limits<R>::digits() + 15) / 16;
  for (int i = 0; i < n; ++i)
    s = s + two_word<R>(std::ldexp(R(chunks[i]), -16 * (i + 1)));
  return s;
}
template <NativeReal R> two_word<R> log(R x)
{
  using P = two_word<R>;
  int e;
  R m = std::frexp(x, &e);
  if (m < R(.75))
  {
    m *= 2;
    --e;
  }
  P t = (P(m) - P(1)) / (P(m) + P(1));
  P t2 = t * t, term = t, s = t;
  for (int k = 1; k < 2 * numeric_limits<R>::digits(); ++k)
  {
    term = term * t2;
    P v = term / P(R(2 * k + 1));
    s = s + v;
    if (std::abs(v.hi) < epsilon<R>() * std::abs(s.hi) / 16) break;
  }
  return P(2) * s + P(R(e)) * log_two<R>();
}
template <NativeReal R> R exp(two_word<R> v)
{
  using L = numeric_limits<R>;
  using P = two_word<R>;
  if (v.hi > R(2 * L::max_exponent)) return L::infinity();
  if (v.hi < R(2 * (L::min_exponent - L::digits()))) return R(0);
  P ln2 = log_two<R>();
  int n = int(std::nearbyint(v.hi / ln2.hi));
  P r = v - P(R(n)) * ln2, s(1), t(1);
  for (int k = 1; k < 2 * L::digits(); ++k)
  {
    t = t * r / P(R(k));
    s = s + t;
    if (std::abs(t.hi) < epsilon<R>() / 16) break;
  }
  // Round subnormal results in integer units, before a native mantissa
  // rounding can erase the deciding low word.
  constexpr int denorm_exp = L::min_exponent - L::digits();
  if (n < L::min_exponent)
  {
    P units = s.scaled(n - denorm_exp);
    R q = std::floor(units.hi);
    P remainder = units - P(q);
    if (remainder.hi < 0)
    {
      q -= 1;
      remainder = remainder + P(1);
    }
    P half = remainder - P(R(.5));
    if (half.hi > 0 || (half.hi == 0 && half.lo > 0) || (half.hi == 0 && half.lo == 0 && std::fmod(q, R(2)) != 0))
      q += 1;
    return std::ldexp(q, denorm_exp);
  }
  return std::ldexp(s.rounded(), n);
}
// Contracted upper-Gamma fraction; positive-input integral bounds justify
// the small-x limiting shortcuts before iteration.
template <NativeReal R> two_word<R> fraction(R av, R xv)
{
  using P = two_word<R>;
  P a(av), x(xv), previous;
  R eps = epsilon<R>();
  if (av < -1 && xv / (-av - 1) < eps / 4) return P(1) / (-a);
  if (av <= R(-.5) && xv < eps * eps / 64) return P(1) / (-a);
  // For b=-a>0, U=integral exp(-x*t)/(1+t)^(b+1) dt. Integrating
  // Taylor's remainder bounds the error after n terms by the next term,
  // x^n/[b*(b-1)*...*(b-n)], provided n<b. This converges quickly near
  // a small-x range boundary where the fraction can converge too slowly.
  if (av < -1 && xv <= 1)
  {
    P term = P(1) / (-a), sum = term;
    for (int n = 1; n < 2 * numeric_limits<R>::digits() && R(n) < -av; ++n)
    {
      P next = -term * x / (-a - P(R(n)));
      if (std::abs(next.hi) <= eps * std::abs(sum.hi) / 16) return sum;
      if (std::abs(next.hi) >= std::abs(term.hi)) break;
      sum = sum + next;
      term = next;
    }
  }
  for (int depth = 8; depth <= 128 * numeric_limits<R>::digits(); depth *= 2)
  {
    P value;
    for (int k = depth; k >= 1; --k)
      value = (a - P(R(k))) * (P(R(k)) / (x - a + P(R(2 * k + 1)) + value));
    value = P(1) / (x - a + P(1) + value);
    if (std::abs((value - previous).hi) <= eps * std::abs(value.hi)) return value;
    previous = value;
  }
  return P(numeric_limits<R>::quiet_NaN());
}
template <NativeReal R> R upper_gamma(R a, R x)
{
  using P = two_word<R>;
  using L = numeric_limits<R>;
  P lx = log(x);
  int shift;
  std::frexp(std::max(R(1), std::abs(lx.hi)), &shift);
  // Keep both terms below the native range before the two-word product:
  // a finite fused difference does not imply that its product is finite.
  // Decide the range while scaled, then restore the bounded exponent.
  P half = P(a / 2) * lx.scaled(-shift) - P(std::ldexp(x / 2, -shift));
  R limit = std::ldexp(R(4 * L::max_exponent), -shift);
  if (half.hi > limit) return L::infinity();
  if (half.hi < -limit) return R(0);
  half = half.scaled(shift);
  P u = fraction(a, x);
  if (std::isnan(u.hi)) return u.hi;
  P lu = log(u.hi) + P(u.lo) / P(u.hi);
  return exp(P(2) * half + lu);
}
template <NativeReal R> two_word<R> log_pair(two_word<R> x)
{
  return log(x.hi) + two_word<R>(x.lo) / two_word<R>(x.hi);
}
template <NativeReal R> two_word<R> log_two_pi()
{
  constexpr unsigned chunks[] = {0xd67f, 0x1c86, 0x4beb, 0x4a69, 0x2979, 0x2002, 0x8832, 0x4047,
                                 0x9f61, 0x1f1a, 0x268b, 0x169b, 0xbd8d, 0x4626, 0x7b54, 0x2aba};
  two_word<R> s(R(1));
  constexpr int n = (2 * numeric_limits<R>::digits() + 15) / 16;
  for (int i = 0; i < n; ++i)
    s = s + two_word<R>(std::ldexp(R(chunks[i]), -16 * (i + 1)));
  return s;
}
template <NativeReal R> two_word<R> log_gamma(R a)
{
  using P = two_word<R>;
  constexpr int p = numeric_limits<R>::digits(),
                cutoff = p <= 24   ? 32
                         : p <= 64 ? 64
                                   : 128,
                terms = p <= 24   ? 4
                        : p <= 53 ? 8
                        : p <= 64 ? 10
                                  : 16;
  constexpr long long nums[] = {1,
                                -1,
                                1,
                                -1,
                                1,
                                -691,
                                1,
                                -3617,
                                43867,
                                -174611,
                                77683,
                                -236364091LL,
                                657931,
                                -3392780147LL,
                                1723168255201LL,
                                -7709321041217LL};
  constexpr long long dens[] = {12,     360,    1260, 1680,    1188, 360360, 156,     122400,
                                244188, 125400, 5796, 1506960, 300,  93960,  2492028, 505920};
  P z(a), prefix;
  while (z.hi < cutoff)
  {
    prefix = prefix + log_pair(z);
    z = z + P(1);
  }
  P inverse = P(1) / z, square = inverse * inverse, c = P(R(nums[terms - 1])) / P(R(dens[terms - 1]));
  for (int k = terms - 2; k >= 0; --k)
    c = c * square + P(R(nums[k])) / P(R(dens[k]));
  return (z - P(R(.5))) * log_pair(z) - z + log_two_pi<R>() / P(2) + inverse * c - prefix;
}
template <NativeReal R> two_word<R> exp_pair(two_word<R> v)
{
  using L = numeric_limits<R>;
  using P = two_word<R>;
  if (v.hi < R(2 * (L::min_exponent - L::digits()))) return P(0);
  P ln2 = log_two<R>();
  int n = int(std::nearbyint(v.hi / ln2.hi));
  P r = v - P(R(n)) * ln2, s(1), t(1);
  for (int k = 1; k < 2 * L::digits(); ++k)
  {
    t = t * r / P(R(k));
    s = s + t;
    if (std::abs(t.hi) < epsilon<R>() / 16) break;
  }
  return s.scaled(n);
}
template <NativeReal R> R lower_complement(R av, R xv)
{
  using P = two_word<R>;
  using L = numeric_limits<R>;
  P a(av), x(xv), g = log_gamma(av), sum(1), term(1);
  for (int k = 1; k <= 128 * L::digits(); ++k)
  {
    term = term * x / (a + P(R(k)));
    sum = sum + term;
    if (std::abs(term.hi) <= epsilon<R>() * std::abs(sum.hi) / 16)
    {
      P lower = exp_pair(a * log(xv) - x - log(av) - g) * sum;
      return exp(g + log_pair(P(1) - lower));
    }
  }
  return L::quiet_NaN();
}
 
} // namespace uni20::detail::native_math::upper_gamma_guard