123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292 |
- ///////////////////////////////////////////////////////////////////////////////
- // Copyright 2018 John Maddock
- // 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_ADDITION_THEOREMS_ON_Z_HPP
- #define BOOST_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP
- #include <boost/math/tools/series.hpp>
- //
- // This file implements the addition theorems for 1F1 on z, specifically
- // each function returns 1F1[a, b, z + k] for some integer k - there's
- // no particular reason why k needs to be an integer, but no reason why
- // it shouldn't be either.
- //
- // The functions are named hypergeometric_1f1_recurrence_on_z_[plus|minus|zero]_[plus|minus|zero]
- // where a "plus" indicates forward recurrence, minus backwards recurrence, and zero no recurrence.
- // So for example hypergeometric_1f1_recurrence_on_z_zero_plus uses forward recurrence on b and
- // hypergeometric_1f1_recurrence_on_z_minus_minus uses backwards recurrence on both a and b.
- //
- // See https://dlmf.nist.gov/13.13
- //
- namespace boost { namespace math { namespace detail {
- //
- // This works moderately well for a < 0, but has some very strange behaviour with
- // strings of values of the same sign followed by a sign switch then another
- // series all the same sign and so on.... doesn't converge smoothly either
- // but rises and falls in wave-like behaviour.... very slow to converge...
- //
- template <class T, class Policy>
- struct hypergeometric_1f1_recurrence_on_z_minus_zero_series
- {
- typedef T result_type;
- hypergeometric_1f1_recurrence_on_z_minus_zero_series(const T& a, const T& b, const T& z, int k_, const Policy& pol)
- : term(1), b_minus_a_plus_n(b - a), a_(a), b_(b), z_(z), n(0), k(k_)
- {
- BOOST_MATH_STD_USING
- int scale1(0), scale2(0);
- M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol, scale1);
- M_next = boost::math::detail::hypergeometric_1F1_imp(T(a - 1), b, z, pol, scale2);
- if (scale1 != scale2)
- M_next *= exp(scale2 - scale1);
- if (M > 1e10f)
- {
- // rescale:
- int rescale = itrunc(log(fabs(M)));
- M *= exp(T(-rescale));
- M_next *= exp(T(-rescale));
- scale1 += rescale;
- }
- scaling = scale1;
- }
- T operator()()
- {
- T result = term * M;
- term *= b_minus_a_plus_n * k / ((z_ + k) * ++n);
- b_minus_a_plus_n += 1;
- T M2 = -((2 * (a_ - n) - b_ + z_) * M_next - (a_ - n) * M) / (b_ - (a_ - n));
- M = M_next;
- M_next = M2;
- return result;
- }
- int scale()const { return scaling; }
- private:
- T term, b_minus_a_plus_n, M, M_next, a_, b_, z_;
- int n, k, scaling;
- };
- template <class T, class Policy>
- T hypergeometric_1f1_recurrence_on_z_minus_zero(const T& a, const T& b, const T& z, int k, const Policy& pol, int& log_scaling)
- {
- BOOST_MATH_STD_USING
- BOOST_ASSERT((z + k) / z > 0.5f);
- hypergeometric_1f1_recurrence_on_z_minus_zero_series<T, Policy> s(a, b, z, k, pol);
- boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
- T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
- log_scaling += s.scale();
- boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol);
- return result * exp(T(k)) * pow(z / (z + k), b - a);
- }
- #if 0
- //
- // These are commented out as they are currently unused, but may find use in the future:
- //
- template <class T, class Policy>
- struct hypergeometric_1f1_recurrence_on_z_plus_plus_series
- {
- typedef T result_type;
- hypergeometric_1f1_recurrence_on_z_plus_plus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol)
- : term(1), a_plus_n(a), b_plus_n(b), z_(z), n(0), k(k_)
- {
- M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol);
- M_next = boost::math::detail::hypergeometric_1F1_imp(a + 1, b + 1, z, pol);
- }
- T operator()()
- {
- T result = term * M;
- term *= a_plus_n * k / (b_plus_n * ++n);
- a_plus_n += 1;
- b_plus_n += 1;
- // The a_plus_n == 0 case below isn't actually correct, but doesn't matter as that term will be zero
- // anyway, we just need to not divde by zero and end up with a NaN in the result.
- T M2 = (a_plus_n == -1) ? 1 : (a_plus_n == 0) ? 0 : (M_next * b_plus_n * (1 - b_plus_n + z_) + b_plus_n * (b_plus_n - 1) * M) / (a_plus_n * z_);
- M = M_next;
- M_next = M2;
- return result;
- }
- T term, a_plus_n, b_plus_n, M, M_next, z_;
- int n, k;
- };
- template <class T, class Policy>
- T hypergeometric_1f1_recurrence_on_z_plus_plus(const T& a, const T& b, const T& z, int k, const Policy& pol)
- {
- hypergeometric_1f1_recurrence_on_z_plus_plus_series<T, Policy> s(a, b, z, k, pol);
- boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
- T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
- boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol);
- return result;
- }
- template <class T, class Policy>
- struct hypergeometric_1f1_recurrence_on_z_zero_minus_series
- {
- typedef T result_type;
- hypergeometric_1f1_recurrence_on_z_zero_minus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol)
- : term(1), b_pochhammer(1 - b), x_k_power(-k_ / z), b_minus_n(b), a_(a), z_(z), b_(b), n(0), k(k_)
- {
- M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol);
- M_next = boost::math::detail::hypergeometric_1F1_imp(a, b - 1, z, pol);
- }
- T operator()()
- {
- BOOST_MATH_STD_USING
- T result = term * M;
- term *= b_pochhammer * x_k_power / ++n;
- b_pochhammer += 1;
- b_minus_n -= 1;
- T M2 = (M_next * b_minus_n * (1 - b_minus_n - z_) + z_ * (b_minus_n - a_) * M) / (-b_minus_n * (b_minus_n - 1));
- M = M_next;
- M_next = M2;
- return result;
- }
- T term, b_pochhammer, x_k_power, M, M_next, b_minus_n, a_, z_, b_;
- int n, k;
- };
- template <class T, class Policy>
- T hypergeometric_1f1_recurrence_on_z_zero_minus(const T& a, const T& b, const T& z, int k, const Policy& pol)
- {
- BOOST_MATH_STD_USING
- BOOST_ASSERT(abs(k) < fabs(z));
- hypergeometric_1f1_recurrence_on_z_zero_minus_series<T, Policy> s(a, b, z, k, pol);
- boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
- T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
- boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol);
- return result * pow((z + k) / z, 1 - b);
- }
- template <class T, class Policy>
- struct hypergeometric_1f1_recurrence_on_z_plus_zero_series
- {
- typedef T result_type;
- hypergeometric_1f1_recurrence_on_z_plus_zero_series(const T& a, const T& b, const T& z, int k_, const Policy& pol)
- : term(1), a_pochhammer(a), z_plus_k(z + k_), b_(b), a_(a), z_(z), n(0), k(k_)
- {
- M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol);
- M_next = boost::math::detail::hypergeometric_1F1_imp(a + 1, b, z, pol);
- }
- T operator()()
- {
- T result = term * M;
- term *= a_pochhammer * k / (++n * z_plus_k);
- a_pochhammer += 1;
- T M2 = (a_pochhammer == -1) ? 1 : (a_pochhammer == 0) ? 0 : (M_next * (2 * a_pochhammer - b_ + z_) + (b_ - a_pochhammer) * M) / a_pochhammer;
- M = M_next;
- M_next = M2;
- return result;
- }
- T term, a_pochhammer, z_plus_k, M, M_next, b_minus_n, a_, b_, z_;
- int n, k;
- };
- template <class T, class Policy>
- T hypergeometric_1f1_recurrence_on_z_plus_zero(const T& a, const T& b, const T& z, int k, const Policy& pol)
- {
- BOOST_MATH_STD_USING
- BOOST_ASSERT(k / z > -0.5f);
- //BOOST_ASSERT(floor(a) != a || a > 0);
- hypergeometric_1f1_recurrence_on_z_plus_zero_series<T, Policy> s(a, b, z, k, pol);
- boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
- T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
- boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol);
- return result * pow(z / (z + k), a);
- }
- template <class T, class Policy>
- struct hypergeometric_1f1_recurrence_on_z_zero_plus_series
- {
- typedef T result_type;
- hypergeometric_1f1_recurrence_on_z_zero_plus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol)
- : term(1), b_minus_a_plus_n(b - a), b_plus_n(b), a_(a), z_(z), n(0), k(k_)
- {
- M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol);
- M_next = boost::math::detail::hypergeometric_1F1_imp(a, b + 1, z, pol);
- }
- T operator()()
- {
- T result = term * M;
- term *= b_minus_a_plus_n * -k / (b_plus_n * ++n);
- b_minus_a_plus_n += 1;
- b_plus_n += 1;
- T M2 = (b_plus_n * (b_plus_n - 1) * M + b_plus_n * (1 - b_plus_n - z_) * M_next) / (-z_ * b_minus_a_plus_n);
- M = M_next;
- M_next = M2;
- return result;
- }
- T term, b_minus_a_plus_n, M, M_next, b_minus_n, a_, b_plus_n, z_;
- int n, k;
- };
- template <class T, class Policy>
- T hypergeometric_1f1_recurrence_on_z_zero_plus(const T& a, const T& b, const T& z, int k, const Policy& pol)
- {
- BOOST_MATH_STD_USING
- hypergeometric_1f1_recurrence_on_z_zero_plus_series<T, Policy> s(a, b, z, k, pol);
- boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
- T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
- boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol);
- return result * exp(T(k));
- }
- //
- // I'm unable to find any situation where this series isn't divergent and therefore
- // is probably quite useless:
- //
- template <class T, class Policy>
- struct hypergeometric_1f1_recurrence_on_z_minus_minus_series
- {
- typedef T result_type;
- hypergeometric_1f1_recurrence_on_z_minus_minus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol)
- : term(1), one_minus_b_plus_n(1 - b), a_(a), b_(b), z_(z), n(0), k(k_)
- {
- M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol);
- M_next = boost::math::detail::hypergeometric_1F1_imp(a - 1, b - 1, z, pol);
- }
- T operator()()
- {
- T result = term * M;
- term *= one_minus_b_plus_n * k / (z_ * ++n);
- one_minus_b_plus_n += 1;
- T M2 = -((b_ - n) * (1 - b_ + n + z_) * M_next - (a_ - n) * z_ * M) / ((b_ - n) * (b_ - n - 1));
- M = M_next;
- M_next = M2;
- return result;
- }
- T term, one_minus_b_plus_n, M, M_next, a_, b_, z_;
- int n, k;
- };
- template <class T, class Policy>
- T hypergeometric_1f1_recurrence_on_z_minus_minus(const T& a, const T& b, const T& z, int k, const Policy& pol)
- {
- BOOST_MATH_STD_USING
- hypergeometric_1f1_recurrence_on_z_minus_minus_series<T, Policy> s(a, b, z, k, pol);
- boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
- T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
- boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol);
- return result * exp(T(k)) * pow((z + k) / z, 1 - b);
- }
- #endif
- } } } // namespaces
- #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP
|