special_functions/detail/hypergeometric_1F1_large_abz.hpp

0001
0002 ///////////////////////////////////////////////////////////////////////////////
0003 //  Copyright 2018 John Maddock
0004 //  Distributed under the Boost
0005 //  Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007
0008 #ifndef BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_
0009 #define BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_
0010
0011 #include <boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp>
0012 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
0013 #include <boost/math/special_functions/gamma.hpp>
0014 #include <boost/math/special_functions/trunc.hpp>
0015
0016   namespace boost { namespace math { namespace detail {
0017
0018      template <class T>
0019      inline bool is_negative_integer(const T& x)
0020      {
0021         using std::floor;
0022         return (x <= 0) && (floor(x) == x);
0023      }
0024
0025
0026      template <class T, class Policy>
0027      struct hypergeometric_1F1_igamma_series
0028      {
0029         enum{ cache_size = 64 };
0030
0031         typedef T result_type;
0032         hypergeometric_1F1_igamma_series(const T& alpha, const T& delta, const T& x, const Policy& pol)
0033            : delta_poch(-delta), alpha_poch(alpha), x(x), k(0), cache_offset(0), pol(pol)
0034         {
0035            BOOST_MATH_STD_USING
0036            T log_term = log(x) * -alpha;
0037            log_scaling = lltrunc(log_term - 3 - boost::math::tools::log_min_value<T>() / 50);
0038            term = exp(log_term - log_scaling);
0039            refill_cache();
0040         }
0041         T operator()()
0042         {
0043            if (k - cache_offset >= cache_size)
0044            {
0045               cache_offset += cache_size;
0046               refill_cache();
0047            }
0048            T result = term * gamma_cache[k - cache_offset];
0049            term *= delta_poch * alpha_poch / (++k * x);
0050            delta_poch += 1;
0051            alpha_poch += 1;
0052            return result;
0053         }
0054         void refill_cache()
0055         {
0056            typedef typename lanczos::lanczos<T, Policy>::type lanczos_type;
0057
0058            gamma_cache[cache_size - 1] = boost::math::gamma_p(alpha_poch + ((int)cache_size - 1), x, pol);
0059            for (int i = cache_size - 1; i > 0; --i)
0060            {
0061               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)));
0062            }
0063         }
0064         T delta_poch, alpha_poch, x, term;
0065         T gamma_cache[cache_size];
0066         int k;
0067         long long log_scaling;
0068         int cache_offset;
0069         Policy pol;
0070      };
0071
0072      template <class T, class Policy>
0073      T hypergeometric_1F1_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, long long& log_scaling)
0074      {
0075         BOOST_MATH_STD_USING
0076         if (b_minus_a == 0)
0077         {
0078            // special case: M(a,a,z) == exp(z)
0079            long long scale = lltrunc(x, pol);
0080            log_scaling += scale;
0081            return exp(x - scale);
0082         }
0083         hypergeometric_1F1_igamma_series<T, Policy> s(b_minus_a, a - 1, x, pol);
0084         log_scaling += s.log_scaling;
0085         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0086         T result = boost::math::tools::sum_series(s, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0087         boost::math::policies::check_series_iterations<T>("boost::math::tgamma<%1%>(%1%,%1%)", max_iter, pol);
0088         T log_prefix = x + boost::math::lgamma(b, pol) - boost::math::lgamma(a, pol);
0089         long long scale = lltrunc(log_prefix);
0090         log_scaling += scale;
0091         return result * exp(log_prefix - scale);
0092      }
0093
0094      template <class T, class Policy>
0095      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, long long& log_scaling)
0096      {
0097         BOOST_MATH_STD_USING
0098         T a = a_local + a_shift;
0099         if (a_shift == 0)
0100            return h;
0101         else if (a_shift > 0)
0102         {
0103            //
0104            // Forward recursion on a is stable as long as 2a-b+z > 0.
0105            // If 2a-b+z < 0 then backwards recursion is stable even though
0106            // the function may be strictly increasing with a.  Potentially
0107            // we may need to split the recurrence in 2 sections - one using
0108            // forward recursion, and one backwards.
0109            //
0110            // We will get the next seed value from the ratio
0111            // on b as that's stable and quick to compute.
0112            //
0113
0114            T crossover_a = (b_local - x) / 2;
0115            int crossover_shift = itrunc(crossover_a - a_local);
0116
0117            if (crossover_shift > 1)
0118            {
0119               //
0120               // Forwards recursion will start off unstable, but may switch to the stable direction later.
0121               // Start in the middle and go in both directions:
0122               //
0123               if (crossover_shift > a_shift)
0124                  crossover_shift = a_shift;
0125               crossover_a = a_local + crossover_shift;
0126               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(crossover_a, b_local, x);
0127               std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0128               T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0129               boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
0130               //
0131               // Convert to a ratio:
0132               //         (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0
0133               //
0134               //  hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio
0135               //
0136               T first = 1;
0137               T second = ((1 + crossover_a - b_local) / crossover_a) + ((b_local - 1) / crossover_a) / b_ratio;
0138               //
0139               // Recurse down to a_local, compare values and re-normalise first and second:
0140               //
0141               boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(crossover_a, b_local, x);
0142               long long backwards_scale = 0;
0143               T comparitor = boost::math::tools::apply_recurrence_relation_backward(a_coef, crossover_shift, second, first, &backwards_scale);
0144               log_scaling -= backwards_scale;
0145               if ((h < 1) && (tools::max_value<T>() * h > comparitor))
0146               {
0147                  // Need to rescale!
0148                  long long scale = lltrunc(log(h), pol) + 1;
0149                  h *= exp(T(-scale));
0150                  log_scaling += scale;
0151               }
0152               comparitor /= h;
0153               first /= comparitor;
0154               second /= comparitor;
0155               //
0156               // Now we can recurse forwards for the rest of the range:
0157               //
0158               if (crossover_shift < a_shift)
0159               {
0160                  boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef_2(crossover_a + 1, b_local, x);
0161                  h = boost::math::tools::apply_recurrence_relation_forward(a_coef_2, a_shift - crossover_shift - 1, first, second, &log_scaling);
0162               }
0163               else
0164                  h = first;
0165            }
0166            else
0167            {
0168               //
0169               // Regular case where forwards iteration is stable right from the start:
0170               //
0171               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a_local, b_local, x);
0172               std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0173               T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0174               boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
0175               //
0176               // Convert to a ratio:
0177               //         (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0
0178               //
0179               //  hence: M(a+1,b,z) = ((1+a-b) / a) M(a,b,z) + ((b-1) / a) M(a,b,z)/b_ratio
0180               //
0181               T second = ((1 + a_local - b_local) / a_local) * h + ((b_local - 1) / a_local) * h / b_ratio;
0182               boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a_local + 1, b_local, x);
0183               h = boost::math::tools::apply_recurrence_relation_forward(a_coef, --a_shift, h, second, &log_scaling);
0184            }
0185         }
0186         else
0187         {
0188            //
0189            // We've calculated h for a larger value of a than we want, and need to recurse down.
0190            // However, only forward iteration is stable, so calculate the ratio, compare values,
0191            // and normalise.  Note that we calculate the ratio on b and convert to a since the
0192            // direction is the minimal solution for N->+INF.
0193            //
0194            // IMPORTANT: this is only currently called for a > b and therefore forwards iteration
0195            // is the only stable direction as we will only iterate down until a ~ b, but we
0196            // will check this with an assert:
0197            //
0198            BOOST_MATH_ASSERT(2 * a - b_local + x > 0);
0199            boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x);
0200            std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0201            T b_ratio = boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0202            boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
0203            //
0204            // Convert to a ratio:
0205            //         (1+a-b)M(a, b, z) - aM(a+1, b, z) + (b-1)M(a, b-1, z) = 0
0206            //
0207            //  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))
0208            //
0209            T first = 1;  // arbitrary value;
0210            T second = ((1 + a - b_local) / a) + ((b_local - 1) / a) * (1 / b_ratio);
0211
0212            if (a_shift == -1)
0213               h = h / second;
0214            else
0215            {
0216               boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> a_coef(a + 1, b_local, x);
0217               T comparitor = boost::math::tools::apply_recurrence_relation_forward(a_coef, -(a_shift + 1), first, second);
0218               if (boost::math::tools::min_value<T>() * comparitor > h)
0219               {
0220                  // Ooops, need to rescale h:
0221                  long long rescale = lltrunc(log(fabs(h)));
0222                  T scale = exp(T(-rescale));
0223                  h *= scale;
0224                  log_scaling += rescale;
0225               }
0226               h = h / comparitor;
0227            }
0228         }
0229         return h;
0230      }
0231
0232      template <class T, class Policy>
0233      T hypergeometric_1F1_shift_on_b(T h, const T& a, const T& b_local, const T& x, int b_shift, const Policy& pol, long long& log_scaling)
0234      {
0235         BOOST_MATH_STD_USING
0236
0237         T b = b_local + b_shift;
0238         if (b_shift == 0)
0239            return h;
0240         else if (b_shift > 0)
0241         {
0242            //
0243            // We get here for b_shift > 0 when b > z.  We can't use forward recursion on b - it's unstable,
0244            // so grab the ratio and work backwards to b - b_shift and normalise.
0245            //
0246            boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b, x);
0247            std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0248
0249            T first = 1;  // arbitrary value;
0250            T second = 1 / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0251            boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
0252            if (b_shift == 1)
0253               h = h / second;
0254            else
0255            {
0256               //
0257               // Reset coefficients and recurse:
0258               //
0259               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b - 1, x);
0260               long long local_scale = 0;
0261               T comparitor = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, --b_shift, first, second, &local_scale);
0262               log_scaling -= local_scale;
0263               if (boost::math::tools::min_value<T>() * comparitor > h)
0264               {
0265                  // Ooops, need to rescale h:
0266                  long long rescale = lltrunc(log(fabs(h)));
0267                  T scale = exp(T(-rescale));
0268                  h *= scale;
0269                  log_scaling += rescale;
0270               }
0271               h = h / comparitor;
0272            }
0273         }
0274         else
0275         {
0276            T second;
0277            if (a == b_local)
0278            {
0279                // recurrence is trivial for a == b and method of ratios fails as the c-term goes to zero:
0280               second = -b_local * (1 - b_local - x) * h / (b_local * (b_local - 1));
0281            }
0282            else
0283            {
0284               BOOST_MATH_ASSERT(!is_negative_integer(b - a));
0285               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef(a, b_local, x);
0286               std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0287               second = h / boost::math::tools::function_ratio_from_backwards_recurrence(b_coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0288               boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_large_abz<%1%>(%1%,%1%,%1%)", max_iter, pol);
0289            }
0290            if (b_shift == -1)
0291               h = second;
0292            else
0293            {
0294               boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> b_coef_2(a, b_local - 1, x);
0295               h = boost::math::tools::apply_recurrence_relation_backward(b_coef_2, -(++b_shift), h, second, &log_scaling);
0296            }
0297         }
0298         return h;
0299      }
0300
0301
0302      template <class T, class Policy>
0303      T hypergeometric_1F1_large_igamma(const T& a, const T& b, const T& x, const T& b_minus_a, const Policy& pol, long long& log_scaling)
0304      {
0305         BOOST_MATH_STD_USING
0306         //
0307         // We need a < b < z in order to ensure there's at least a chance of convergence,
0308         // we can use recurrence relations to shift forwards on a+b or just a to achieve this,
0309         // for decent accuracy, try to keep 2b - 1 < a < 2b < z
0310         //
0311         int b_shift = b * 2 < x ? 0 : itrunc(b - x / 2);
0312         int a_shift = a > b - b_shift ? -itrunc(b - b_shift - a - 1) : -itrunc(b - b_shift - a);
0313
0314         if (a_shift < 0)
0315         {
0316            // Might as well do all the shifting on b as scale a downwards:
0317            b_shift -= a_shift;
0318            a_shift = 0;
0319         }
0320
0321         T a_local = a - a_shift;
0322         T b_local = b - b_shift;
0323         T b_minus_a_local = (a_shift == 0) && (b_shift == 0) ? b_minus_a : b_local - a_local;
0324         long long local_scaling = 0;
0325         T h = hypergeometric_1F1_igamma(a_local, b_local, x, b_minus_a_local, pol, local_scaling);
0326         log_scaling += local_scaling;
0327
0328         //
0329         // Apply shifts on a and b as required:
0330         //
0331         h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, x, a_shift, pol, log_scaling);
0332         h = hypergeometric_1F1_shift_on_b(h, a, b_local, x, b_shift, pol, log_scaling);
0333
0334         return h;
0335      }
0336
0337      template <class T, class Policy>
0338      T hypergeometric_1F1_large_series(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0339      {
0340         BOOST_MATH_STD_USING
0341         //
0342         // We make a small, and b > z:
0343         //
0344         int a_shift(0), b_shift(0);
0345         if (a * z > b)
0346         {
0347            a_shift = itrunc(a) - 5;
0348            b_shift = b < z ? itrunc(b - z - 1) : 0;
0349         }
0350         //
0351         // If a_shift is trivially small, there's really not much point in losing
0352         // accuracy to save a couple of iterations:
0353         //
0354         if (a_shift < 5)
0355            a_shift = 0;
0356         T a_local = a - a_shift;
0357         T b_local = b - b_shift;
0358         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)");
0359         //
0360         // Apply shifts on a and b as required:
0361         //
0362         if (a_shift && (a_local == 0))
0363         {
0364            //
0365            // Shifting on a via method of ratios in hypergeometric_1F1_shift_on_a fails when
0366            // a_local == 0.  However, the value of h calculated was trivial (unity), so
0367            // calculate a second 1F1 for a == 1 and recurse as normal:
0368            //
0369            long long scale = 0;
0370            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)");
0371            if (scale != log_scaling)
0372            {
0373               h2 *= exp(T(scale - log_scaling));
0374            }
0375            boost::math::detail::hypergeometric_1F1_recurrence_a_coefficients<T> coef(a_local + 1, b_local, z);
0376            h = boost::math::tools::apply_recurrence_relation_forward(coef, a_shift - 1, h, h2, &log_scaling);
0377            h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling);
0378         }
0379         else
0380         {
0381            h = hypergeometric_1F1_shift_on_a(h, a_local, b_local, z, a_shift, pol, log_scaling);
0382            h = hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling);
0383         }
0384         return h;
0385      }
0386
0387      template <class T, class Policy>
0388      T hypergeometric_1F1_large_13_3_6_series(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0389      {
0390         BOOST_MATH_STD_USING
0391         //
0392         // A&S 13.3.6 is good only when a ~ b, but isn't too fussy on the size of z.
0393         // So shift b to match a (b shifting seems to be more stable via method of ratios).
0394         //
0395         int b_shift = itrunc(b - a);
0396         if ((b_shift < 0) && (b - b_shift != a))
0397            b_shift -= 1;
0398         T b_local = b - b_shift;
0399         if ((b_local - a - 0.5 <= 0) && (b_local != a))
0400         {
0401            // Make sure b_local - a - 0.5 > 0
0402            b_shift -= 1;
0403            b_local += 1;
0404         }
0405         T h = boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b_local, z, T(b_local - a), pol, log_scaling);
0406         return hypergeometric_1F1_shift_on_b(h, a, b_local, z, b_shift, pol, log_scaling);
0407      }
0408
0409      template <class T, class Policy>
0410      T hypergeometric_1F1_large_abz(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0411      {
0412         BOOST_MATH_STD_USING
0413         //
0414         // This is the selection logic to pick the "best" method.
0415         // We have a,b,z >> 0 and need to compute the approximate cost of each method
0416         // and then select whichever wins out.
0417         //
0418         enum method
0419         {
0420            method_series = 0,
0421            method_shifted_series,
0422            method_gamma,
0423            method_bessel
0424         };
0425         //
0426         // Cost of direct series, is the approx number of terms required until we hit convergence:
0427         //
0428         T current_cost = (sqrt(16 * z * (3 * a + z) + 9 * b * b - 24 * b * z) - 3 * b + 4 * z) / 6;
0429         method current_method = method_series;
0430         //
0431         // Cost of shifted series, is the number of recurrences required to move to a zone where
0432         // the series is convergent right from the start.
0433         // Note that recurrence relations fail for very small b, and too many recurrences on a
0434         // will completely destroy all our digits.
0435         // Also note that the method fails when b-a is a negative integer unless b is already
0436         // larger than z and thus does not need shifting.
0437         //
0438         T cost = a + ((b < z) ? T(z - b) : T(0));
0439         if((b > 1) && (cost < current_cost) && ((b > z) || !is_negative_integer(b-a)))
0440         {
0441            current_method = method_shifted_series;
0442            current_cost = cost;
0443         }
0444         //
0445         // Cost for gamma function method is the number of recurrences required to move it
0446         // into a convergent zone, note that recurrence relations fail for very small b.
0447         // Also add on a fudge factor to account for the fact that this method is both
0448         // more expensive to compute (requires gamma functions), and less accurate than the
0449         // methods above:
0450         //
0451         T b_shift = fabs(b * 2 < z ? T(0) : T(b - z / 2));
0452         T a_shift = fabs(a > b - b_shift ? T(-(b - b_shift - a - 1)) : T(-(b - b_shift - a)));
0453         cost = 1000 + b_shift + a_shift;
0454         if((b > 1) && (cost <= current_cost))
0455         {
0456            current_method = method_gamma;
0457            current_cost = cost;
0458         }
0459         //
0460         // Cost for bessel approximation is the number of recurrences required to make a ~ b,
0461         // Note that recurrence relations fail for very small b.  We also have issue with large
0462         // z: either overflow/numeric instability or else the series goes divergent.  We seem to be
0463         // OK for z smaller than log_max_value<Quad> though, maybe we can stretch a little further
0464         // but that's not clear...
0465         // Also need to add on a fudge factor to the cost to account for the fact that we need
0466         // to calculate the Bessel functions, this is not quite as high as the gamma function
0467         // method above as this is generally more accurate and so preferred if the methods are close:
0468         //
0469         cost = 50 + fabs(b - a);
0470         if((b > 1) && (cost <= current_cost) && (z < tools::log_max_value<T>()) && (z < 11356) && (b - a != 0.5f))
0471         {
0472            current_method = method_bessel;
0473            current_cost = cost;
0474         }
0475
0476         switch (current_method)
0477         {
0478         case method_series:
0479            return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, "hypergeometric_1f1_large_abz<%1%>(a,b,z)");
0480         case method_shifted_series:
0481            return detail::hypergeometric_1F1_large_series(a, b, z, pol, log_scaling);
0482         case method_gamma:
0483            return detail::hypergeometric_1F1_large_igamma(a, b, z, T(b - a), pol, log_scaling);
0484         case method_bessel:
0485            return detail::hypergeometric_1F1_large_13_3_6_series(a, b, z, pol, log_scaling);
0486         }
0487         return 0; // We don't get here!
0488      }
0489
0490   } } } // namespaces
0491
0492 #endif // BOOST_HYPERGEOMETRIC_1F1_LARGE_ABZ_HPP_