native_beta_math.hpp

Go to the documentation of this file.

#pragma once
 
#include "native_gamma_math.hpp"
#include <uni20/core/math_constants.hpp>
 
namespace uni20::detail::native_math
{
 
template <NativeReal R> R portable_beta(R a, R b)
{
  using L = numeric_limits<R>;
  if (!(a > 0) || !(b > 0)) return L::quiet_NaN();
  if (std::isinf(a) || std::isinf(b)) return R(0);
  if (a > b) std::swap(a, b);
  if (a == 1) return R(1) / b;
  if (b == 1) return R(1) / a;
  if (a < L::min())
  {
    // Here B(a,b) differs from 1/a+1/b by at most a modest logarithmic
    // correction, far below one result spacing. Ordinary reciprocal addition
    // can still overflow prematurely. Sum half-reciprocals with their FMA
    // division residuals before restoring the scale.
    R qa = R(.5) / a, qb = R(.5) / b, sum = qa + qb;
    if (std::isinf(sum)) return L::infinity();
    R low = qb - (sum - qa);
    low += std::fma(-qa, a, R(.5)) / a + std::fma(-qb, b, R(.5)) / b;
    return R(2) * (sum + low);
  }
  constexpr int bits = L::digits();
  constexpr int cutoff = bits <= 24 ? 8 : bits <= 64 ? 16 : 32;
  // Reduce modest arguments before Gamma evaluation: adding rounded arguments
  // near the cutoff otherwise amplifies their lost low bits through Gamma.
  if (b < cutoff)
  {
    R factor = 1;
    while (b > 2)
    {
      b -= 1;
      factor /= R(1) + a / b;
    }
    while (a > 2)
    {
      a -= 1;
      factor /= R(1) + b / a;
    }
    R ratio = (std::tgamma(R(1) + a) / std::tgamma(R(1) + (a + b))) * std::tgamma(R(1) + b);
    return (R(1) / a + R(1) / b) * (ratio * factor);
  }
  // Reciprocal of the sum avoids overflowing a+b.
  R r = a / b, l = std::log1p(r), invb = R(1) / b, invsum = invb / (R(1) + r);
  R logbeta;
  if (a < cutoff)
  {
    R lgamma = std::log(std::tgamma(R(1) + a));
    // a-b*log1p(a/b) is O(a*a/b); below epsilon its contribution is negligible.
    R remainder = r <= L::epsilon() ? R(0) : std::fma(-b, l, a);
    logbeta = std::fma(-a, std::log(b), lgamma) + remainder - (a - R(.5)) * l +
              native_math::stirling_correction(invb) - native_math::stirling_correction(invsum);
    // Keep 1/a out of the exponent near the final overflow boundary.
    // Here b>=cutoff and a<1, so exp(logbeta) itself remains representable.
    if (a < 1) return std::exp(logbeta) / a;
    logbeta -= std::log(a);
  }
  else
  {
    // Entropy form: never subtract separately rounded, enormous log-Gamma values.
    logbeta = a * std::log(r) - a * l - b * l + (std::log(R(2) * pi<R>) - std::log(a) + l) / 2 +
              native_math::stirling_correction(R(1) / a) + native_math::stirling_correction(invb) -
              native_math::stirling_correction(invsum);
  }
  return std::exp(logbeta);
}
 
} // namespace uni20::detail::native_math