factorial.hpp

Go to the documentation of this file.

#pragma once
 
#include "detail/native_math.hpp"
#include "runtime_precision.hpp"
#if UNI20_HAS_MPREAL
#include "mpreal.hpp"
#endif
 
namespace uni20
{
template <detail::native_math::NativeReal R, detail::MathInteger I> R factorial(I n)
{
  auto order = detail::math_integer(n);
  if (order < 0) throw std::domain_error("factorial requires a nonnegative order");
  using L = numeric_limits<R>;
  if (order > L::max_exponent) return L::infinity(); // n! >= 2^(n-1).
  R high = 1, low = 0;
  for (std::int64_t k = 2; k <= order; ++k)
  {
    R factor = R(k);
    R product = high * factor;
    if (std::isinf(product)) return product;
    // Recover the high product's rounding residual and propagate the low word.
    // FastTwoSum renormalizes the much smaller correction into a new low word.
    R error = std::fma(high, factor, -product) + low * factor;
    R sum = product + error;
    low = error - (sum - product);
    high = sum;
  }
  return high;
}
 
template <detail::native_math::NativeReal R, detail::MathInteger I> R factorial(I n, Precision p)
{
  require_precision<R>(p);
  return uni20::factorial<R>(n);
}
 
#if UNI20_HAS_MPREAL
template <class R, detail::MathInteger I>
  requires std::same_as<R, mpreal>
R factorial(I n, Precision p)
{
  return uni20::factorial(n, p);
}
#endif
} // namespace uni20