hypergeometric_3F2.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_PRIM_FUN_HYPERGEOMETRIC_3F2_HPP
#define STAN_MATH_PRIM_FUN_HYPERGEOMETRIC_3F2_HPP
 
#include <stan/math/prim/meta.hpp>
#include <stan/math/prim/err.hpp>
#include <stan/math/prim/fun/append_row.hpp>
#include <stan/math/prim/fun/as_array_or_scalar.hpp>
#include <stan/math/prim/fun/to_vector.hpp>
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/prim/fun/fabs.hpp>
#include <stan/math/prim/fun/hypergeometric_pFq.hpp>
#include <stan/math/prim/fun/sum.hpp>
#include <stan/math/prim/fun/sign.hpp>
#include <stan/math/prim/fun/value_of_rec.hpp>
 
namespace stan {
namespace math {
namespace internal {
template <typename Ta, typename Tb, typename Tz,
          require_all_vector_t<Ta, Tb>* = nullptr,
          require_stan_scalar_t<Tz>* = nullptr>
inline return_type_t<Ta, Tb, Tz> hypergeometric_3F2_infsum(
    const Ta& a, const Tb& b, const Tz& z, double precision = 1e-6,
    int max_steps = 1e5) {
  using T_return = return_type_t<Ta, Tb, Tz>;
  Eigen::Array<scalar_type_t<Ta>, 3, 1> a_array = as_array_or_scalar(a);
  Eigen::Array<scalar_type_t<Tb>, 3, 1> b_array
      = append_row(as_array_or_scalar(b), 1.0);
  check_3F2_converges("hypergeometric_3F2", a_array[0], a_array[1], a_array[2],
                      b_array[0], b_array[1], z);
 
  T_return t_acc = 1.0;
  T_return log_t = 0.0;
  auto log_z = log(fabs(z));
  Eigen::Array<int, 3, 1> a_signs = sign(value_of_rec(a_array));
  Eigen::Array<int, 3, 1> b_signs = sign(value_of_rec(b_array));
  int z_sign = sign(value_of_rec(z));
  int t_sign = z_sign * a_signs.prod() * b_signs.prod();
 
  int k = 0;
  const double log_precision = log(precision);
  while (k <= max_steps && log_t >= log_precision) {
    // Replace zero values with 1 prior to taking the log so that we accumulate
    // 0.0 rather than -inf
    const auto& abs_apk = math::fabs((a_array == 0).select(1.0, a_array));
    const auto& abs_bpk = math::fabs((b_array == 0).select(1.0, b_array));
    auto p = sum(log(abs_apk)) - sum(log(abs_bpk));
    if (p == NEGATIVE_INFTY) {
      return t_acc;
    }
 
    log_t += p + log_z;
    t_acc += t_sign * exp(log_t);
 
    if (is_inf(t_acc)) {
      throw_domain_error("hypergeometric_3F2", "sum (output)", t_acc,
                         "overflow hypergeometric function did not converge.");
    }
    k++;
    a_array += 1.0;
    b_array += 1.0;
    a_signs = sign(value_of_rec(a_array));
    b_signs = sign(value_of_rec(b_array));
    t_sign = a_signs.prod() * b_signs.prod() * t_sign;
  }
  if (k == max_steps) {
    throw_domain_error("hypergeometric_3F2", "k (internal counter)", max_steps,
                       "exceeded  iterations, hypergeometric function did not ",
                       "converge.");
  }
  return t_acc;
}
}  // namespace internal
 
template <typename Ta, typename Tb, typename Tz,
          require_all_vector_t<Ta, Tb>* = nullptr,
          require_stan_scalar_t<Tz>* = nullptr>
inline auto hypergeometric_3F2(const Ta& a, const Tb& b, const Tz& z) {
  check_3F2_converges("hypergeometric_3F2", a[0], a[1], a[2], b[0], b[1], z);
  // Boost's pFq throws convergence errors in some cases, fallback to naive
  // infinite-sum approach (tests pass for these)
  if (z == 1.0 && (sum(b) - sum(a)) < 0.0) {
    return internal::hypergeometric_3F2_infsum(a, b, z);
  }
  return hypergeometric_pFq(to_vector(a), to_vector(b), z);
}
 
template <typename Ta, typename Tb, typename Tz,
          require_all_stan_scalar_t<Ta, Tb, Tz>* = nullptr>
inline auto hypergeometric_3F2(const std::initializer_list<Ta>& a,
                               const std::initializer_list<Tb>& b,
                               const Tz& z) {
  return hypergeometric_3F2(std::vector<Ta>(a), std::vector<Tb>(b), z);
}
 
}  // namespace math
}  // namespace stan
#endif