special_functions/detail/hypergeometric_1F1_bessel.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2014 Anton Bikineev
0003 //  Copyright 2014 Christopher Kormanyos
0004 //  Copyright 2014 John Maddock
0005 //  Copyright 2014 Paul Bristow
0006 //  Distributed under the Boost
0007 //  Software License, Version 1.0. (See accompanying file
0008 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0009 //
0010 #ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP
0011 #define BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP
0012
0013 #include <boost/math/tools/series.hpp>
0014 #include <boost/math/special_functions/bessel.hpp>
0015 #include <boost/math/special_functions/laguerre.hpp>
0016 #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
0017 #include <boost/math/special_functions/bessel_iterators.hpp>
0018
0019
0020   namespace boost { namespace math { namespace detail {
0021
0022      template <class T, class Policy>
0023      T hypergeometric_1F1_divergent_fallback(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling);
0024
0025      template <class T>
0026      bool hypergeometric_1F1_is_tricomi_viable_positive_b(const T& a, const T& b, const T& z)
0027      {
0028         BOOST_MATH_STD_USING
0029            if ((z < b) && (a > -50))
0030               return false;  // might as well fall through to recursion
0031         if (b <= 100)
0032            return true;
0033         // Even though we're in a reasonable domain for Tricomi's approximation,
0034         // the arguments to the Bessel functions may be so large that we can't
0035         // actually evaluate them:
0036         T x = sqrt(fabs(2 * z * b - 4 * a * z));
0037         T v = b - 1;
0038         return log(boost::math::constants::e<T>() * x / (2 * v)) * v > tools::log_min_value<T>();
0039      }
0040
0041      //
0042      // Returns an arbitrarily small value compared to "target" for use as a seed
0043      // value for Bessel recurrences.  Note that we'd better not make it too small
0044      // or underflow may occur resulting in either one of the terms in the
0045      // recurrence being zero, or else the result being zero.  Using 1/epsilon
0046      // as a safety factor ensures that if we do underflow to zero, all of the digits
0047      // will have been cancelled out anyway:
0048      //
0049      template <class T>
0050      T arbitrary_small_value(const T& target)
0051      {
0052         using std::fabs;
0053         return (tools::min_value<T>() / tools::epsilon<T>()) * (fabs(target) > 1 ? target : 1);
0054      }
0055
0056
0057      template <class T, class Policy>
0058      struct hypergeometric_1F1_AS_13_3_7_tricomi_series
0059      {
0060         typedef T result_type;
0061
0062         enum { cache_size = 64 };
0063
0064         hypergeometric_1F1_AS_13_3_7_tricomi_series(const T& a, const T& b, const T& z, const Policy& pol_)
0065            : A_minus_2(1), A_minus_1(0), A(b / 2), mult(z / 2), term(1), b_minus_1_plus_n(b - 1),
0066             bessel_arg((b / 2 - a) * z),
0067            two_a_minus_b(2 * a - b), pol(pol_), n(2)
0068         {
0069            BOOST_MATH_STD_USING
0070            term /= pow(fabs(bessel_arg), b_minus_1_plus_n / 2);
0071            mult /= sqrt(fabs(bessel_arg));
0072            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);
0073            if (fabs(bessel_cache[cache_size - 1]) < tools::min_value<T>() / tools::epsilon<T>())
0074            {
0075               // We get very limited precision due to rapid denormalisation/underflow of the Bessel values, raise an exception and try something else:
0076               policies::raise_evaluation_error("hypergeometric_1F1_AS_13_3_7_tricomi_series<%1%>", "Underflow in Bessel functions", bessel_cache[cache_size - 1], pol);
0077            }
0078            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>())))
0079            {
0080               term = -log(fabs(bessel_arg)) * b_minus_1_plus_n / 2;
0081               log_scale = lltrunc(term);
0082               term -= log_scale;
0083               term = exp(term);
0084            }
0085            else
0086               log_scale = 0;
0087 #ifndef BOOST_NO_CXX17_IF_CONSTEXPR
0088            if constexpr (std::numeric_limits<T>::has_infinity)
0089            {
0090               if (!(boost::math::isfinite)(bessel_cache[cache_size - 1]))
0091                  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);
0092            }
0093            else
0094               if ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))
0095                  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);
0096 #else
0097            if ((std::numeric_limits<T>::has_infinity && !(boost::math::isfinite)(bessel_cache[cache_size - 1]))
0098               || (!std::numeric_limits<T>::has_infinity && ((boost::math::isnan)(bessel_cache[cache_size - 1]) || (fabs(bessel_cache[cache_size - 1]) >= tools::max_value<T>()))))
0099               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);
0100 #endif
0101            cache_offset = -cache_size;
0102            refill_cache();
0103         }
0104         T operator()()
0105         {
0106            //
0107            // We return the n-2 term, and do 2 terms at once as every other term can be
0108            // very small (or zero) when b == 2a:
0109            //
0110            BOOST_MATH_STD_USING
0111            if(n - 2 - cache_offset >= cache_size)
0112               refill_cache();
0113            T result = A_minus_2 * term * bessel_cache[n - 2 - cache_offset];
0114            term *= mult;
0115            ++n;
0116            T A_next = ((b_minus_1_plus_n + 2) * A_minus_1 + two_a_minus_b * A_minus_2) / n;
0117            b_minus_1_plus_n += 1;
0118            A_minus_2 = A_minus_1;
0119            A_minus_1 = A;
0120            A = A_next;
0121
0122            if (A_minus_2 != 0)
0123            {
0124               if (n - 2 - cache_offset >= cache_size)
0125                  refill_cache();
0126               result += A_minus_2 * term * bessel_cache[n - 2 - cache_offset];
0127            }
0128            term *= mult;
0129            ++n;
0130            A_next = ((b_minus_1_plus_n + 2) * A_minus_1 + two_a_minus_b * A_minus_2) / n;
0131            b_minus_1_plus_n += 1;
0132            A_minus_2 = A_minus_1;
0133            A_minus_1 = A;
0134            A = A_next;
0135
0136            return result;
0137         }
0138
0139         long long scale()const
0140         {
0141            return log_scale;
0142         }
0143
0144      private:
0145         T A_minus_2, A_minus_1, A, mult, term, b_minus_1_plus_n, bessel_arg, two_a_minus_b;
0146         std::array<T, cache_size> bessel_cache;
0147         const Policy& pol;
0148         int n, cache_offset;
0149         long long log_scale;
0150
0151         hypergeometric_1F1_AS_13_3_7_tricomi_series operator=(const hypergeometric_1F1_AS_13_3_7_tricomi_series&) = delete;
0152
0153         void refill_cache()
0154         {
0155            BOOST_MATH_STD_USING
0156            //
0157            // We don't calculate a new bessel I/J value: instead start our iterator off
0158            // with an arbitrary small value, then when we get back to the last value in the previous cache
0159            // calculate the ratio and use it to renormalise all the new values.  This is more efficient, but
0160            // also avoids problems with J_v(x) or I_v(x) underflowing to zero.
0161            //
0162            cache_offset += cache_size;
0163            T last_value = bessel_cache.back();
0164            T ratio;
0165            if (bessel_arg > 0)
0166            {
0167               //
0168               // We will be calculating Bessel J.
0169               // We need a different recurrence strategy for positive and negative orders:
0170               //
0171               if (b_minus_1_plus_n > 0)
0172               {
0173                  bessel_j_backwards_iterator<T, Policy> i(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(bessel_arg), arbitrary_small_value(last_value));
0174
0175                  for (int j = cache_size - 1; j >= 0; --j, ++i)
0176                  {
0177                     bessel_cache[j] = *i;
0178                     //
0179                     // Depending on the value of bessel_arg, the values stored in the cache can grow so
0180                     // large as to overflow, if that looks likely then we need to rescale all the
0181                     // existing terms (most of which will then underflow to zero).  In this situation
0182                     // it's likely that our series will only need 1 or 2 terms of the series but we
0183                     // can't be sure of that:
0184                     //
0185                     if ((j < cache_size - 2) && (tools::max_value<T>() / fabs(64 * bessel_cache[j] / bessel_cache[j + 1]) < fabs(bessel_cache[j])))
0186                     {
0187                        T rescale = static_cast<T>(pow(fabs(bessel_cache[j] / bessel_cache[j + 1]), T(j + 1)) * 2);
0188                        if (!((boost::math::isfinite)(rescale)))
0189                           rescale = tools::max_value<T>();
0190                        for (int k = j; k < cache_size; ++k)
0191                           bessel_cache[k] /= rescale;
0192                        bessel_j_backwards_iterator<T, Policy> ti(b_minus_1_plus_n + j, 2 * sqrt(bessel_arg), bessel_cache[j + 1], bessel_cache[j]);
0193                        i = ti;
0194                     }
0195                  }
0196                  ratio = last_value / *i;
0197               }
0198               else
0199               {
0200                  //
0201                  // Negative order is difficult: the J_v(x) recurrence relations are unstable
0202                  // *in both directions* for v < 0, except as v -> -INF when we have
0203                  // J_-v(x)  ~= -sin(pi.v)Y_v(x).
0204                  // For small v what we can do is compute every other Bessel function and
0205                  // then fill in the gaps using the recurrence relation.  This *is* stable
0206                  // provided that v is not so negative that the above approximation holds.
0207                  //
0208                  bessel_cache[1] = cyl_bessel_j(b_minus_1_plus_n + 1, 2 * sqrt(bessel_arg), pol);
0209                  bessel_cache[0] = (last_value + bessel_cache[1]) / (b_minus_1_plus_n / sqrt(bessel_arg));
0210                  int pos = 2;
0211                  while ((pos < cache_size - 1) && (b_minus_1_plus_n + pos < 0))
0212                  {
0213                     bessel_cache[pos + 1] = cyl_bessel_j(b_minus_1_plus_n + pos + 1, 2 * sqrt(bessel_arg), pol);
0214                     bessel_cache[pos] = (bessel_cache[pos-1] + bessel_cache[pos+1]) / ((b_minus_1_plus_n + pos) / sqrt(bessel_arg));
0215                     pos += 2;
0216                  }
0217                  if (pos < cache_size)
0218                  {
0219                     //
0220                     // We have crossed over into the region where backward recursion is the stable direction
0221                     // start from arbitrary value and recurse down to "pos" and normalise:
0222                     //
0223                     bessel_j_backwards_iterator<T, Policy> i2(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(bessel_arg), arbitrary_small_value(bessel_cache[pos-1]));
0224                     for (int loc = cache_size - 1; loc >= pos; --loc)
0225                        bessel_cache[loc] = *i2++;
0226                     ratio = bessel_cache[pos - 1] / *i2;
0227                     //
0228                     // Sanity check, if we normalised to an unusually small value then it was likely
0229                     // to be near a root and the calculated ratio is garbage, if so perform one
0230                     // more J_v(x) evaluation at position and normalise again:
0231                     //
0232                     if (fabs(bessel_cache[pos] * ratio / bessel_cache[pos - 1]) > 5)
0233                        ratio = cyl_bessel_j(b_minus_1_plus_n + pos, 2 * sqrt(bessel_arg), pol) / bessel_cache[pos];
0234                     while (pos < cache_size)
0235                        bessel_cache[pos++] *= ratio;
0236                  }
0237                  ratio = 1;
0238               }
0239            }
0240            else
0241            {
0242               //
0243               // Bessel I.
0244               // We need a different recurrence strategy for positive and negative orders:
0245               //
0246               if (b_minus_1_plus_n > 0)
0247               {
0248                  bessel_i_backwards_iterator<T, Policy> i(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(-bessel_arg), arbitrary_small_value(last_value));
0249
0250                  for (int j = cache_size - 1; j >= 0; --j, ++i)
0251                  {
0252                     bessel_cache[j] = *i;
0253                     //
0254                     // Depending on the value of bessel_arg, the values stored in the cache can grow so
0255                     // large as to overflow, if that looks likely then we need to rescale all the
0256                     // existing terms (most of which will then underflow to zero).  In this situation
0257                     // it's likely that our series will only need 1 or 2 terms of the series but we
0258                     // can't be sure of that:
0259                     //
0260                     if ((j < cache_size - 2) && (tools::max_value<T>() / fabs(64 * bessel_cache[j] / bessel_cache[j + 1]) < fabs(bessel_cache[j])))
0261                     {
0262                        T rescale = static_cast<T>(pow(fabs(bessel_cache[j] / bessel_cache[j + 1]), T(j + 1)) * 2);
0263                        if (!((boost::math::isfinite)(rescale)))
0264                           rescale = tools::max_value<T>();
0265                        for (int k = j; k < cache_size; ++k)
0266                           bessel_cache[k] /= rescale;
0267                        i = bessel_i_backwards_iterator<T, Policy>(b_minus_1_plus_n + j, 2 * sqrt(-bessel_arg), bessel_cache[j + 1], bessel_cache[j]);
0268                     }
0269                  }
0270                  ratio = last_value / *i;
0271               }
0272               else
0273               {
0274                  //
0275                  // Forwards iteration is stable:
0276                  //
0277                  bessel_i_forwards_iterator<T, Policy> i(b_minus_1_plus_n, 2 * sqrt(-bessel_arg));
0278                  int pos = 0;
0279                  while ((pos < cache_size) && (b_minus_1_plus_n + pos < 0.5))
0280                  {
0281                     bessel_cache[pos++] = *i++;
0282                  }
0283                  if (pos < cache_size)
0284                  {
0285                     //
0286                     // We have crossed over into the region where backward recursion is the stable direction
0287                     // start from arbitrary value and recurse down to "pos" and normalise:
0288                     //
0289                     bessel_i_backwards_iterator<T, Policy> i2(b_minus_1_plus_n + (int)cache_size - 1, 2 * sqrt(-bessel_arg), arbitrary_small_value(last_value));
0290                     for (int loc = cache_size - 1; loc >= pos; --loc)
0291                        bessel_cache[loc] = *i2++;
0292                     ratio = bessel_cache[pos - 1] / *i2;
0293                     while (pos < cache_size)
0294                        bessel_cache[pos++] *= ratio;
0295                  }
0296                  ratio = 1;
0297               }
0298            }
0299            if(ratio != 1)
0300               for (auto j = bessel_cache.begin(); j != bessel_cache.end(); ++j)
0301                  *j *= ratio;
0302            //
0303            // Very occasionally our normalisation fails because the normalisztion value
0304            // is sitting right on top of a root (or very close to it).  When that happens
0305            // best to calculate a fresh Bessel evaluation and normalise again.
0306            //
0307            if (fabs(bessel_cache[0] / last_value) > 5)
0308            {
0309               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];
0310               if (ratio != 1)
0311                  for (auto j = bessel_cache.begin(); j != bessel_cache.end(); ++j)
0312                     *j *= ratio;
0313            }
0314         }
0315      };
0316
0317      template <class T, class Policy>
0318      T hypergeometric_1F1_AS_13_3_7_tricomi(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scale)
0319      {
0320         BOOST_MATH_STD_USING
0321         //
0322         // Works for a < 0, b < 0, z > 0.
0323         //
0324         // For convergence we require A * term to be converging otherwise we get
0325         // a divergent alternating series.  It's actually really hard to analyse this
0326         // and the best purely heuristic policy we've found is
0327         // z < fabs((2 * a - b) / (sqrt(fabs(a)))) ; b > 0  or:
0328         // z < fabs((2 * a - b) / (sqrt(fabs(ab)))) ; b < 0
0329         //
0330         T prefix(0);
0331         int prefix_sgn(0);
0332         bool use_logs = false;
0333         long long scale = 0;
0334         //
0335         // We can actually support the b == 2a case within here, but we haven't
0336         // as we appear never to get here in practice.  Which means this get out
0337         // clause is a bit of defensive programming....
0338         //
0339         if(b == 2 * a)
0340            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
0341
0342 #ifndef BOOST_NO_EXCEPTIONS
0343         try
0344 #endif
0345         {
0346            prefix = boost::math::tgamma(b, pol);
0347            prefix *= exp(z / 2);
0348         }
0349 #ifndef BOOST_NO_EXCEPTIONS
0350         catch (const std::runtime_error&)
0351         {
0352            use_logs = true;
0353         }
0354 #endif
0355         if (use_logs || (prefix == 0) || !(boost::math::isfinite)(prefix) || (!std::numeric_limits<T>::has_infinity && (fabs(prefix) >= tools::max_value<T>())))
0356         {
0357            use_logs = true;
0358            prefix = boost::math::lgamma(b, &prefix_sgn, pol) + z / 2;
0359            scale = lltrunc(prefix);
0360            log_scale += scale;
0361            prefix -= scale;
0362         }
0363         T result(0);
0364         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0365         bool retry = false;
0366         long long series_scale = 0;
0367 #ifndef BOOST_NO_EXCEPTIONS
0368         try
0369 #endif
0370         {
0371            hypergeometric_1F1_AS_13_3_7_tricomi_series<T, Policy> s(a, b, z, pol);
0372            series_scale = s.scale();
0373            log_scale += s.scale();
0374 #ifndef BOOST_NO_EXCEPTIONS
0375            try
0376 #endif
0377            {
0378               T norm = 0;
0379               result = 0;
0380               if((a < 0) && (b < 0))
0381                  result = boost::math::tools::checked_sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, result, norm);
0382               else
0383                  result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter, result);
0384               if (!(boost::math::isfinite)(result) || (result == 0) || (!std::numeric_limits<T>::has_infinity && (fabs(result) >= tools::max_value<T>())))
0385                  retry = true;
0386               if (norm / fabs(result) > 1 / boost::math::tools::root_epsilon<T>())
0387                  retry = true;  // fatal cancellation
0388            }
0389 #ifndef BOOST_NO_EXCEPTIONS
0390            catch (const std::overflow_error&)
0391            {
0392               retry = true;
0393            }
0394            catch (const boost::math::evaluation_error&)
0395            {
0396               retry = true;
0397            }
0398 #endif
0399         }
0400 #ifndef BOOST_NO_EXCEPTIONS
0401         catch (const std::overflow_error&)
0402         {
0403            log_scale -= scale;
0404            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
0405         }
0406         catch (const boost::math::evaluation_error&)
0407         {
0408            log_scale -= scale;
0409            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
0410         }
0411 #endif
0412         if (retry)
0413         {
0414            log_scale -= scale;
0415            log_scale -= series_scale;
0416            return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scale);
0417         }
0418         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_7<%1%>(%1%,%1%,%1%)", max_iter, pol);
0419         if (use_logs)
0420         {
0421            int sgn = boost::math::sign(result);
0422            prefix += log(fabs(result));
0423            result = sgn * prefix_sgn * exp(prefix);
0424         }
0425         else
0426         {
0427            if ((fabs(result) > 1) && (fabs(prefix) > 1) && (tools::max_value<T>() / fabs(result) < fabs(prefix)))
0428            {
0429               // Overflow:
0430               scale = lltrunc(tools::log_max_value<T>()) - 10;
0431               log_scale += scale;
0432               result /= exp(T(scale));
0433            }
0434            result *= prefix;
0435         }
0436         return result;
0437      }
0438
0439
0440      template <class T>
0441      struct cyl_bessel_i_large_x_sum
0442      {
0443         typedef T result_type;
0444
0445         cyl_bessel_i_large_x_sum(const T& v, const T& x) : v(v), z(x), term(1), k(0) {}
0446
0447         T operator()()
0448         {
0449            T result = term;
0450            ++k;
0451            term *= -(4 * v * v - (2 * k - 1) * (2 * k - 1)) / (8 * k * z);
0452            return result;
0453         }
0454         T v, z, term;
0455         int k;
0456      };
0457
0458      template <class T, class Policy>
0459      T cyl_bessel_i_large_x_scaled(const T& v, const T& x, long long& log_scaling, const Policy& pol)
0460      {
0461         BOOST_MATH_STD_USING
0462            cyl_bessel_i_large_x_sum<T> s(v, x);
0463         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0464         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0465         boost::math::policies::check_series_iterations<T>("boost::math::cyl_bessel_i_large_x<%1%>(%1%,%1%)", max_iter, pol);
0466         long long scale = boost::math::lltrunc(x);
0467         log_scaling += scale;
0468         return result * exp(x - scale) / sqrt(boost::math::constants::two_pi<T>() * x);
0469      }
0470
0471
0472
0473      template <class T, class Policy>
0474      struct hypergeometric_1F1_AS_13_3_6_series
0475      {
0476         typedef T result_type;
0477
0478         enum { cache_size = 64 };
0479         //
0480         // This series is only convergent/useful for a and b approximately equal
0481         // (ideally |a-b| < 1).  The series can also go divergent after a while
0482         // when b < 0, which limits precision to around that of double.  In that
0483         // situation we return 0 to terminate the series as otherwise the divergent
0484         // terms will destroy all the bits in our result before they do eventually
0485         // converge again.  One important use case for this series is for z < 0
0486         // and |a| << |b| so that either b-a == b or at least most of the digits in a
0487         // are lost in the subtraction.  Note that while you can easily convince yourself
0488         // that the result should be unity when b-a == b, in fact this is not (quite)
0489         // the case for large z.
0490         //
0491         hypergeometric_1F1_AS_13_3_6_series(const T& a, const T& b, const T& z, const T& b_minus_a, const Policy& pol_)
0492            : 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_)
0493         {
0494            bessel_i_cache[cache_size - 1] = half_z > tools::log_max_value<T>() ?
0495               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);
0496            refill_cache();
0497         }
0498         T operator()()
0499         {
0500            BOOST_MATH_STD_USING
0501            if(n - cache_offset >= cache_size)
0502               refill_cache();
0503
0504            T result = term * (b_minus_a - 0.5f + n) * sign * bessel_i_cache[n - cache_offset];
0505            ++n;
0506            term *= poch_1;
0507            poch_1 = (n == 1) ? T(2 * b_minus_a) : T(poch_1 + 1);
0508            term *= poch_2;
0509            poch_2 += 1;
0510            term /= n;
0511            term /= b_poch;
0512            b_poch += 1;
0513            sign = -sign;
0514
0515            if ((fabs(result) > fabs(last_result)) && (n > 100))
0516               return 0;  // series has gone divergent!
0517
0518            last_result = result;
0519
0520            return result;
0521         }
0522
0523         long long scaling()const
0524         {
0525            return scale;
0526         }
0527
0528      private:
0529         T b_minus_a, half_z, poch_1, poch_2, b_poch, term, last_result;
0530         int sign;
0531         int n, cache_offset;
0532         long long scale;
0533         const Policy& pol;
0534         std::array<T, cache_size> bessel_i_cache;
0535
0536         void refill_cache()
0537         {
0538            BOOST_MATH_STD_USING
0539            //
0540            // We don't calculate a new bessel I value: instead start our iterator off
0541            // with an arbitrary small value, then when we get back to the last value in the previous cache
0542            // calculate the ratio and use it to renormalise all the values.  This is more efficient, but
0543            // also avoids problems with I_v(x) underflowing to zero.
0544            //
0545            cache_offset += cache_size;
0546            T last_value = bessel_i_cache.back();
0547            bessel_i_backwards_iterator<T, Policy> 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>());
0548
0549            for (int j = cache_size - 1; j >= 0; --j, ++i)
0550            {
0551               bessel_i_cache[j] = *i;
0552               //
0553               // Depending on the value of half_z, the values stored in the cache can grow so
0554               // large as to overflow, if that looks likely then we need to rescale all the
0555               // existing terms (most of which will then underflow to zero).  In this situation
0556               // it's likely that our series will only need 1 or 2 terms of the series but we
0557               // can't be sure of that:
0558               //
0559               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])))
0560               {
0561                  T rescale = static_cast<T>(pow(fabs(bessel_i_cache[j] / bessel_i_cache[j + 1]), T(j + 1)) * 2);
0562                  if (rescale > tools::max_value<T>())
0563                     rescale = tools::max_value<T>();
0564                  for (int k = j; k < cache_size; ++k)
0565                     bessel_i_cache[k] /= rescale;
0566                  i = bessel_i_backwards_iterator<T, Policy>(b_minus_a -0.5f + cache_offset + j, half_z, bessel_i_cache[j + 1], bessel_i_cache[j]);
0567               }
0568            }
0569            T ratio = last_value / *i;
0570            for (auto j = bessel_i_cache.begin(); j != bessel_i_cache.end(); ++j)
0571               *j *= ratio;
0572         }
0573
0574         hypergeometric_1F1_AS_13_3_6_series() = delete;
0575         hypergeometric_1F1_AS_13_3_6_series(const hypergeometric_1F1_AS_13_3_6_series&) = delete;
0576         hypergeometric_1F1_AS_13_3_6_series& operator=(const hypergeometric_1F1_AS_13_3_6_series&) = delete;
0577      };
0578
0579      template <class T, class Policy>
0580      T hypergeometric_1F1_AS_13_3_6(const T& a, const T& b, const T& z, const T& b_minus_a, const Policy& pol, long long& log_scaling)
0581      {
0582         BOOST_MATH_STD_USING
0583         if(b_minus_a == 0)
0584         {
0585            // special case: M(a,a,z) == exp(z)
0586            long long scale = lltrunc(z, pol);
0587            log_scaling += scale;
0588            return exp(z - scale);
0589         }
0590         hypergeometric_1F1_AS_13_3_6_series<T, Policy> s(a, b, z, b_minus_a, pol);
0591         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0592         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0593         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_6<%1%>(%1%,%1%,%1%)", max_iter, pol);
0594         result *= boost::math::tgamma(b_minus_a - 0.5f, pol) * pow(z / 4, -b_minus_a + T(0.5f));
0595         long long scale = lltrunc(z / 2);
0596         log_scaling += scale;
0597         log_scaling += s.scaling();
0598         result *= exp(z / 2 - scale);
0599         return result;
0600      }
0601
0602      /****************************************************************************************************************/
0603      //
0604      // The following are not used at present and are commented out for that reason:
0605      //
0606      /****************************************************************************************************************/
0607
0608 #if 0
0609
0610      template <class T, class Policy>
0611      struct hypergeometric_1F1_AS_13_3_8_series
0612      {
0613         //
0614         // TODO: store and cache Bessel function evaluations via backwards recurrence.
0615         //
0616         // The C term grows by at least an order of magnitude with each iteration, and
0617         // rate of growth is largely independent of the arguments.  Free parameter h
0618         // seems to give accurate results for small values (almost zero) or h=1.
0619         // Convergence and accuracy, only when -a/z > 100, this appears to have no
0620         // or little benefit over 13.3.7 as it generally requires more iterations?
0621         //
0622         hypergeometric_1F1_AS_13_3_8_series(const T& a, const T& b, const T& z, const T& h, const Policy& pol_)
0623            : C_minus_2(1), C_minus_1(-b * h), C(b * (b + 1) * h * h / 2 - (2 * h - 1) * a / 2),
0624            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)),
0625            a_(a), b_(b), z_(z), h_(h), n(2), pol(pol_)
0626         {
0627         }
0628         T operator()()
0629         {
0630            // we actually return the n-2 term:
0631            T result = C_minus_2 * power_term * boost::math::cyl_bessel_j(bessel_order, bessel_arg, pol);
0632            bessel_order += 1;
0633            power_term *= mult;
0634            ++n;
0635            T C_next = ((1 - 2 * h_) * (n - 1) - b_ * h_) * C
0636               + ((1 - 2 * h_) * a_ - h_ * (h_ - 1) *(b_ + n - 2)) * C_minus_1
0637               - h_ * (h_ - 1) * a_ * C_minus_2;
0638            C_next /= n;
0639            C_minus_2 = C_minus_1;
0640            C_minus_1 = C;
0641            C = C_next;
0642            return result;
0643         }
0644         T C, C_minus_1, C_minus_2, bessel_arg, bessel_order, power_term, mult, a_, b_, z_, h_;
0645         const Policy& pol;
0646         int n;
0647
0648         typedef T result_type;
0649      };
0650
0651      template <class T, class Policy>
0652      T hypergeometric_1F1_AS_13_3_8(const T& a, const T& b, const T& z, const T& h, const Policy& pol)
0653      {
0654         BOOST_MATH_STD_USING
0655         T prefix = exp(h * z) * boost::math::tgamma(b);
0656         hypergeometric_1F1_AS_13_3_8_series<T, Policy> s(a, b, z, h, pol);
0657         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0658         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0659         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_AS_13_3_8<%1%>(%1%,%1%,%1%)", max_iter, pol);
0660         result *= prefix;
0661         return result;
0662      }
0663
0664      //
0665      // This is the series from https://dlmf.nist.gov/13.11
0666      // It appears to be unusable for a,z < 0, and for
0667      // b < 0 appears to never be better than the defining series
0668      // for 1F1.
0669      //
0670      template <class T, class Policy>
0671      struct hypergeometric_1f1_13_11_1_series
0672      {
0673         typedef T result_type;
0674
0675         hypergeometric_1f1_13_11_1_series(const T& a, const T& b, const T& z, const Policy& pol_)
0676            : 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_)
0677         {
0678         }
0679         T operator()()
0680         {
0681            T result = term * a_minus_half_plus_s * boost::math::cyl_bessel_i(a_minus_half_plus_s, half_z, pol);
0682
0683            term *= two_a_minus_1_plus_s * two_a_minus_b_plus_s / (b_plus_s * ++s);
0684            two_a_minus_1_plus_s += 1;
0685            a_minus_half_plus_s += 1;
0686            two_a_minus_b_plus_s += 1;
0687            b_plus_s += 1;
0688
0689            return result;
0690         }
0691         T term, two_a_minus_1_plus_s, two_a_minus_b_plus_s, b_plus_s, a_minus_half_plus_s, half_z;
0692         long long s;
0693         const Policy& pol;
0694      };
0695
0696      template <class T, class Policy>
0697      T hypergeometric_1f1_13_11_1(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0698      {
0699         BOOST_MATH_STD_USING
0700            bool use_logs = false;
0701         T prefix;
0702         int prefix_sgn = 1;
0703         if (true/*(a < boost::math::max_factorial<T>::value) && (a > 0)*/)
0704            prefix = boost::math::tgamma(a - 0.5f, pol);
0705         else
0706         {
0707            prefix = boost::math::lgamma(a - 0.5f, &prefix_sgn, pol);
0708            use_logs = true;
0709         }
0710         if (use_logs)
0711         {
0712            prefix += z / 2;
0713            prefix += log(z / 4) * (0.5f - a);
0714         }
0715         else if (z > 0)
0716         {
0717            prefix *= pow(z / 4, 0.5f - a);
0718            prefix *= exp(z / 2);
0719         }
0720         else
0721         {
0722            prefix *= exp(z / 2);
0723            prefix *= pow(z / 4, 0.5f - a);
0724         }
0725
0726         hypergeometric_1f1_13_11_1_series<T, Policy> s(a, b, z, pol);
0727         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0728         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0729         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1f1_13_11_1<%1%>(%1%,%1%,%1%)", max_iter, pol);
0730         if (use_logs)
0731         {
0732            long long scaling = lltrunc(prefix);
0733            log_scaling += scaling;
0734            prefix -= scaling;
0735            result *= exp(prefix) * prefix_sgn;
0736         }
0737         else
0738            result *= prefix;
0739
0740         return result;
0741      }
0742
0743 #endif
0744
0745   } } } // namespaces
0746
0747 #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_BESSEL_HPP