boost/math/special_functions/detail/hypergeometric_1F1_addition_theorems_on_z.hpp

   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
  40         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         }
  59         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         }
  70         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>
  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 divde 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