From 47fa45bee4c6daa6e948efd99812cbce688af5c1 Mon Sep 17 00:00:00 2001 From: jzmaddock Date: Sun, 29 Oct 2017 19:33:39 +0000 Subject: [PATCH] Import SOC headers, Test and fix up 1F0 and 0F1. --- .../detail/hypergeometric_0f1_bessel.hpp | 45 +++ .../detail/hypergeometric_1f1_bessel.hpp | 312 ++++++++++++++++++ .../detail/hypergeometric_1f1_recurrence.hpp | 200 +++++++++++ .../detail/hypergeometric_asym.hpp | 87 +++++ .../detail/hypergeometric_cf.hpp | 228 +++++++++++++ .../detail/hypergeometric_pade.hpp | 131 ++++++++ .../detail/hypergeometric_rational.hpp | 165 +++++++++ .../hypergeometric_separated_series.hpp | 51 +++ .../detail/hypergeometric_series.hpp | 242 ++++++++++++++ .../special_functions/hypergeometric_0F1.hpp | 110 ++++++ .../special_functions/hypergeometric_1F0.hpp | 73 ++++ test/test_0F1.cpp | 85 +++++ test/test_0F1.hpp | 197 +++++++++++ test/test_1F0.cpp | 60 ++++ test/test_1F0.hpp | 60 ++++ 15 files changed, 2046 insertions(+) create mode 100644 include/boost/math/special_functions/detail/hypergeometric_0f1_bessel.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_1f1_recurrence.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_asym.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_cf.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_pade.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_rational.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_separated_series.hpp create mode 100644 include/boost/math/special_functions/detail/hypergeometric_series.hpp create mode 100644 include/boost/math/special_functions/hypergeometric_0F1.hpp create mode 100644 include/boost/math/special_functions/hypergeometric_1F0.hpp create mode 100644 test/test_0F1.cpp create mode 100644 test/test_0F1.hpp create mode 100644 test/test_1F0.cpp create mode 100644 test/test_1F0.hpp diff --git a/include/boost/math/special_functions/detail/hypergeometric_0f1_bessel.hpp b/include/boost/math/special_functions/detail/hypergeometric_0f1_bessel.hpp new file mode 100644 index 000000000..04a523045 --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_0f1_bessel.hpp @@ -0,0 +1,45 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP + #define BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP + + namespace boost { namespace math { namespace detail { + + template + inline T hypergeometric_0f1_bessel(const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING + + const bool is_z_nonpositive = z <= 0; + + const T sqrt_z = is_z_nonpositive ? T(sqrt(-z)) : T(sqrt(z)); + const T bessel_mult = is_z_nonpositive ? + boost::math::cyl_bessel_j(b - 1, 2 * sqrt_z, pol) : + boost::math::cyl_bessel_i(b - 1, 2 * sqrt_z, pol) ; + + if (b > boost::math::max_factorial::value) + { + const T lsqrt_z = log(sqrt_z); + const T lsqrt_z_pow_b = (b - 1) * lsqrt_z; + T lg = (boost::math::lgamma(b, pol) - lsqrt_z_pow_b); + lg = exp(lg); + return lg * bessel_mult; + } + else + { + const T sqrt_z_pow_b = pow(sqrt_z, b - 1); + return (boost::math::tgamma(b, pol) / sqrt_z_pow_b) * bessel_mult; + } + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_0F1_BESSEL_HPP diff --git a/include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp b/include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp new file mode 100644 index 000000000..ad661595d --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_1f1_bessel.hpp @@ -0,0 +1,312 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP + #define BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP + + #include + #include + #include + + namespace boost { namespace math { namespace detail { + + // declarations of helpers for 13_3_7 and 13_3_8: + // bessel_j recurrence relation based on finite continued fractions; + // it is stable forward recurrence algorithm, see William J. Lentz + // "Generating Bessel functions in Mie scattering calculations + // using continued fractions" for details + template + inline T hypergeometric_bessel_j_recurrence_next(const T& jvm1, const T& v, const T& z); + + // swap values of Bessel functions: + template + inline void hypergeometric_bessel_j_recurrence_iterate(T& jvm1, T& jv, const T& v, const T& z); + + // next coefficient for 13_3_7 + template + inline T hypergeometric_13_3_7_coefficient_next(const T& anm3, const T& anm2, const T& a, const T& b, const unsigned n); + + // next coefficient for 13_3_8 + template + inline T hypergeometric_13_3_8_coefficient_next(const T& cnm3, const T& cnm2, const T& cnm1, const T& a, const T& b, const unsigned n); + + // swap values of coefficients: + template + inline void hypergeometric_coefficient_13_3_7_iterate(T& anm3, T& anm2, T& anm1, T& an, const T& a, const T& b, const unsigned n); + + template + inline void hypergeometric_coefficient_13_3_8_iterate(T& cnm3, T& cnm2, T& cnm1, T& cn, const T& a, const T& b, const unsigned n); + + // term class of Abramowitz & Stegun 13_3_7 formula + template + struct hypergeometric_1f1_13_3_7_series_term + { + typedef T result_type; + + hypergeometric_1f1_13_3_7_series_term(const T& a, const T& b, const T& z): + a(a), b(b), z(z), n(0u) + { + BOOST_MATH_STD_USING + + sqrt_z_pow_n = sqrt_z = sqrt(z); + sqrt_2b_minus_4a_pow_n = sqrt_2b_minus_4a = sqrt(2 * (b - (2 * a))); + sqrt_2zb_minus_4za = sqrt_z * sqrt_2b_minus_4a; + + v_current = b - 1; + + anm3 = 1; + anm2 = 0; + anm1 = b / 2; + an = detail::hypergeometric_13_3_7_coefficient_next(anm3, anm2, a, b, 3u); + + jvm1 = boost::math::cyl_bessel_j(v_current, sqrt_2zb_minus_4za); + v_current++; + jv = detail::hypergeometric_bessel_j_recurrence_next(jvm1, v_current, sqrt_2zb_minus_4za); + + term = jvm1; + ++n; + } + + T operator()() + { + const T result = term; + + iterate(); + term = ((anm2 * sqrt_z_pow_n) / sqrt_2b_minus_4a_pow_n) * jv; + + return result; + } + + private: + void iterate() + { + ++n; ++v_current; + sqrt_z_pow_n *= sqrt_z; + sqrt_2b_minus_4a_pow_n *= sqrt_2b_minus_4a; + + detail::hypergeometric_coefficient_13_3_7_iterate(anm3, anm2, anm1, an, a, b, n + 2); + detail::hypergeometric_bessel_j_recurrence_iterate(jvm1, jv, v_current, sqrt_2zb_minus_4za); + } + + const T a, b, z; + unsigned n; + T sqrt_z, sqrt_2b_minus_4a, sqrt_2zb_minus_4za; + T sqrt_z_pow_n, sqrt_2b_minus_4a_pow_n; + T v_current; + T anm3, anm2, anm1, an; + T jvm1, jv; + T term; + }; + + // function for 13_3_7 evaluation + template + inline T hypergeometric_1f1_13_3_7_series(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING + + const T sqrt_bz_div_2_minus_az = sqrt(((b * z) / 2) - (a * z)); + const T prefix = ((boost::math::tgamma(b, pol) * sqrt_bz_div_2_minus_az) / + pow(sqrt_bz_div_2_minus_az, b)) * exp(z / 2); + + detail::hypergeometric_1f1_13_3_7_series_term s(a, b, z); + boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations(); +#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) + T zero = 0; + T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter, zero); +#else + T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter); +#endif + boost::math::policies::check_series_iterations("boost::math::hypergeometric_1f1_13_3_7_series<%1%>(%1%,%1%,%1%)", max_iter, pol); + return prefix * result; + } + + // term class of Abramowitz & Stegun 13_3_8 formula + template + struct hypergeometric_1f1_13_3_8_series_term + { + typedef T result_type; + + static const T h; + + hypergeometric_1f1_13_3_8_series_term(const T& a, const T& b, const T& z): + a(a), b(b), z(z), n(0u) + { + BOOST_MATH_STD_USING + + sqrt_minus_az = sqrt(-a * z); + const T double_sqrt_minus_az = 2 * sqrt_minus_az; + + v_current = b - 1; + z_pow_n = z; sqrt_minus_az_pow_n = sqrt_minus_az; + + cnm3 = 1; + cnm2 = -b * h; + cnm1 = (((1 - (2 * h)) * a) + ((b * (b + 1)) * (h * h))) / 2; + cn = detail::hypergeometric_13_3_8_coefficient_next(cnm3, cnm2, cnm1, a, b, 3u); + + jvm1 = boost::math::cyl_bessel_j(v_current, double_sqrt_minus_az); + ++v_current; + jv = detail::hypergeometric_bessel_j_recurrence_next(jvm1, b, double_sqrt_minus_az); + } + + T operator()() + { + ++n; + switch (n - 1) + { + case 0u: + return jvm1; + case 1u: + detail::hypergeometric_bessel_j_recurrence_iterate(jvm1, jv, v_current, T(2 * sqrt_minus_az)); + return ((cnm2 * z) / sqrt_minus_az) * jvm1; + } + + ++v_current; + z_pow_n *= z; + sqrt_minus_az_pow_n *= sqrt_minus_az; + const T result = ((cnm1 * z_pow_n) / sqrt_minus_az_pow_n) * jv; + + detail::hypergeometric_coefficient_13_3_8_iterate(cnm3, cnm2, cnm1, cn, a, b, n + 1); + detail::hypergeometric_bessel_j_recurrence_iterate(jvm1, jv, v_current, T(2 * sqrt_minus_az)); + return result; + } + + private: + const T a, b, z; + unsigned n; + T sqrt_minus_az; + T z_pow_n, sqrt_minus_az_pow_n; + T v_current; + T cnm3, cnm2, cnm1, cn; + T jvm1, jv; + }; + + template + const T hypergeometric_1f1_13_3_8_series_term::h = -boost::math::constants::pi() / T(10.); + + // function for 13_3_8 evaluation + template + inline T hypergeometric_1f1_13_3_8_series(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING + + const T sqrt_minus_az_pow_b_minus_one = pow(sqrt(-a * z), b - 1); + const T prefix = (boost::math::tgamma(b, pol) / sqrt_minus_az_pow_b_minus_one) * + exp(detail::hypergeometric_1f1_13_3_8_series_term::h * z); + + detail::hypergeometric_1f1_13_3_8_series_term s(a, b, z); + boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations(); +#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) + T zero = 0; + T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter, zero); +#else + T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon(), max_iter); +#endif + boost::math::policies::check_series_iterations("boost::math::hypergeometric_1f1_13_3_8_series<%1%>(%1%,%1%,%1%)", max_iter, pol); + return prefix * result; + } + + // definitions of helpers: + + template + inline T hypergeometric_bessel_j_recurrence_next(const T& jvm1, const T& v, const T& z) + { + BOOST_MATH_STD_USING + + // TODO: we can estimate a possible capacity for these vectors + std::vector finite_cf_nums; + std::vector finite_cf_denoms; + + const T a_coef_when_i_is_one = (2 * v) / z; + const T a_coef_when_i_is_two = (-2 * (v + 1)) / z; + + finite_cf_nums.push_back(a_coef_when_i_is_one); + finite_cf_nums.push_back(a_coef_when_i_is_two + (1. / a_coef_when_i_is_one)); + + finite_cf_denoms.push_back(a_coef_when_i_is_two); + + unsigned i = 3u; + + do { + const T temp = (2 * (v + (i - 1))) / z; + const T next_coef = i & 1 ? temp : T(-temp); + + finite_cf_nums.push_back(next_coef + (1. / finite_cf_nums.back())); + finite_cf_denoms.push_back(next_coef + (1. / finite_cf_denoms.back())); + + ++i; + } while (fabs(finite_cf_nums.back() - finite_cf_denoms.back()) > boost::math::tools::epsilon()); + + T ratio = finite_cf_nums.front(); + + for (i = 0u; i < finite_cf_denoms.size(); ++i) + ratio *= finite_cf_nums[i+1] / finite_cf_denoms[i]; + + return jvm1 / ratio; + } + + template + inline void hypergeometric_bessel_j_recurrence_iterate(T& jvm1, T& jv, const T& v, const T& z) + { + using std::swap; + + swap(jvm1, jv); + jv = detail::hypergeometric_bessel_j_recurrence_next(jvm1, v, z); + } + + template + inline T hypergeometric_13_3_7_coefficient_next(const T& anm3, const T& anm2, const T& a, const T& b, const unsigned n) + { + const T term_anm2 = ((n + b) - 2) * anm2; + const T term_anm3 = ((2 * a) - b) * anm3; + + return (term_anm2 + term_anm3) / n; + + } + + template + inline T hypergeometric_13_3_8_coefficient_next(const T& cnm3, const T& cnm2, const T& cnm1, const T& a, const T& b, const unsigned n) + { + static const T& h = detail::hypergeometric_1f1_13_3_8_series_term::h; + const T one_minus_two_h = 1 - (2 * h); + const T h_minus_one = h - 1; + + const T term_cnm1 = ((one_minus_two_h * n) - (b * h)) * cnm1; + const T term_cnm2 = ((one_minus_two_h * a) - ((h * h_minus_one) * (b + (n - 1)))) * cnm2; + const T term_cnm3 = ((-h * h_minus_one) * a) * cnm3; + + return ((term_cnm1 + term_cnm2) + term_cnm3) / (n + 1); + } + + template + inline void hypergeometric_coefficient_13_3_7_iterate(T& anm3, T& anm2, T& anm1, T& an, const T& a, const T& b, const unsigned n) + { + using std::swap; + + swap(anm3, anm2); + swap(anm2, anm1); + swap(anm1, an); + an = detail::hypergeometric_13_3_7_coefficient_next(anm3, anm2, a, b, n); + } + + template + inline void hypergeometric_coefficient_13_3_8_iterate(T& cnm3, T& cnm2, T& cnm1, T& cn, const T& a, const T& b, const unsigned n) + { + using std::swap; + + swap(cnm3, cnm2); + swap(cnm2, cnm1); + swap(cnm1, cn); + cn = detail::hypergeometric_13_3_8_coefficient_next(cnm3, cnm2, cnm1, a, b, n); + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP diff --git a/include/boost/math/special_functions/detail/hypergeometric_1f1_recurrence.hpp b/include/boost/math/special_functions/detail/hypergeometric_1f1_recurrence.hpp new file mode 100644 index 000000000..7482c3cae --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_1f1_recurrence.hpp @@ -0,0 +1,200 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_ + #define BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_ + + #include + #include + + #include + + namespace boost { namespace math { namespace detail { + + template + struct hypergeometric_1f1_recurrence_a_coefficients + { + typedef boost::math::tuple result_type; + + hypergeometric_1f1_recurrence_a_coefficients(const T& a, const T& b, const T& z): + a(a), b(b), z(z) + { + } + + result_type operator()(boost::intmax_t i) const + { + const T ai = a + i; + + const T an = -ai; + const T bn = (b - (2 * ai)) - z; + const T cn = b - ai; + + return boost::math::make_tuple(an, bn, cn); + } + + private: + const T a, b, z; + }; + + template + struct hypergeometric_1f1_recurrence_b_coefficients + { + typedef boost::math::tuple result_type; + + hypergeometric_1f1_recurrence_b_coefficients(const T& a, const T& b, const T& z): + a(a), b(b), z(z) + { + } + + result_type operator()(boost::intmax_t i) const + { + const T bi = b + i; + + const T an = z * (bi - a); + const T bn = bi * ((z + bi) - 1); + const T cn = bi * (bi - 1); + + return boost::math::make_tuple(an, bn, cn); + } + + private: + const T a, b, z; + }; + + template + struct hypergeometric_1f1_recurrence_a_and_b_coefficients + { + typedef boost::math::tuple result_type; + + hypergeometric_1f1_recurrence_a_and_b_coefficients(const T& a, const T& b, const T& z): + a(a), b(b), z(z) + { + } + + result_type operator()(boost::intmax_t i) const + { + const T ai = a + i; + const T bi = b + i; + + const T an = ai * z; + const T bn = bi * ((1 - bi) + z); + const T cn = bi * (1 - bi); + + return boost::math::make_tuple(an, bn, cn); + } + + private: + const T a, b, z; + }; + + // forward declaration for initial values + template + inline T hypergeometric_1f1_imp(const T& a, const T& b, const T& z, const Policy& pol); + + template + inline T hypergeometric_1f1_backward_recurrence_for_negative_a(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING // modf, frexp, fabs, pow + + boost::intmax_t integer_part = 0; + const T bk = modf(b, &integer_part); + T ak = modf(a, &integer_part); + + int exp_of_a = 0; frexp(a, &exp_of_a); + int exp_of_b = 0; frexp(b, &exp_of_b); + + const bool are_fractional_parts_close_enough = + fabs(boost::math::float_distance(ak, bk)) <= pow(2, (std::max)(exp_of_a, exp_of_b)); + + if ((a < b) && (b < 0) && (are_fractional_parts_close_enough)) // TODO: has to be researched deeper + { + ak = b - 1; + integer_part -= (boost::math::lltrunc(ceil(b)) - 1); + } + + T first = detail::hypergeometric_1f1_imp(ak, b, z, pol); + --ak; + T second = detail::hypergeometric_1f1_imp(ak, b, z, pol); + + detail::hypergeometric_1f1_recurrence_a_coefficients s(ak, b, z); + + return tools::solve_recurrence_relation_backward(s, + static_cast(std::abs(integer_part)), + first, + second); + } + + template + inline T hypergeometric_1f1_forward_recurrence_for_positive_a(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING // modf, fabs + + boost::intmax_t integer_part = 0; + T ak = modf(a, &integer_part); + + T first = detail::hypergeometric_1f1_imp(ak, b, z, pol); + ++ak; + T second = detail::hypergeometric_1f1_imp(ak, b, z, pol); + + detail::hypergeometric_1f1_recurrence_a_coefficients s(ak, b, z); + + return tools::solve_recurrence_relation_forward(s, integer_part, first, second); + } + + template + inline T hypergeometric_1f1_backward_recurrence_for_negative_b(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING // modf, fabs + + boost::intmax_t integer_part = 0; + T bk = modf(b, &integer_part); + + T first = detail::hypergeometric_1f1_imp(a, bk, z, pol); + --bk; + T second = detail::hypergeometric_1f1_imp(a, bk, z, pol); + + detail::hypergeometric_1f1_recurrence_b_coefficients s(a, bk, z); + + return tools::solve_recurrence_relation_backward(s, + static_cast(std::abs(integer_part)), + first, + second); + } + + // this method works provided that integer part of a is the same as integer part of b + // we don't make this check inside it. + template + inline T hypergeometric_1f1_backward_recurrence_for_negative_a_and_b(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING // modf, fabs + + boost::intmax_t integer_part = 0; + T ak = modf(a, &integer_part); + T bk = modf(b, &integer_part); + + T first = detail::hypergeometric_1f1_imp(ak, bk, z, pol); + --ak; --bk; + T second = detail::hypergeometric_1f1_imp(ak, bk, z, pol); + + detail::hypergeometric_1f1_recurrence_a_and_b_coefficients s(ak, bk, z); + + return tools::solve_recurrence_relation_backward(s, fabs(integer_part), first, second); + } + + // ranges + template + inline bool hypergeometric_1f1_is_a_small_enough(const T& a) + { + return a < -10; // TODO: make dependent on precision + } + + } } } // namespaces + +#endif // BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_ diff --git a/include/boost/math/special_functions/detail/hypergeometric_asym.hpp b/include/boost/math/special_functions/detail/hypergeometric_asym.hpp new file mode 100644 index 000000000..e515a89a6 --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_asym.hpp @@ -0,0 +1,87 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP + #define BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP + + #include + + namespace boost { namespace math { + + // forward declaration of used function 2f0 + template + inline typename tools::promote_args::type hypergeometric_2f0(T1 a1, T2 a2, T3 z, const Policy& /* pol */); + + namespace detail { + + // assumes a and b are not non-positive integers + template + inline T hypergeometric_1f1_asym_positive_series(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING + static const char* const function = "boost::math::hypergeometric_1f1_asym_positive_series<%1%>(%1%,%1%,%1%)"; + + const bool is_a_integer = (a == floor(a)); + const bool is_b_integer = (b == floor(b)); + + BOOST_ASSERT(!(is_a_integer && a <= 0) && !(is_b_integer && b <= 0)); + + const T gamma_ratio = (b > 0 && a > 0) ? + boost::math::tgamma_ratio(b, a, pol) : + T(boost::math::tgamma(b, pol) / boost::math::tgamma(a, pol)); + const T prefix_a = (exp(z) * gamma_ratio) * pow(z, (a - b)); + + return prefix_a * boost::math::hypergeometric_2f0(b - a, 1 - a, 1 / z, pol); + } + + // assumes b and (b - a) are not non-positive integers + template + inline T hypergeometric_1f1_asym_negative_series(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING + static const char* const function = "boost::math::hypergeometric_1f1_asym_negative_series<%1%>(%1%,%1%,%1%)"; + + const T b_minus_a = b - a; + + const bool is_a_integer = (a == floor(a)); + const bool is_b_minus_a_integer = (b_minus_a == floor(b_minus_a)); + + BOOST_ASSERT(!(is_a_integer && a <= 0) && !(is_b_minus_a_integer && b_minus_a <= 0)); + + const T gamma_ratio = (b > 0 && b_minus_a > 0) ? + boost::math::tgamma_ratio(b, b_minus_a, pol) : + boost::math::tgamma(b) / boost::math::tgamma((b_minus_a), pol); + const T prefix_b = gamma_ratio / pow(-z, a); + + return prefix_b * boost::math::hypergeometric_2f0(a, 1 - b_minus_a, -1 / z, pol); + } + + // experimental range + template + inline bool hypergeometric_1f1_asym_region(const T& a, const T& b, const T& z) + { + BOOST_MATH_STD_USING + + const T the_max_of_one_and_b_minus_a ((std::max)(T(1), fabs(b - a))); + const T the_max_of_one_and_one_minus_a((std::max)(T(1), fabs(1 - a))); + + const T the_product_of_these_maxima(the_max_of_one_and_b_minus_a * the_max_of_one_and_one_minus_a); + + const T abs_of_z(fabs(z)); + + if ((abs_of_z > 100) && (the_product_of_these_maxima < (abs_of_z / 2))) // TODO: not a good way + return true; + + return false; + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP diff --git a/include/boost/math/special_functions/detail/hypergeometric_cf.hpp b/include/boost/math/special_functions/detail/hypergeometric_cf.hpp new file mode 100644 index 000000000..bac89f1f4 --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_cf.hpp @@ -0,0 +1,228 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_CF_HPP + #define BOOST_MATH_HYPERGEOMETRIC_CF_HPP + + namespace boost { namespace math { namespace detail { + + // primary template for term of continued fraction + template + struct hypergeometric_pfq_cf_term; + + // partial specialization for 0f1 + template + struct hypergeometric_pfq_cf_term + { + typedef std::pair result_type; + + hypergeometric_pfq_cf_term(const T& b, const T& z): + n(1), b(b), z(z), + term(std::make_pair(T(0), T(1))) + { + } + + result_type operator()() + { + const result_type result = term; + ++b; ++n; + numer = -(z / (b * n)); + term = std::make_pair(numer, 1 - numer); + return result; + } + + private: + unsigned n; + T b; + const T z; + T numer; + result_type term; + }; + + // partial specialization for 1f0 + template + struct hypergeometric_pfq_cf_term + { + typedef std::pair result_type; + + hypergeometric_pfq_cf_term(const T& a, const T& z): + n(1), a(a), z(z), + term(std::make_pair(T(0), T(1))) + { + } + + result_type operator()() + { + const result_type result = term; + ++a; ++n; + numer = -((a * z) / n); + term = std::make_pair(numer, 1 - numer); + return result; + } + + private: + unsigned n; + T a; + const T z; + T numer; + result_type term; + }; + + // partial specialization for 1f1 + template + struct hypergeometric_pfq_cf_term + { + typedef std::pair result_type; + + hypergeometric_pfq_cf_term(const T& a, const T& b, const T& z): + n(1), a(a), b(b), z(z), + term(std::make_pair(T(0), T(1))) + { + } + + result_type operator()() + { + const result_type result = term; + ++a; ++b; ++n; + numer = -((a * z) / (b * n)); + term = std::make_pair(numer, 1 - numer); + return result; + } + + private: + unsigned n; + T a, b; + const T z; + T numer; + result_type term; + }; + + // partial specialization for 1f2 + template + struct hypergeometric_pfq_cf_term + { + typedef std::pair result_type; + + hypergeometric_pfq_cf_term(const T& a, const T& b, const T& c, const T& z): + n(1), a(a), b(b), c(c), z(z), + term(std::make_pair(T(0), T(1))) + { + } + + result_type operator()() + { + const result_type result = term; + ++a; ++b; ++c; ++n; + numer = -((a * z) / ((b * c) * n)); + term = std::make_pair(numer, 1 - numer); + return result; + } + + private: + unsigned n; + T a, b, c; + const T z; + T numer; + result_type term; + }; + + // partial specialization for 2f1 + template + struct hypergeometric_pfq_cf_term + { + typedef std::pair result_type; + + hypergeometric_pfq_cf_term(const T& a, const T& b, const T& c, const T& z): + n(1), a(a), b(b), c(c), z(z), + term(std::make_pair(T(0), T(1))) + { + } + + result_type operator()() + { + const result_type result = term; + ++a; ++b; ++c; ++n; + numer = -(((a * b) * z) / (c * n)); + term = std::make_pair(numer, 1 - numer); + return result; + } + + private: + unsigned n; + T a, b, c; + const T z; + T numer; + result_type term; + }; + + template + inline T compute_cf_pfq(detail::hypergeometric_pfq_cf_term& term, const Policy& pol) + { + BOOST_MATH_STD_USING + boost::uintmax_t max_iter = policies::get_max_series_iterations(); + const T result = tools::continued_fraction_b( + term, + boost::math::policies::get_epsilon(), + max_iter); + boost::math::policies::check_series_iterations( + "boost::math::hypergeometric_pfq_cf<%1%>(%1%,%1%,%1%)", + max_iter, + pol); + return result; + } + + template + inline T hypergeometric_0f1_cf(const T& b, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_cf_term f(b, z); + T result = detail::compute_cf_pfq(f, pol); + result = ((z / b) / result) + 1; + return result; + } + + template + inline T hypergeometric_1f0_cf(const T& a, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_cf_term f(a, z); + T result = detail::compute_cf_pfq(f, pol); + result = ((a * z) / result) + 1; + return result; + } + + template + inline T hypergeometric_1f1_cf(const T& a, const T& b, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_cf_term f(a, b, z); + T result = detail::compute_cf_pfq(f, pol); + result = (((a * z) / b) / result) + 1; + return result; + } + + template + inline T hypergeometric_1f2_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_cf_term f(a, b, c, z); + T result = detail::compute_cf_pfq(f, pol); + result = (((a * z) / (b * c)) / result) + 1; + return result; + } + + template + inline T hypergeometric_2f1_cf(const T& a, const T& b, const T& c, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_cf_term f(a, b, c, z); + T result = detail::compute_cf_pfq(f, pol); + result = ((((a * b) * z) / c) / result) + 1; + return result; + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_CF_HPP diff --git a/include/boost/math/special_functions/detail/hypergeometric_pade.hpp b/include/boost/math/special_functions/detail/hypergeometric_pade.hpp new file mode 100644 index 000000000..328a35c4b --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_pade.hpp @@ -0,0 +1,131 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_PADE_HPP + #define BOOST_MATH_HYPERGEOMETRIC_PADE_HPP + + namespace boost{ namespace math{ namespace detail{ + + // Luke: C ---------- SUBROUTINE R1F1P(CP, Z, A, B, N) ---------- + // Luke: C ----- PADE APPROXIMATION OF 1F1( 1 ; CP ; -Z ) ------- + template + inline T hypergeometric_1f1_pade(const T& cp, const T& zp, const Policy& pol) + { + BOOST_MATH_STD_USING + + static const T one = T(1); + + // Luke: C ------------- INITIALIZATION ------------- + const T z = -zp; + const T zz = z * z; + T b0 = one; + T a0 = one; + T xi1 = one; + T ct1 = cp + one; + T cp1 = cp - one; + + T b1 = one + (z / ct1); + T a1 = b1 - (z / cp); + + const unsigned max_iterations = boost::math::policies::get_max_series_iterations(); + + T b2 = T(0), a2 = T(0); + T result = T(0), prev_result = a1 / b1; + + for (unsigned k = 1; k < max_iterations; ++k) + { + // Luke: C ----- CALCULATION OF THE MULTIPLIERS ----- + // Luke: C ----------- FOR THE RECURSION ------------ + const T ct2 = ct1 * ct1; + const T g1 = one + ((cp1 / (ct2 + ct1 + ct1)) * z); + const T g2 = ((xi1 / (ct2 - one)) * ((xi1 + cp1) / ct2)) * zz; + + // Luke: C ------- THE RECURRENCE RELATIONS --------- + // Luke: C ------------ ARE AS FOLLOWS -------------- + b2 = (g1 * b1) + (g2 * b0); + a2 = (g1 * a1) + (g2 * a0); + + prev_result = result; + result = a2 / b2; + + // condition for interruption + if ((fabs(result) * boost::math::tools::epsilon()) > fabs(result - prev_result)) + break; + + b0 = b1; b1 = b2; + a0 = a1; a1 = a2; + + ct1 += 2; + ++xi1; + } + + return a2 / b2; + } + + // Luke: C -------- SUBROUTINE R2F1P(BP, CP, Z, A, B, N) -------- + // Luke: C ---- PADE APPROXIMATION OF 2F1( 1 , BP; CP ; -Z ) ---- + template + inline T hypergeometric_2f1_pade(const T& bp, const T& cp, const T& zp, const Policy& pol) + { + BOOST_MATH_STD_USING + + static const T one = T(1); + + // Luke: C ---------- INITIALIZATION ----------- + const T z = -zp; + const T zz = z * z; + T b0 = one; + T a0 = one; + T xi1 = one; + T ct1 = cp; + const T b1c1 = (cp - one) * (bp - one); + + T b1 = one + ((z / (cp + one)) * (bp + one)); + T a1 = b1 - ((bp / cp) * z); + + const unsigned max_iterations = boost::math::policies::get_max_series_iterations(); + + T b2 = T(0), a2 = T(0); + T result = T(0), prev_result = a1 / b1; + + for (unsigned k = 1; k < max_iterations; ++k) + { + // Luke: C ----- CALCULATION OF THE MULTIPLIERS ----- + // Luke: C ----------- FOR THE RECURSION ------------ + const T ct2 = ct1 + xi1; + const T ct3 = ct2 * ct2; + const T g2 = (((((ct1 / ct3) * (bp - ct1)) / (ct3 - one)) * xi1) * (bp + xi1)) * zz; + ++xi1; + const T g1 = one + (((((xi1 + xi1) * ct1) + b1c1) / (ct3 + ct2 + ct2)) * z); + + // Luke: C ------- THE RECURRENCE RELATIONS --------- + // Luke: C ------------ ARE AS FOLLOWS -------------- + b2 = (g1 * b1) + (g2 * b0); + a2 = (g1 * a1) + (g2 * a0); + + prev_result = result; + result = a2 / b2; + + // condition for interruption + if ((fabs(result) * boost::math::tools::epsilon()) > fabs(result - prev_result)) + break; + + b0 = b1; b1 = b2; + a0 = a1; a1 = a2; + + ++ct1; + } + + return a2 / b2; + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_PADE_HPP diff --git a/include/boost/math/special_functions/detail/hypergeometric_rational.hpp b/include/boost/math/special_functions/detail/hypergeometric_rational.hpp new file mode 100644 index 000000000..b9261a912 --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_rational.hpp @@ -0,0 +1,165 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP + #define BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP + + #include + + namespace boost{ namespace math{ namespace detail{ + + // Luke: C ------- SUBROUTINE R1F1P(AP, CP, Z, A, B, N) --------- + // Luke: C --- RATIONAL APPROXIMATION OF 1F1( AP ; CP ; -Z ) ---- + template + inline T hypergeometric_1f1_rational(const T& ap, const T& cp, const T& zp, const Policy& pol) + { + BOOST_MATH_STD_USING + + static const T zero = T(0), one = T(1), two = T(2), three = T(3); + + // Luke: C ------------- INITIALIZATION ------------- + const T z = -zp; + const T z2 = z / two; + + T ct1 = ap * (z / cp); + T ct2 = z2 / (one + cp); + T xn3 = zero; + T xn2 = one; + T xn1 = two; + T xn0 = three; + + T b1 = one; + T a1 = one; + T b2 = one + ((one + ap) * (z2 / cp)); + T a2 = b2 - ct1; + T b3 = one + ((two + b2) * (((two + ap) / three) * ct2)); + T a3 = b3 - ((one + ct2) * ct1); + ct1 = three; + + const unsigned max_iterations = boost::math::policies::get_max_series_iterations(); + + T a4 = T(0), b4 = T(0); + T result = T(0), prev_result = a3 / b3; + + for (unsigned k = 2; k < max_iterations; ++k) + { + // Luke: C ----- CALCULATION OF THE MULTIPLIERS ----- + // Luke: C ----------- FOR THE RECURSION ------------ + ct2 = (z2 / ct1) / (cp + xn1); + const T g1 = one + (ct2 * (xn2 - ap)); + ct2 *= ((ap + xn1) / (cp + xn2)); + const T g2 = ct2 * ((cp - xn1) + (((ap + xn0) / (ct1 + two)) * z2)); + const T g3 = ((ct2 * z2) * (((z2 / ct1) / (ct1 - two)) * ((ap + xn2)) / (cp + xn3))) * (ap - xn2); + + // Luke: C ------- THE RECURRENCE RELATIONS --------- + // Luke: C ------------ ARE AS FOLLOWS -------------- + b4 = (g1 * b3) + (g2 * b2) + (g3 * b1); + a4 = (g1 * a3) + (g2 * a2) + (g3 * a1); + + prev_result = result; + result = a4 / b4; + + // condition for interruption + if ((fabs(result) * boost::math::tools::epsilon()) > fabs(result - prev_result)) + break; + + b1 = b2; b2 = b3; b3 = b4; + a1 = a2; a2 = a3; a3 = a4; + + xn3 = xn2; xn2 = xn1; xn1 = xn0; ++xn0; + ct1 += two; + } + + return result; + } + + // Luke: C ----- SUBROUTINE R2F1P(AB, BP, CP, Z, A, B, N) ------- + // Luke: C -- RATIONAL APPROXIMATION OF 2F1( AB , BP; CP ; -Z ) - + template + inline T hypergeometric_2f1_rational(const T& ap, const T& bp, const T& cp, const T& zp, const unsigned n, const Policy& pol) + { + BOOST_MATH_STD_USING + + static const T one = T(1), two = T(2), three = T(3), four = T(4), + six = T(6), half_7 = T(3.5), half_3 = T(1.5), forth_3 = T(0.75); + + // Luke: C ------------- INITIALIZATION ------------- + const T z = -zp; + const T z2 = z / two; + + T sabz = (ap + bp) * z; + const T ab = ap * bp; + const T abz = ab * z; + const T abz1 = z + (abz + sabz); + const T abz2 = abz1 + (sabz + (three * z)); + const T cp1 = cp + one; + const T ct1 = cp1 + cp1; + + T b1 = one; + T a1 = one; + T b2 = one + (abz1 / (cp + cp)); + T a2 = b2 - (abz / cp); + T b3 = one + ((abz2 / ct1) * (one + (abz1 / ((-six) + (three * ct1))))); + T a3 = b3 - ((abz / cp) * (one + ((abz2 - abz1) / ct1))); + sabz /= four; + + const T abz1_div_4 = abz1 / four; + const T cp1_inc = cp1 + one; + const T cp1_mul_cp1_inc = cp1 * cp1_inc; + + boost::array d = {{ + ((half_7 - ab) * z2) - sabz, + abz1_div_4, + abz1_div_4 - (two * sabz), + cp1_inc, + cp1_mul_cp1_inc, + cp * cp1_mul_cp1_inc, + half_3, + forth_3, + forth_3 * z + }}; + + T xi = three; + T a4 = T(0), b4 = T(0); + for (unsigned k = 2; k < n; ++k) + { + // Luke: C ----- CALCULATION OF THE MULTIPLIERS ----- + // Luke: C ----------- FOR THE RECURSION ------------ + T g3 = (d[2] / d[7]) * (d[1] / d[5]); + d[1] += d[8] + sabz; + d[2] += d[8] - sabz; + g3 *= d[1] / d[6]; + T g1 = one + (((d[1] + d[0]) / d[6]) / d[3]); + T g2 = (d[1] / d[4]) / d[6]; + d[7] += two * d[6]; + ++d[6]; + g2 *= cp1 - (xi + ((d[2] + d[0]) / d[6])); + + // Luke: C ------- THE RECURRENCE RELATIONS --------- + // Luke: C ------------ ARE AS FOLLOWS -------------- + b4 = (g1 * b3) + (g2 * b2) + (g3 * b1); + a4 = (g1 * a3) + (g2 * a2) + (g3 * a1); + b1 = b2; b2 = b3; b3 = b4; + a1 = a2; a2 = a3; a3 = a4; + + d[8] += z2; + d[0] += two * d[8]; + d[5] += three * d[4]; + d[4] += two * d[3]; + ++d[3]; + ++xi; + } + + return a4 / b4; + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_RATIONAL_HPP diff --git a/include/boost/math/special_functions/detail/hypergeometric_separated_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_separated_series.hpp new file mode 100644 index 000000000..dd749e572 --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_separated_series.hpp @@ -0,0 +1,51 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_SEPARATED_SERIES_HPP + #define BOOST_MATH_HYPERGEOMETRIC_SEPARATED_SERIES_HPP + + namespace boost { namespace math { namespace detail { + + template + inline T hypergeometric_1f1_separated_series(const T& a, const T& b, const T& z, const Policy& pol) + { + BOOST_MATH_STD_USING + + boost::uintmax_t max_iter = policies::get_max_series_iterations(); + const T factor = policies::get_epsilon(); + + T denom = 1, numer = 1; + T intermediate_result = 1, result = 1; + T a_pochhammer = a, z_pow = z; + unsigned N = 0; + while (--max_iter) + { + ++N; + const T mult = (((b + N) - 1) * N); + denom *= mult; numer *= mult; + numer += a_pochhammer * z_pow; + + result = numer / denom; + + if (fabs(factor * result) > fabs(result - intermediate_result)) + break; + + intermediate_result = result; + + a_pochhammer *= (a + N); + z_pow *= z; + } + + return result; + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_SEPARATED_SERIES_HPP diff --git a/include/boost/math/special_functions/detail/hypergeometric_series.hpp b/include/boost/math/special_functions/detail/hypergeometric_series.hpp new file mode 100644 index 000000000..36767e11d --- /dev/null +++ b/include/boost/math/special_functions/detail/hypergeometric_series.hpp @@ -0,0 +1,242 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) +// +#ifndef BOOST_MATH_HYPERGEOMETRIC_SERIES_HPP + #define BOOST_MATH_HYPERGEOMETRIC_SERIES_HPP + +#include +#include + + namespace boost { namespace math { namespace detail { + + // primary template for term of Taylor series + template + struct hypergeometric_pfq_generic_series_term; + + // partial specialization for 0F1 + template + struct hypergeometric_pfq_generic_series_term + { + typedef T result_type; + + hypergeometric_pfq_generic_series_term(const T& b, const T& z) + : n(0), term(1), b(b), z(z) + { + } + + T operator()() + { + BOOST_MATH_STD_USING + const T r = term; + term *= ((1 / ((b + n) * (n + 1))) * z); + ++n; + return r; + } + + private: + unsigned n; + T term; + const T b, z; + }; + + // partial specialization for 1F0 + template + struct hypergeometric_pfq_generic_series_term + { + typedef T result_type; + + hypergeometric_pfq_generic_series_term(const T& a, const T& z) + : n(0), term(1), a(a), z(z) + { + } + + T operator()() + { + BOOST_MATH_STD_USING + const T r = term; + term *= (((a + n) / (n + 1)) * z); + ++n; + return r; + } + + private: + unsigned n; + T term; + const T a, z; + }; + + // partial specialization for 1F1 + template + struct hypergeometric_pfq_generic_series_term + { + typedef T result_type; + + hypergeometric_pfq_generic_series_term(const T& a, const T& b, const T& z) + : n(0), term(1), a(a), b(b), z(z) + { + } + + T operator()() + { + BOOST_MATH_STD_USING + const T r = term; + term *= (((a + n) / ((b + n) * (n + 1))) * z); + ++n; + return r; + } + + private: + unsigned n; + T term; + const T a, b, z; + }; + + // partial specialization for 1F2 + template + struct hypergeometric_pfq_generic_series_term + { + typedef T result_type; + + hypergeometric_pfq_generic_series_term(const T& a, const T& b1, const T& b2, const T& z) + : n(0), term(1), a(a), b1(b1), b2(b2), z(z) + { + } + + T operator()() + { + BOOST_MATH_STD_USING + const T r = term; + term *= (((a + n) / ((b1 + n) * (b2 + n) * (n + 1))) * z); + ++n; + return r; + } + + private: + unsigned n; + T term; + const T a, b1, b2, z; + }; + + // partial specialization for 2F0 + template + struct hypergeometric_pfq_generic_series_term + { + typedef T result_type; + + hypergeometric_pfq_generic_series_term(const T& a1, const T& a2, const T& z) + : n(0), term(1), a1(a1), a2(a2), z(z) + { + } + + T operator()() + { + BOOST_MATH_STD_USING + const T r = term; + term *= (((a1 + n) * (a2 + n) / (n + 1)) * z); + ++n; + return r; + } + + private: + unsigned n; + T term; + const T a1, a2, z; + }; + + // partial specialization for 2F1 + template + struct hypergeometric_pfq_generic_series_term + { + typedef T result_type; + + hypergeometric_pfq_generic_series_term(const T& a1, const T& a2, const T& b, const T& z) + : n(0), term(1), a1(a1), a2(a2), b(b), z(z) + { + } + + T operator()() + { + BOOST_MATH_STD_USING + const T r = term; + term *= (((a1 + n) * (a2 + n) / ((b + n) * (n + 1))) * z); + ++n; + return r; + } + + private: + unsigned n; + T term; + const T a1, a2, b, z; + }; + + // we don't need to define extra check and make a polinom from + // series, when p(i) and q(i) are negative integers and p(i) >= q(i) + // as described in functions.wolfram.alpha, because we always + // stop summation when result (in this case numerator) is zero. + template + inline T sum_pfq_series(detail::hypergeometric_pfq_generic_series_term& term, const Policy& pol) + { + BOOST_MATH_STD_USING + boost::uintmax_t max_iter = policies::get_max_series_iterations(); +#if BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582)) + const T zero = 0; + const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon(), max_iter, zero); +#else + const T result = boost::math::tools::sum_series(term, boost::math::policies::get_epsilon(), max_iter); +#endif + policies::check_series_iterations("boost::math::hypergeometric_pfq_generic_series<%1%>(%1%,%1%,%1%)", max_iter, pol); + return result; + } + + template + inline T hypergeometric_0f1_generic_series(const T& b, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_generic_series_term s(b, z); + return detail::sum_pfq_series(s, pol); + } + + template + inline T hypergeometric_1f0_generic_series(const T& a, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_generic_series_term s(a, z); + return detail::sum_pfq_series(s, pol); + } + + template + inline T hypergeometric_1f1_generic_series(const T& a, const T& b, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_generic_series_term s(a, b, z); + return detail::sum_pfq_series(s, pol); + } + + template + inline T hypergeometric_1f2_generic_series(const T& a, const T& b1, const T& b2, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_generic_series_term s(a, b1, b2, z); + return detail::sum_pfq_series(s, pol); + } + + template + inline T hypergeometric_2f0_generic_series(const T& a1, const T& a2, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_generic_series_term s(a1, a2, z); + return detail::sum_pfq_series(s, pol); + } + + template + inline T hypergeometric_2f1_generic_series(const T& a1, const T& a2, const T& b, const T& z, const Policy& pol) + { + detail::hypergeometric_pfq_generic_series_term s(a1, a2, b, z); + return detail::sum_pfq_series(s, pol); + } + + } } } // namespaces + +#endif // BOOST_MATH_HYPERGEOMETRIC_SERIES_HPP diff --git a/include/boost/math/special_functions/hypergeometric_0F1.hpp b/include/boost/math/special_functions/hypergeometric_0F1.hpp new file mode 100644 index 000000000..fd9c64fba --- /dev/null +++ b/include/boost/math/special_functions/hypergeometric_0F1.hpp @@ -0,0 +1,110 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef _BOOST_HYPERGEOMETRIC_0F1_HPP_ + #define _BOOST_HYPERGEOMETRIC_0F1_HPP_ + +#include +#include +#include +#include + +namespace boost { namespace math { namespace detail { + + + template + struct hypergeometric_0F1_cf + { + T b, z; + int k; + hypergeometric_0F1_cf(T b_, T z_) : b(b_), z(z_), k(0) {} + typedef std::pair result_type; + + result_type operator()() + { + ++k; + return std::make_pair(-z / ((k + 1) * (b + k)), 1 + z / ((k + 1) * (b + k))); + } + }; + + template + T hypergeometric_0F1_cf_imp(T b, T z, const Policy& pol, const char* function) + { + hypergeometric_0F1_cf evaluator(b, z); + boost::uintmax_t max_iter = policies::get_max_series_iterations(); + T cf = tools::continued_fraction_a(evaluator, policies::get_epsilon(), max_iter); + policies::check_series_iterations(function, max_iter, pol); + return 1 + z / (b * (1 + cf)); + } + + + template + inline T hypergeometric_0f1_imp(const T& b, const T& z, const Policy& pol) + { + const char* function = "boost::math::hypergeometric_0f1<%1%,%1%>(%1%, %1%)"; + BOOST_MATH_STD_USING + + // some special cases + if (z == 0) + return T(1); + + if ((b <= 0) && (b == floor(b))) + return policies::raise_pole_error( + function, + "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol); + + if (z < -5 && b > -5) + { + // Series is alternating and divergent, need to do something else here, + // Bessel function relation is much more accurate, unless |b| is similarly + // large to |z|, otherwise the CF formula suffers from cancellation when + // the result would be very small. + if (fabs(z / b) > 4) + return hypergeometric_0f1_bessel(b, z, pol); + return hypergeometric_0F1_cf_imp(b, z, pol, function); + } + // evaluation through Taylor series looks + // more precisious than Bessel relation: + // detail::hypergeometric_0f1_bessel(b, z, pol); + return detail::hypergeometric_0f1_generic_series(b, z, pol); + } + +} // namespace detail + +template +inline typename tools::promote_args::type hypergeometric_0F1(T1 b, T2 z, const Policy& /* pol */) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast( + detail::hypergeometric_0f1_imp( + static_cast(b), + static_cast(z), + forwarding_policy()), + "boost::math::hypergeometric_0f1<%1%>(%1%,%1%)"); +} + +template +inline typename tools::promote_args::type hypergeometric_0F1(T1 b, T2 z) +{ + return hypergeometric_0F1(b, z, policies::policy<>()); +} + + +} } // namespace boost::math + +#endif // _BOOST_HYPERGEOMETRIC_2014_04_07_HPP_ diff --git a/include/boost/math/special_functions/hypergeometric_1F0.hpp b/include/boost/math/special_functions/hypergeometric_1F0.hpp new file mode 100644 index 000000000..5a6f4c54d --- /dev/null +++ b/include/boost/math/special_functions/hypergeometric_1F0.hpp @@ -0,0 +1,73 @@ + +/////////////////////////////////////////////////////////////////////////////// +// Copyright 2014 Anton Bikineev +// Copyright 2014 Christopher Kormanyos +// Copyright 2014 John Maddock +// Copyright 2014 Paul Bristow +// Distributed under the Boost +// Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#ifndef _BOOST_HYPERGEOMETRIC_1F0_HPP_ + #define _BOOST_HYPERGEOMETRIC_1F0_HPP_ + +#include +#include + + +namespace boost { namespace math { namespace detail { + +template +inline T hypergeometric_1f0_imp(const T& a, const T& z, const Policy& pol) +{ + static const char* function = "boost::math::hypergeometric_1f0<%1%,%1%>(%1%, %1%)"; + BOOST_MATH_STD_USING // pow + + if (z == 1) + return policies::raise_pole_error( + function, + "Evaluation of 1F0 with z = %1%.", + z, + pol); + if (1 - z < 0) + { + if (floor(a) != a) + return policies::raise_domain_error(function, + "Result is complex when a is non-integral and z > 1, but got z = %1%", z, pol); + } + // more naive and convergent method than series + return pow(T(1 - z), T(-a)); +} + +} // namespace detail + +template +inline typename tools::promote_args::type hypergeometric_1f0(T1 a, T2 z, const Policy&) +{ + BOOST_FPU_EXCEPTION_GUARD + typedef typename tools::promote_args::type result_type; + typedef typename policies::evaluation::type value_type; + typedef typename policies::normalise< + Policy, + policies::promote_float, + policies::promote_double, + policies::discrete_quantile<>, + policies::assert_undefined<> >::type forwarding_policy; + return policies::checked_narrowing_cast( + detail::hypergeometric_1f0_imp( + static_cast(a), + static_cast(z), + forwarding_policy()), + "boost::math::hypergeometric_1f0<%1%>(%1%,%1%)"); +} + +template +inline typename tools::promote_args::type hypergeometric_1f0(T1 a, T2 z) +{ + return hypergeometric_1f0(a, z, policies::policy<>()); +} + + + } } // namespace boost::math + +#endif // _BOOST_HYPERGEOMETRIC_2014_04_07_HPP_ diff --git a/test/test_0F1.cpp b/test/test_0F1.cpp new file mode 100644 index 000000000..f52a1a0a2 --- /dev/null +++ b/test/test_0F1.cpp @@ -0,0 +1,85 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include "test_0F1.hpp" + +#include +#include + +void expected_results() +{ + // + // Define the max and mean errors expected for + // various compilers and platforms. + // + const char* largest_type; +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + if(boost::math::policies::digits >() == boost::math::policies::digits >()) + { + largest_type = "(long\\s+)?double|real_concept|cpp_bin_float_quad|dec_40"; + } + else + { + largest_type = "long double|real_concept"; + } +#else + largest_type = "(long\\s+)?double"; +#endif + + add_expected_result( + ".*", // compiler + ".*", // stdlib + ".*", // platform + largest_type, // test type(s) + "Integer.*", // test data group + ".*", 300, 100); // test function + add_expected_result( + ".*", // compiler + ".*", // stdlib + ".*", // platform + largest_type, // test type(s) + "Real.*", // test data group + ".*", 2000, 1000); // test function + + add_expected_result( + ".*", // compiler + ".*", // stdlib + ".*", // platform + largest_type, // test type(s) + "Large.*", // test data group + ".*", 400, 100); // test function + + + // + // Finish off by printing out the compiler/stdlib/platform names, + // we do this to make it easier to mark up expected error rates. + // + std::cout << "Tests run with " << BOOST_COMPILER << ", " + << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl; +} + +BOOST_AUTO_TEST_CASE( test_main ) +{ + typedef boost::multiprecision::number > dec_40; + expected_results(); + BOOST_MATH_CONTROL_FP; + +#ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS + test_spots(0.0F, "float"); +#endif + test_spots(0.0, "double"); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_spots(0.0L, "long double"); +#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS + test_spots(boost::math::concepts::real_concept(0.1), "real_concept"); +#endif +#endif + test_spots(boost::multiprecision::cpp_bin_float_quad(), "cpp_bin_float_quad"); + test_spots(dec_40(), "dec_40"); +} + + + diff --git a/test/test_0F1.hpp b/test/test_0F1.hpp new file mode 100644 index 000000000..bdcc97f61 --- /dev/null +++ b/test/test_0F1.hpp @@ -0,0 +1,197 @@ +// Copyright John Maddock 2006. +// Copyright Paul A. Bristow 2007, 2009 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#define BOOST_TEST_MAIN +#include +#include +#include +#include +#include +#include +#include +#include +#include "functor.hpp" + +#include "handle_test_result.hpp" +#include "table_type.hpp" + +#include + +template +void do_test_0F1(const T& data, const char* type_name, const char* test_name) +{ + typedef Real value_type; + + typedef value_type(*pg)(value_type, value_type); +#if defined(BOOST_MATH_NO_DEDUCED_FUNCTION_POINTERS) + pg funcp = boost::math::hypergeometric_0F1; +#else + pg funcp = boost::math::hypergeometric_0F1; +#endif + + boost::math::tools::test_result result; + + std::cout << "Testing " << test_name << " with type " << type_name + << "\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n"; + + // + // test hypergeometric_0F1 against data: + // + result = boost::math::tools::test_hetero( + data, + bind_func(funcp, 0, 1), + extract_result(2)); + handle_test_result(result, data[result.worst()], result.worst(), type_name, "hypergeometric_0F1", test_name); + std::cout << std::endl; +} + +template +void test_spots(T, const char* type_name) +{ +#ifndef SC_ +#define SC_(x) BOOST_MATH_BIG_CONSTANT(T, 1000000, x) +#endif + + // + // basic sanity checks, tolerance is 10 epsilon expressed as a percentage: + // + T tolerance = boost::math::tools::epsilon() * 1000; + + BOOST_CHECK_EQUAL(boost::math::hypergeometric_0F1(T(3), T(0)), 1); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_0F1(T(-3), T(0)), 1); + BOOST_CHECK_EQUAL(boost::math::hypergeometric_0F1(T(0), T(0)), 1); + + BOOST_CHECK_THROW(boost::math::hypergeometric_0F1(T(0), T(-1)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_0F1(T(-1), T(-1)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_0F1(T(-10), T(-5)), std::domain_error); + + static const boost::array, 50> hypergeometric_0F1_integer_data = { { + { SC_(4.0), SC_(-20.0), SC_(-0.012889714201783047561923257996127233830940165138385) }, + { SC_(8.0), SC_(-20.0), SC_(0.046498609282365144223175012935939437508273248399881) }, + { SC_(12.0), SC_(-20.0), SC_(0.16608847431869756642136191351311569335145459224622) }, + { SC_(16.0), SC_(-20.0), SC_(0.27230484709157170329168048388841880599105216477631) }, + { SC_(20.0), SC_(-20.0), SC_(0.35865872656868844615709101792040025253126126604266) }, + { SC_(4.0), SC_(-16.0), SC_(-0.027293644412433023379286103818840667403690937153604) }, + { SC_(8.0), SC_(-16.0), SC_(0.098618710511372349330666801041676087431136532039702) }, + { SC_(12.0), SC_(-16.0), SC_(0.24360114226383905073379763460037817885919457531523) }, + { SC_(16.0), SC_(-16.0), SC_(0.35635186318802906043824855864337727878754460163525) }, + { SC_(20.0), SC_(-16.0), SC_(0.44218381382689101428948260613085371477815110358789) }, + { SC_(4.0), SC_(-12.0), SC_(-0.021743572290699436419371120781513860006290363262907) }, + { SC_(8.0), SC_(-12.0), SC_(0.19025625754362006866949730683824627505504067855043) }, + { SC_(12.0), SC_(-12.0), SC_(0.35251228238278927379621049815222218665165551016489) }, + { SC_(16.0), SC_(-12.0), SC_(0.46415411486674623230458980010115972932474705884865) }, + { SC_(20.0), SC_(-12.0), SC_(0.54394918325286018927327004362535051310016558628741) }, + { SC_(4.0), SC_(-8.0), SC_(0.056818744289274872033266550620647787396712125304880) }, + { SC_(8.0), SC_(-8.0), SC_(0.34487371876996263249797701802458885718691612997456) }, + { SC_(12.0), SC_(-8.0), SC_(0.50411654015891701804499796523449656998841355305043) }, + { SC_(16.0), SC_(-8.0), SC_(0.60191459981670594041254437708158847428118361245442) }, + { SC_(20.0), SC_(-8.0), SC_(0.66770752550930138035694866478078941681114294465418) }, + { SC_(4.0), SC_(-4.0), SC_(0.32262860540671645526863760914000166725449779629143) }, + { SC_(8.0), SC_(-4.0), SC_(0.59755773349355150397404772151441126513126998265958) }, + { SC_(12.0), SC_(-4.0), SC_(0.71337465206009117934071859694314971137807212605147) }, + { SC_(16.0), SC_(-4.0), SC_(0.77734333649378860739496954157535257278092349684783) }, + { SC_(20.0), SC_(-4.0), SC_(0.81794177985447769150469288350369205683856312760890) }, + + { SC_(4.0), SC_(4.0), SC_(2.5029568338152582758923890008139391395035041790831) }, + { SC_(8.0), SC_(4.0), SC_(1.6273673128576761227855719910743734060605725722129) }, + { SC_(12.0), SC_(4.0), SC_(1.3898419290864057799739567227851793491657442624207) }, + { SC_(16.0), SC_(4.0), SC_(1.2817098157957427946677711269410726972209834860612) }, + { SC_(20.0), SC_(4.0), SC_(1.2202539302152377230940386181201477276788392792437) }, + { SC_(4.0), SC_(8.0), SC_(5.5616961007411965409200003309686924059253894118586) }, + { SC_(8.0), SC_(8.0), SC_(2.5877053985451664722152913482683136948296873738479) }, + { SC_(12.0), SC_(8.0), SC_(1.9166410733572697158003086323981583993970490592046) }, + { SC_(16.0), SC_(8.0), SC_(1.6370675016890669952237854163997946987362497613701) }, + { SC_(20.0), SC_(8.0), SC_(1.4862852701827990444915220582410007454379891584086) }, + { SC_(4.0), SC_(12.0), SC_(11.419268276211177842169936131590385979116019595164) }, + { SC_(8.0), SC_(12.0), SC_(4.0347215359576567066789638314925802225312840819037) }, + { SC_(12.0), SC_(12.0), SC_(2.6242497527837800417573064942486918368886996538285) }, + { SC_(16.0), SC_(12.0), SC_(2.0840468784170876805932772732753387258909164486511) }, + { SC_(20.0), SC_(12.0), SC_(1.8071042457762091748544382847762106786633952487005) }, + { SC_(4.0), SC_(16.0), SC_(22.132051970576036053853444648907108439504682530918) }, + { SC_(8.0), SC_(16.0), SC_(6.1850485247748975008808779795786699492711191898792) }, + { SC_(12.0), SC_(16.0), SC_(3.5694322843488018916484224923627864928705138154372) }, + { SC_(16.0), SC_(16.0), SC_(2.6447371137201451261118187672029372265909501355722) }, + { SC_(20.0), SC_(16.0), SC_(2.1934058398888071720297525592515838555602675797235) }, + { SC_(4.0), SC_(20.0), SC_(41.021743268279206331672552645354782698296383424328) }, + { SC_(8.0), SC_(20.0), SC_(9.3414225299809886395081381945971250426599939097753) }, + { SC_(12.0), SC_(20.0), SC_(4.8253866205826406499959001774187695527272168375992) }, + { SC_(16.0), SC_(20.0), SC_(3.3462305133519485784864062004430532216764447939942) }, + { SC_(20.0), SC_(20.0), SC_(2.6578698872220394617444624241257799193518140676691) }, + } }; + + static const boost::array, 121> hypergeometric_0F1_real_data = { { + { SC_(-20.25), SC_(-20.25), SC_(2.7957722939837585922721916572828265337379885810463) }, + {SC_(-20.25), SC_(-16.25), SC_(2.2711152061687223613453898226719770018304527545367) }, + {SC_(-20.25), SC_(-12.25), SC_(1.8494593608212709260834907188011226625504792564008) }, + {SC_(-20.25), SC_(-8.25), SC_(1.5096123122209879279439571483106546579344083224475) }, + {SC_(-20.25), SC_(-4.25), SC_(1.2349600261245108147217156468676391406107607937783) }, + {SC_(-20.25), SC_(-0.25), SC_(1.0124262131477737489760461958823145562583738548481) }, + {SC_(-20.25), SC_(3.75), SC_(0.83168105368847042690250207452503923498261925399290) }, + {SC_(-20.25), SC_(7.75), SC_(0.68453616446705348287989858944821561345027038288523) }, + {SC_(-20.25), SC_(11.75), SC_(0.56447959084730590158981806575147730464939958811641) }, + {SC_(-20.25), SC_(15.75), SC_(0.46631672324832840040099031042536652056451149930724) }, + {SC_(-20.25), SC_(19.75), SC_(0.38589178960453729031753371638509489720353970202557) }, + {SC_(-16.25), SC_(-20.25), SC_(3.6853033678389815793858749638872118256381350051933) }, + {SC_(-16.25), SC_(-16.25), SC_(2.8189932976213221794412912826715210157252407675492) }, + {SC_(-16.25), SC_(-12.25), SC_(2.1683206814848382415951069178541485022190389032861) }, + {SC_(-16.25), SC_(-8.25), SC_(1.6762781419304165453124288212144469387497248544980) }, + {SC_(-16.25), SC_(-4.25), SC_(1.3019174596541545423250019596925072929338360502457) }, + {SC_(-16.25), SC_(-0.25), SC_(1.0155114597289924802082858072926838795825629097345) }, + {SC_(-16.25), SC_(3.75), SC_(0.79528104007050457443200675888120572482230255556703) }, + {SC_(-16.25), SC_(7.75), SC_(0.62514113457518737301424222399549180276283312092061) }, + {SC_(-16.25), SC_(11.75), SC_(0.49312611939645449055396232495817215814052621005108) }, + {SC_(-16.25), SC_(15.75), SC_(0.39027640231656357551572959324213784623094177949066) }, + {SC_(-16.25), SC_(19.75), SC_(0.30983458381312295039033218076191326310645958859520) }, + {SC_(-12.25), SC_(-20.25), SC_(6.1704666666278938283835490432795642056267750516111) }, + {SC_(-12.25), SC_(-16.25), SC_(4.1552272888330207621302106989548328150305062793812) }, + {SC_(-12.25), SC_(-12.25), SC_(2.8629405204600569631773103187639297768454467434670) }, + {SC_(-12.25), SC_(-8.25), SC_(2.0050800521989931571210090020740770238121125995302) }, + {SC_(-12.25), SC_(-4.25), SC_(1.4226882175494078094616273681832235056660520864074) }, + {SC_(-12.25), SC_(-0.25), SC_(1.0206367767231596623643982724086720984175407326239) }, + {SC_(-12.25), SC_(3.75), SC_(0.73925379130945509368577275378254380866098045869975) }, + {SC_(-12.25), SC_(7.75), SC_(0.54000881441780813642748748511053584681016992673867) }, + {SC_(-12.25), SC_(11.75), SC_(0.39734252271724847449760568202775839997979417321096) }, + {SC_(-12.25), SC_(15.75), SC_(0.28584659061590258036270241201657559564217479119707) }, + {SC_(-12.25), SC_(19.75), SC_(-0.0081173945118881955152340569022581720505805271716412) }, + {SC_(-8.25), SC_(-20.25), SC_(30.214776746264667028399437825207186990952563144565) }, + {SC_(-8.25), SC_(-16.25), SC_(13.206923736113233408873312298765667487728366521206) }, + {SC_(-8.25), SC_(-12.25), SC_(5.8765290507182112008139505719133766922110698377617) }, + {SC_(-8.25), SC_(-8.25), SC_(2.9943833929883303847073818502561706078874721587174) }, + {SC_(-8.25), SC_(-4.25), SC_(1.7091345745360580812285240901667243909385176010261) }, + {SC_(-8.25), SC_(-0.25), SC_(1.0308325464760906146441689208198464885112314788329) }, + {SC_(-8.25), SC_(3.75), SC_(0.64304308306518963175886701666419407366549702609687) }, + {SC_(-8.25), SC_(7.75), SC_(0.37629752350234203158515342524983705294147434733076) }, + {SC_(-8.25), SC_(11.75), SC_(-2.1055871696987101269944140698868945037062228729050) }, + {SC_(-8.25), SC_(15.75), SC_(-50.562697627448074052585932806347970082382096647974) }, + {SC_(-8.25), SC_(19.75), SC_(-579.94915881084897988862450916571551198001894291852) }, + { SC_(-4.25), SC_(-20.25), SC_(-66.582207191755226871227185886830083823068161088090) }, {SC_(-4.25), SC_(-16.25), SC_(0.99585574579280575697804972744074637223473969232632)}, {SC_(-4.25), SC_(-12.25), SC_(23.510774658631484713205463401664873191546111331469)}, {SC_(-4.25), SC_(-8.25), SC_(14.506347207276003438062329264569459812455033677057)}, {SC_(-4.25), SC_(-4.25), SC_(3.8047736693317086794526961473100619698786965398531)}, {SC_(-4.25), SC_(-0.25), SC_(1.0611747490912035255122776417926285437513018011867)}, {SC_(-4.25), SC_(3.75), SC_(-0.80390254520675152395475412235031817082929562419835)}, {SC_(-4.25), SC_(7.75), SC_(-100.43339815718209139439327680390551018715836820761)}, {SC_(-4.25), SC_(11.75), SC_(-1533.8025906900325906352282963205257001995564543560)}, {SC_(-4.25), SC_(15.75), SC_(-11903.890928701912286627232258091180714843387637568)}, {SC_(-4.25), SC_(19.75), SC_(-63573.720858843295034857393096900086420282088096223)}, + {SC_(-0.25), SC_(-20.25), SC_(5.8154983455887482996117997327072965277711303658190)}, {SC_(-0.25), SC_(-16.25), SC_(7.8645427110771775078914759561882836266181443022557)}, {SC_(-0.25), SC_(-12.25), SC_(2.8859855860040835605817890532890777322459344278127)}, {SC_(-0.25), SC_(-8.25), SC_(-4.4978057721216329075474975538674330456610136767402)}, {SC_(-0.25), SC_(-4.25), SC_(-3.2843951886889091457664917868727091433634846305430)}, {SC_(-0.25), SC_(-0.25), SC_(1.8410918501257885334774939532672407660260932497146)}, {SC_(-0.25), SC_(3.75), SC_(-89.507839434346979238664123270851349643390027303630)}, {SC_(-0.25), SC_(7.75), SC_(-684.94531486572802102565352112548840783460445744423)}, {SC_(-0.25), SC_(11.75), SC_(-2980.9405282185801785979381426998934753615183894674)}, {SC_(-0.25), SC_(15.75), SC_(-9963.0984480949236929126495173451788015488533057312)}, {SC_(-0.25), SC_(19.75), SC_(-28354.285610765325603366060327671144508742992535586)}, + {SC_(3.75), SC_(-20.25), SC_(-0.0076845585942792341317694895749400053976192347461355)}, {SC_(3.75), SC_(-16.25), SC_(-0.026676473833363805218474733943073126852098252692238)}, {SC_(3.75), SC_(-12.25), SC_(-0.033150612620140155332437908902071743835476014605824)}, {SC_(3.75), SC_(-8.25), SC_(0.024862128358835876745073743789586165594233386818764)}, {SC_(3.75), SC_(-4.25), SC_(0.26643741037095553396026886662906805496932883791011)}, {SC_(3.75), SC_(-0.25), SC_(0.93506252732788015337201598317744583689077922422313)}, {SC_(3.75), SC_(3.75), SC_(2.4937080902782278272865074460355462658482459750891)}, {SC_(3.75), SC_(7.75), SC_(5.7782165977734310384947791703219081203794593104898)}, {SC_(3.75), SC_(11.75), SC_(12.239149970819448162027824359733914867194838641207)}, {SC_(3.75), SC_(15.75), SC_(24.315766317657856526802495882902104213004243807834)}, {SC_(3.75), SC_(19.75), SC_(46.004706133409223855903017790328079349160741103989)}, + {SC_(7.75), SC_(-20.25), SC_(0.037235565108204273745081008045192031428127589963378)}, {SC_(7.75), SC_(-16.25), SC_(0.084729597975458928169783385558211341184710436894527)}, {SC_(7.75), SC_(-12.25), SC_(0.17080635394402470513520780978929502652414093665998)}, {SC_(7.75), SC_(-8.25), SC_(0.31950343789501896907243089625679164966361926427901)}, {SC_(7.75), SC_(-4.25), SC_(0.56722331717592488544254097481360697514802761113579)}, {SC_(7.75), SC_(-0.25), SC_(0.96819884906571311048909491350402784950555589172778)}, {SC_(7.75), SC_(3.75), SC_(1.6020878195055145404559679357988911544832570195929)}, {SC_(7.75), SC_(7.75), SC_(2.5844402449306647890994627013439918401353618660618)}, {SC_(7.75), SC_(11.75), SC_(4.0810182332680718846036289812422612468750062913766)}, {SC_(7.75), SC_(15.75), SC_(6.3272529834696435764802870864047851742114357592781)}, {SC_(7.75), SC_(19.75), SC_(9.6545155398613375064364794042284225393502577025757)}, + {SC_(11.75), SC_(-20.25), SC_(0.15486394311085924113988037667235580632108514049961)}, {SC_(11.75), SC_(-16.25), SC_(0.22985383829087750746997781544418442650948306731355)}, {SC_(11.75), SC_(-12.25), SC_(0.33615798547981440387258985310070073263439571574352)}, {SC_(11.75), SC_(-8.25), SC_(0.48533791924865732402484841770387435054807176843922)}, {SC_(11.75), SC_(-4.25), SC_(0.69279476096089396016927378719976382816612949045723)}, {SC_(11.75), SC_(-0.25), SC_(0.97893073946457038719047235276954159274687951133239)}, {SC_(11.75), SC_(3.75), SC_(1.3706338991326610410785896448412400026741719976537)}, {SC_(11.75), SC_(7.75), SC_(1.9031693475249062854591923959807180024296317864633)}, {SC_(11.75), SC_(11.75), SC_(2.6225804467534900775963810848242641845248345641566)}, {SC_(11.75), SC_(15.75), SC_(3.5887279276489256865660553406508246843289687389280)}, {SC_(11.75), SC_(19.75), SC_(4.8791249650630001302695046693426702384634788655376)}, + {SC_(15.75), SC_(-20.25), SC_(0.26169993773941504714021934731773052404003435302548)}, {SC_(15.75), SC_(-16.25), SC_(0.34428767856756521281521994646213632566799672170823)}, {SC_(15.75), SC_(-12.25), SC_(0.45068705303931291468643647219845752956581924151602)}, {SC_(15.75), SC_(-8.25), SC_(0.58722904807165396939912732843015800226017006383889)}, {SC_(15.75), SC_(-4.25), SC_(0.76180845119510741657129800934863890343557924548267)}, {SC_(15.75), SC_(-0.25), SC_(0.98424488519002985986673227194811800795924545598747)}, {SC_(15.75), SC_(3.75), SC_(1.2667221438640806060367418270328561473306297946124)}, {SC_(15.75), SC_(7.75), SC_(1.6243218504073229104605986746109545497016155557456)}, {SC_(15.75), SC_(11.75), SC_(2.0756705639798239579947551789739745396066053149652)}, {SC_(15.75), SC_(15.75), SC_(2.6437231357748258767712950042187028922110335308912)}, {SC_(15.75), SC_(19.75), SC_(3.3567094629489291868172072628630988653584236011026)}, + {SC_(19.75), SC_(-20.25), SC_(0.34910186209886263988734004369646266800659689838760)}, {SC_(19.75), SC_(-16.25), SC_(0.43171851514846723176830271990344447636246814031258)}, {SC_(19.75), SC_(-12.25), SC_(0.53264760691226214051563276361223610040193575209722)}, {SC_(19.75), SC_(-8.25), SC_(0.65570932018638949205413006445398352201545799769694)}, {SC_(19.75), SC_(-4.25), SC_(0.80547695536244789152425915957330496255335464733020)}, {SC_(19.75), SC_(-0.25), SC_(0.98741773517508964436868231216434995093107999788450)}, {SC_(19.75), SC_(3.75), SC_(1.2080586524527299846909832323308446155063281314444)}, {SC_(19.75), SC_(7.75), SC_(1.4751816095060559707610879513647971547078203326785)}, {SC_(19.75), SC_(11.75), SC_(1.7980527870029919079278807774802021351478928537606)}, {SC_(19.75), SC_(15.75), SC_(2.1876919765295108137252526717389200156554872202073)}, {SC_(19.75), SC_(19.75), SC_(2.6571885305131654484163938203293169727636168525176)} + } }; + + static const boost::array, 49> hypergeometric_0F1_large_data = { { + {SC_(-3000.25), SC_(-3000.25), SC_(2.7187352281773126038196142076996225709478475910049)}, {SC_(-3000.25), SC_(-2000.25), SC_(1.9479325178881941977364868969702984561194629278005)}, {SC_(-3000.25), SC_(-1000.25), SC_(1.3957158200680510708579318399559946925035296441969)}, {SC_(-3000.25), SC_(-0.25), SC_(1.0000833298623651641409421309527751141654844918029)}, {SC_(-3000.25), SC_(999.75), SC_(0.71662418615385142079900540026845986168175943851435)}, {SC_(-3000.25), SC_(1999.75), SC_(0.51352644595687380807192108174566627121840655509659)}, {SC_(-3000.25), SC_(2999.75), SC_(0.36800205062066223494767361354167341477112311421338)}, + {SC_(-2000.25), SC_(-3000.25), SC_(4.4839337747691491664160512023395421809359987226033)}, {SC_(-2000.25), SC_(-2000.25), SC_(2.7189621930215506574513067999560632490479415009312)}, {SC_(-2000.25), SC_(-1000.25), SC_(1.6489274597179783576990334595113842446837612984637)}, {SC_(-2000.25), SC_(-0.25), SC_(1.0001249921917327582862611787131655498347374570534)}, {SC_(-2000.25), SC_(999.75), SC_(0.60668227026545158239160396829655637629150925375338)}, {SC_(-2000.25), SC_(1999.75), SC_(0.36806334257221693918034782383378883527025156371561)}, {SC_(-2000.25), SC_(2999.75), SC_(0.22332534501129284494390126255673652609551748168190)}, + {SC_(-1000.25), SC_(-3000.25), SC_(20.166446124477887217071731764646075153043381039788)}, {SC_(-1000.25), SC_(-2000.25), SC_(7.4020458503369278169047297486603438977667911685454)}, {SC_(-1000.25), SC_(-1000.25), SC_(2.7196441509079794494949614854599837302121989306579)}, {SC_(-1000.25), SC_(-0.25), SC_(1.0002499687838699774154682651484339556685615759685)}, {SC_(-1000.25), SC_(999.75), SC_(0.36824716733727347121444226394598870471042459648225)}, {SC_(-1000.25), SC_(1999.75), SC_(0.13570722006231629831063556421918175042108431159794)}, {SC_(-1000.25), SC_(2999.75), SC_(0.050060760590698310218753944980019339485378387824874)}, + {SC_(-0.25), SC_(-3000.25), SC_(39.685967392358901409909771219162599134220400862017)}, {SC_(-0.25), SC_(-2000.25), SC_(42.589937906362216198404510603403791771675699582116)}, {SC_(-0.25), SC_(-1000.25), SC_(1.4484609059158563254575589143006796017980497019204)}, {SC_(-0.25), SC_(-0.25), SC_(1.8410918501257885334774939532672407660260932497146)}, {SC_(-0.25), SC_(999.75), { SC_(-5.3078213765443269585708737506907854937173735893125e28) }}, {SC_(-0.25), SC_(1999.75), { SC_(-1.6498771472158460728992067699031965281087177690027e40) }}, {SC_(-0.25), SC_(2999.75), { SC_(-1.0342523123207467753461742345710027886431938539942e49) }}, + {SC_(999.75), SC_(-3000.25), SC_(0.049513102314849089089531272645546209356947364503693)}, {SC_(999.75), SC_(-2000.25), SC_(0.13496287555692361766808821455578584612232691847368)}, {SC_(999.75), SC_(-1000.25), SC_(0.36751140843960205420082946295266027163679747815171)}, {SC_(999.75), SC_(-0.25), SC_(0.99974996871616159816731662896550362438335384158579)}, {SC_(999.75), SC_(999.75), SC_(2.7169258487486599795213704968075693651612629011137)}, {SC_(999.75), SC_(1999.75), SC_(7.3761834991423536961023917193398034973414374920724)}, {SC_(999.75), SC_(2999.75), SC_(20.005752597249211333476751024563565243417014841294)}, + {SC_(1999.75), SC_(-3000.25), SC_(0.22293486720691878909733143630236248153349047769471)}, {SC_(1999.75), SC_(-2000.25), SC_(0.36769546312912945115503617374640562125026416585911)}, {SC_(1999.75), SC_(-1000.25), SC_(0.60637900493362087991135014744274427763909334895653)}, {SC_(1999.75), SC_(-0.25), SC_(0.99987499218326921518475889542511447173965188103381)}, {SC_(1999.75), SC_(999.75), SC_(1.6485152791856256326123376967789594359576461105786)}, {SC_(1999.75), SC_(1999.75), SC_(2.7176030495659989803962821602572269144434394541710)}, {SC_(1999.75), SC_(2999.75), SC_(4.4794520688544123057192034298995619570726142156707)}, + {SC_(2999.75), SC_(-3000.25), SC_(0.36775679765931355729803647936785321939707796855752)}, {SC_(2999.75), SC_(-2000.25), SC_(0.51330776829498225175667208437735196999610552588713)}, {SC_(2999.75), SC_(-1000.25), SC_(0.71643841877594009499597675075273584291266593195763)}, {SC_(2999.75), SC_(-0.25), SC_(0.99991666319319078123594326189204686183292938405581)}, {SC_(2999.75), SC_(999.75), SC_(1.3955090626492889292016427682906884819283986108381)}, {SC_(2999.75), SC_(1999.75), SC_(1.9475357570910758403975860806320078849962664462281)}, {SC_(2999.75), SC_(2999.75), SC_(2.7178291334815105222784044221313459883434486746788)}, + } }; + + do_test_0F1(hypergeometric_0F1_integer_data, type_name, "Integer Arguments"); + do_test_0F1(hypergeometric_0F1_real_data, type_name, "Real Arguments"); + do_test_0F1(hypergeometric_0F1_large_data, type_name, "Large Arguments"); +} + diff --git a/test/test_1F0.cpp b/test/test_1F0.cpp new file mode 100644 index 000000000..1a971cb33 --- /dev/null +++ b/test/test_1F0.cpp @@ -0,0 +1,60 @@ +// (C) Copyright John Maddock 2006. +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#include +#include "test_1F0.hpp" + +#include +#include + +void expected_results() +{ + // + // Define the max and mean errors expected for + // various compilers and platforms. + // + const char* largest_type; +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + if(boost::math::policies::digits >() == boost::math::policies::digits >()) + { + largest_type = "(long\\s+)?double"; + } + else + { + largest_type = "long double"; + } +#else + largest_type = "(long\\s+)?double"; +#endif + + // + // Finish off by printing out the compiler/stdlib/platform names, + // we do this to make it easier to mark up expected error rates. + // + std::cout << "Tests run with " << BOOST_COMPILER << ", " + << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl; +} + +BOOST_AUTO_TEST_CASE( test_main ) +{ + expected_results(); + BOOST_MATH_CONTROL_FP; + +#ifndef BOOST_MATH_BUGGY_LARGE_FLOAT_CONSTANTS + test_spots(0.0F); +#endif + test_spots(0.0); +#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS + test_spots(0.0L); +#ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS + test_spots(boost::math::concepts::real_concept(0.1)); +#endif +#endif + test_spots(boost::multiprecision::cpp_bin_float_quad()); + test_spots(boost::multiprecision::cpp_dec_float_50()); +} + + + diff --git a/test/test_1F0.hpp b/test/test_1F0.hpp new file mode 100644 index 000000000..03c0c2dff --- /dev/null +++ b/test/test_1F0.hpp @@ -0,0 +1,60 @@ +// Copyright John Maddock 2006. +// Copyright Paul A. Bristow 2007, 2009 +// Use, modification and distribution are subject to the +// Boost Software License, Version 1.0. (See accompanying file +// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) + +#define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error + +#include +#include +#define BOOST_TEST_MAIN +#include +#include +#include +#include +#include +#include +#include +#include "functor.hpp" + +#include "handle_test_result.hpp" +#include "table_type.hpp" + +#include + +#ifndef SC_ +#define SC_(x) static_cast::type>(BOOST_JOIN(x, L)) +#endif + + +template +void test_spots(T) +{ + using std::pow; + // + // basic sanity checks, tolerance is 10 epsilon expressed as a percentage: + // + T tolerance = boost::math::tools::epsilon() * 1000; + + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(-3), T(2)), T(-1), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(-3), T(4)), T(-27), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(-3), T(0.5)), T(0.125), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(3), T(0.5)), T(8), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(3), T(2)), T(-1), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(3), T(4)), T(T(-1) / 27), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(3), T(-0.5)), pow(T(1.5), -3), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(3), T(-2)), T(1 / T(27)), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(3), T(-4)), T(T(1) / 125), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(-3), T(-0.5)), pow(T(1.5), 3), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(-3), T(-2)), T(27), tolerance); + BOOST_CHECK_CLOSE(boost::math::hypergeometric_1f0(T(-3), T(-4)), T(125), tolerance); + + BOOST_CHECK_THROW(boost::math::hypergeometric_1f0(T(3), T(1)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1f0(T(-3), T(1)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1f0(T(3.25), T(1)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1f0(T(-3.25), T(1)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1f0(T(3.25), T(2)), std::domain_error); + BOOST_CHECK_THROW(boost::math::hypergeometric_1f0(T(-3.25), T(2)), std::domain_error); +} +