1 2 /////////////////////////////////////////////////////////////////////////////// 3 // Copyright 2018 John Maddock 4 // Distributed under the Boost 5 // Software License, Version 1.0. (See accompanying file 6 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 7 // 8 #ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP 9 #define BOOST_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP 10 11 #include <boost/math/tools/series.hpp> 12 13 // 14 // This file implements the addition theorems for 1F1 on z, specifically 15 // each function returns 1F1[a, b, z + k] for some integer k - there's 16 // no particular reason why k needs to be an integer, but no reason why 17 // it shouldn't be either. 18 // 19 // The functions are named hypergeometric_1f1_recurrence_on_z_[plus|minus|zero]_[plus|minus|zero] 20 // where a "plus" indicates forward recurrence, minus backwards recurrence, and zero no recurrence. 21 // So for example hypergeometric_1f1_recurrence_on_z_zero_plus uses forward recurrence on b and 22 // hypergeometric_1f1_recurrence_on_z_minus_minus uses backwards recurrence on both a and b. 23 // 24 // See https://dlmf.nist.gov/13.13 25 // 26 27 namespace boost { namespace math { namespace detail { 28 29 // 30 // This works moderately well for a < 0, but has some very strange behaviour with 31 // strings of values of the same sign followed by a sign switch then another 32 // series all the same sign and so on.... doesn't converge smoothly either 33 // but rises and falls in wave-like behaviour.... very slow to converge... 34 // 35 template <class T, class Policy> 36 struct hypergeometric_1f1_recurrence_on_z_minus_zero_series 37 { 38 typedef T result_type; 39 hypergeometric_1f1_recurrence_on_z_minus_zero_seriesboost::math::detail::hypergeometric_1f1_recurrence_on_z_minus_zero_series40 hypergeometric_1f1_recurrence_on_z_minus_zero_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) 41 : term(1), b_minus_a_plus_n(b - a), a_(a), b_(b), z_(z), n(0), k(k_) 42 { 43 BOOST_MATH_STD_USING 44 int scale1(0), scale2(0); 45 M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol, scale1); 46 M_next = boost::math::detail::hypergeometric_1F1_imp(T(a - 1), b, z, pol, scale2); 47 if (scale1 != scale2) 48 M_next *= exp(scale2 - scale1); 49 if (M > 1e10f) 50 { 51 // rescale: 52 int rescale = itrunc(log(fabs(M))); 53 M *= exp(T(-rescale)); 54 M_next *= exp(T(-rescale)); 55 scale1 += rescale; 56 } 57 scaling = scale1; 58 } operator ()boost::math::detail::hypergeometric_1f1_recurrence_on_z_minus_zero_series59 T operator()() 60 { 61 T result = term * M; 62 term *= b_minus_a_plus_n * k / ((z_ + k) * ++n); 63 b_minus_a_plus_n += 1; 64 T M2 = -((2 * (a_ - n) - b_ + z_) * M_next - (a_ - n) * M) / (b_ - (a_ - n)); 65 M = M_next; 66 M_next = M2; 67 68 return result; 69 } scaleboost::math::detail::hypergeometric_1f1_recurrence_on_z_minus_zero_series70 int scale()const { return scaling; } 71 private: 72 T term, b_minus_a_plus_n, M, M_next, a_, b_, z_; 73 int n, k, scaling; 74 }; 75 76 template <class T, class Policy> hypergeometric_1f1_recurrence_on_z_minus_zero(const T & a,const T & b,const T & z,int k,const Policy & pol,int & log_scaling)77 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) 78 { 79 BOOST_MATH_STD_USING 80 BOOST_ASSERT((z + k) / z > 0.5f); 81 hypergeometric_1f1_recurrence_on_z_minus_zero_series<T, Policy> s(a, b, z, k, pol); 82 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 83 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 84 log_scaling += s.scale(); 85 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); 86 return result * exp(T(k)) * pow(z / (z + k), b - a); 87 } 88 89 #if 0 90 91 // 92 // These are commented out as they are currently unused, but may find use in the future: 93 // 94 95 template <class T, class Policy> 96 struct hypergeometric_1f1_recurrence_on_z_plus_plus_series 97 { 98 typedef T result_type; 99 100 hypergeometric_1f1_recurrence_on_z_plus_plus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) 101 : term(1), a_plus_n(a), b_plus_n(b), z_(z), n(0), k(k_) 102 { 103 M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); 104 M_next = boost::math::detail::hypergeometric_1F1_imp(a + 1, b + 1, z, pol); 105 } 106 T operator()() 107 { 108 T result = term * M; 109 term *= a_plus_n * k / (b_plus_n * ++n); 110 a_plus_n += 1; 111 b_plus_n += 1; 112 // The a_plus_n == 0 case below isn't actually correct, but doesn't matter as that term will be zero 113 // anyway, we just need to not divide by zero and end up with a NaN in the result. 114 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_); 115 M = M_next; 116 M_next = M2; 117 return result; 118 } 119 T term, a_plus_n, b_plus_n, M, M_next, z_; 120 int n, k; 121 }; 122 123 template <class T, class Policy> 124 T hypergeometric_1f1_recurrence_on_z_plus_plus(const T& a, const T& b, const T& z, int k, const Policy& pol) 125 { 126 hypergeometric_1f1_recurrence_on_z_plus_plus_series<T, Policy> s(a, b, z, k, pol); 127 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 128 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 129 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); 130 return result; 131 } 132 133 template <class T, class Policy> 134 struct hypergeometric_1f1_recurrence_on_z_zero_minus_series 135 { 136 typedef T result_type; 137 138 hypergeometric_1f1_recurrence_on_z_zero_minus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) 139 : 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_) 140 { 141 M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); 142 M_next = boost::math::detail::hypergeometric_1F1_imp(a, b - 1, z, pol); 143 } 144 T operator()() 145 { 146 BOOST_MATH_STD_USING 147 T result = term * M; 148 term *= b_pochhammer * x_k_power / ++n; 149 b_pochhammer += 1; 150 b_minus_n -= 1; 151 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)); 152 M = M_next; 153 M_next = M2; 154 return result; 155 } 156 T term, b_pochhammer, x_k_power, M, M_next, b_minus_n, a_, z_, b_; 157 int n, k; 158 }; 159 160 template <class T, class Policy> 161 T hypergeometric_1f1_recurrence_on_z_zero_minus(const T& a, const T& b, const T& z, int k, const Policy& pol) 162 { 163 BOOST_MATH_STD_USING 164 BOOST_ASSERT(abs(k) < fabs(z)); 165 hypergeometric_1f1_recurrence_on_z_zero_minus_series<T, Policy> s(a, b, z, k, pol); 166 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 167 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 168 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); 169 return result * pow((z + k) / z, 1 - b); 170 } 171 172 template <class T, class Policy> 173 struct hypergeometric_1f1_recurrence_on_z_plus_zero_series 174 { 175 typedef T result_type; 176 177 hypergeometric_1f1_recurrence_on_z_plus_zero_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) 178 : term(1), a_pochhammer(a), z_plus_k(z + k_), b_(b), a_(a), z_(z), n(0), k(k_) 179 { 180 M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); 181 M_next = boost::math::detail::hypergeometric_1F1_imp(a + 1, b, z, pol); 182 } 183 T operator()() 184 { 185 T result = term * M; 186 term *= a_pochhammer * k / (++n * z_plus_k); 187 a_pochhammer += 1; 188 T M2 = (a_pochhammer == -1) ? 1 : (a_pochhammer == 0) ? 0 : (M_next * (2 * a_pochhammer - b_ + z_) + (b_ - a_pochhammer) * M) / a_pochhammer; 189 M = M_next; 190 M_next = M2; 191 192 return result; 193 } 194 T term, a_pochhammer, z_plus_k, M, M_next, b_minus_n, a_, b_, z_; 195 int n, k; 196 }; 197 198 template <class T, class Policy> 199 T hypergeometric_1f1_recurrence_on_z_plus_zero(const T& a, const T& b, const T& z, int k, const Policy& pol) 200 { 201 BOOST_MATH_STD_USING 202 BOOST_ASSERT(k / z > -0.5f); 203 //BOOST_ASSERT(floor(a) != a || a > 0); 204 hypergeometric_1f1_recurrence_on_z_plus_zero_series<T, Policy> s(a, b, z, k, pol); 205 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 206 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 207 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); 208 return result * pow(z / (z + k), a); 209 } 210 211 template <class T, class Policy> 212 struct hypergeometric_1f1_recurrence_on_z_zero_plus_series 213 { 214 typedef T result_type; 215 216 hypergeometric_1f1_recurrence_on_z_zero_plus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) 217 : term(1), b_minus_a_plus_n(b - a), b_plus_n(b), a_(a), z_(z), n(0), k(k_) 218 { 219 M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); 220 M_next = boost::math::detail::hypergeometric_1F1_imp(a, b + 1, z, pol); 221 } 222 T operator()() 223 { 224 T result = term * M; 225 term *= b_minus_a_plus_n * -k / (b_plus_n * ++n); 226 b_minus_a_plus_n += 1; 227 b_plus_n += 1; 228 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); 229 M = M_next; 230 M_next = M2; 231 232 return result; 233 } 234 T term, b_minus_a_plus_n, M, M_next, b_minus_n, a_, b_plus_n, z_; 235 int n, k; 236 }; 237 238 template <class T, class Policy> 239 T hypergeometric_1f1_recurrence_on_z_zero_plus(const T& a, const T& b, const T& z, int k, const Policy& pol) 240 { 241 BOOST_MATH_STD_USING 242 hypergeometric_1f1_recurrence_on_z_zero_plus_series<T, Policy> s(a, b, z, k, pol); 243 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 244 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 245 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); 246 return result * exp(T(k)); 247 } 248 // 249 // I'm unable to find any situation where this series isn't divergent and therefore 250 // is probably quite useless: 251 // 252 template <class T, class Policy> 253 struct hypergeometric_1f1_recurrence_on_z_minus_minus_series 254 { 255 typedef T result_type; 256 257 hypergeometric_1f1_recurrence_on_z_minus_minus_series(const T& a, const T& b, const T& z, int k_, const Policy& pol) 258 : term(1), one_minus_b_plus_n(1 - b), a_(a), b_(b), z_(z), n(0), k(k_) 259 { 260 M = boost::math::detail::hypergeometric_1F1_imp(a, b, z, pol); 261 M_next = boost::math::detail::hypergeometric_1F1_imp(a - 1, b - 1, z, pol); 262 } 263 T operator()() 264 { 265 T result = term * M; 266 term *= one_minus_b_plus_n * k / (z_ * ++n); 267 one_minus_b_plus_n += 1; 268 T M2 = -((b_ - n) * (1 - b_ + n + z_) * M_next - (a_ - n) * z_ * M) / ((b_ - n) * (b_ - n - 1)); 269 M = M_next; 270 M_next = M2; 271 272 return result; 273 } 274 T term, one_minus_b_plus_n, M, M_next, a_, b_, z_; 275 int n, k; 276 }; 277 278 template <class T, class Policy> 279 T hypergeometric_1f1_recurrence_on_z_minus_minus(const T& a, const T& b, const T& z, int k, const Policy& pol) 280 { 281 BOOST_MATH_STD_USING 282 hypergeometric_1f1_recurrence_on_z_minus_minus_series<T, Policy> s(a, b, z, k, pol); 283 boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>(); 284 T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter); 285 boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_recurrence_on_z_plus_plus<%1%>(%1%,%1%,%1%)", max_iter, pol); 286 return result * exp(T(k)) * pow((z + k) / z, 1 - b); 287 } 288 #endif 289 290 } } } // namespaces 291 292 #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_ADDITION_THEOREMS_ON_Z_HPP 293