boost/math/special_functions/detail/hypergeometric_1F1_large_abz.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_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_
   9 #define BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_
  10
  11 #include <boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp>
  12 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
  13 #include <boost/math/special_functions/gamma.hpp>
  14
  15
  16   namespace boost { namespace math { namespace detail {
  17
  18      template <class T>
  19      inline bool is_negative_integer(const T& x)
  20      {
  21         using std::floor;
  22         return (x <= 0) && (floor(x) == x);
  23      }
  24
  25
  26      template <class T, class Policy>
  27      struct hypergeometric_1F1_igamma_series
  28      {
  29         enum{ cache_size = 64 };
  30
  31         typedef T result_type;
  32         hypergeometric_1F1_igamma_series(const T& alpha, const T& delta, const T& x, const Policy& pol)
  33            : delta_poch(-delta), alpha_poch(alpha), x(x), k(0), cache_offset(0), pol(pol)
  34         {
  35            BOOST_MATH_STD_USING
  36            T log_term = log(x) * -alpha;
  37            log_scaling = itrunc(log_term - 3 - boost::math::tools::log_min_value<T>() / 50);
  38            term = exp(log_term - log_scaling);
  39            refill_cache();
  40         }
  41         T operator()()
  42         {
  43            if (k - cache_offset >= cache_size)
  44            {
  45               cache_offset += cache_size;
  46               refill_cache();
  47            }
  48            T result = term * gamma_cache[k - cache_offset];
  49            term *= delta_poch * alpha_poch / (++k * x);
  50            delta_poch += 1;
  51            alpha_poch += 1;
  52            return result;
  53         }
  54         void refill_cache()
  55         {
  56            typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
  57
  58            gamma_cache[cache_size - 1] = boost::math::gamma_p(alpha_poch + ((int)cache_size - 1), x, pol);
  59            for (int i = cache_size - 1; i > 0; --i)
  60            {
  61               gamma_cache[i - 1] = gamma_cache[i] >= 1 ? T(1) : T(gamma_cache[i] + regularised_gamma_prefix(T(alpha_poch + (i - 1)), x, pol, lanczos_type()) / (alpha_poch + (i - 1)));
  62            }
  63         }
  64         T delta_poch, alpha_poch, x, term;
  65         T gamma_cache[cache_size];
  66         int k;
  67         int log_scaling;
  68         int cache_offset;
  69         Policy pol;
  70      };
  71
  72      template <class T, class Policy>
  73      T hypergeometric_1F1_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, int& log_scaling)
  74      {
  75         BOOST_MATH_STD_USING
  76         if (b_minus_a == 0)
  77         {
  78            // special case: M(a,a,z) == exp(z)
  79            int scale = itrunc(x, pol);
  80            log_scaling += scale;
  81            return exp(x - scale);
  82         }
  83         hypergeometric_1F1_igamma_series<T, Policy> s(b_minus_a, a - 1, x, pol);
  84         log_scaling += s.log_scaling;
  85         boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
  86         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
  87         boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol);
  88         T log_prefix = x + boost::math::lgamma(b, pol) - boost::math::lgamma(a, pol);
  89         int scale = itrunc(log_prefix);
  90         log_scaling += scale;
  91         return result * exp(log_prefix - scale);
  92      }
  93
  94      template <class T, class Policy>
  95      T hypergeometric_1F1_shift_on_a(T h, const T& a_local, const T& b_local, const T& x, int a_shift, const Policy& pol, int& log_scaling)
  96      {
  97         BOOST_MATH_STD_USING
  98         T a = a_local + a_shift;
  99         if (a_shift == 0)
 100            return h;
 101         else if (a_shift > 0)
 102         {
 103            //
 104            // Forward recursion on a is stable as long as 2a-b+z > 0.
 105            // If 2a-b+z < 0 then backwards recursion is stable even though
 106            // the function may be strictly increasing with a.  Potentially
 107            // we may need to split the recurrnce in 2 sections - one using
 108            // forward recursion, and one backwards.
 109            //
 110            // We will get the next seed value from the ratio
 111            // on b as that's stable and quick to compute.
 112            //
 113
 114            T crossover_a = (b_local - x) / 2;
 115            int crossover_shift = itrunc(crossover_a - a_local);
 116
 117            if (crossover_shift > 1)
 118            {
 119               //
 120               // Forwards recursion will start off unstable, but may switch to the stable direction later.
 121               // Start in the middle and go in both directions:
 122               //
 123               if (crossover_shift > a_shift)
 124                  crossover_shift = a_shift;
 125               crossover_a = a_local + crossover_shift;
 126               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(crossover_a, b_local, x);
 127               boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 128               T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 129               boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
 130               //
 131               // Convert to a ratio:
 132               //         (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0
 133               //
 134               //  hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio
 135               //
 136               T first = 1;
 137               T second = ((1 + crossover_a - b_local) / crossover_a) + ((b_local - 1) / crossover_a) / b_ratio;
 138               //
 139               // Recurse down to a_local, compare values and re-nomralise first and second:
 140               //
 141               boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(crossover_a, b_local, x);
 142               int backwards_scale = 0;
 143               T comparitor = boost::math::tools::apply_recurrence_relation_backward(a_coef, crossover_shift, second, first, &backwards_scale);
 144               log_scaling -= backwards_scale;
 145               if ((h < 1) && (tools::max_value<T>() * h > comparitor))
 146               {
 147                  // Need to rescale!
 148                  int scale = itrunc(log(h), pol) + 1;
 149                  h *= exp(T(-scale));
 150                  log_scaling += scale;
 151               }
 152               comparitor /= h;
 153               first /= comparitor;
 154               second /= comparitor;
 155               //
 156               // Now we can recurse forwards for the rest of the range:
 157               //
 158               if (crossover_shift < a_shift)
 159               {
 160                  boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef_2(crossover_a + 1, b_local, x);
 161                  h = boost::math::tools::apply_recurrence_relation_forward(a_coef_2, a_shift - crossover_shift - 1, first, second, &log_scaling);
 162               }
 163               else
 164                  h = first;
 165            }
 166            else
 167            {
 168               //
 169               // Regular case where forwards iteration is stable right from the start:
 170               //
 171               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a_local, b_local, x);
 172               boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 173               T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 174               boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
 175               //
 176               // Convert to a ratio:
 177               //         (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0
 178               //
 179               //  hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio
 180               //
 181               T second = ((1 + a_local - b_local) / a_local) * h + ((b_local - 1) / a_local) * h / b_ratio;
 182               boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a_local + 1, b_local, x);
 183               h = boost::math::tools::apply_recurrence_relation_forward(a_coef, --a_shift, h, second, &log_scaling);
 184            }
 185         }
 186         else
 187         {
 188            //
 189            // We've calculated h for a larger value of a than we want, and need to recurse down.
 190            // However, only forward iteration is stable, so calculate the ratio, compare values,
 191            // and normalise.  Note that we calculate the ratio on b and convert to a since the
 192            // direction is the minimal solution for N->+INF.
 193            //
 194            // IMPORTANT: this is only currently called for a > b and therefore forwards iteration
 195            // is the only stable direction as we will only iterate down until a ~ b, but we
 196            // will check this with an assert:
 197            //
 198            BOOST_ASSERT(2 * a - b_local + x > 0);
 199            boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x);
 200            boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 201            T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 202            boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
 203            //
 204            // Convert to a ratio:
 205            //         (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0
 206            //
 207            //  hence: M(a+1,b,z) = (1+a-b) / a M(a,b,z) + (b-1) / a M(a,b,z)/ (M(a,b,z)/M(a,b-1,z))
 208            //
 209            T first = 1;  // arbitrary value;
 210            T second = ((1 + a - b_local) / a) + ((b_local - 1) / a) * (1 / b_ratio);
 211
 212            if (a_shift == -1)
 213               h = h / second;
 214            else
 215            {
 216               boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a + 1, b_local, x);
 217               T comparitor = boost::math::tools::apply_recurrence_relation_forward(a_coef, -(a_shift + 1), first, second);
 218               if (boost::math::tools::min_value<T>() * comparitor > h)
 219               {
 220                  // Ooops, need to rescale h:
 221                  int rescale = itrunc(log(fabs(h)));
 222                  T scale = exp(T(-rescale));
 223                  h *= scale;
 224                  log_scaling += rescale;
 225               }
 226               h = h / comparitor;
 227            }
 228         }
 229         return h;
 230      }
 231
 232      template <class T, class Policy>
 233      T hypergeometric_1F1_shift_on_b(T h, const T& a, const T& b_local, const T& x, int b_shift, const Policy& pol, int& log_scaling)
 234      {
 235         BOOST_MATH_STD_USING
 236
 237         T b = b_local + b_shift;
 238         if (b_shift == 0)
 239            return h;
 240         else if (b_shift > 0)
 241         {
 242            //
 243            // We get here for b_shift > 0 when b > z.  We can't use forward recursion on b - it's unstable,
 244            // so grab the ratio and work backwards to b - b_shift and normalise.
 245            //
 246            boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b, x);
 247            boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 248
 249            T first = 1;  // arbitrary value;
 250            T second = 1 / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 251            boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
 252            if (b_shift == 1)
 253               h = h / second;
 254            else
 255            {
 256               //
 257               // Reset coefficients and recurse:
 258               //
 259               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b - 1, x);
 260               int local_scale = 0;
 261               T comparitor = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, --b_shift, first, second, &local_scale);
 262               log_scaling -= local_scale;
 263               if (boost::math::tools::min_value<T>() * comparitor > h)
 264               {
 265                  // Ooops, need to rescale h:
 266                  int rescale = itrunc(log(fabs(h)));
 267                  T scale = exp(T(-rescale));
 268                  h *= scale;
 269                  log_scaling += rescale;
 270               }
 271               h = h / comparitor;
 272            }
 273         }
 274         else
 275         {
 276            T second;
 277            if (a == b_local)
 278            {
 279                // recurrence is trivial for a == b and method of ratios fails as the c-term goes to zero:
 280               second = -b_local * (1 - b_local - x) * h / (b_local * (b_local - 1));
 281            }
 282            else
 283            {
 284               BOOST_ASSERT(!is_negative_integer(b - a));
 285               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x);
 286               boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 287               second = h / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 288               boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
 289            }
 290            if (b_shift == -1)
 291               h = second;
 292            else
 293            {
 294               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b_local - 1, x);
 295               h = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, -(++b_shift), h, second, &log_scaling);
 296            }
 297         }
 298         return h;
 299      }
 300
 301
 302      template <class T, class Policy>
 303      T hypergeometric_1F1_large_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, int& log_scaling)
 304      {
 305         BOOST_MATH_STD_USING
 306         //
 307         // We need a < b < z in order to ensure there's at least a chance of convergence,
 308         // we can use recurrence relations to shift forwards on a+b or just a to achieve this,
 309         // for decent accuracy, try to keep 2b - 1 < a < 2b < z
 310         //
 311         int b_shift = b * 2 < x ? 0 : itrunc(b - x / 2);
 312         int a_shift = a > b - b_shift ? -itrunc(b - b_shift - a - 1) : -itrunc(b - b_shift - a);
 313
 314         if (a_shift < 0)
 315         {
 316            // Might as well do all the shifting on b as scale a downwards:
 317            b_shift -= a_shift;
 318            a_shift = 0;
 319         }
 320
 321         T a_local = a - a_shift;
 322         T b_local = b - b_shift;
 323         T b_minus_a_local = (a_shift == 0) && (b_shift == 0) ? b_minus_a : b_local - a_local;
 324         int local_scaling = 0;
 325         T h = hypergeometric_1F1_igamma(a_local, b_local, x, b_minus_a_local, pol, local_scaling);
 326         log_scaling += local_scaling;
 327
 328         //
 329         // Apply shifts on a and b as required:
 330         //
 331         h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, x, a_shift, pol, log_scaling);
 332         h = hypergeometric_1F1_shift_on_b(h, a, b_local, x, b_shift, pol, log_scaling);
 333
 334         return h;
 335      }
 336
 337      template <class T, class Policy>
 338      T hypergeometric_1F1_large_series(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
 339      {
 340         BOOST_MATH_STD_USING
 341         //
 342         // We make a small, and b > z:
 343         //
 344         int a_shift(0), b_shift(0);
 345         if (a * z > b)
 346         {
 347            a_shift = itrunc(a) - 5;
 348            b_shift = b < z ? itrunc(b - z - 1) : 0;
 349         }
 350         //
 351         // If a_shift is trivially small, there's really not much point in losing
 352         // accuracy to save a couple of iterations:
 353         //
 354         if (a_shift < 5)
 355            a_shift = 0;
 356         T a_local = a - a_shift;
 357         T b_local = b - b_shift;
 358         T h = boost::math::detail::hypergeometric_1F1_generic_series(a_local, b_local, z, pol, log_scaling, "hypergeometric_1F1_large_series<%1%>(a,b,z)");
 359         //
 360         // Apply shifts on a and b as required:
 361         //
 362         if (a_shift && (a_local == 0))
 363         {
 364            //
 365            // Shifting on a via method of ratios in hypergeometric_1F1_shift_on_a fails when
 366            // a_local == 0.  However, the value of h calculated was trivial (unity), so
 367            // calculate a second 1F1 for a == 1 and recurse as normal:
 368            //
 369            int scale = 0;
 370            T h2 = boost::math::detail::hypergeometric_1F1_generic_series(T(a_local + 1), b_local, z, pol, scale, "hypergeometric_1F1_large_series<%1%>(a,b,z)");
 371            if (scale != log_scaling)
 372            {
 373               h2 *= exp(T(scale - log_scaling));
 374            }
 375            boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> coef(a_local + 1, b_local, z);
 376            h = boost::math::tools::apply_recurrence_relation_forward(coef, a_shift - 1, h, h2, &log_scaling);
 377            h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling);
 378         }
 379         else
 380         {
 381            h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, z, a_shift, pol, log_scaling);
 382            h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling);
 383         }
 384         return h;
 385      }
 386
 387      template <class T, class Policy>
 388      T hypergeometric_1F1_large_13_3_6_series(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
 389      {
 390         BOOST_MATH_STD_USING
 391         //
 392         // A&S 13.3.6 is good only when a ~ b, but isn't too fussy on the size of z.
 393         // So shift b to match a (b shifting seems to be more stable via method of ratios).
 394         //
 395         int b_shift = itrunc(b - a);
 396         T b_local = b - b_shift;
 397         T h = boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b_local, z, T(b_local - a), pol, log_scaling);
 398         return hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling);
 399      }
 400
 401      template <class T, class Policy>
 402      T hypergeometric_1F1_large_abz(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
 403      {
 404         BOOST_MATH_STD_USING
 405         //
 406         // This is the selection logic to pick the "best" method.
 407         // We have a,b,z >> 0 and need to comute the approximate cost of each method
 408         // and then select whichever wins out.
 409         //
 410         enum method
 411         {
 412            method_series = 0,
 413            method_shifted_series,
 414            method_gamma,
 415            method_bessel
 416         };
 417         //
 418         // Cost of direct series, is the approx number of terms required until we hit convergence:
 419         //
 420         T current_cost = (sqrt(16 * z * (3 * a + z) + 9 * b * b - 24 * b * z) - 3 * b + 4 * z) / 6;
 421         method current_method = method_series;
 422         //
 423         // Cost of shifted series, is the number of recurrences required to move to a zone where
 424         // the series is convergent right from the start.
 425         // Note that recurrence relations fail for very small b, and too many recurrences on a
 426         // will completely destroy all our digits.
 427         // Also note that the method fails when b-a is a negative integer unless b is already
 428         // larger than z and thus does not need shifting.
 429         //
 430         T cost = a + ((b < z) ? T(z - b) : T(0));
 431         if((b > 1) && (cost < current_cost) && ((b > z) || !is_negative_integer(b-a)))
 432         {
 433            current_method = method_shifted_series;
 434            current_cost = cost;
 435         }
 436         //
 437         // Cost for gamma function method is the number of recurrences required to move it
 438         // into a convergent zone, note that recurrence relations fail for very small b.
 439         // Also add on a fudge factor to account for the fact that this method is both
 440         // more expensive to compute (requires gamma functions), and less accurate than the
 441         // methods above:
 442         //
 443         T b_shift = fabs(b * 2 < z ? T(0) : T(b - z / 2));
 444         T a_shift = fabs(a > b - b_shift ? T(-(b - b_shift - a - 1)) : T(-(b - b_shift - a)));
 445         cost = 1000 + b_shift + a_shift;
 446         if((b > 1) && (cost <= current_cost))
 447         {
 448            current_method = method_gamma;
 449            current_cost = cost;
 450         }
 451         //
 452         // Cost for bessel approximation is the number of recurrences required to make a ~ b,
 453         // Note that recurrence relations fail for very small b.  We also have issue with large
 454         // z: either overflow/numeric instability or else the series goes divergent.  We seem to be
 455         // OK for z smaller than log_max_value<Quad> though, maybe we can stretch a little further
 456         // but that's not clear...
 457         // Also need to add on a fudge factor to the cost to account for the fact that we need
 458         // to calculate the Bessel functions, this is not quite as high as the gamma function
 459         // method above as this is generally more accurate and so prefered if the methods are close:
 460         //
 461         cost = 50 + fabs(b - a);
 462         if((b > 1) && (cost <= current_cost) && (z < tools::log_max_value<T>()) && (z < 11356) && (b - a != 0.5f))
 463         {
 464            current_method = method_bessel;
 465            current_cost = cost;
 466         }
 467
 468         switch (current_method)
 469         {
 470         case method_series:
 471            return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, "hypergeometric_1f1_large_abz<%1%>(a,b,z)");
 472         case method_shifted_series:
 473            return detail::hypergeometric_1F1_large_series(a, b, z, pol, log_scaling);
 474         case method_gamma:
 475            return detail::hypergeometric_1F1_large_igamma(a, b, z, T(b - a), pol, log_scaling);
 476         case method_bessel:
 477            return detail::hypergeometric_1F1_large_13_3_6_series(a, b, z, pol, log_scaling);
 478         }
 479         return 0; // We don't get here!
 480      }
 481
 482   } } } // namespaces
 483
 484 #endif // BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_