mpreal_special.hpp

Go to the documentation of this file.

#pragma once
 
#include "detail/exact_bernoulli.hpp"
#include "detail/math_integer.hpp"
#include "math_results.hpp"
#include "mpreal.hpp"
#include <vector>
 
namespace uni20
{
namespace detail
{
inline bool rational_integer(exact_rational const& x) { return mpz_cmp_ui(mpq_denref(x.native_handle()), 1) == 0; }
inline unsigned long exact_nonnegative_order(exact_rational const& x)
{
  auto q = x.native_handle();
  if (!rational_integer(x) || mpq_sgn(q) < 0)
    throw std::domain_error("special function requires a nonnegative integer order");
  if (!mpz_fits_slong_p(mpq_numref(q))) throw std::out_of_range("special function order does not fit int64_t");
  return mpz_get_ui(mpq_numref(q));
}
inline exact_rational exact_factorial(unsigned long n)
{
  return mpreal_access::rational_result([&](mpq_ptr out) { mpz_fac_ui(mpq_numref(out), n); });
}
inline mpreal exact_gamma(exact_rational const& x)
{
  if (x > 0 && rational_integer(x)) return mpreal(exact_factorial(exact_nonnegative_order(x - exact_rational{1})));
  throw std::logic_error("Gamma requires finite precision unless its input is a positive integer");
}
inline mpreal exact_zeta(exact_rational const& x)
{
  if (x == 0) return mpreal(exact_rational{-1} / exact_rational{2});
  if (x < 0 && rational_integer(x))
  {
    if (mpz_even_p(mpq_numref(x.native_handle()))) return mpreal{};
    // Zeta(-m) = -B_(m+1)/(m+1) for odd positive m. This exact
    // Bernoulli recurrence has no floating-point approximation or ambient precision.
    auto order = exact_nonnegative_order(-x) + 1;
    return mpreal(-exact_bernoulli(order) / exact_rational{order});
  }
  throw std::logic_error("Zeta requires finite precision except at nonpositive integers");
}
 
template <auto Operation> mpreal special_unary(mpreal const& x)
{
  auto p = x.precision();
  if (p.is_exact())
  {
    auto const& q = x.exact_value();
    if constexpr (same_function_v<Operation, mpfr_gamma>)
      return exact_gamma(q);
    else if constexpr (same_function_v<Operation, mpfr_zeta>)
      return exact_zeta(q);
    else if constexpr (same_function_v<Operation, mpfr_erf> || same_function_v<Operation, mpfr_li2>)
    {
      if (q == 0) return mpreal{};
    }
    else if constexpr (same_function_v<Operation, mpfr_erfc>)
    {
      if (q == 0) return mpreal{1};
    }
    throw std::logic_error("special function requires finite working precision");
  }
  return mpreal_access::finite_result(p, [&](mpfr_ptr out) { Operation(out, x.native_handle(), MPFR_RNDN); });
}
template <auto Operation> mpreal special_unary(mpreal const& x, Precision p)
{
  (void)p.bit_count();
  return special_unary<Operation>(x.at(p));
}
 
inline mpreal exact_beta(exact_rational const& a, exact_rational const& b)
{
  if (a > 0 && b > 0 && (rational_integer(a) || rational_integer(b)))
  {
    // If both are integers, use the shorter recurrence. In particular B(a,1)
    // needs no large-order conversion even when a is an arbitrary-size integer.
    bool const use_a = rational_integer(a) && (!rational_integer(b) || a <= b);
    auto n = exact_nonnegative_order(use_a ? a : b);
    auto x = use_a ? b : a;
    // B(x,n) = (n-1)! / (x (x+1) ... (x+n-1)). The recurrence
    // avoids computing three enormous factorials for small rational arguments.
    exact_rational result = exact_rational{1} / x;
    for (unsigned long k = 1; k < n; ++k)
      result = result * exact_rational{k} / (x + exact_rational{k});
    return mpreal(result);
  }
  throw std::logic_error("Beta requires finite precision except for positive rational/integer pairs");
}
template <auto Operation> mpreal special_binary(mpreal const& a, mpreal const& b)
{
  auto p = common_precision(a.precision(), b.precision());
  if (p.is_exact())
  {
    if constexpr (same_function_v<Operation, mpfr_beta>)
      return exact_beta(a.exact_value(), b.exact_value());
    else if constexpr (same_function_v<Operation, mpfr_gamma_inc>)
    {
      if (b == 0) return exact_gamma(a.exact_value());
    }
    else if constexpr (same_function_v<Operation, mpfr_agm>)
    {
      if (a >= 0 && b >= 0)
      {
        if (a == b) return a;
        if (a == 0 || b == 0) return mpreal{};
      }
    }
    throw std::logic_error("special function requires finite working precision");
  }
  auto av = a.at(p), bv = b.at(p);
  return mpreal_access::finite_result(
      p, [&](mpfr_ptr out) { Operation(out, av.native_handle(), bv.native_handle(), MPFR_RNDN); });
}
template <auto Operation> mpreal special_binary(mpreal const& a, mpreal const& b, Precision p)
{
  (void)p.bit_count();
  return special_binary<Operation>(a.at(p), b.at(p));
}
} // namespace detail
 
inline mpreal tgamma(mpreal const& x) { return detail::special_unary<mpfr_gamma>(x); }
inline mpreal tgamma(mpreal const& x, Precision p) { return detail::special_unary<mpfr_gamma>(x, p); }
 
inline mpreal digamma(mpreal const& x) { return detail::special_unary<mpfr_digamma>(x); }
inline mpreal digamma(mpreal const& x, Precision p) { return detail::special_unary<mpfr_digamma>(x, p); }
 
inline mpreal erf(mpreal const& x) { return detail::special_unary<mpfr_erf>(x); }
inline mpreal erf(mpreal const& x, Precision p) { return detail::special_unary<mpfr_erf>(x, p); }
 
inline mpreal erfc(mpreal const& x) { return detail::special_unary<mpfr_erfc>(x); }
inline mpreal erfc(mpreal const& x, Precision p) { return detail::special_unary<mpfr_erfc>(x, p); }
 
inline mpreal zeta(mpreal const& x) { return detail::special_unary<mpfr_zeta>(x); }
inline mpreal zeta(mpreal const& x, Precision p) { return detail::special_unary<mpfr_zeta>(x, p); }
 
inline mpreal expint(mpreal const& x) { return detail::special_unary<mpfr_eint>(x); }
inline mpreal expint(mpreal const& x, Precision p) { return detail::special_unary<mpfr_eint>(x, p); }
 
inline mpreal dilog_real(mpreal const& x) { return detail::special_unary<mpfr_li2>(x); }
inline mpreal dilog_real(mpreal const& x, Precision p) { return detail::special_unary<mpfr_li2>(x, p); }
 
inline mpreal airy_ai(mpreal const& x) { return detail::special_unary<mpfr_ai>(x); }
inline mpreal airy_ai(mpreal const& x, Precision p) { return detail::special_unary<mpfr_ai>(x, p); }
 
inline mpreal beta(mpreal const& a, mpreal const& x) { return detail::special_binary<mpfr_beta>(a, x); }
inline mpreal beta(mpreal const& a, mpreal const& x, Precision p) { return detail::special_binary<mpfr_beta>(a, x, p); }
 
inline mpreal upper_gamma(mpreal const& a, mpreal const& x) { return detail::special_binary<mpfr_gamma_inc>(a, x); }
inline mpreal upper_gamma(mpreal const& a, mpreal const& x, Precision p)
{
  return detail::special_binary<mpfr_gamma_inc>(a, x, p);
}
 
inline mpreal agm(mpreal const& a, mpreal const& x) { return detail::special_binary<mpfr_agm>(a, x); }
inline mpreal agm(mpreal const& a, mpreal const& x, Precision p) { return detail::special_binary<mpfr_agm>(a, x, p); }
 
inline lgamma_result<mpreal> lgamma_sign(mpreal const& x)
{
  auto p = x.precision();
  if (p.is_exact())
  {
    if (x == 1 || x == 2) return {mpreal{}, 1};
    throw std::logic_error("log-Gamma requires finite precision except at one and two");
  }
  int sign = 0;
  auto value = detail::mpreal_access::finite_result(
      p, [&](mpfr_ptr out) { mpfr_lgamma(out, &sign, x.native_handle(), MPFR_RNDN); });
  if (isnan(x) || (x < 0 && (isinf(x) || mpfr_integer_p(x.native_handle())))) sign = 0;
  return {std::move(value), sign};
}
inline lgamma_result<mpreal> lgamma_sign(mpreal const& x, Precision p)
{
  (void)p.bit_count();
  return lgamma_sign(x.at(p));
}
inline mpreal lgamma(mpreal const& x) { return lgamma_sign(x).value; }
inline mpreal lgamma(mpreal const& x, Precision p) { return lgamma_sign(x, p).value; }
 
namespace detail
{
template <auto Operation> mpreal integer_bessel(mpreal const& x, std::int64_t order)
{
  auto p = x.precision();
  if (p.is_exact())
  {
    if constexpr (same_function_v<Operation, mpfr_jn>)
      if (x == 0) return mpreal{order == 0 ? 1 : 0};
    throw std::logic_error("Bessel function requires finite working precision");
  }
  static_assert(sizeof(long) >= sizeof(std::int64_t));
  return mpreal_access::finite_result(
      p, [&](mpfr_ptr out) { Operation(out, static_cast<long>(order), x.native_handle(), MPFR_RNDN); });
}
} // namespace detail
 
template <detail::MathInteger I> inline mpreal bessel_j(I n, mpreal const& x)
{
  return detail::integer_bessel<mpfr_jn>(x, detail::math_integer(n));
}
template <detail::MathInteger I> inline mpreal bessel_j(I n, mpreal const& x, Precision p)
{
  (void)p.bit_count();
  return bessel_j(n, x.at(p));
}
 
template <detail::MathInteger I> inline mpreal bessel_y(I n, mpreal const& x)
{
  return detail::integer_bessel<mpfr_yn>(x, detail::math_integer(n));
}
template <detail::MathInteger I> inline mpreal bessel_y(I n, mpreal const& x, Precision p)
{
  (void)p.bit_count();
  return bessel_y(n, x.at(p));
}
 
template <detail::MathInteger I> inline mpreal factorial(I n, Precision p)
{
  auto order = detail::math_integer(n);
  if (order < 0) throw std::domain_error("factorial requires a nonnegative order");
  if (p.is_exact()) return mpreal(detail::exact_factorial(static_cast<unsigned long>(order)));
  return detail::mpreal_access::finite_result(
      p, [&](mpfr_ptr out) { mpfr_fac_ui(out, static_cast<unsigned long>(order), MPFR_RNDN); });
}
} // namespace uni20