]>
Commit | Line | Data |
---|---|---|
92f5a8d4 TL |
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_2F0_HPP | |
11 | #define BOOST_MATH_HYPERGEOMETRIC_2F0_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/laguerre.hpp> | |
17 | #include <boost/math/special_functions/hermite.hpp> | |
18 | #include <boost/math/tools/fraction.hpp> | |
19 | ||
20 | namespace boost { namespace math { namespace detail { | |
21 | ||
22 | template <class T> | |
23 | struct hypergeometric_2F0_cf | |
24 | { | |
25 | // | |
26 | // We start this continued fraction at b on index -1 | |
27 | // and treat the -1 and 0 cases as special cases. | |
28 | // We do this to avoid adding the continued fraction result | |
29 | // to 1 so that we can accurately evaluate for small results | |
30 | // as well as large ones. See http://functions.wolfram.com/07.31.10.0002.01 | |
31 | // | |
32 | T a1, a2, z; | |
33 | int k; | |
34 | hypergeometric_2F0_cf(T a1_, T a2_, T z_) : a1(a1_), a2(a2_), z(z_), k(-2) {} | |
35 | typedef std::pair<T, T> result_type; | |
36 | ||
37 | result_type operator()() | |
38 | { | |
39 | ++k; | |
40 | if (k <= 0) | |
41 | return std::make_pair(z * a1 * a2, 1); | |
42 | return std::make_pair(-z * (a1 + k) * (a2 + k) / (k + 1), 1 + z * (a1 + k) * (a2 + k) / (k + 1)); | |
43 | } | |
44 | }; | |
45 | ||
46 | template <class T, class Policy> | |
47 | T hypergeometric_2F0_cf_imp(T a1, T a2, T z, const Policy& pol, const char* function) | |
48 | { | |
49 | using namespace boost::math; | |
50 | hypergeometric_2F0_cf<T> evaluator(a1, a2, z); | |
51 | boost::uintmax_t max_iter = policies::get_max_series_iterations<Policy>(); | |
52 | T cf = tools::continued_fraction_b(evaluator, policies::get_epsilon<T, Policy>(), max_iter); | |
53 | policies::check_series_iterations<T>(function, max_iter, pol); | |
54 | return cf; | |
55 | } | |
56 | ||
57 | ||
58 | template <class T, class Policy> | |
59 | inline T hypergeometric_2F0_imp(T a1, T a2, const T& z, const Policy& pol, bool asymptotic = false) | |
60 | { | |
61 | // | |
62 | // The terms in this series go to infinity unless one of a1 and a2 is a negative integer. | |
63 | // | |
64 | using std::swap; | |
65 | BOOST_MATH_STD_USING | |
66 | ||
67 | static const char* const function = "boost::math::hypergeometric_2F0<%1%,%1%,%1%>(%1%,%1%,%1%)"; | |
68 | ||
69 | if (z == 0) | |
70 | return 1; | |
71 | ||
72 | bool is_a1_integer = (a1 == floor(a1)); | |
73 | bool is_a2_integer = (a2 == floor(a2)); | |
74 | ||
75 | if (!asymptotic && !is_a1_integer && !is_a2_integer) | |
76 | return boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
77 | if (!is_a1_integer || (a1 > 0)) | |
78 | { | |
79 | swap(a1, a2); | |
80 | swap(is_a1_integer, is_a2_integer); | |
81 | } | |
82 | // | |
83 | // At this point a1 must be a negative integer: | |
84 | // | |
85 | if(!asymptotic && (!is_a1_integer || (a1 > 0))) | |
86 | return boost::math::policies::raise_overflow_error<T>(function, 0, pol); | |
87 | // | |
88 | // Special cases first: | |
89 | // | |
90 | if (a1 == 0) | |
91 | return 1; | |
92 | if ((a1 == a2 - 0.5f) && (z < 0)) | |
93 | { | |
94 | // http://functions.wolfram.com/07.31.03.0083.01 | |
95 | int n = static_cast<int>(static_cast<boost::uintmax_t>(boost::math::lltrunc(-2 * a1))); | |
96 | T smz = sqrt(-z); | |
20effc67 | 97 | return pow(2 / smz, -n) * boost::math::hermite(n, 1 / smz, pol); |
92f5a8d4 TL |
98 | } |
99 | ||
100 | if (is_a1_integer && is_a2_integer) | |
101 | { | |
102 | if ((a1 < 1) && (a2 <= a1)) | |
103 | { | |
104 | const unsigned int n = static_cast<unsigned int>(static_cast<boost::uintmax_t>(boost::math::lltrunc(-a1))); | |
105 | const unsigned int m = static_cast<unsigned int>(static_cast<boost::uintmax_t>(boost::math::lltrunc(-a2 - n))); | |
106 | ||
107 | return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) * | |
108 | boost::math::laguerre(n, m, -(1 / z), pol); | |
109 | } | |
110 | else if ((a2 < 1) && (a1 <= a2)) | |
111 | { | |
112 | // function is symmetric for a1 and a2 | |
113 | const unsigned int n = static_cast<unsigned int>(static_cast<boost::uintmax_t>(boost::math::lltrunc(-a2))); | |
114 | const unsigned int m = static_cast<unsigned int>(static_cast<boost::uintmax_t>(boost::math::lltrunc(-a1 - n))); | |
115 | ||
116 | return (pow(z, T(n)) * boost::math::factorial<T>(n, pol)) * | |
117 | boost::math::laguerre(n, m, -(1 / z), pol); | |
118 | } | |
119 | } | |
120 | ||
121 | if ((a1 * a2 * z < 0) && (a2 < -5) && (fabs(a1 * a2 * z) > 0.5)) | |
122 | { | |
123 | // Series is alternating and maybe divergent at least for the first few terms | |
124 | // (until a2 goes positive), try the continued fraction: | |
125 | return hypergeometric_2F0_cf_imp(a1, a2, z, pol, function); | |
126 | } | |
127 | ||
128 | return detail::hypergeometric_2F0_generic_series(a1, a2, z, pol); | |
129 | } | |
130 | ||
131 | } // namespace detail | |
132 | ||
133 | template <class T1, class T2, class T3, class Policy> | |
134 | inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z, const Policy& /* pol */) | |
135 | { | |
136 | BOOST_FPU_EXCEPTION_GUARD | |
137 | typedef typename tools::promote_args<T1, T2, T3>::type result_type; | |
138 | typedef typename policies::evaluation<result_type, Policy>::type value_type; | |
139 | typedef typename policies::normalise< | |
140 | Policy, | |
141 | policies::promote_float<false>, | |
142 | policies::promote_double<false>, | |
143 | policies::discrete_quantile<>, | |
144 | policies::assert_undefined<> >::type forwarding_policy; | |
145 | return policies::checked_narrowing_cast<result_type, Policy>( | |
146 | detail::hypergeometric_2F0_imp<value_type>( | |
147 | static_cast<value_type>(a1), | |
148 | static_cast<value_type>(a2), | |
149 | static_cast<value_type>(z), | |
150 | forwarding_policy()), | |
151 | "boost::math::hypergeometric_2F0<%1%>(%1%,%1%,%1%)"); | |
152 | } | |
153 | ||
154 | template <class T1, class T2, class T3> | |
155 | inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_2F0(T1 a1, T2 a2, T3 z) | |
156 | { | |
157 | return hypergeometric_2F0(a1, a2, z, policies::policy<>()); | |
158 | } | |
159 | ||
160 | ||
161 | } } // namespace boost::math | |
162 | ||
163 | #endif // BOOST_MATH_HYPERGEOMETRIC_HPP |