math/special_functions/hypergeometric_1F1.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_HPP
0011 #define BOOST_MATH_HYPERGEOMETRIC_1F1_HPP
0012
0013 #include <boost/math/tools/config.hpp>
0014 #include <boost/math/policies/policy.hpp>
0015 #include <boost/math/policies/error_handling.hpp>
0016 #include <boost/math/special_functions/detail/hypergeometric_series.hpp>
0017 #include <boost/math/special_functions/detail/hypergeometric_asym.hpp>
0018 #include <boost/math/special_functions/detail/hypergeometric_rational.hpp>
0019 #include <boost/math/special_functions/detail/hypergeometric_1F1_recurrence.hpp>
0020 #include <boost/math/special_functions/detail/hypergeometric_1F1_by_ratios.hpp>
0021 #include <boost/math/special_functions/detail/hypergeometric_pade.hpp>
0022 #include <boost/math/special_functions/detail/hypergeometric_1F1_bessel.hpp>
0023 #include <boost/math/special_functions/detail/hypergeometric_1F1_scaled_series.hpp>
0024 #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
0025 #include <boost/math/special_functions/detail/hypergeometric_1F1_addition_theorems_on_z.hpp>
0026 #include <boost/math/special_functions/detail/hypergeometric_1F1_large_abz.hpp>
0027 #include <boost/math/special_functions/detail/hypergeometric_1F1_small_a_negative_b_by_ratio.hpp>
0028 #include <boost/math/special_functions/detail/hypergeometric_1F1_negative_b_regions.hpp>
0029
0030 namespace boost { namespace math { namespace detail {
0031
0032    // check when 1F1 series can't decay to polynom
0033    template <class T>
0034    inline bool check_hypergeometric_1F1_parameters(const T& a, const T& b)
0035    {
0036       BOOST_MATH_STD_USING
0037
0038          if ((b <= 0) && (b == floor(b)))
0039          {
0040             if ((a >= 0) || (a < b) || (a != floor(a)))
0041                return false;
0042          }
0043
0044       return true;
0045    }
0046
0047    template <class T, class Policy>
0048    T hypergeometric_1F1_divergent_fallback(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0049    {
0050       BOOST_MATH_STD_USING
0051       const char* function = "hypergeometric_1F1_divergent_fallback<%1%>(%1%,%1%,%1%)";
0052       //
0053       // We get here if either:
0054       // 1) We decide up front that Tricomi's method won't work, or:
0055       // 2) We've called Tricomi's method and it's failed.
0056       //
0057       if (b > 0)
0058       {
0059          // Commented out since recurrence seems to always be better?
0060 #if 0
0061          if ((z < b) && (a > -50))
0062             // Might as well use a recurrence in preference to z-recurrence:
0063             return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);
0064          T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));
0065          int k = 1 + itrunc(z - z_limit);
0066          // If k is too large we destroy all the digits in the result:
0067          T convergence_at_50 = (b - a + 50) * k / (z * 50);
0068          if ((k > 0) && (k < 50) && (fabs(convergence_at_50) < 1) && (z > z_limit))
0069          {
0070             return boost::math::detail::hypergeometric_1f1_recurrence_on_z_minus_zero(a, b, T(z - k), k, pol, log_scaling);
0071          }
0072 #endif
0073          if (z < b)
0074             return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);
0075          else
0076             return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);
0077       }
0078       else  // b < 0
0079       {
0080          if (a < 0)
0081          {
0082             if ((b < a) && (z < -b / 4))
0083                return hypergeometric_1F1_from_function_ratio_negative_ab(a, b, z, pol, log_scaling);
0084             else
0085             {
0086                //
0087                // Solve (a+n)z/((b+n)n) == 1 for n, the number of iterations till the series starts to converge.
0088                // If this is well away from the origin then it's probably better to use the series to evaluate this.
0089                // Note that if sqr is negative then we have no solution, so assign an arbitrarily large value to the
0090                // number of iterations.
0091                //
0092                bool can_use_recursion = (z - b + 100 < boost::math::policies::get_max_series_iterations<Policy>()) && (100 - a < boost::math::policies::get_max_series_iterations<Policy>());
0093                T sqr = 4 * a * z + b * b - 2 * b * z + z * z;
0094                T iterations_to_convergence = sqr > 0 ? T(0.5f * (-sqrt(sqr) - b + z)) : T(-a - b);
0095                if(can_use_recursion && ((std::max)(a, b) + iterations_to_convergence > -300))
0096                   return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);
0097                //
0098                // When a < b and if we fall through to the series, then we get divergent behaviour when b crosses the origin
0099                // so ideally we would pick another method.  Otherwise the terms immediately after b crosses the origin may
0100                // suffer catastrophic cancellation....
0101                //
0102                if((a < b) && can_use_recursion)
0103                   return hypergeometric_1F1_backwards_recursion_on_b_for_negative_a(a, b, z, pol, function, log_scaling);
0104             }
0105          }
0106          else
0107          {
0108             //
0109             // Start by getting the domain of the recurrence relations, we get either:
0110             //   -1     Backwards recursion is stable and the CF will converge to double precision.
0111             //   +1     Forwards recursion is stable and the CF will converge to double precision.
0112             //    0     No man's land, we're not far enough away from the crossover point to get double precision from either CF.
0113             //
0114             // At higher than double precision we need to be further away from the crossover location to
0115             // get full converge, but it's not clear how much further - indeed at quad precision it's
0116             // basically impossible to ever get forwards iteration to work.  Backwards seems to work
0117             // OK as long as a > 1 whatever the precision though.
0118             //
0119             int domain = hypergeometric_1F1_negative_b_recurrence_region(a, b, z);
0120             if ((domain < 0) && ((a > 1) || (boost::math::policies::digits<T, Policy>() <= 64)))
0121                return hypergeometric_1F1_from_function_ratio_negative_b(a, b, z, pol, log_scaling);
0122             else if (domain > 0)
0123             {
0124                if (boost::math::policies::digits<T, Policy>() <= 64)
0125                   return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling);
0126 #ifndef BOOST_NO_EXCEPTIONS
0127                try
0128 #endif
0129                {
0130                   return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
0131                }
0132 #ifndef BOOST_NO_EXCEPTIONS
0133                catch (const evaluation_error&)
0134                {
0135                   //
0136                   // The series failed, try the recursions instead and hope we get at least double precision:
0137                   //
0138                   return hypergeometric_1F1_from_function_ratio_negative_b_forwards(a, b, z, pol, log_scaling);
0139                }
0140 #endif
0141             }
0142             //
0143             // We could fall back to Tricomi's approximation if we're in the transition zone
0144             // between the above two regions.  However, I've been unable to find any examples
0145             // where this is better than the series, and there are many cases where it leads to
0146             // quite grievous errors.
0147             /*
0148             else if (allow_tricomi)
0149             {
0150                T aa = a < 1 ? T(1) : a;
0151                if (z < fabs((2 * aa - b) / (sqrt(fabs(aa * b)))))
0152                   return hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
0153             }
0154             */
0155          }
0156       }
0157
0158       // If we get here, then we've run out of methods to try, use the checked series which will
0159       // raise an error if the result is garbage:
0160       return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
0161    }
0162
0163    template <class T>
0164    bool is_convergent_negative_z_series(const T& a, const T& b, const T& z, const T& b_minus_a)
0165    {
0166       BOOST_MATH_STD_USING
0167       //
0168       // Filter out some cases we don't want first:
0169       //
0170       if((b_minus_a > 0) && (b > 0))
0171       {
0172          if (a < 0)
0173             return false;
0174       }
0175       //
0176       // Generic check: we have small initial divergence and are convergent after 10 terms:
0177       //
0178       if ((fabs(z * a / b) < 2) && (fabs(z * (a + 10) / ((b + 10) * 10)) < 1))
0179       {
0180          // Double check for divergence when we cross the origin on a and b:
0181          if (a < 0)
0182          {
0183             T n = 300 - floor(a);
0184             if (fabs((a + n) * z / ((b + n) * n)) < 1)
0185             {
0186                if (b < 0)
0187                {
0188                   T m = 3 - floor(b);
0189                   if (fabs((a + m) * z / ((b + m) * m)) < 1)
0190                      return true;
0191                }
0192                else
0193                   return true;
0194             }
0195          }
0196          else if (b < 0)
0197          {
0198             T n = 3 - floor(b);
0199             if (fabs((a + n) * z / ((b + n) * n)) < 1)
0200                return true;
0201          }
0202       }
0203       if ((b > 0) && (a < 0))
0204       {
0205          //
0206          // For a and z both negative, we're OK with some initial divergence as long as
0207          // it occurs before we hit the origin, as to start with all the terms have the
0208          // same sign.
0209          //
0210          // https://www.wolframalpha.com/input/?i=solve+(a%2Bn)z+%2F+((b%2Bn)n)+%3D%3D+1+for+n
0211          //
0212          T sqr = 4 * a * z + b * b - 2 * b * z + z * z;
0213          T iterations_to_convergence = sqr > 0 ? T(0.5f * (-sqrt(sqr) - b + z)) : T(-a + b);
0214          if (iterations_to_convergence < 0)
0215             iterations_to_convergence = 0.5f * (sqrt(sqr) - b + z);
0216          if (a + iterations_to_convergence < -50)
0217          {
0218             // Need to check for divergence when we cross the origin on a:
0219             if (a > -1)
0220                return true;
0221             T n = 300 - floor(a);
0222             if(fabs((a + n) * z / ((b + n) * n)) < 1)
0223                return true;
0224          }
0225       }
0226       return false;
0227    }
0228
0229    template <class T>
0230    inline T cyl_bessel_i_shrinkage_rate(const T& z)
0231    {
0232       // Approximately the ratio I_10.5(z/2) / I_9.5(z/2), this gives us an idea of how quickly
0233       // the Bessel terms in A&S 13.6.4 are converging:
0234       if (z < -160)
0235          return 1;
0236       if (z < -40)
0237          return 0.75f;
0238       if (z < -20)
0239          return 0.5f;
0240       if (z < -7)
0241          return 0.25f;
0242       if (z < -2)
0243          return 0.1f;
0244       return 0.05f;
0245    }
0246
0247    template <class T>
0248    inline bool hypergeometric_1F1_is_13_3_6_region(const T& a, const T& b, const T& z)
0249    {
0250       BOOST_MATH_STD_USING
0251       if(fabs(a) == 0.5)
0252          return false;
0253       if ((z < 0) && (fabs(10 * a / b) < 1) && (fabs(a) < 50))
0254       {
0255          T shrinkage = cyl_bessel_i_shrinkage_rate(z);
0256          // We want the first term not too divergent, and convergence by term 10:
0257          if ((fabs((2 * a - 1) * (2 * a - b) / b) < 2) && (fabs(shrinkage * (2 * a + 9) * (2 * a - b + 10) / (10 * (b + 10))) < 0.75))
0258             return true;
0259       }
0260       return false;
0261    }
0262
0263    template <class T>
0264    inline bool hypergeometric_1F1_need_kummer_reflection(const T& a, const T& b, const T& z)
0265    {
0266       BOOST_MATH_STD_USING
0267       //
0268       // Check to see if we should apply Kummer's relation or not:
0269       //
0270       if (z > 0)
0271          return false;
0272       if (z < -1)
0273          return true;
0274       //
0275       // When z is small and negative, things get more complex.
0276       // More often than not we do not need apply Kummer's relation and the
0277       // series is convergent as is, but we do need to check:
0278       //
0279       if (a > 0)
0280       {
0281          if (b > 0)
0282          {
0283             return fabs((a + 10) * z / (10 * (b + 10))) < 1;  // Is the 10'th term convergent?
0284          }
0285          else
0286          {
0287             return true;  // Likely to be divergent as b crosses the origin
0288          }
0289       }
0290       else // a < 0
0291       {
0292          if (b > 0)
0293          {
0294             return false;  // Terms start off all positive and then by the time a crosses the origin we *must* be convergent.
0295          }
0296          else
0297          {
0298             return true;  // Likely to be divergent as b crosses the origin, but hard to rationalise about!
0299          }
0300       }
0301    }
0302
0303
0304    template <class T, class Policy>
0305    T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0306    {
0307       BOOST_MATH_STD_USING // exp, fabs, sqrt
0308
0309       static const char* const function = "boost::math::hypergeometric_1F1<%1%,%1%,%1%>(%1%,%1%,%1%)";
0310
0311       if ((z == 0) || (a == 0))
0312          return T(1);
0313
0314       // undefined result:
0315       if (!detail::check_hypergeometric_1F1_parameters(a, b))
0316          return policies::raise_domain_error<T>(
0317             function,
0318             "Function is indeterminate for negative integer b = %1%.",
0319             b,
0320             pol);
0321
0322       // other checks:
0323       if (a == -1)
0324       {
0325          T r = 1 - (z / b);
0326          if (fabs(r) < 0.5)
0327             r = (b - z) / b;
0328          return r;
0329       }
0330
0331       const T b_minus_a = b - a;
0332
0333       // 0f0 a == b case;
0334       if (b_minus_a == 0)
0335       {
0336          if ((a < 0) && (floor(a) == a))
0337          {
0338             // Special case, use the truncated series to match what Mathematica does.
0339             if ((a < -20) && (z > 0) && (z < 1))
0340             {
0341                // https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/03/01/04/02/0002/
0342                return exp(z) * boost::math::gamma_q(1 - a, z, pol);
0343             }
0344             // https://functions.wolfram.com/HypergeometricFunctions/Hypergeometric1F1/03/01/04/02/0003/
0345             return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
0346          }
0347          long long scale = lltrunc(z, pol);
0348          log_scaling += scale;
0349          return exp(z - scale);
0350       }
0351       // Special case for b-a = -1, we don't use for small a as it throws the digits of a away and leads to large errors:
0352       if ((b_minus_a == -1) && (fabs(a) > 0.5))
0353       {
0354          // for negative small integer a it is reasonable to use truncated series - polynomial
0355          if ((a < 0) && (a == ceil(a)) && (a > -50))
0356             return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);
0357
0358          log_scaling = lltrunc(floor(z));
0359          T local_z = z - log_scaling;
0360          return (b + z) * exp(local_z) / b;
0361       }
0362
0363       if ((a == 1) && (b == 2))
0364          return boost::math::expm1(z, pol) / z;
0365
0366       if ((b - a == b) && (fabs(z / b) < policies::get_epsilon<T, Policy>()))
0367          return 1;
0368       //
0369       // Special case for A&S 13.3.6:
0370       //
0371       if (z < 0)
0372       {
0373          if (hypergeometric_1F1_is_13_3_6_region(a, b, z))
0374          {
0375             // a is tiny compared to b, and z < 0
0376             // 13.3.6 appears to be the most efficient and often the most accurate method.
0377             T r = boost::math::detail::hypergeometric_1F1_AS_13_3_6(b_minus_a, b, T(-z), a, pol, log_scaling);
0378             long long scale = lltrunc(z, pol);
0379             log_scaling += scale;
0380             return r * exp(z - scale);
0381          }
0382          if ((b < 0) && (fabs(a) < 1e-2))
0383          {
0384             //
0385             // This is a tricky area, potentially we have no good method at all:
0386             //
0387             if (b - ceil(b) == a)
0388             {
0389                // Fractional parts of a and b are genuinely equal, we might as well
0390                // apply Kummer's relation and get a truncated series:
0391                long long scaling = lltrunc(z);
0392                T r = exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);
0393                log_scaling += scaling;
0394                return r;
0395             }
0396             if ((b < -1) && (max_b_for_1F1_small_a_negative_b_by_ratio(z) < b))
0397                return hypergeometric_1F1_small_a_negative_b_by_ratio(a, b, z, pol, log_scaling);
0398             if ((b > -1) && (b < -0.5f))
0399             {
0400                // Recursion is meta-stable:
0401                T first = hypergeometric_1F1_imp(a, T(b + 2), z, pol);
0402                T second = hypergeometric_1F1_imp(a, T(b + 1), z, pol);
0403                return tools::apply_recurrence_relation_backward(hypergeometric_1F1_recurrence_small_b_coefficients<T>(a, b, z, 1), 1, first, second);
0404             }
0405             //
0406             // We've got nothing left but 13.3.6, even though it may be initially divergent:
0407             //
0408             T r = boost::math::detail::hypergeometric_1F1_AS_13_3_6(b_minus_a, b, T(-z), a, pol, log_scaling);
0409             long long scale = lltrunc(z, pol);
0410             log_scaling += scale;
0411             return r * exp(z - scale);
0412          }
0413       }
0414       //
0415       // Asymptotic expansion for large z
0416       // TODO: check region for higher precision types.
0417       // Use recurrence relations to move to this region when a and b are also large.
0418       //
0419       if (detail::hypergeometric_1F1_asym_region(a, b, z, pol))
0420       {
0421          long long saved_scale = log_scaling;
0422 #ifndef BOOST_NO_EXCEPTIONS
0423          try
0424 #endif
0425          {
0426             return hypergeometric_1F1_asym_large_z_series(a, b, z, pol, log_scaling);
0427          }
0428 #ifndef BOOST_NO_EXCEPTIONS
0429          catch (const evaluation_error&)
0430          {
0431          }
0432 #endif
0433          //
0434          // Very occasionally our convergence criteria don't quite go to full precision
0435          // and we have to try another method:
0436          //
0437          log_scaling = saved_scale;
0438       }
0439
0440       if ((fabs(a * z / b) < 3.5) && (fabs(z * 100) < fabs(b)) && ((fabs(a) > 1e-2) || (b < -5)))
0441          return detail::hypergeometric_1F1_rational(a, b, z, pol);
0442
0443       if (hypergeometric_1F1_need_kummer_reflection(a, b, z))
0444       {
0445          if (a == 1)
0446             return detail::hypergeometric_1F1_pade(b, z, pol);
0447          if (is_convergent_negative_z_series(a, b, z, b_minus_a))
0448          {
0449             if ((boost::math::sign(b_minus_a) == boost::math::sign(b)) && ((b > 0) || (b < -200)))
0450             {
0451                // Series is close enough to convergent that we should be OK,
0452                // In this domain b - a ~ b and since 1F1[a, a, z] = e^z 1F1[b-a, b, -z]
0453                // and 1F1[a, a, -z] = e^-z the result must necessarily be somewhere near unity.
0454                // We have to rule out b small and negative because if b crosses the origin early
0455                // in the series (before we're pretty much converged) then all bets are off.
0456                // Note that this can go badly wrong when b and z are both large and negative,
0457                // in that situation the series goes in waves of large and small values which
0458                // may or may not cancel out.  Likewise the initial part of the series may or may
0459                // not converge, and even if it does may or may not give a correct answer!
0460                // For example 1F1[-small, -1252.5, -1043.7] can loose up to ~800 digits due to
0461                // cancellation and is basically incalculable via this method.
0462                return hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
0463             }
0464          }
0465          if ((b < 0) && (floor(b) == b))
0466          {
0467             // Negative integer b, so a must be a negative integer too.
0468             // Kummer's transformation fails here!
0469             if(a > -50)
0470                return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);
0471             // Is there anything better than this??
0472             return hypergeometric_1F1_imp(a, float_next(b), z, pol, log_scaling);
0473          }
0474          else
0475          {
0476             // Let's otherwise make z positive (almost always)
0477             // by Kummer's transformation
0478             // (we also don't transform if z belongs to [-1,0])
0479             // Also note that Kummer's transformation fails when b is
0480             // a negative integer, although this seems to be unmentioned
0481             // in the literature...
0482             long long scaling = lltrunc(z);
0483             T r = exp(z - scaling) * detail::hypergeometric_1F1_imp<T>(b_minus_a, b, -z, pol, log_scaling);
0484             log_scaling += scaling;
0485             return r;
0486          }
0487       }
0488       //
0489       // Check for initial divergence:
0490       //
0491       bool series_is_divergent = (a + 1) * z / (b + 1) < -1;
0492       if (series_is_divergent && (a < 0) && (b < 0) && (a > -1))
0493          series_is_divergent = false;   // Best off taking the series in this situation
0494       //
0495       // If series starts off non-divergent, and becomes divergent later
0496       // then it's because both a and b are negative, so check for later
0497       // divergence as well:
0498       //
0499       if (!series_is_divergent && (a < 0) && (b < 0) && (b > a))
0500       {
0501          //
0502          // We need to exclude situations where we're over the initial "hump"
0503          // in the series terms (ie series has already converged by the time
0504          // b crosses the origin:
0505          //
0506          //T fa = fabs(a);
0507          //T fb = fabs(b);
0508          T convergence_point = sqrt((a - 1) * (a - b)) - a;
0509          if (-b < convergence_point)
0510          {
0511             T n = -floor(b);
0512             series_is_divergent = (a + n) * z / ((b + n) * n) < -1;
0513          }
0514       }
0515       else if (!series_is_divergent && (b < 0) && (a > 0))
0516       {
0517          // Series almost always become divergent as b crosses the origin:
0518          series_is_divergent = true;
0519       }
0520       if (series_is_divergent && (b < -1) && (b > -5) && (a > b))
0521          series_is_divergent = false;  // don't bother with divergence, series will be OK
0522
0523       //
0524       // Test for alternating series due to negative a,
0525       // in particular, see if the series is initially divergent
0526       // If so use the recurrence relation on a:
0527       //
0528       if (series_is_divergent)
0529       {
0530          if((a < 0) && (floor(a) == a) && (-a < policies::get_max_series_iterations<Policy>()))
0531             // This works amazingly well for negative integer a:
0532             return hypergeometric_1F1_backward_recurrence_for_negative_a(a, b, z, pol, function, log_scaling);
0533          //
0534          // In what follows we have to set limits on how large z can be otherwise
0535          // the Bessel series become large and divergent and all the digits cancel out.
0536          // The criteria are distinctly empiracle rather than based on a firm analysis
0537          // of the terms in the series.
0538          //
0539          if (b > 0)
0540          {
0541             T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));
0542             if ((z < z_limit) && hypergeometric_1F1_is_tricomi_viable_positive_b(a, b, z))
0543                return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
0544          }
0545          else  // b < 0
0546          {
0547             if (a < 0)
0548             {
0549                T z_limit = fabs((2 * a - b) / (sqrt(fabs(a))));
0550                //
0551                // I hate these hard limits, but they're about the best we can do to try and avoid
0552                // Bessel function internal failures: these will be caught and handled
0553                // but up the expense of this function call:
0554                //
0555                if (((z < z_limit) || (a > -500)) && ((b > -500) || (b - 2 * a > 0)) && (z < -a))
0556                {
0557                   //
0558                   // Outside this domain we will probably get better accuracy from the recursive methods.
0559                   //
0560                   if(!(((a < b) && (z > -b)) || (z > z_limit)))
0561                      return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
0562                   //
0563                   // When b and z are both very small, we get large errors from the recurrence methods
0564                   // in the fallbacks.  Tricomi seems to work well here, as does direct series evaluation
0565                   // at least some of the time.  Picking the right method is not easy, and sometimes this
0566                   // is much worse than the fallback.  Overall though, it's a reasonable choice that keeps
0567                   // the very worst errors under control.
0568                   //
0569                   if(b > -1)
0570                      return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
0571                }
0572             }
0573             //
0574             // We previously used Tricomi here, but it appears to be worse than
0575             // the recurrence-based algorithms in hypergeometric_1F1_divergent_fallback.
0576             /*
0577             else
0578             {
0579                T aa = a < 1 ? T(1) : a;
0580                if (z < fabs((2 * aa - b) / (sqrt(fabs(aa * b)))))
0581                   return detail::hypergeometric_1F1_AS_13_3_7_tricomi(a, b, z, pol, log_scaling);
0582             }*/
0583          }
0584
0585          return hypergeometric_1F1_divergent_fallback(a, b, z, pol, log_scaling);
0586       }
0587
0588       if (hypergeometric_1F1_is_13_3_6_region(b_minus_a, b, T(-z)))
0589       {
0590          // b_minus_a is tiny compared to b, and -z < 0
0591          // 13.3.6 appears to be the most efficient and often the most accurate method.
0592          return boost::math::detail::hypergeometric_1F1_AS_13_3_6(a, b, z, b_minus_a, pol, log_scaling);
0593       }
0594 #if 0
0595       if ((a > 0) && (b > 0) && (a * z / b > 2))
0596       {
0597          //
0598          // Series is initially divergent and slow to converge, see if applying
0599          // Kummer's relation can improve things:
0600          //
0601          if (is_convergent_negative_z_series(b_minus_a, b, T(-z), b_minus_a))
0602          {
0603             long long scaling = lltrunc(z);
0604             T r = exp(z - scaling) * detail::hypergeometric_1F1_checked_series_impl(b_minus_a, b, T(-z), pol, log_scaling);
0605             log_scaling += scaling;
0606             return r;
0607          }
0608
0609       }
0610 #endif
0611       if ((a > 0) && (b > 0) && (a * z > 50))
0612          return detail::hypergeometric_1F1_large_abz(a, b, z, pol, log_scaling);
0613
0614       if (b < 0)
0615          return detail::hypergeometric_1F1_checked_series_impl(a, b, z, pol, log_scaling);
0616
0617       return detail::hypergeometric_1F1_generic_series(a, b, z, pol, log_scaling, function);
0618    }
0619
0620    template <class T, class Policy>
0621    inline T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol)
0622    {
0623       BOOST_MATH_STD_USING // exp, fabs, sqrt
0624       long long log_scaling = 0;
0625       T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);
0626       //
0627       // Actual result will be result * e^log_scaling.
0628       //
0629       static const thread_local long long max_scaling = lltrunc(boost::math::tools::log_max_value<T>()) - 2;
0630       static const thread_local T max_scale_factor = exp(T(max_scaling));
0631
0632       while (log_scaling > max_scaling)
0633       {
0634          result *= max_scale_factor;
0635          log_scaling -= max_scaling;
0636       }
0637       while (log_scaling < -max_scaling)
0638       {
0639          result /= max_scale_factor;
0640          log_scaling += max_scaling;
0641       }
0642       if (log_scaling)
0643          result *= exp(T(log_scaling));
0644       return result;
0645    }
0646
0647    template <class T, class Policy>
0648    inline T log_hypergeometric_1F1_imp(const T& a, const T& b, const T& z, int* sign, const Policy& pol)
0649    {
0650       BOOST_MATH_STD_USING // exp, fabs, sqrt
0651       long long log_scaling = 0;
0652       T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);
0653       if (sign)
0654       *sign = result < 0 ? -1 : 1;
0655      result = log(fabs(result)) + log_scaling;
0656       return result;
0657    }
0658
0659    template <class T, class Policy>
0660    inline T hypergeometric_1F1_regularized_imp(const T& a, const T& b, const T& z, const Policy& pol)
0661    {
0662       BOOST_MATH_STD_USING // exp, fabs, sqrt
0663       long long log_scaling = 0;
0664       T result = hypergeometric_1F1_imp(a, b, z, pol, log_scaling);
0665       //
0666       // Actual result will be result * e^log_scaling / tgamma(b).
0667       //
0668       int result_sign = 1;
0669       T scale = log_scaling - boost::math::lgamma(b, &result_sign, pol);
0670
0671       static const thread_local T max_scaling = boost::math::tools::log_max_value<T>() - 2;
0672       static const thread_local T max_scale_factor = exp(max_scaling);
0673
0674       while (scale > max_scaling)
0675       {
0676          result *= max_scale_factor;
0677          scale -= max_scaling;
0678       }
0679       while (scale < -max_scaling)
0680       {
0681          result /= max_scale_factor;
0682      scale += max_scaling;
0683       }
0684       if (scale != 0)
0685          result *= exp(scale);
0686       return result * result_sign;
0687    }
0688
0689 } // namespace detail
0690
0691 template <class T1, class T2, class T3, class Policy>
0692 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1(T1 a, T2 b, T3 z, const Policy& /* pol */)
0693 {
0694    BOOST_FPU_EXCEPTION_GUARD
0695       typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0696    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0697    typedef typename policies::normalise<
0698       Policy,
0699       policies::promote_float<false>,
0700       policies::promote_double<false>,
0701       policies::discrete_quantile<>,
0702       policies::assert_undefined<> >::type forwarding_policy;
0703    return policies::checked_narrowing_cast<result_type, Policy>(
0704       detail::hypergeometric_1F1_imp<value_type>(
0705          static_cast<value_type>(a),
0706          static_cast<value_type>(b),
0707          static_cast<value_type>(z),
0708          forwarding_policy()),
0709       "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
0710 }
0711
0712 template <class T1, class T2, class T3>
0713 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1(T1 a, T2 b, T3 z)
0714 {
0715    return hypergeometric_1F1(a, b, z, policies::policy<>());
0716 }
0717
0718 template <class T1, class T2, class T3, class Policy>
0719 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1_regularized(T1 a, T2 b, T3 z, const Policy& /* pol */)
0720 {
0721    BOOST_FPU_EXCEPTION_GUARD
0722       typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0723    typedef typename policies::evaluation<result_type, Policy>::type value_type;
0724    typedef typename policies::normalise<
0725       Policy,
0726       policies::promote_float<false>,
0727       policies::promote_double<false>,
0728       policies::discrete_quantile<>,
0729       policies::assert_undefined<> >::type forwarding_policy;
0730    return policies::checked_narrowing_cast<result_type, Policy>(
0731       detail::hypergeometric_1F1_regularized_imp<value_type>(
0732          static_cast<value_type>(a),
0733          static_cast<value_type>(b),
0734          static_cast<value_type>(z),
0735          forwarding_policy()),
0736       "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
0737 }
0738
0739 template <class T1, class T2, class T3>
0740 inline typename tools::promote_args<T1, T2, T3>::type hypergeometric_1F1_regularized(T1 a, T2 b, T3 z)
0741 {
0742    return hypergeometric_1F1_regularized(a, b, z, policies::policy<>());
0743 }
0744
0745 template <class T1, class T2, class T3, class Policy>
0746 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, const Policy& /* pol */)
0747 {
0748   BOOST_FPU_EXCEPTION_GUARD
0749     typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0750   typedef typename policies::evaluation<result_type, Policy>::type value_type;
0751   typedef typename policies::normalise<
0752     Policy,
0753     policies::promote_float<false>,
0754     policies::promote_double<false>,
0755     policies::discrete_quantile<>,
0756     policies::assert_undefined<> >::type forwarding_policy;
0757   return policies::checked_narrowing_cast<result_type, Policy>(
0758     detail::log_hypergeometric_1F1_imp<value_type>(
0759       static_cast<value_type>(a),
0760       static_cast<value_type>(b),
0761       static_cast<value_type>(z),
0762       0,
0763       forwarding_policy()),
0764     "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
0765 }
0766
0767 template <class T1, class T2, class T3>
0768 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z)
0769 {
0770   return log_hypergeometric_1F1(a, b, z, policies::policy<>());
0771 }
0772
0773 template <class T1, class T2, class T3, class Policy>
0774 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign, const Policy& /* pol */)
0775 {
0776   BOOST_FPU_EXCEPTION_GUARD
0777     typedef typename tools::promote_args<T1, T2, T3>::type result_type;
0778   typedef typename policies::evaluation<result_type, Policy>::type value_type;
0779   typedef typename policies::normalise<
0780     Policy,
0781     policies::promote_float<false>,
0782     policies::promote_double<false>,
0783     policies::discrete_quantile<>,
0784     policies::assert_undefined<> >::type forwarding_policy;
0785   return policies::checked_narrowing_cast<result_type, Policy>(
0786     detail::log_hypergeometric_1F1_imp<value_type>(
0787       static_cast<value_type>(a),
0788       static_cast<value_type>(b),
0789       static_cast<value_type>(z),
0790       sign,
0791       forwarding_policy()),
0792     "boost::math::hypergeometric_1F1<%1%>(%1%,%1%,%1%)");
0793 }
0794
0795 template <class T1, class T2, class T3>
0796 inline typename tools::promote_args<T1, T2, T3>::type log_hypergeometric_1F1(T1 a, T2 b, T3 z, int* sign)
0797 {
0798   return log_hypergeometric_1F1(a, b, z, sign, policies::policy<>());
0799 }
0800
0801
0802   } } // namespace boost::math
0803
0804 #endif // BOOST_MATH_HYPERGEOMETRIC_HPP