boost/math/special_functions/detail/hypergeometric_1F1_bessel.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_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP
  12 #define BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP
  13
  14 #include <boost/math/tools/series.hpp>
  15 #include <boost/math/special_functions/bessel.hpp>
  16 #include <boost/math/special_functions/laguerre.hpp>
  17 #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
  18 #include <boost/math/special_functions/bessel_iterators.hpp>
  19
  20
  21   namespace boost { namespace math { namespace detail {
  22
  23      template <class T, class Policy>
  24      T hypergeometric_1F1_divergent_fallback(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling);
  25
  26      template <class T>
  27      bool hypergeometric_1F1_is_tricomi_viable_positive_b(const T& a, const T& b, const T& z)
  28      {
  29         BOOST_MATH_STD_USING
  30            if ((z < b) && (a > -50))
  31               return false;  // might as well fall through to recursion
  32         if (b <= 100)
  33            return true;
  34         // Even though we're in a reasonable domain for Tricomi's approximation,
  35         // the arguments to the Bessel functions may be so large that we can't
  36         // actually evaluate them:
  37         T x = sqrt(fabs(2 * z * b - 4 * a * z));
  38         T v = b - 1;
  39         return log(boost::math::constants::e<T>() * x / (2 * v)) * v > tools::log_min_value<T>();
  40      }
  41
  42      //
  43      // Returns an arbitrarily small value compared to "target" for use as a seed
  44      // value for Bessel recurrences.  Note that we'd better not make it too small
  45      // or underflow may occur resulting in either one of the terms in the
  46      // recurrence being zero, or else the result being zero.  Using 1/epsilon
  47      // as a safety factor ensures that if we do underflow to zero, all of the digits
  48      // will have been cancelled out anyway:
  49      //
  50      template <class T>
  51      T arbitrary_small_value(const T& target)
  52      {
  53         using std::fabs;
  54         return (tools::min_value<T>() / tools::epsilon<T>()) * (fabs(target) > 1 ? target : 1);
  55      }
  56
  57
  58      template <class T, class Policy>
  59      struct hypergeometric_1F1_AS_13_3_7_tricomi_series
  60      {
  61         typedef T result_type;
  62
  63         enum { cache_size = 64 };
  64
  65         hypergeometric_1F1_AS_13_3_7_tricomi_series(const T& a, const T& b, const T& z, const Policy& pol_)
  66            : A_minus_2(1), A_minus_1(0), A(b / 2), mult(z / 2), term(1), b_minus_1_plus_n(b - 1),
  67             bessel_arg((b / 2 - a) * z),
  68            two_a_minus_b(2 * a - b), pol(pol_), n(2)
  69         {
  70            BOOST_MATH_STD_USING
  71            term /= pow(fabs(bessel_arg), b_minus_1_plus_n / 2);
  72            mult /= sqrt(fabs(bessel_arg));
  73            bessel_cache[cache_size - 1] = bessel_arg > 0 ? boost::math::cyl_bessel_j(b_minus_1_plus_n - 1, 2 * sqrt(bessel_arg), pol) : boost::math::cyl_bessel_i(b_minus_1_plus_n - 1, 2 * sqrt(-bessel_arg), pol);
  74            if (fabs(bessel_cache[cache_size - 1]) < tools::min_value<T>() / tools::epsilon<T>())
  75            {
  76               // We get very limited precision due to rapid denormalisation/underflow of the Bessel values, raise an exception and try something else:
  77               policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Underflow in Bessel functions", bessel_cache[cache_size - 1], pol);
  78            }
  79            if ((term * bessel_cache[cache_size - 1] < tools::min_value<T>() / (tools::epsilon<T>() * tools::epsilon<T>())) || !(boost::math::isfinite)(term) || (!std::numeric_limits<T>::has_infinity && (fabs(term) > tools::max_value<T>())))
  80            {
  81               term = -log(fabs(bessel_arg)) * b_minus_1_plus_n / 2;
  82               log_scale = itrunc(term);
  83               term -= log_scale;
  84               term = exp(term);
  85            }
  86            else
  87               log_scale = 0;
  88 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
  89            if constexpr (std::numeric_limits<T>::has_infinity)
  90            {
  91               if (!(boost::math::isfinite)(bessel_cache[cache_size - 1]))
  92                  policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol);
  93            }
  94            else
  95               if ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))
  96                  policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol);
  97 #else
  98            if ((std::numeric_limits<T>::has_infinity && !(boost::math::isfinite)(bessel_cache[cache_size - 1]))
  99               || (!std::numeric_limits<T>::has_infinity && ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))))
 100               policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Expected finite Bessel function result but got %1%", bessel_cache[cache_size - 1], pol);
 101 #endif
 102            cache_offset = -cache_size;
 103            refill_cache();
 104         }
 105         T operator()()
 106         {
 107            //
 108            // We return the n-2 term, and do 2 terms at once as every other term can be
 109            // very small (or zero) when b == 2a:
 110            //
 111            BOOST_MATH_STD_USING
 112            if(n - 2 - cache_offset >= cache_size)
 113               refill_cache();
 114            T result = A_minus_2 * term * bessel_cache[n - 2 - cache_offset];
 115            term *= mult;
 116            ++n;
 117            T A_next = ((b_minus_1_plus_n + 2) * A_minus_1 + two_a_minus_b * A_minus_2) / n;
 118            b_minus_1_plus_n += 1;
 119            A_minus_2 = A_minus_1;
 120            A_minus_1 = A;
 121            A = A_next;
 122
 123            if (A_minus_2 != 0)
 124            {
 125               if (n - 2 - cache_offset >= cache_size)
 126                  refill_cache();
 127               result += A_minus_2 * term * bessel_cache[n - 2 - cache_offset];
 128            }
 129            term *= mult;
 130            ++n;
 131            A_next = ((b_minus_1_plus_n + 2) * A_minus_1 + two_a_minus_b * A_minus_2) / n;
 132            b_minus_1_plus_n += 1;
 133            A_minus_2 = A_minus_1;
 134            A_minus_1 = A;
 135            A = A_next;
 136
 137            return result;
 138         }
 139
 140         int scale()const
 141         {
 142            return log_scale;
 143         }
 144
 145      private:
 146         T A_minus_2, A_minus_1, A, mult, term, b_minus_1_plus_n, bessel_arg, two_a_minus_b;
 147         std::array<T, cache_size> bessel_cache;
 148         const Policy& pol;
 149         int n, log_scale, cache_offset;
 150
 151         hypergeometric_1F1_AS_13_3_7_tricomi_series operator=(const hypergeometric_1F1_AS_13_3_7_tricomi_series&);
 152
 153         void refill_cache()
 154         {
 155            BOOST_MATH_STD_USING
 156            //
 157            // We don't calculate a new bessel I/J value: instead start our iterator off
 158            // with an arbitrary small value, then when we get back to the last value in the previous cache
 159            // calculate the ratio and use it to renormalise all the new values.  This is more efficient, but
 160            // also avoids problems with J_v(x) or I_v(x) underflowing to zero.
 161            //
 162            cache_offset += cache_size;
 163            T last_value = bessel_cache.back();
 164            T ratio;
 165            if (bessel_arg > 0)
 166            {
 167               //
 168               // We will be calculating Bessel J.
 169               // We need a different recurrence strategy for positive and negative orders:
 170               //
 171               if (b_minus_1_plus_n > 0)
 172               {
 173                  bessel_j_backwards_iterator<T> i(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(bessel_arg), arbitrary_small_value(last_value));
 174
 175                  for (int j = cache_size - 1; j >= 0; --j, ++i)
 176                  {
 177                     bessel_cache[j] = *i;
 178                     //
 179                     // Depending on the value of bessel_arg, the values stored in the cache can grow so
 180                     // large as to overflow, if that looks likely then we need to rescale all the
 181                     // existing terms (most of which will then underflow to zero).  In this situation
 182                     // it's likely that our series will only need 1 or 2 terms of the series but we
 183                     // can't be sure of that:
 184                     //
 185                     if ((j < cache_size - 2) && (tools::max_value<T>() / fabs(64 * bessel_cache[j] / bessel_cache[j + 1]) < fabs(bessel_cache[j])))
 186                     {
 187                        T rescale = pow(fabs(bessel_cache[j] / bessel_cache[j + 1]), j + 1) * 2;
 188                        if (!((boost::math::isfinite)(rescale)))
 189                           rescale = tools::max_value<T>();
 190                        for (int k = j; k < cache_size; ++k)
 191                           bessel_cache[k] /= rescale;
 192                        bessel_j_backwards_iterator<T> ti(b_minus_1_plus_n + j, 2 * sqrt(bessel_arg), bessel_cache[j + 1], bessel_cache[j]);
 193                        i = ti;
 194                     }
 195                  }
 196                  ratio = last_value / *i;
 197               }
 198               else
 199               {
 200                  //
 201                  // Negative order is difficult: the J_v(x) recurrence relations are unstable
 202                  // *in both directions* for v < 0, except as v -> -INF when we have
 203                  // J_-v(x)  ~= -sin(pi.v)Y_v(x).
 204                  // For small v what we can do is compute every other Bessel function and
 205                  // then fill in the gaps using the recurrence relation.  This *is* stable
 206                  // provided that v is not so negative that the above approximation holds.
 207                  //
 208                  bessel_cache[1] = cyl_bessel_j(b_minus_1_plus_n + 1, 2 * sqrt(bessel_arg), pol);
 209                  bessel_cache[0] = (last_value + bessel_cache[1]) / (b_minus_1_plus_n / sqrt(bessel_arg));
 210                  int pos = 2;
 211                  while ((pos < cache_size - 1) && (b_minus_1_plus_n + pos < 0))
 212                  {
 213                     bessel_cache[pos + 1] = cyl_bessel_j(b_minus_1_plus_n + pos + 1, 2 * sqrt(bessel_arg), pol);
 214                     bessel_cache[pos] = (bessel_cache[pos-1] + bessel_cache[pos+1]) / ((b_minus_1_plus_n + pos) / sqrt(bessel_arg));
 215                     pos += 2;
 216                  }
 217                  if (pos < cache_size)
 218                  {
 219                     //
 220                     // We have crossed over into the region where backward recursion is the stable direction
 221                     // start from arbitrary value and recurse down to "pos" and normalise:
 222                     //
 223                     bessel_j_backwards_iterator<T> i2(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(bessel_arg), arbitrary_small_value(bessel_cache[pos-1]));
 224                     for (int loc = cache_size - 1; loc >= pos; --loc)
 225                        bessel_cache[loc] = *i2++;
 226                     ratio = bessel_cache[pos - 1] / *i2;
 227                     //
 228                     // Sanity check, if we normalised to an unusually small value then it was likely
 229                     // to be near a root and the calculated ratio is garbage, if so perform one
 230                     // more J_v(x) evaluation at position and normalise again:
 231                     //
 232                     if (fabs(bessel_cache[pos] * ratio / bessel_cache[pos - 1]) > 5)
 233                        ratio = cyl_bessel_j(b_minus_1_plus_n + pos, 2 * sqrt(bessel_arg), pol) / bessel_cache[pos];
 234                     while (pos < cache_size)
 235                        bessel_cache[pos++] *= ratio;
 236                  }
 237                  ratio = 1;
 238               }
 239            }
 240            else
 241            {
 242               //
 243               // Bessel I.
 244               // We need a different recurrence strategy for positive and negative orders:
 245               //
 246               if (b_minus_1_plus_n > 0)
 247               {
 248                  bessel_i_backwards_iterator<T> i(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(-bessel_arg), arbitrary_small_value(last_value));
 249
 250                  for (int j = cache_size - 1; j >= 0; --j, ++i)
 251                  {
 252                     bessel_cache[j] = *i;
 253                     //
 254                     // Depending on the value of bessel_arg, the values stored in the cache can grow so
 255                     // large as to overflow, if that looks likely then we need to rescale all the
 256                     // existing terms (most of which will then underflow to zero).  In this situation
 257                     // it's likely that our series will only need 1 or 2 terms of the series but we
 258                     // can't be sure of that:
 259                     //
 260                     if ((j < cache_size - 2) && (tools::max_value<T>() / fabs(64 * bessel_cache[j] / bessel_cache[j + 1]) < fabs(bessel_cache[j])))
 261                     {
 262                        T rescale = pow(fabs(bessel_cache[j] / bessel_cache[j + 1]), j + 1) * 2;
 263                        if (!((boost::math::isfinite)(rescale)))
 264                           rescale = tools::max_value<T>();
 265                        for (int k = j; k < cache_size; ++k)
 266                           bessel_cache[k] /= rescale;
 267                        i = bessel_i_backwards_iterator<T>(b_minus_1_plus_n + j, 2 * sqrt(-bessel_arg), bessel_cache[j + 1], bessel_cache[j]);
 268                     }
 269                  }
 270                  ratio = last_value / *i;
 271               }
 272               else
 273               {
 274                  //
 275                  // Forwards iteration is stable:
 276                  //
 277                  bessel_i_forwards_iterator<T> i(b_minus_1_plus_n, 2 * sqrt(-bessel_arg));
 278                  int pos = 0;
 279                  while ((pos < cache_size) && (b_minus_1_plus_n + pos < 0.5))
 280                  {
 281                     bessel_cache[pos++] = *i++;
 282                  }
 283                  if (pos < cache_size)
 284                  {
 285                     //
 286                     // We have crossed over into the region where backward recursion is the stable direction
 287                     // start from arbitrary value and recurse down to "pos" and normalise:
 288                     //
 289                     bessel_i_backwards_iterator<T> i2(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(-bessel_arg), arbitrary_small_value(last_value));
 290                     for (int loc = cache_size - 1; loc >= pos; --loc)
 291                        bessel_cache[loc] = *i2++;
 292                     ratio = bessel_cache[pos - 1] / *i2;
 293                     while (pos < cache_size)
 294                        bessel_cache[pos++] *= ratio;
 295                  }
 296                  ratio = 1;
 297               }
 298            }
 299            if(ratio != 1)
 300               for (auto j = bessel_cache.begin(); j != bessel_cache.end(); ++j)
 301                  *j *= ratio;
 302            //
 303            // Very occationally our normalisation fails because the normalisztion value
 304            // is sitting right on top of a root (or very close to it).  When that happens
 305            // best to calculate a fresh Bessel evaluation and normalise again.
 306            //
 307            if (fabs(bessel_cache[0] / last_value) > 5)
 308            {
 309               ratio = (bessel_arg < 0 ? cyl_bessel_i(b_minus_1_plus_n, 2 * sqrt(-bessel_arg), pol) : cyl_bessel_j(b_minus_1_plus_n, 2 * sqrt(bessel_arg), pol)) / bessel_cache[0];
 310               if (ratio != 1)
 311                  for (auto j = bessel_cache.begin(); j != bessel_cache.end(); ++j)
 312                     *j *= ratio;
 313            }
 314         }
 315      };
 316
 317      template <class T, class Policy>
 318      T hypergeometric_1F1_AS_13_3_7_tricomi(const T& a, const T& b, const T& z, const Policy& pol, int& log_scale)
 319      {
 320         BOOST_MATH_STD_USING
 321         //
 322         // Works for a < 0, b < 0, z > 0.
 323         //
 324         // For convergence we require A * term to be converging otherwise we get
 325         // a divergent alternating series.  It's actually really hard to analyse this
 326         // and the best purely heuristic policy we've found is
 327         // z < fabs((2 * a - b) / (sqrt(fabs(a)))) ; b > 0  or:
 328         // z < fabs((2 * a - b) / (sqrt(fabs(ab)))) ; b < 0
 329         //
 330         T prefix(0);
 331         int prefix_sgn(0);
 332         bool use_logs = false;
 333         int scale = 0;
 334         //
 335         // We can actually support the b == 2a case within here, but we haven't
 336         // as we appear never to get here in practice.  Which means this get out
 337         // clause is a bit of defensive programming....
 338         //
 339         if(b == 2 * a)
 340            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
 341
 342         try
 343         {
 344            prefix = boost::math::tgamma(b, pol);
 345            prefix *= exp(z / 2);
 346         }
 347         catch (const std::runtime_error&)
 348         {
 349            use_logs = true;
 350         }
 351         if (use_logs || (prefix == 0) || !(boost::math::isfinite)(prefix) || (!std::numeric_limits<T>::has_infinity && (fabs(prefix) >= tools::max_value<T>())))
 352         {
 353            use_logs = true;
 354            prefix = boost::math::lgamma(b, &prefix_sgn, pol) + z / 2;
 355            scale = itrunc(prefix);
 356            log_scale += scale;
 357            prefix -= scale;
 358         }
 359         T result(0);
 360         boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 361         bool retry = false;
 362         int series_scale = 0;
 363         try
 364         {
 365            hypergeometric_1F1_AS_13_3_7_tricomi_series<T, Policy> s(a, b, z, pol);
 366            series_scale = s.scale();
 367            log_scale += s.scale();
 368            try
 369            {
 370               T norm = 0;
 371               result = 0;
 372               if((a < 0) && (b < 0))
 373                  result = boost::math::tools::checked_sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, result, norm);
 374               else
 375                  result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, result);
 376               if (!(boost::math::isfinite)(result) || (result == 0) || (!std::numeric_limits<T>::has_infinity && (fabs(result) >= tools::max_value<T>())))
 377                  retry = true;
 378               if (norm / fabs(result) > 1 / boost::math::tools::root_epsilon<T>())
 379                  retry = true;  // fatal cancellation
 380            }
 381            catch (const std::overflow_error&)
 382            {
 383               retry = true;
 384            }
 385            catch (const boost::math::evaluation_error&)
 386            {
 387               retry = true;
 388            }
 389         }
 390         catch (const std::overflow_error&)
 391         {
 392            log_scale -= scale;
 393            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
 394         }
 395         catch (const boost::math::evaluation_error&)
 396         {
 397            log_scale -= scale;
 398            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
 399         }
 400         if (retry)
 401         {
 402            log_scale -= scale;
 403            log_scale -= series_scale;
 404            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
 405         }
 406         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_7<%1%>(%1%,%1%,%1%)", max_iter, pol);
 407         if (use_logs)
 408         {
 409            int sgn = boost::math::sign(result);
 410            prefix += log(fabs(result));
 411            result = sgn * prefix_sgn * exp(prefix);
 412         }
 413         else
 414         {
 415            if ((fabs(result) > 1) && (fabs(prefix) > 1) && (tools::max_value<T>() / fabs(result) < fabs(prefix)))
 416            {
 417               // Overflow:
 418               scale = itrunc(tools::log_max_value<T>()) - 10;
 419               log_scale += scale;
 420               result /= exp(T(scale));
 421            }
 422            result *= prefix;
 423         }
 424         return result;
 425      }
 426
 427
 428      template <class T>
 429      struct cyl_bessel_i_large_x_sum
 430      {
 431         typedef T result_type;
 432
 433         cyl_bessel_i_large_x_sum(const T& v, const T& x) : v(v), z(x), term(1), k(0) {}
 434
 435         T operator()()
 436         {
 437            T result = term;
 438            ++k;
 439            term *= -(4 * v * v - (2 * k - 1) * (2 * k - 1)) / (8 * k * z);
 440            return result;
 441         }
 442         T v, z, term;
 443         int k;
 444      };
 445
 446      template <class T, class Policy>
 447      T cyl_bessel_i_large_x_scaled(const T& v, const T& x, int& log_scaling, const Policy& pol)
 448      {
 449         BOOST_MATH_STD_USING
 450            cyl_bessel_i_large_x_sum<T> s(v, x);
 451         boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 452         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 453         boost::math::policies::check_series_iterations<T>("boost::math::cyl_bessel_i_large_x<%1%>(%1%,%1%)", max_iter, pol);
 454         int scale = boost::math::itrunc(x);
 455         log_scaling += scale;
 456         return result * exp(x - scale) / sqrt(boost::math::constants::two_pi<T>() * x);
 457      }
 458
 459
 460
 461      template <class T, class Policy>
 462      struct hypergeometric_1F1_AS_13_3_6_series
 463      {
 464         typedef T result_type;
 465
 466         enum { cache_size = 64 };
 467         //
 468         // This series is only convergent/useful for a and b approximately equal
 469         // (ideally |a-b| < 1).  The series can also go divergent after a while
 470         // when b < 0, which limits precision to around that of double.  In that
 471         // situation we return 0 to terminate the series as otherwise the divergent
 472         // terms will destroy all the bits in our result before they do eventually
 473         // converge again.  One important use case for this series is for z < 0
 474         // and |a| << |b| so that either b-a == b or at least most of the digits in a
 475         // are lost in the subtraction.  Note that while you can easily convince yourself
 476         // that the result should be unity when b-a == b, in fact this is not (quite)
 477         // the case for large z.
 478         //
 479         hypergeometric_1F1_AS_13_3_6_series(const T& a, const T& b, const T& z, const T& b_minus_a, const Policy& pol_)
 480            : b_minus_a(b_minus_a), half_z(z / 2), poch_1(2 * b_minus_a - 1), poch_2(b_minus_a - a), b_poch(b), term(1), last_result(1), sign(1), n(0), cache_offset(-cache_size), scale(0), pol(pol_)
 481         {
 482            bessel_i_cache[cache_size - 1] = half_z > tools::log_max_value<T>() ?
 483               cyl_bessel_i_large_x_scaled(T(b_minus_a - 1.5f), half_z, scale, pol) : boost::math::cyl_bessel_i(b_minus_a - 1.5f, half_z, pol);
 484            refill_cache();
 485         }
 486         T operator()()
 487         {
 488            BOOST_MATH_STD_USING
 489            if(n - cache_offset >= cache_size)
 490               refill_cache();
 491
 492            T result = term * (b_minus_a - 0.5f + n) * sign * bessel_i_cache[n - cache_offset];
 493            ++n;
 494            term *= poch_1;
 495            poch_1 = (n == 1) ? T(2 * b_minus_a) : T(poch_1 + 1);
 496            term *= poch_2;
 497            poch_2 += 1;
 498            term /= n;
 499            term /= b_poch;
 500            b_poch += 1;
 501            sign = -sign;
 502
 503            if ((fabs(result) > fabs(last_result)) && (n > 100))
 504               return 0;  // series has gone divergent!
 505
 506            last_result = result;
 507
 508            return result;
 509         }
 510
 511         int scaling()const
 512         {
 513            return scale;
 514         }
 515
 516      private:
 517         T b_minus_a, half_z, poch_1, poch_2, b_poch, term, last_result;
 518         int sign;
 519         int n, cache_offset, scale;
 520         const Policy& pol;
 521         std::array<T, cache_size> bessel_i_cache;
 522
 523         void refill_cache()
 524         {
 525            BOOST_MATH_STD_USING
 526            //
 527            // We don't calculate a new bessel I value: instead start our iterator off
 528            // with an arbitrary small value, then when we get back to the last value in the previous cache
 529            // calculate the ratio and use it to renormalise all the values.  This is more efficient, but
 530            // also avoids problems with I_v(x) underflowing to zero.
 531            //
 532            cache_offset += cache_size;
 533            T last_value = bessel_i_cache.back();
 534            bessel_i_backwards_iterator<T> i(b_minus_a + cache_offset + (int)cache_size - 1.5f, half_z, tools::min_value<T>() * (fabs(last_value) > 1 ? last_value : 1) / tools::epsilon<T>());
 535
 536            for (int j = cache_size - 1; j >= 0; --j, ++i)
 537            {
 538               bessel_i_cache[j] = *i;
 539               //
 540               // Depending on the value of half_z, the values stored in the cache can grow so
 541               // large as to overflow, if that looks likely then we need to rescale all the
 542               // existing terms (most of which will then underflow to zero).  In this situation
 543               // it's likely that our series will only need 1 or 2 terms of the series but we
 544               // can't be sure of that:
 545               //
 546               if((j < cache_size - 2) && (bessel_i_cache[j + 1] != 0) && (tools::max_value<T>() / fabs(64 * bessel_i_cache[j] / bessel_i_cache[j + 1]) < fabs(bessel_i_cache[j])))
 547               {
 548                  T rescale = pow(fabs(bessel_i_cache[j] / bessel_i_cache[j + 1]), j + 1) * 2;
 549                  if (rescale > tools::max_value<T>())
 550                     rescale = tools::max_value<T>();
 551                  for (int k = j; k < cache_size; ++k)
 552                     bessel_i_cache[k] /= rescale;
 553                  i = bessel_i_backwards_iterator<T>(b_minus_a -0.5f + cache_offset + j, half_z, bessel_i_cache[j + 1], bessel_i_cache[j]);
 554               }
 555            }
 556            T ratio = last_value / *i;
 557            for (auto j = bessel_i_cache.begin(); j != bessel_i_cache.end(); ++j)
 558               *j *= ratio;
 559         }
 560
 561         hypergeometric_1F1_AS_13_3_6_series();
 562         hypergeometric_1F1_AS_13_3_6_series(const hypergeometric_1F1_AS_13_3_6_series&);
 563         hypergeometric_1F1_AS_13_3_6_series& operator=(const hypergeometric_1F1_AS_13_3_6_series&);
 564      };
 565
 566      template <class T, class Policy>
 567      T hypergeometric_1F1_AS_13_3_6(const T& a, const T& b, const T& z, const T& b_minus_a, const Policy& pol, int& log_scaling)
 568      {
 569         BOOST_MATH_STD_USING
 570         if(b_minus_a == 0)
 571         {
 572            // special case: M(a,a,z) == exp(z)
 573            int scale = itrunc(z, pol);
 574            log_scaling += scale;
 575            return exp(z - scale);
 576         }
 577         hypergeometric_1F1_AS_13_3_6_series<T, Policy> s(a, b, z, b_minus_a, pol);
 578         boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 579         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 580         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_6<%1%>(%1%,%1%,%1%)", max_iter, pol);
 581         result *= boost::math::tgamma(b_minus_a - 0.5f) * pow(z / 4, -b_minus_a + 0.5f);
 582         int scale = itrunc(z / 2);
 583         log_scaling += scale;
 584         log_scaling += s.scaling();
 585         result *= exp(z / 2 - scale);
 586         return result;
 587      }
 588
 589      /****************************************************************************************************************/
 590      //
 591      // The following are not used at present and are commented out for that reason:
 592      //
 593      /****************************************************************************************************************/
 594
 595 #if 0
 596
 597      template <class T, class Policy>
 598      struct hypergeometric_1F1_AS_13_3_8_series
 599      {
 600         //
 601         // TODO: store and cache Bessel function evaluations via backwards recurrence.
 602         //
 603         // The C term grows by at least an order of magnitude with each iteration, and
 604         // rate of growth is largely independent of the arguments.  Free parameter h
 605         // seems to give accurate results for small values (almost zero) or h=1.
 606         // Convergence and accuracy, only when -a/z > 100, this appears to have no
 607         // or little benefit over 13.3.7 as it generally requires more iterations?
 608         //
 609         hypergeometric_1F1_AS_13_3_8_series(const T& a, const T& b, const T& z, const T& h, const Policy& pol_)
 610            : C_minus_2(1), C_minus_1(-b * h), C(b * (b + 1) * h * h / 2 - (2 * h - 1) * a / 2),
 611            bessel_arg(2 * sqrt(-a * z)), bessel_order(b - 1), power_term(std::pow(-a * z, (1 - b) / 2)), mult(z / std::sqrt(-a * z)),
 612            a_(a), b_(b), z_(z), h_(h), n(2), pol(pol_)
 613         {
 614         }
 615         T operator()()
 616         {
 617            // we actually return the n-2 term:
 618            T result = C_minus_2 * power_term * boost::math::cyl_bessel_j(bessel_order, bessel_arg, pol);
 619            bessel_order += 1;
 620            power_term *= mult;
 621            ++n;
 622            T C_next = ((1 - 2 * h_) * (n - 1) - b_ * h_) * C
 623               + ((1 - 2 * h_) * a_ - h_ * (h_ - 1) *(b_ + n - 2)) * C_minus_1
 624               - h_ * (h_ - 1) * a_ * C_minus_2;
 625            C_next /= n;
 626            C_minus_2 = C_minus_1;
 627            C_minus_1 = C;
 628            C = C_next;
 629            return result;
 630         }
 631         T C, C_minus_1, C_minus_2, bessel_arg, bessel_order, power_term, mult, a_, b_, z_, h_;
 632         const Policy& pol;
 633         int n;
 634
 635         typedef T result_type;
 636      };
 637
 638      template <class T, class Policy>
 639      T hypergeometric_1F1_AS_13_3_8(const T& a, const T& b, const T& z, const T& h, const Policy& pol)
 640      {
 641         BOOST_MATH_STD_USING
 642         T prefix = exp(h * z) * boost::math::tgamma(b);
 643         hypergeometric_1F1_AS_13_3_8_series<T, Policy> s(a, b, z, h, pol);
 644         boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 645         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 646         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_8<%1%>(%1%,%1%,%1%)", max_iter, pol);
 647         result *= prefix;
 648         return result;
 649      }
 650
 651      //
 652      // This is the series from https://dlmf.nist.gov/13.11
 653      // It appears to be unusable for a,z < 0, and for
 654      // b < 0 appears to never be better than the defining series
 655      // for 1F1.
 656      //
 657      template <class T, class Policy>
 658      struct hypergeometric_1f1_13_11_1_series
 659      {
 660         typedef T result_type;
 661
 662         hypergeometric_1f1_13_11_1_series(const T& a, const T& b, const T& z, const Policy& pol_)
 663            : term(1), two_a_minus_1_plus_s(2 * a - 1), two_a_minus_b_plus_s(2 * a - b), b_plus_s(b), a_minus_half_plus_s(a - 0.5f), half_z(z / 2), s(0), pol(pol_)
 664         {
 665         }
 666         T operator()()
 667         {
 668            T result = term * a_minus_half_plus_s * boost::math::cyl_bessel_i(a_minus_half_plus_s, half_z, pol);
 669
 670            term *= two_a_minus_1_plus_s * two_a_minus_b_plus_s / (b_plus_s * ++s);
 671            two_a_minus_1_plus_s += 1;
 672            a_minus_half_plus_s += 1;
 673            two_a_minus_b_plus_s += 1;
 674            b_plus_s += 1;
 675
 676            return result;
 677         }
 678         T term, two_a_minus_1_plus_s, two_a_minus_b_plus_s, b_plus_s, a_minus_half_plus_s, half_z;
 679         int s;
 680         const Policy& pol;
 681      };
 682
 683      template <class T, class Policy>
 684      T hypergeometric_1f1_13_11_1(const T& a, const T& b, const T& z, const Policy& pol, int& log_scaling)
 685      {
 686         BOOST_MATH_STD_USING
 687            bool use_logs = false;
 688         T prefix;
 689         int prefix_sgn = 1;
 690         if (true/*(a < boost::math::max_factorial<T>::value) && (a > 0)*/)
 691            prefix = boost::math::tgamma(a - 0.5f, pol);
 692         else
 693         {
 694            prefix = boost::math::lgamma(a - 0.5f, &prefix_sgn, pol);
 695            use_logs = true;
 696         }
 697         if (use_logs)
 698         {
 699            prefix += z / 2;
 700            prefix += log(z / 4) * (0.5f - a);
 701         }
 702         else if (z > 0)
 703         {
 704            prefix *= pow(z / 4, 0.5f - a);
 705            prefix *= exp(z / 2);
 706         }
 707         else
 708         {
 709            prefix *= exp(z / 2);
 710            prefix *= pow(z / 4, 0.5f - a);
 711         }
 712
 713         hypergeometric_1f1_13_11_1_series<T, Policy> s(a, b, z, pol);
 714         boost::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
 715         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
 716         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_13_11_1<%1%>(%1%,%1%,%1%)", max_iter, pol);
 717         if (use_logs)
 718         {
 719            int scaling = itrunc(prefix);
 720            log_scaling += scaling;
 721            prefix -= scaling;
 722            result *= exp(prefix) * prefix_sgn;
 723         }
 724         else
 725            result *= prefix;
 726
 727         return result;
 728      }
 729
 730 #endif
 731
 732   } } } // namespaces
 733
 734 #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP