boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp

   1
   2 ///////////////////////////////////////////////////////////////////////////////
   3 //  Copyright 2014 Anton Bikineev
   4 //  Copyright 2014 Christopher Kormanyos
   5 //  Copyright 2014 John Maddock
   6 //  Copyright 2014 Paul Bristow
   7 //  Distributed under the Boost
   8 //  Software License, Version 1.0. (See accompanying file
   9 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
  10
  11 #ifndef BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_
  12 #define BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_
  13
  14 #include <boost/math/special_functions/modf.hpp>
  15 #include <boost/math/special_functions/next.hpp>
  16
  17 #include <boost/math/tools/recurrence.hpp>
  18 #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
  19
  20   namespace boost { namespace math { namespace detail {
  21
  22   // forward declaration for initial values
  23   template <class T, class Policy>
  24   inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol);
  25
  26   template <class T, class Policy>
  27   inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling);
  28
  29   template <class T>
  30   struct hypergeometric_1F1_recurrence_a_coefficients
  31   {
  32     typedef boost::math::tuple<T, T, T> result_type;
  33
  34     hypergeometric_1F1_recurrence_a_coefficients(const T& a, const T& b, const T& z):
  35     a(a), b(b), z(z)
  36     {
  37     }
  38
  39     result_type operator()(boost::intmax_t i) const
  40     {
  41       const T ai = a + i;
  42
  43       const T an = b - ai;
  44       const T bn = (2 * ai - b + z);
  45       const T cn = -ai;
  46
  47       return boost::math::make_tuple(an, bn, cn);
  48     }
  49
  50   private:
  51     const T a, b, z;
  52     hypergeometric_1F1_recurrence_a_coefficients operator=(const hypergeometric_1F1_recurrence_a_coefficients&);
  53   };
  54
  55   template <class T>
  56   struct hypergeometric_1F1_recurrence_b_coefficients
  57   {
  58     typedef boost::math::tuple<T, T, T> result_type;
  59
  60     hypergeometric_1F1_recurrence_b_coefficients(const T& a, const T& b, const T& z):
  61     a(a), b(b), z(z)
  62     {
  63     }
  64
  65     result_type operator()(boost::intmax_t i) const
  66     {
  67       const T bi = b + i;
  68
  69       const T an = bi * (bi - 1);
  70       const T bn = bi * (1 - bi - z);
  71       const T cn = z * (bi - a);
  72
  73       return boost::math::make_tuple(an, bn, cn);
  74     }
  75
  76   private:
  77     const T a, b, z;
  78     hypergeometric_1F1_recurrence_b_coefficients& operator=(const hypergeometric_1F1_recurrence_b_coefficients&);
  79   };
  80   //
  81   // for use when we're recursing to a small b:
  82   //
  83   template <class T>
  84   struct hypergeometric_1F1_recurrence_small_b_coefficients
  85   {
  86      typedef boost::math::tuple<T, T, T> result_type;
  87
  88      hypergeometric_1F1_recurrence_small_b_coefficients(const T& a, const T& b, const T& z, int N) :
  89         a(a), b(b), z(z), N(N)
  90      {
  91      }
  92
  93      result_type operator()(boost::intmax_t i) const
  94      {
  95         const T bi = b + (i + N);
  96         const T bi_minus_1 = b + (i + N - 1);
  97
  98         const T an = bi * bi_minus_1;
  99         const T bn = bi * (-bi_minus_1 - z);
 100         const T cn = z * (bi - a);
 101
 102         return boost::math::make_tuple(an, bn, cn);
 103      }
 104
 105   private:
 106      hypergeometric_1F1_recurrence_small_b_coefficients operator=(const hypergeometric_1F1_recurrence_small_b_coefficients&);
 107      const T a, b, z;
 108      int N;
 109   };
 110
 111   template <class T>
 112   struct hypergeometric_1F1_recurrence_a_and_b_coefficients
 113   {
 114     typedef boost::math::tuple<T, T, T> result_type;
 115
 116     hypergeometric_1F1_recurrence_a_and_b_coefficients(const T& a, const T& b, const T& z, int offset = 0):
 117     a(a), b(b), z(z), offset(offset)
 118     {
 119     }
 120
 121     result_type operator()(boost::intmax_t i) const
 122     {
 123       const T ai = a + (offset + i);
 124       const T bi = b + (offset + i);
 125
 126       const T an = bi * (b + (offset + i - 1));
 127       const T bn = bi * (z - (b + (offset + i - 1)));
 128       const T cn = -ai * z;
 129
 130       return boost::math::make_tuple(an, bn, cn);
 131     }
 132
 133   private:
 134     const T a, b, z;
 135     int offset;
 136     hypergeometric_1F1_recurrence_a_and_b_coefficients operator=(const hypergeometric_1F1_recurrence_a_and_b_coefficients&);
 137   };
 138 #if 0
 139   //
 140   // These next few recurrence relations are archived for future refernece, some of them are novel, though all
 141   // are trivially derived from the existing well known relations:
 142   //
 143   // Recurrence relation for double-stepping on both a and b:
 144   // - b(b-1)(b-2) / (2-b+z) M(a-2,b-2,z) + [b(a-1)z / (2-b+z) + b(1-b+z) + abz(b+1) /(b+1)(z-b)] M(a,b,z) - a(a+1)z^2 / (b+1)(z-b) M(a+2,b+2,z)
 145   //
 146   template <class T>
 147   struct hypergeometric_1F1_recurrence_2a_and_2b_coefficients
 148   {
 149      typedef boost::math::tuple<T, T, T> result_type;
 150
 151      hypergeometric_1F1_recurrence_2a_and_2b_coefficients(const T& a, const T& b, const T& z, int offset = 0) :
 152         a(a), b(b), z(z), offset(offset)
 153      {
 154      }
 155
 156      result_type operator()(boost::intmax_t i) const
 157      {
 158         i *= 2;
 159         const T ai = a + (offset + i);
 160         const T bi = b + (offset + i);
 161
 162         const T an = -bi * (b + (offset + i - 1)) * (b + (offset + i - 2)) / (-(b + (offset + i - 2)) + z);
 163         const T bn = bi * (a + (offset + i - 1)) * z / (z - (b + (offset + i - 2)))
 164            + bi * (z - (b + (offset + i - 1)))
 165            + ai * bi * z * (b + (offset + i + 1)) / ((b + (offset + i + 1)) * (z - bi));
 166         const T cn = -ai * (a + (offset + i + 1)) * z * z / ((b + (offset + i + 1)) * (z - bi));
 167
 168         return boost::math::make_tuple(an, bn, cn);
 169      }
 170
 171   private:
 172      const T a, b, z;
 173      int offset;
 174      hypergeometric_1F1_recurrence_2a_and_2b_coefficients operator=(const hypergeometric_1F1_recurrence_2a_and_2b_coefficients&);
 175   };
 176
 177   //
 178   // Recurrence relation for double-stepping on a:
 179   // -(b-a)(1 + b - a)/(2a-2-b+z)M(a-2,b,z)  + [(b-a)(a-1)/(2a-2-b+z) + (2a-b+z) + a(b-a-1)/(2a+2-b+z)]M(a,b,z)   -a(a+1)/(2a+2-b+z)M(a+2,b,z)
 180   //
 181   template <class T>
 182   struct hypergeometric_1F1_recurrence_2a_coefficients
 183   {
 184      typedef boost::math::tuple<T, T, T> result_type;
 185
 186      hypergeometric_1F1_recurrence_2a_coefficients(const T& a, const T& b, const T& z, int offset = 0) :
 187         a(a), b(b), z(z), offset(offset)
 188      {
 189      }
 190
 191      result_type operator()(boost::intmax_t i) const
 192      {
 193         i *= 2;
 194         const T ai = a + (offset + i);
 195         // -(b-a)(1 + b - a)/(2a-2-b+z)
 196         const T an = -(b - ai) * (b - (a + (offset + i - 1))) / (2 * (a + (offset + i - 1)) - b + z);
 197         const T bn = (b - ai) * (a + (offset + i - 1)) / (2 * (a + (offset + i - 1)) - b + z) + (2 * ai - b + z) + ai * (b - (a + (offset + i + 1))) / (2 * (a + (offset + i + 1)) - b + z);
 198         const T cn = -ai * (a + (offset + i + 1)) / (2 * (a + (offset + i + 1)) - b + z);
 199
 200         return boost::math::make_tuple(an, bn, cn);
 201      }
 202
 203   private:
 204      const T a, b, z;
 205      int offset;
 206      hypergeometric_1F1_recurrence_2a_coefficients operator=(const hypergeometric_1F1_recurrence_2a_coefficients&);
 207   };
 208
 209   //
 210   // Recurrence relation for double-stepping on b:
 211   // b(b-1)^2(b-2)/((1-b)(2-b-z)) M(a,b-2,z)  + [zb(b-1)(b-1-a)/((1-b)(2-b-z)) + b(1-b-z) + z(b-a)(b+1)b/((b+1)(b+z)) ] M(a,b,z) + z^2(b-a)(b+1-a)/((b+1)(b+z)) M(a,b+2,z)
 212   //
 213   template <class T>
 214   struct hypergeometric_1F1_recurrence_2b_coefficients
 215   {
 216      typedef boost::math::tuple<T, T, T> result_type;
 217
 218      hypergeometric_1F1_recurrence_2b_coefficients(const T& a, const T& b, const T& z, int offset = 0) :
 219         a(a), b(b), z(z), offset(offset)
 220      {
 221      }
 222
 223      result_type operator()(boost::intmax_t i) const
 224      {
 225         i *= 2;
 226         const T bi = b + (offset + i);
 227         const T bi_m1 = b + (offset + i - 1);
 228         const T bi_p1 = b + (offset + i + 1);
 229         const T bi_m2 = b + (offset + i - 2);
 230
 231         const T an = bi * (bi_m1) * (bi_m1) * (bi_m2) / (-bi_m1 * (-bi_m2 - z));
 232         const T bn = z * bi * bi_m1 * (bi_m1 - a) / (-bi_m1 * (-bi_m2 - z)) + bi * (-bi_m1 - z) + z * (bi - a) * bi_p1 * bi / (bi_p1 * (bi + z));
 233         const T cn = z * z * (bi - a) * (bi_p1 - a) / (bi_p1 * (bi + z));
 234
 235         return boost::math::make_tuple(an, bn, cn);
 236      }
 237
 238   private:
 239      const T a, b, z;
 240      int offset;
 241      hypergeometric_1F1_recurrence_2b_coefficients operator=(const hypergeometric_1F1_recurrence_2b_coefficients&);
 242   };
 243
 244   //
 245   // Recurrence relation for a+ b-:
 246   // -z(b-a)(a-1-b)/(b(a-1+z)) M(a-1,b+1,z) + [(b-a)(a-1)b/(b(a-1+z)) + (2a-b+z) + a(b-a-1)/(a+z)] M(a,b,z) + a(1-b)/(a+z) M(a+1,b-1,z)
 247   //
 248   // This is potentially the most useful of these novel recurrences.
 249   //              -                                      -                  +        -                           +
 250   template <class T>
 251   struct hypergeometric_1F1_recurrence_a_plus_b_minus_coefficients
 252   {
 253      typedef boost::math::tuple<T, T, T> result_type;
 254
 255      hypergeometric_1F1_recurrence_a_plus_b_minus_coefficients(const T& a, const T& b, const T& z, int offset = 0) :
 256         a(a), b(b), z(z), offset(offset)
 257      {
 258      }
 259
 260      result_type operator()(boost::intmax_t i) const
 261      {
 262         const T ai = a + (offset + i);
 263         const T bi = b - (offset + i);
 264
 265         const T an = -z * (bi - ai) * (ai - 1 - bi) / (bi * (ai - 1 + z));
 266         const T bn = z * ((-1 / (ai + z) - 1 / (ai + z - 1)) * (bi + z - 1) + 3) + bi - 1;
 267         const T cn = ai * (1 - bi) / (ai + z);
 268
 269         return boost::math::make_tuple(an, bn, cn);
 270      }
 271
 272   private:
 273      const T a, b, z;
 274      int offset;
 275      hypergeometric_1F1_recurrence_a_plus_b_minus_coefficients operator=(const hypergeometric_1F1_recurrence_a_plus_b_minus_coefficients&);
 276   };
 277 #endif
 278
 279   template <class T, class Policy>
 280   inline T hypergeometric_1F1_backward_recurrence_for_negative_a(const T& a, const T& b, const T& z, const Policy& pol, const char* function, int& log_scaling)
 281   {
 282     BOOST_MATH_STD_USING // modf, frexp, fabs, pow
 283
 284     boost::intmax_t integer_part = 0;
 285     T ak = modf(a, &integer_part);
 286     //
 287     // We need ak-1 positive to avoid infinite recursion below:
 288     //
 289     if (0 != ak)
 290     {
 291        ak += 2;
 292        integer_part -= 2;
 293     }
 294
 295     if (-integer_part > static_cast<boost::intmax_t>(policies::get_max_series_iterations<Policy>()))
 296        return policies::raise_evaluation_error<T>(function, "1F1 arguments sit in a range with a so negative that we have no evaluation method, got a = %1%", std::numeric_limits<T>::quiet_NaN(), pol);
 297
 298     T first, second;
 299     if(ak == 0)
 300     {
 301        first = 1;
 302        ak -= 1;
 303        second = 1 - z / b;
 304     }
 305     else
 306     {
 307        int scaling1(0), scaling2(0);
 308        first = detail::hypergeometric_1F1_imp(ak, b, z, pol, scaling1);
 309        ak -= 1;
 310        second = detail::hypergeometric_1F1_imp(ak, b, z, pol, scaling2);
 311        if (scaling1 != scaling2)
 312        {
 313           second *= exp(T(scaling2 - scaling1));
 314        }
 315        log_scaling += scaling1;
 316     }
 317     ++integer_part;
 318
 319     detail::hypergeometric_1F1_recurrence_a_coefficients<T> s(ak, b, z);
 320
 321     return tools::apply_recurrence_relation_backward(s,
 322                                                      static_cast<unsigned int>(std::abs(integer_part)),
 323                                                      first,
 324                                                      second, &log_scaling);
 325   }
 326
 327
 328   template <class T, class Policy>
 329   T hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(const T& a, const T& b, const T& z, const Policy& pol, const char*, int& log_scaling)
 330   {
 331      using std::swap;
 332      BOOST_MATH_STD_USING // modf, frexp, fabs, pow
 333      //
 334      // We compute
 335      //
 336      // M[a + a_shift, b + b_shift; z]
 337      //
 338      // and recurse backwards on a and b down to
 339      //
 340      // M[a, b, z]
 341      //
 342      // With a + a_shift > 1 and b + b_shift > z
 343      //
 344      // There are 3 distinct regions to ensure stability during the recursions:
 345      //
 346      // a > 0         :  stable for backwards on a
 347      // a < 0, b > 0  :  stable for backwards on a and b
 348      // a < 0, b < 0  :  stable for backwards on b (as long as |b| is small).
 349      //
 350      // We could simplify things by ignoring the middle region, but it's more efficient
 351      // to recurse on a and b together when we can.
 352      //
 353
 354      BOOST_ASSERT(a < -1); // Not tested nor taken for -1 < a < 0
 355
 356      int b_shift = itrunc(z - b) + 2;
 357
 358      int a_shift = itrunc(-a);
 359      if (a + a_shift != 0)
 360      {
 361         a_shift += 2;
 362      }
 363      //
 364      // If the shifts are so large that we would throw an evaluation_error, try the series instead,
 365      // even though this will almost certainly throw as well:
 366      //
 367      if (b_shift > static_cast<boost::intmax_t>(boost::math::policies::get_max_series_iterations<Policy>()))
 368         return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
 369
 370      if (a_shift > static_cast<boost::intmax_t>(boost::math::policies::get_max_series_iterations<Policy>()))
 371         return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
 372
 373      int a_b_shift = b < 0 ? itrunc(b + b_shift) : b_shift;   // The max we can shift on a and b together
 374      int leading_a_shift = (std::min)(3, a_shift);        // Just enough to make a negative
 375      if (a_b_shift > a_shift - 3)
 376      {
 377         a_b_shift = a_shift < 3 ? 0 : a_shift - 3;
 378      }
 379      else
 380      {
 381         // Need to ensure that leading_a_shift is large enough that a will reach it's target
 382         // after the first 2 phases (-,0) and (-,-) are over:
 383         leading_a_shift = a_shift - a_b_shift;
 384      }
 385      int trailing_b_shift = b_shift - a_b_shift;
 386      if (a_b_shift < 5)
 387      {
 388         // Might as well do things in two steps rather than 3:
 389         if (a_b_shift > 0)
 390         {
 391            leading_a_shift += a_b_shift;
 392            trailing_b_shift += a_b_shift;
 393         }
 394         a_b_shift = 0;
 395         --leading_a_shift;
 396      }
 397
 398      BOOST_ASSERT(leading_a_shift > 1);
 399      BOOST_ASSERT(a_b_shift + leading_a_shift + (a_b_shift == 0 ? 1 : 0) == a_shift);
 400      BOOST_ASSERT(a_b_shift + trailing_b_shift == b_shift);
 401
 402      if ((trailing_b_shift == 0) && (fabs(b) < 0.5) && a_b_shift)
 403      {
 404         // Better to have the final recursion on b alone, otherwise we lose precision when b is very small:
 405         int diff = (std::min)(a_b_shift, 3);
 406         a_b_shift -= diff;
 407         leading_a_shift += diff;
 408         trailing_b_shift += diff;
 409      }
 410
 411      T first, second;
 412      int scale1(0), scale2(0);
 413      first = boost::math::detail::hypergeometric_1F1_imp(T(a + a_shift), T(b + b_shift), z, pol, scale1);
 414      //
 415      // It would be good to compute "second" from first and the ratio - unfortunately we are right on the cusp
 416      // recursion on a switching from stable backwards to stable forwards behaviour and so this is not possible here.
 417      //
 418      second = boost::math::detail::hypergeometric_1F1_imp(T(a + a_shift - 1), T(b + b_shift), z, pol, scale2);
 419      if (scale1 != scale2)
 420         second *= exp(T(scale2 - scale1));
 421      log_scaling += scale1;
 422
 423      //
 424      // Now we have [a + a_shift, b + b_shift, z] and [a + a_shift - 1, b + b_shift, z]
 425      // and want to recurse until [a + a_shift - leading_a_shift, b + b_shift, z] and [a + a_shift - leadng_a_shift - 1, b + b_shift, z]
 426      // which is leading_a_shift -1 steps.
 427      //
 428      second = boost::math::tools::apply_recurrence_relation_backward(
 429         hypergeometric_1F1_recurrence_a_coefficients<T>(a + a_shift - 1, b + b_shift, z),
 430         leading_a_shift, first, second, &log_scaling, &first);
 431
 432      if (a_b_shift)
 433      {
 434         //
 435         // Now we need to switch to an a+b shift so that we have:
 436         // [a + a_shift - leading_a_shift, b + b_shift, z] and [a + a_shift - leadng_a_shift - 1, b + b_shift - 1, z]
 437         // A&S 13.4.3 gives us what we need:
 438         //
 439         {
 440            // local a's and b's:
 441            T la = a + a_shift - leading_a_shift - 1;
 442            T lb = b + b_shift;
 443            second = ((1 + la - lb) * second - la * first) / (1 - lb);
 444         }
 445         //
 446         // Now apply a_b_shift - 1 recursions to get down to
 447         // [a + 1, b + trailing_b_shift + 1, z] and [a, b + trailing_b_shift, z]
 448         //
 449         second = boost::math::tools::apply_recurrence_relation_backward(
 450            hypergeometric_1F1_recurrence_a_and_b_coefficients<T>(a, b + b_shift - a_b_shift, z, a_b_shift - 1),
 451            a_b_shift - 1, first, second, &log_scaling, &first);
 452         //
 453         // Now we need to switch to a b shift, a different application of A&S 13.4.3
 454         // will get us there, we leave "second" where it is, and move "first" sideways:
 455         //
 456         {
 457            T lb = b + trailing_b_shift + 1;
 458            first = (second * (lb - 1) - a * first) / -(1 + a - lb);
 459         }
 460      }
 461      else
 462      {
 463         //
 464         // We have M[a+1, b+b_shift, z] and M[a, b+b_shift, z] and need M[a, b+b_shift-1, z] for
 465         // recursion on b: A&S 13.4.3 gives us what we need.
 466         //
 467         T third = -(second * (1 + a - b - b_shift) - first * a) / (b + b_shift - 1);
 468         swap(first, second);
 469         swap(second, third);
 470         --trailing_b_shift;
 471      }
 472      //
 473      // Finish off by applying trailing_b_shift recursions:
 474      //
 475      if (trailing_b_shift)
 476      {
 477         second = boost::math::tools::apply_recurrence_relation_backward(
 478            hypergeometric_1F1_recurrence_small_b_coefficients<T>(a, b, z, trailing_b_shift),
 479            trailing_b_shift, first, second, &log_scaling);
 480      }
 481      return second;
 482   }
 483
 484
 485
 486   } } } // namespaces
 487
 488 #endif // BOOST_HYPERGEOMETRIC_1F1_RECURRENCE_HPP_