1 ///////////////////////////////////////////////////////////////////////////////
2 // Copyright 2014 Anton Bikineev
3 // Copyright 2014 Christopher Kormanyos
4 // Copyright 2014 John Maddock
5 // Copyright 2014 Paul Bristow
6 // Distributed under the Boost
7 // Software License, Version 1.0. (See accompanying file
8 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
9
10 #ifndef BOOST_MATH_HYPERGEOMETRIC_0F1_HPP
11 #define BOOST_MATH_HYPERGEOMETRIC_0F1_HPP
12
13 #include <boost/math/policies/policy.hpp>
14 #include <boost/math/policies/error_handling.hpp>
15 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
16 #include <boost/math/special_functions/detail/hypergeometric_0F1_bessel.hpp>
17
18 namespace boost { namespace math { namespace detail {
19
20
21 template <class T>
22 struct hypergeometric_0F1_cf
23 {
24 //
25 // We start this continued fraction at b on index -1
26 // and treat the -1 and 0 cases as special cases.
27 // We do this to avoid adding the continued fraction result
28 // to 1 so that we can accurately evaluate for small results
29 // as well as large ones. See http://functions.wolfram.com/07.17.10.0002.01
30 //
31 T b, z;
32 int k;
hypergeometric_0F1_cfboost::math::detail::hypergeometric_0F1_cf33 hypergeometric_0F1_cf(T b_, T z_) : b(b_), z(z_), k(-2) {}
34 typedef std::pair<T, T> result_type;
35
operator ()boost::math::detail::hypergeometric_0F1_cf36 result_type operator()()
37 {
38 ++k;
39 if (k <= 0)
40 return std::make_pair(z / b, 1);
41 return std::make_pair(-z / ((k + 1) * (b + k)), 1 + z / ((k + 1) * (b + k)));
42 }
43 };
44
45 template <class T, class Policy>
46 T hypergeometric_0F1_cf_imp(T b, T z, const Policy& pol, const char* function)
47 {
48 hypergeometric_0F1_cf<T> evaluator(b, z);
49 std::uintmax_t max_iter = policies::get_max_series_iterations<Policy>();
50 T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter);
51 policies::check_series_iterations<T>(function, max_iter, pol);
52 return cf;
53 }
54
55
56 template <class T, class Policy>
hypergeometric_0F1_imp(const T & b,const T & z,const Policy & pol)57 inline T hypergeometric_0F1_imp(const T& b, const T& z, const Policy& pol)
58 {
59 const char* function = "boost::math::hypergeometric_0f1<%1%,%1%>(%1%, %1%)";
60 BOOST_MATH_STD_USING
61
62 // some special cases
63 if (z == 0)
64 return T(1);
65
66 if ((b <= 0) && (b == floor(b)))
67 return policies::raise_pole_error<T>(
68 function,
69 "Evaluation of 0f1 with nonpositive integer b = %1%.", b, pol);
70
71 if (z < -5 && b > -5)
72 {
73 // Series is alternating and divergent, need to do something else here,
74 // Bessel function relation is much more accurate, unless |b| is similarly
75 // large to |z|, otherwise the CF formula suffers from cancellation when
76 // the result would be very small.
77 if (fabs(z / b) > 4)
78 return hypergeometric_0F1_bessel(b, z, pol);
79 return hypergeometric_0F1_cf_imp(b, z, pol, function);
80 }
81 // evaluation through Taylor series looks
82 // more precisious than Bessel relation:
83 // detail::hypergeometric_0f1_bessel(b, z, pol);
84 return detail::hypergeometric_0F1_generic_series(b, z, pol);
85 }
86
87 } // namespace detail
88
89 template <class T1, class T2, class Policy>
hypergeometric_0F1(T1 b,T2 z,const Policy &)90 inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z, const Policy& /* pol */)
91 {
92 BOOST_FPU_EXCEPTION_GUARD
93 typedef typename tools::promote_args<T1, T2>::type result_type;
94 typedef typename policies::evaluation<result_type, Policy>::type value_type;
95 typedef typename policies::normalise<
96 Policy,
97 policies::promote_float<false>,
98 policies::promote_double<false>,
99 policies::discrete_quantile<>,
100 policies::assert_undefined<> >::type forwarding_policy;
101 return policies::checked_narrowing_cast<result_type, Policy>(
102 detail::hypergeometric_0F1_imp<value_type>(
103 static_cast<value_type>(b),
104 static_cast<value_type>(z),
105 forwarding_policy()),
106 "boost::math::hypergeometric_0F1<%1%>(%1%,%1%)");
107 }
108
109 template <class T1, class T2>
hypergeometric_0F1(T1 b,T2 z)110 inline typename tools::promote_args<T1, T2>::type hypergeometric_0F1(T1 b, T2 z)
111 {
112 return hypergeometric_0F1(b, z, policies::policy<>());
113 }
114
115
116 } } // namespace boost::math
117
118 #endif // BOOST_MATH_HYPERGEOMETRIC_HPP
119