special_functions/detail/hypergeometric_1F1_by_ratios.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_BY_RATIOS_HPP_
0009 #define BOOST_HYPERGEOMETRIC_1F1_BY_RATIOS_HPP_
0010
0011 #include <boost/math/tools/recurrence.hpp>
0012 #include <boost/math/policies/error_handling.hpp>
0013
0014   namespace boost { namespace math { namespace detail {
0015
0016      template <class T, class Policy>
0017      T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling);
0018
0019      /*
0020       Evaluation by method of ratios for domain b < 0 < a,z
0021
0022       We first convert the recurrence relation into a ratio
0023       of M(a+1, b+1, z) / M(a, b, z) using Shintan's equivalence
0024       between solving a recurrence relation using Miller's method
0025       and continued fractions.  The continued fraction is VERY rapid
0026       to converge (typically < 10 terms), but may converge to entirely
0027       the wrong value if we're in a bad part of the domain.  Strangely
0028       it seems to matter not whether we use recurrence on a, b or a and b
0029       they all work or not work about the same, so we might as well make
0030       life easy for ourselves and use the a and b recurrence to avoid
0031       having to apply one extra recurrence to convert from an a or b
0032       recurrence to an a+b one.
0033
0034       See: H. Shintan, Note on Miller's recurrence algorithm, J. Sci. Hiroshima Univ. Ser. A-I Math., 29 (1965), pp. 121-133.
0035       Also: Computational Aspects of Three Term Recurrence Relations, SIAM Review, January 1967.
0036
0037       The following table lists by experiment, how large z can be in order to
0038       ensure the continued fraction converges to the correct value:
0039
0040           a         b    max  z
0041          13,      -130,      22
0042          13,     -1300,     335
0043          13,    -13000,    3585
0044         130,      -130,       8
0045         130,     -1300,     263
0046         130,   - 13000,    3420
0047        1300,      -130,       1
0048        1300,     -1300,      90
0049        1300,    -13000,    2650
0050       13000,       -13,       -
0051       13000,      -130,       -
0052       13000,     -1300,      13
0053       13000,    -13000,     750
0054
0055       So try z_limit = -b / (4 - 5 * sqrt(log(a)) * a / b);
0056       or     z_limit = -b / (4 - 5 * (log(a)) * a / b)  for a < 100
0057
0058       This still isn't quite right for both a and b small, but we'll be using a Bessel approximation
0059       in that region anyway.
0060
0061       Normalization using wronksian {1,2} from A&S 13.1.20 (also 13.1.12, 13.1.13):
0062
0063       W{ M(a,b,z), z^(1-b)M(1+a-b, 2-b, z) } = (1-b)z^-b e^z
0064
0065        = M(a,b,z) M2'(a,b,z) - M'(a,b,z) M2(a,b,z)
0066        = M(a,b,z) [(a-b+1)z^(1-b)/(2-b) M2(a+1,b+1,z) + (1-b)z^-b M2(a,b,z)] - a/b M(a+1,b+1,z) z^(1-b)M2(a,b,z)
0067        = M(a,b,z) [(a-b+1)z^(1-b)/(2-b) M2(a+1,b+1,z) + (1-b)z^-b M2(a,b,z)] - a/b R(a,b,z) M(a,b,z) z^(1-b)M2(a,b,z)
0068        = M(a,b,z) [(a-b+1)z^(1-b)/(2-b) M2(a+1,b+1,z) + (1-b)z^-b M2(a,b,z) - a/b R(a,b,z) z^(1-b)M2(a,b,z) ]
0069        so:
0070        (1-b)e^z = M(a,b,z) [(a-b+1)z/(2-b) M2(a+1,b+1,z) + (1-b) M2(a,b,z) - a/b z R(a,b,z) M2(a,b,z) ]
0071
0072       */
0073
0074      template <class T>
0075      inline bool is_in_hypergeometric_1F1_from_function_ratio_negative_b_region(const T& a, const T& b, const T& z)
0076      {
0077         BOOST_MATH_STD_USING
0078         if (a < 100)
0079            return z < -b / (4 - 5 * (log(a)) * a / b);
0080         else
0081            return z < -b / (4 - 5 * sqrt(log(a)) * a / b);
0082      }
0083
0084      template <class T, class Policy>
0085      T hypergeometric_1F1_from_function_ratio_negative_b(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling, const T& ratio)
0086      {
0087         BOOST_MATH_STD_USING
0088         //
0089         // Let M2 = M(1+a-b, 2-b, z)
0090         // This is going to be a mighty big number:
0091         //
0092         long long local_scaling = 0;
0093         T M2 = boost::math::detail::hypergeometric_1F1_imp(T(1 + a - b), T(2 - b), z, pol, local_scaling);
0094         log_scaling -= local_scaling; // all the M2 terms are in the denominator
0095         //
0096         // Since a, b and z are all likely to be large we need the Wronksian
0097         // calculation below to not overflow, so scale everything right down:
0098         //
0099         if (fabs(M2) > 1)
0100         {
0101            long long s = lltrunc(log(fabs(M2)));
0102            log_scaling -= s;  // M2 will be in the denominator, so subtract the scaling!
0103            local_scaling += s;
0104            M2 *= exp(T(-s));
0105         }
0106         //
0107         // Let M3 = M(1+a-b + 1, 2-b + 1, z)
0108         // we can get to this from the ratio which is cheaper to calculate:
0109         //
0110         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0111         boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<T> coef2(1 + a - b + 1, 2 - b + 1, z);
0112         T M3 = boost::math::tools::function_ratio_from_backwards_recurrence(coef2, boost::math::policies::get_epsilon<T, Policy>(), max_iter) * M2;
0113         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_from_function_ratio_negative_b_positive_a<%1%>(%1%,%1%,%1%)", max_iter, pol);
0114         //
0115         // Get the RHS of the Wronksian:
0116         //
0117         long long fz = lltrunc(z);
0118         log_scaling += fz;
0119         T rhs = (1 - b) * exp(z - fz);
0120         //
0121         // We need to divide that by:
0122         // [(a-b+1)z/(2-b) M2(a+1,b+1,z) + (1-b) M2(a,b,z) - a/b z^b R(a,b,z) M2(a,b,z) ]
0123         // Note that at this stage, both M3 and M2 are scaled by exp(local_scaling).
0124         //
0125         T lhs = (a - b + 1) * z * M3 / (2 - b) + (1 - b) * M2 - a * z * ratio * M2 / b;
0126
0127         return rhs / lhs;
0128      }
0129
0130      template <class T, class Policy>
0131      T hypergeometric_1F1_from_function_ratio_negative_b(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0132      {
0133         BOOST_MATH_STD_USING
0134         //
0135         // Get the function ratio, M(a+1, b+1, z)/M(a,b,z):
0136         //
0137         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0138         boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<T> coef(a + 1, b + 1, z);
0139         T ratio = boost::math::tools::function_ratio_from_backwards_recurrence(coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0140         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_from_function_ratio_negative_b_positive_a<%1%>(%1%,%1%,%1%)", max_iter, pol);
0141         return hypergeometric_1F1_from_function_ratio_negative_b(a, b, z, pol, log_scaling, ratio);
0142      }
0143      //
0144      // And over again, this time via forwards recurrence when z is large enough:
0145      //
0146      template <class T>
0147      bool hypergeometric_1F1_is_in_forwards_recurence_for_negative_b_region(const T& a, const T& b, const T& z)
0148      {
0149         //
0150         // There's no easy relation between a, b and z that tells us whether we're in the region
0151         // where forwards recursion is stable, so use a lookup table, note that the minimum
0152         // permissible z-value is decreasing with a, and increasing with |b|:
0153         //
0154         static const float data[][3] = {
0155            {7.500e+00f, -7.500e+00f, 8.409e+00f },
0156            {7.500e+00f, -1.125e+01f, 8.409e+00f },
0157            {7.500e+00f, -1.688e+01f, 9.250e+00f },
0158            {7.500e+00f, -2.531e+01f, 1.119e+01f },
0159            {7.500e+00f, -3.797e+01f, 1.354e+01f },
0160            {7.500e+00f, -5.695e+01f, 1.983e+01f },
0161            {7.500e+00f, -8.543e+01f, 2.639e+01f },
0162            {7.500e+00f, -1.281e+02f, 3.864e+01f },
0163            {7.500e+00f, -1.922e+02f, 5.657e+01f },
0164            {7.500e+00f, -2.883e+02f, 8.283e+01f },
0165            {7.500e+00f, -4.325e+02f, 1.213e+02f },
0166            {7.500e+00f, -6.487e+02f, 1.953e+02f },
0167            {7.500e+00f, -9.731e+02f, 2.860e+02f },
0168            {7.500e+00f, -1.460e+03f, 4.187e+02f },
0169            {7.500e+00f, -2.189e+03f, 6.130e+02f },
0170            {7.500e+00f, -3.284e+03f, 9.872e+02f },
0171            {7.500e+00f, -4.926e+03f, 1.445e+03f },
0172            {7.500e+00f, -7.389e+03f, 2.116e+03f },
0173            {7.500e+00f, -1.108e+04f, 3.098e+03f },
0174            {7.500e+00f, -1.663e+04f, 4.990e+03f },
0175            {1.125e+01f, -7.500e+00f, 6.318e+00f },
0176            {1.125e+01f, -1.125e+01f, 6.950e+00f },
0177            {1.125e+01f, -1.688e+01f, 7.645e+00f },
0178            {1.125e+01f, -2.531e+01f, 9.250e+00f },
0179            {1.125e+01f, -3.797e+01f, 1.231e+01f },
0180            {1.125e+01f, -5.695e+01f, 1.639e+01f },
0181            {1.125e+01f, -8.543e+01f, 2.399e+01f },
0182            {1.125e+01f, -1.281e+02f, 3.513e+01f },
0183            {1.125e+01f, -1.922e+02f, 5.657e+01f },
0184            {1.125e+01f, -2.883e+02f, 8.283e+01f },
0185            {1.125e+01f, -4.325e+02f, 1.213e+02f },
0186            {1.125e+01f, -6.487e+02f, 1.776e+02f },
0187            {1.125e+01f, -9.731e+02f, 2.860e+02f },
0188            {1.125e+01f, -1.460e+03f, 4.187e+02f },
0189            {1.125e+01f, -2.189e+03f, 6.130e+02f },
0190            {1.125e+01f, -3.284e+03f, 9.872e+02f },
0191            {1.125e+01f, -4.926e+03f, 1.445e+03f },
0192            {1.125e+01f, -7.389e+03f, 2.116e+03f },
0193            {1.125e+01f, -1.108e+04f, 3.098e+03f },
0194            {1.125e+01f, -1.663e+04f, 4.990e+03f },
0195            {1.688e+01f, -7.500e+00f, 4.747e+00f },
0196            {1.688e+01f, -1.125e+01f, 5.222e+00f },
0197            {1.688e+01f, -1.688e+01f, 5.744e+00f },
0198            {1.688e+01f, -2.531e+01f, 7.645e+00f },
0199            {1.688e+01f, -3.797e+01f, 1.018e+01f },
0200            {1.688e+01f, -5.695e+01f, 1.490e+01f },
0201            {1.688e+01f, -8.543e+01f, 2.181e+01f },
0202            {1.688e+01f, -1.281e+02f, 3.193e+01f },
0203            {1.688e+01f, -1.922e+02f, 5.143e+01f },
0204            {1.688e+01f, -2.883e+02f, 7.530e+01f },
0205            {1.688e+01f, -4.325e+02f, 1.213e+02f },
0206            {1.688e+01f, -6.487e+02f, 1.776e+02f },
0207            {1.688e+01f, -9.731e+02f, 2.600e+02f },
0208            {1.688e+01f, -1.460e+03f, 4.187e+02f },
0209            {1.688e+01f, -2.189e+03f, 6.130e+02f },
0210            {1.688e+01f, -3.284e+03f, 9.872e+02f },
0211            {1.688e+01f, -4.926e+03f, 1.445e+03f },
0212            {1.688e+01f, -7.389e+03f, 2.116e+03f },
0213            {1.688e+01f, -1.108e+04f, 3.098e+03f },
0214            {1.688e+01f, -1.663e+04f, 4.990e+03f },
0215            {2.531e+01f, -7.500e+00f, 3.242e+00f },
0216            {2.531e+01f, -1.125e+01f, 3.566e+00f },
0217            {2.531e+01f, -1.688e+01f, 4.315e+00f },
0218            {2.531e+01f, -2.531e+01f, 5.744e+00f },
0219            {2.531e+01f, -3.797e+01f, 7.645e+00f },
0220            {2.531e+01f, -5.695e+01f, 1.231e+01f },
0221            {2.531e+01f, -8.543e+01f, 1.803e+01f },
0222            {2.531e+01f, -1.281e+02f, 2.903e+01f },
0223            {2.531e+01f, -1.922e+02f, 4.676e+01f },
0224            {2.531e+01f, -2.883e+02f, 6.845e+01f },
0225            {2.531e+01f, -4.325e+02f, 1.102e+02f },
0226            {2.531e+01f, -6.487e+02f, 1.776e+02f },
0227            {2.531e+01f, -9.731e+02f, 2.600e+02f },
0228            {2.531e+01f, -1.460e+03f, 4.187e+02f },
0229            {2.531e+01f, -2.189e+03f, 6.130e+02f },
0230            {2.531e+01f, -3.284e+03f, 8.974e+02f },
0231            {2.531e+01f, -4.926e+03f, 1.445e+03f },
0232            {2.531e+01f, -7.389e+03f, 2.116e+03f },
0233            {2.531e+01f, -1.108e+04f, 3.098e+03f },
0234            {2.531e+01f, -1.663e+04f, 4.990e+03f },
0235            {3.797e+01f, -7.500e+00f, 2.214e+00f },
0236            {3.797e+01f, -1.125e+01f, 2.679e+00f },
0237            {3.797e+01f, -1.688e+01f, 3.242e+00f },
0238            {3.797e+01f, -2.531e+01f, 4.315e+00f },
0239            {3.797e+01f, -3.797e+01f, 6.318e+00f },
0240            {3.797e+01f, -5.695e+01f, 9.250e+00f },
0241            {3.797e+01f, -8.543e+01f, 1.490e+01f },
0242            {3.797e+01f, -1.281e+02f, 2.399e+01f },
0243            {3.797e+01f, -1.922e+02f, 3.864e+01f },
0244            {3.797e+01f, -2.883e+02f, 6.223e+01f },
0245            {3.797e+01f, -4.325e+02f, 1.002e+02f },
0246            {3.797e+01f, -6.487e+02f, 1.614e+02f },
0247            {3.797e+01f, -9.731e+02f, 2.600e+02f },
0248            {3.797e+01f, -1.460e+03f, 3.806e+02f },
0249            {3.797e+01f, -2.189e+03f, 6.130e+02f },
0250            {3.797e+01f, -3.284e+03f, 8.974e+02f },
0251            {3.797e+01f, -4.926e+03f, 1.445e+03f },
0252            {3.797e+01f, -7.389e+03f, 2.116e+03f },
0253            {3.797e+01f, -1.108e+04f, 3.098e+03f },
0254            {3.797e+01f, -1.663e+04f, 4.990e+03f },
0255            {5.695e+01f, -7.500e+00f, 1.513e+00f },
0256            {5.695e+01f, -1.125e+01f, 1.830e+00f },
0257            {5.695e+01f, -1.688e+01f, 2.214e+00f },
0258            {5.695e+01f, -2.531e+01f, 3.242e+00f },
0259            {5.695e+01f, -3.797e+01f, 4.315e+00f },
0260            {5.695e+01f, -5.695e+01f, 7.645e+00f },
0261            {5.695e+01f, -8.543e+01f, 1.231e+01f },
0262            {5.695e+01f, -1.281e+02f, 1.983e+01f },
0263            {5.695e+01f, -1.922e+02f, 3.513e+01f },
0264            {5.695e+01f, -2.883e+02f, 5.657e+01f },
0265            {5.695e+01f, -4.325e+02f, 9.111e+01f },
0266            {5.695e+01f, -6.487e+02f, 1.467e+02f },
0267            {5.695e+01f, -9.731e+02f, 2.363e+02f },
0268            {5.695e+01f, -1.460e+03f, 3.806e+02f },
0269            {5.695e+01f, -2.189e+03f, 5.572e+02f },
0270            {5.695e+01f, -3.284e+03f, 8.974e+02f },
0271            {5.695e+01f, -4.926e+03f, 1.314e+03f },
0272            {5.695e+01f, -7.389e+03f, 2.116e+03f },
0273            {5.695e+01f, -1.108e+04f, 3.098e+03f },
0274            {5.695e+01f, -1.663e+04f, 4.990e+03f },
0275            {8.543e+01f, -7.500e+00f, 1.250e+00f },
0276            {8.543e+01f, -1.125e+01f, 1.250e+00f },
0277            {8.543e+01f, -1.688e+01f, 1.513e+00f },
0278            {8.543e+01f, -2.531e+01f, 2.214e+00f },
0279            {8.543e+01f, -3.797e+01f, 3.242e+00f },
0280            {8.543e+01f, -5.695e+01f, 5.222e+00f },
0281            {8.543e+01f, -8.543e+01f, 9.250e+00f },
0282            {8.543e+01f, -1.281e+02f, 1.639e+01f },
0283            {8.543e+01f, -1.922e+02f, 2.903e+01f },
0284            {8.543e+01f, -2.883e+02f, 5.143e+01f },
0285            {8.543e+01f, -4.325e+02f, 8.283e+01f },
0286            {8.543e+01f, -6.487e+02f, 1.334e+02f },
0287            {8.543e+01f, -9.731e+02f, 2.148e+02f },
0288            {8.543e+01f, -1.460e+03f, 3.460e+02f },
0289            {8.543e+01f, -2.189e+03f, 5.572e+02f },
0290            {8.543e+01f, -3.284e+03f, 8.974e+02f },
0291            {8.543e+01f, -4.926e+03f, 1.314e+03f },
0292            {8.543e+01f, -7.389e+03f, 2.116e+03f },
0293            {8.543e+01f, -1.108e+04f, 3.098e+03f },
0294            {8.543e+01f, -1.663e+04f, 4.536e+03f },
0295            {1.281e+02f, -7.500e+00f, 1.250e+00f },
0296            {1.281e+02f, -1.125e+01f, 1.250e+00f },
0297            {1.281e+02f, -1.688e+01f, 1.250e+00f },
0298            {1.281e+02f, -2.531e+01f, 1.513e+00f },
0299            {1.281e+02f, -3.797e+01f, 2.214e+00f },
0300            {1.281e+02f, -5.695e+01f, 3.923e+00f },
0301            {1.281e+02f, -8.543e+01f, 6.950e+00f },
0302            {1.281e+02f, -1.281e+02f, 1.231e+01f },
0303            {1.281e+02f, -1.922e+02f, 2.181e+01f },
0304            {1.281e+02f, -2.883e+02f, 4.250e+01f },
0305            {1.281e+02f, -4.325e+02f, 6.845e+01f },
0306            {1.281e+02f, -6.487e+02f, 1.213e+02f },
0307            {1.281e+02f, -9.731e+02f, 1.953e+02f },
0308            {1.281e+02f, -1.460e+03f, 3.460e+02f },
0309            {1.281e+02f, -2.189e+03f, 5.572e+02f },
0310            {1.281e+02f, -3.284e+03f, 8.159e+02f },
0311            {1.281e+02f, -4.926e+03f, 1.314e+03f },
0312            {1.281e+02f, -7.389e+03f, 1.924e+03f },
0313            {1.281e+02f, -1.108e+04f, 3.098e+03f },
0314            {1.281e+02f, -1.663e+04f, 4.536e+03f },
0315            {1.922e+02f, -7.500e+00f, 1.250e+00f },
0316            {1.922e+02f, -1.125e+01f, 1.250e+00f },
0317            {1.922e+02f, -1.688e+01f, 1.250e+00f },
0318            {1.922e+02f, -2.531e+01f, 1.250e+00f },
0319            {1.922e+02f, -3.797e+01f, 1.664e+00f },
0320            {1.922e+02f, -5.695e+01f, 2.679e+00f },
0321            {1.922e+02f, -8.543e+01f, 5.222e+00f },
0322            {1.922e+02f, -1.281e+02f, 9.250e+00f },
0323            {1.922e+02f, -1.922e+02f, 1.803e+01f },
0324            {1.922e+02f, -2.883e+02f, 3.193e+01f },
0325            {1.922e+02f, -4.325e+02f, 5.657e+01f },
0326            {1.922e+02f, -6.487e+02f, 1.002e+02f },
0327            {1.922e+02f, -9.731e+02f, 1.776e+02f },
0328            {1.922e+02f, -1.460e+03f, 3.145e+02f },
0329            {1.922e+02f, -2.189e+03f, 5.066e+02f },
0330            {1.922e+02f, -3.284e+03f, 8.159e+02f },
0331            {1.922e+02f, -4.926e+03f, 1.194e+03f },
0332            {1.922e+02f, -7.389e+03f, 1.924e+03f },
0333            {1.922e+02f, -1.108e+04f, 3.098e+03f },
0334            {1.922e+02f, -1.663e+04f, 4.536e+03f },
0335            {2.883e+02f, -7.500e+00f, 1.250e+00f },
0336            {2.883e+02f, -1.125e+01f, 1.250e+00f },
0337            {2.883e+02f, -1.688e+01f, 1.250e+00f },
0338            {2.883e+02f, -2.531e+01f, 1.250e+00f },
0339            {2.883e+02f, -3.797e+01f, 1.250e+00f },
0340            {2.883e+02f, -5.695e+01f, 2.013e+00f },
0341            {2.883e+02f, -8.543e+01f, 3.566e+00f },
0342            {2.883e+02f, -1.281e+02f, 6.950e+00f },
0343            {2.883e+02f, -1.922e+02f, 1.354e+01f },
0344            {2.883e+02f, -2.883e+02f, 2.399e+01f },
0345            {2.883e+02f, -4.325e+02f, 4.676e+01f },
0346            {2.883e+02f, -6.487e+02f, 8.283e+01f },
0347            {2.883e+02f, -9.731e+02f, 1.614e+02f },
0348            {2.883e+02f, -1.460e+03f, 2.600e+02f },
0349            {2.883e+02f, -2.189e+03f, 4.605e+02f },
0350            {2.883e+02f, -3.284e+03f, 7.417e+02f },
0351            {2.883e+02f, -4.926e+03f, 1.194e+03f },
0352            {2.883e+02f, -7.389e+03f, 1.924e+03f },
0353            {2.883e+02f, -1.108e+04f, 2.817e+03f },
0354            {2.883e+02f, -1.663e+04f, 4.536e+03f },
0355            {4.325e+02f, -7.500e+00f, 1.250e+00f },
0356            {4.325e+02f, -1.125e+01f, 1.250e+00f },
0357            {4.325e+02f, -1.688e+01f, 1.250e+00f },
0358            {4.325e+02f, -2.531e+01f, 1.250e+00f },
0359            {4.325e+02f, -3.797e+01f, 1.250e+00f },
0360            {4.325e+02f, -5.695e+01f, 1.375e+00f },
0361            {4.325e+02f, -8.543e+01f, 2.436e+00f },
0362            {4.325e+02f, -1.281e+02f, 4.747e+00f },
0363            {4.325e+02f, -1.922e+02f, 9.250e+00f },
0364            {4.325e+02f, -2.883e+02f, 1.803e+01f },
0365            {4.325e+02f, -4.325e+02f, 3.513e+01f },
0366            {4.325e+02f, -6.487e+02f, 6.845e+01f },
0367            {4.325e+02f, -9.731e+02f, 1.334e+02f },
0368            {4.325e+02f, -1.460e+03f, 2.363e+02f },
0369            {4.325e+02f, -2.189e+03f, 3.806e+02f },
0370            {4.325e+02f, -3.284e+03f, 6.743e+02f },
0371            {4.325e+02f, -4.926e+03f, 1.086e+03f },
0372            {4.325e+02f, -7.389e+03f, 1.749e+03f },
0373            {4.325e+02f, -1.108e+04f, 2.817e+03f },
0374            {4.325e+02f, -1.663e+04f, 4.536e+03f },
0375            {6.487e+02f, -7.500e+00f, 1.250e+00f },
0376            {6.487e+02f, -1.125e+01f, 1.250e+00f },
0377            {6.487e+02f, -1.688e+01f, 1.250e+00f },
0378            {6.487e+02f, -2.531e+01f, 1.250e+00f },
0379            {6.487e+02f, -3.797e+01f, 1.250e+00f },
0380            {6.487e+02f, -5.695e+01f, 1.250e+00f },
0381            {6.487e+02f, -8.543e+01f, 1.664e+00f },
0382            {6.487e+02f, -1.281e+02f, 3.242e+00f },
0383            {6.487e+02f, -1.922e+02f, 6.950e+00f },
0384            {6.487e+02f, -2.883e+02f, 1.354e+01f },
0385            {6.487e+02f, -4.325e+02f, 2.639e+01f },
0386            {6.487e+02f, -6.487e+02f, 5.143e+01f },
0387            {6.487e+02f, -9.731e+02f, 1.002e+02f },
0388            {6.487e+02f, -1.460e+03f, 1.953e+02f },
0389            {6.487e+02f, -2.189e+03f, 3.460e+02f },
0390            {6.487e+02f, -3.284e+03f, 6.130e+02f },
0391            {6.487e+02f, -4.926e+03f, 9.872e+02f },
0392            {6.487e+02f, -7.389e+03f, 1.590e+03f },
0393            {6.487e+02f, -1.108e+04f, 2.561e+03f },
0394            {6.487e+02f, -1.663e+04f, 4.124e+03f },
0395            {9.731e+02f, -7.500e+00f, 1.250e+00f },
0396            {9.731e+02f, -1.125e+01f, 1.250e+00f },
0397            {9.731e+02f, -1.688e+01f, 1.250e+00f },
0398            {9.731e+02f, -2.531e+01f, 1.250e+00f },
0399            {9.731e+02f, -3.797e+01f, 1.250e+00f },
0400            {9.731e+02f, -5.695e+01f, 1.250e+00f },
0401            {9.731e+02f, -8.543e+01f, 1.250e+00f },
0402            {9.731e+02f, -1.281e+02f, 2.214e+00f },
0403            {9.731e+02f, -1.922e+02f, 4.747e+00f },
0404            {9.731e+02f, -2.883e+02f, 9.250e+00f },
0405            {9.731e+02f, -4.325e+02f, 1.983e+01f },
0406            {9.731e+02f, -6.487e+02f, 3.864e+01f },
0407            {9.731e+02f, -9.731e+02f, 7.530e+01f },
0408            {9.731e+02f, -1.460e+03f, 1.467e+02f },
0409            {9.731e+02f, -2.189e+03f, 2.860e+02f },
0410            {9.731e+02f, -3.284e+03f, 5.066e+02f },
0411            {9.731e+02f, -4.926e+03f, 8.974e+02f },
0412            {9.731e+02f, -7.389e+03f, 1.445e+03f },
0413            {9.731e+02f, -1.108e+04f, 2.561e+03f },
0414            {9.731e+02f, -1.663e+04f, 4.124e+03f },
0415            {1.460e+03f, -7.500e+00f, 1.250e+00f },
0416            {1.460e+03f, -1.125e+01f, 1.250e+00f },
0417            {1.460e+03f, -1.688e+01f, 1.250e+00f },
0418            {1.460e+03f, -2.531e+01f, 1.250e+00f },
0419            {1.460e+03f, -3.797e+01f, 1.250e+00f },
0420            {1.460e+03f, -5.695e+01f, 1.250e+00f },
0421            {1.460e+03f, -8.543e+01f, 1.250e+00f },
0422            {1.460e+03f, -1.281e+02f, 1.513e+00f },
0423            {1.460e+03f, -1.922e+02f, 3.242e+00f },
0424            {1.460e+03f, -2.883e+02f, 6.950e+00f },
0425            {1.460e+03f, -4.325e+02f, 1.354e+01f },
0426            {1.460e+03f, -6.487e+02f, 2.903e+01f },
0427            {1.460e+03f, -9.731e+02f, 5.657e+01f },
0428            {1.460e+03f, -1.460e+03f, 1.213e+02f },
0429            {1.460e+03f, -2.189e+03f, 2.148e+02f },
0430            {1.460e+03f, -3.284e+03f, 4.187e+02f },
0431            {1.460e+03f, -4.926e+03f, 7.417e+02f },
0432            {1.460e+03f, -7.389e+03f, 1.314e+03f },
0433            {1.460e+03f, -1.108e+04f, 2.328e+03f },
0434            {1.460e+03f, -1.663e+04f, 3.749e+03f },
0435            {2.189e+03f, -7.500e+00f, 1.250e+00f },
0436            {2.189e+03f, -1.125e+01f, 1.250e+00f },
0437            {2.189e+03f, -1.688e+01f, 1.250e+00f },
0438            {2.189e+03f, -2.531e+01f, 1.250e+00f },
0439            {2.189e+03f, -3.797e+01f, 1.250e+00f },
0440            {2.189e+03f, -5.695e+01f, 1.250e+00f },
0441            {2.189e+03f, -8.543e+01f, 1.250e+00f },
0442            {2.189e+03f, -1.281e+02f, 1.250e+00f },
0443            {2.189e+03f, -1.922e+02f, 2.214e+00f },
0444            {2.189e+03f, -2.883e+02f, 4.747e+00f },
0445            {2.189e+03f, -4.325e+02f, 9.250e+00f },
0446            {2.189e+03f, -6.487e+02f, 1.983e+01f },
0447            {2.189e+03f, -9.731e+02f, 4.250e+01f },
0448            {2.189e+03f, -1.460e+03f, 8.283e+01f },
0449            {2.189e+03f, -2.189e+03f, 1.776e+02f },
0450            {2.189e+03f, -3.284e+03f, 3.460e+02f },
0451            {2.189e+03f, -4.926e+03f, 6.130e+02f },
0452            {2.189e+03f, -7.389e+03f, 1.086e+03f },
0453            {2.189e+03f, -1.108e+04f, 1.924e+03f },
0454            {2.189e+03f, -1.663e+04f, 3.408e+03f },
0455            {3.284e+03f, -7.500e+00f, 1.250e+00f },
0456            {3.284e+03f, -1.125e+01f, 1.250e+00f },
0457            {3.284e+03f, -1.688e+01f, 1.250e+00f },
0458            {3.284e+03f, -2.531e+01f, 1.250e+00f },
0459            {3.284e+03f, -3.797e+01f, 1.250e+00f },
0460            {3.284e+03f, -5.695e+01f, 1.250e+00f },
0461            {3.284e+03f, -8.543e+01f, 1.250e+00f },
0462            {3.284e+03f, -1.281e+02f, 1.250e+00f },
0463            {3.284e+03f, -1.922e+02f, 1.513e+00f },
0464            {3.284e+03f, -2.883e+02f, 3.242e+00f },
0465            {3.284e+03f, -4.325e+02f, 6.318e+00f },
0466            {3.284e+03f, -6.487e+02f, 1.354e+01f },
0467            {3.284e+03f, -9.731e+02f, 2.903e+01f },
0468            {3.284e+03f, -1.460e+03f, 6.223e+01f },
0469            {3.284e+03f, -2.189e+03f, 1.334e+02f },
0470            {3.284e+03f, -3.284e+03f, 2.600e+02f },
0471            {3.284e+03f, -4.926e+03f, 5.066e+02f },
0472            {3.284e+03f, -7.389e+03f, 9.872e+02f },
0473            {3.284e+03f, -1.108e+04f, 1.749e+03f },
0474            {3.284e+03f, -1.663e+04f, 3.098e+03f },
0475            {4.926e+03f, -7.500e+00f, 1.250e+00f },
0476            {4.926e+03f, -1.125e+01f, 1.250e+00f },
0477            {4.926e+03f, -1.688e+01f, 1.250e+00f },
0478            {4.926e+03f, -2.531e+01f, 1.250e+00f },
0479            {4.926e+03f, -3.797e+01f, 1.250e+00f },
0480            {4.926e+03f, -5.695e+01f, 1.250e+00f },
0481            {4.926e+03f, -8.543e+01f, 1.250e+00f },
0482            {4.926e+03f, -1.281e+02f, 1.250e+00f },
0483            {4.926e+03f, -1.922e+02f, 1.250e+00f },
0484            {4.926e+03f, -2.883e+02f, 2.013e+00f },
0485            {4.926e+03f, -4.325e+02f, 4.315e+00f },
0486            {4.926e+03f, -6.487e+02f, 9.250e+00f },
0487            {4.926e+03f, -9.731e+02f, 2.181e+01f },
0488            {4.926e+03f, -1.460e+03f, 4.250e+01f },
0489            {4.926e+03f, -2.189e+03f, 9.111e+01f },
0490            {4.926e+03f, -3.284e+03f, 1.953e+02f },
0491            {4.926e+03f, -4.926e+03f, 3.806e+02f },
0492            {4.926e+03f, -7.389e+03f, 7.417e+02f },
0493            {4.926e+03f, -1.108e+04f, 1.445e+03f },
0494            {4.926e+03f, -1.663e+04f, 2.561e+03f },
0495            {7.389e+03f, -7.500e+00f, 1.250e+00f },
0496            {7.389e+03f, -1.125e+01f, 1.250e+00f },
0497            {7.389e+03f, -1.688e+01f, 1.250e+00f },
0498            {7.389e+03f, -2.531e+01f, 1.250e+00f },
0499            {7.389e+03f, -3.797e+01f, 1.250e+00f },
0500            {7.389e+03f, -5.695e+01f, 1.250e+00f },
0501            {7.389e+03f, -8.543e+01f, 1.250e+00f },
0502            {7.389e+03f, -1.281e+02f, 1.250e+00f },
0503            {7.389e+03f, -1.922e+02f, 1.250e+00f },
0504            {7.389e+03f, -2.883e+02f, 1.375e+00f },
0505            {7.389e+03f, -4.325e+02f, 2.947e+00f },
0506            {7.389e+03f, -6.487e+02f, 6.318e+00f },
0507            {7.389e+03f, -9.731e+02f, 1.490e+01f },
0508            {7.389e+03f, -1.460e+03f, 3.193e+01f },
0509            {7.389e+03f, -2.189e+03f, 6.845e+01f },
0510            {7.389e+03f, -3.284e+03f, 1.334e+02f },
0511            {7.389e+03f, -4.926e+03f, 2.860e+02f },
0512            {7.389e+03f, -7.389e+03f, 5.572e+02f },
0513            {7.389e+03f, -1.108e+04f, 1.086e+03f },
0514            {7.389e+03f, -1.663e+04f, 2.116e+03f },
0515            {1.108e+04f, -7.500e+00f, 1.250e+00f },
0516            {1.108e+04f, -1.125e+01f, 1.250e+00f },
0517            {1.108e+04f, -1.688e+01f, 1.250e+00f },
0518            {1.108e+04f, -2.531e+01f, 1.250e+00f },
0519            {1.108e+04f, -3.797e+01f, 1.250e+00f },
0520            {1.108e+04f, -5.695e+01f, 1.250e+00f },
0521            {1.108e+04f, -8.543e+01f, 1.250e+00f },
0522            {1.108e+04f, -1.281e+02f, 1.250e+00f },
0523            {1.108e+04f, -1.922e+02f, 1.250e+00f },
0524            {1.108e+04f, -2.883e+02f, 1.250e+00f },
0525            {1.108e+04f, -4.325e+02f, 2.013e+00f },
0526            {1.108e+04f, -6.487e+02f, 4.315e+00f },
0527            {1.108e+04f, -9.731e+02f, 1.018e+01f },
0528            {1.108e+04f, -1.460e+03f, 2.181e+01f },
0529            {1.108e+04f, -2.189e+03f, 4.676e+01f },
0530            {1.108e+04f, -3.284e+03f, 1.002e+02f },
0531            {1.108e+04f, -4.926e+03f, 2.148e+02f },
0532            {1.108e+04f, -7.389e+03f, 4.187e+02f },
0533            {1.108e+04f, -1.108e+04f, 8.974e+02f },
0534            {1.108e+04f, -1.663e+04f, 1.749e+03f },
0535            {1.663e+04f, -7.500e+00f, 1.250e+00f },
0536            {1.663e+04f, -1.125e+01f, 1.250e+00f },
0537            {1.663e+04f, -1.688e+01f, 1.250e+00f },
0538            {1.663e+04f, -2.531e+01f, 1.250e+00f },
0539            {1.663e+04f, -3.797e+01f, 1.250e+00f },
0540            {1.663e+04f, -5.695e+01f, 1.250e+00f },
0541            {1.663e+04f, -8.543e+01f, 1.250e+00f },
0542            {1.663e+04f, -1.281e+02f, 1.250e+00f },
0543            {1.663e+04f, -1.922e+02f, 1.250e+00f },
0544            {1.663e+04f, -2.883e+02f, 1.250e+00f },
0545            {1.663e+04f, -4.325e+02f, 1.375e+00f },
0546            {1.663e+04f, -6.487e+02f, 2.947e+00f },
0547            {1.663e+04f, -9.731e+02f, 6.318e+00f },
0548            {1.663e+04f, -1.460e+03f, 1.490e+01f },
0549            {1.663e+04f, -2.189e+03f, 3.193e+01f },
0550            {1.663e+04f, -3.284e+03f, 6.845e+01f },
0551            {1.663e+04f, -4.926e+03f, 1.467e+02f },
0552            {1.663e+04f, -7.389e+03f, 3.145e+02f },
0553            {1.663e+04f, -1.108e+04f, 6.743e+02f },
0554            {1.663e+04f, -1.663e+04f, 1.314e+03f },
0555         };
0556         if ((a > 1.663e+04) || (-b > 1.663e+04))
0557            return z > -b;  // Way overly conservative?
0558         if (a < data[0][0])
0559            return false;
0560         int index = 0;
0561         while (data[index][0] < a)
0562            ++index;
0563         if(a != data[index][0])
0564            --index;
0565         while ((data[index][1] < b) && (data[index][2] > 1.25))
0566            --index;
0567         ++index;
0568         BOOST_MATH_ASSERT(a > data[index][0]);
0569         BOOST_MATH_ASSERT(-b < -data[index][1]);
0570         return z > data[index][2];
0571      }
0572      template <class T, class Policy>
0573      T hypergeometric_1F1_from_function_ratio_negative_b_forwards(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0574      {
0575         BOOST_MATH_STD_USING
0576         //
0577         // Get the function ratio, M(a+1, b+1, z)/M(a,b,z):
0578         //
0579         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0580         boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<T> coef(a, b, z);
0581         T ratio = 1 / boost::math::tools::function_ratio_from_forwards_recurrence(coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0582         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_from_function_ratio_negative_b_positive_a<%1%>(%1%,%1%,%1%)", max_iter, pol);
0583         //
0584         // We can't normalise via the Wronksian as the subtraction in the Wronksian will suffer an exquisite amount of cancellation -
0585         // potentially many hundreds of digits in this region.  However, if forwards iteration is stable at this point
0586         // it will also be stable for M(a, b+1, z) etc all the way up to the origin, and hopefully one step beyond.  So
0587         // use a reference value just over the origin to normalise:
0588         //
0589         long long scale = 0;
0590         int steps = itrunc(ceil(-b));
0591         T reference_value = hypergeometric_1F1_imp(T(a + steps), T(b + steps), z, pol, log_scaling);
0592         T found = boost::math::tools::apply_recurrence_relation_forward(boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<T>(a + 1, b + 1, z), steps - 1, T(1), ratio, &scale);
0593         log_scaling -= scale;
0594         if ((fabs(reference_value) < 1) && (fabs(reference_value) < tools::min_value<T>() * fabs(found)))
0595         {
0596            // Possible underflow, rescale
0597            long long s = lltrunc(tools::log_max_value<T>());
0598            log_scaling -= s;
0599            reference_value *= exp(T(s));
0600         }
0601         else if ((fabs(found) < 1) && (fabs(reference_value) > tools::max_value<T>() * fabs(found)))
0602         {
0603            // Overflow, rescale:
0604            long long s = lltrunc(tools::log_max_value<T>());
0605            log_scaling += s;
0606            reference_value /= exp(T(s));
0607         }
0608         return reference_value / found;
0609      }
0610
0611
0612
0613      //
0614      // This next version is largely the same as above, but calculates the ratio for the b recurrence relation
0615      // which has a larger area of stability than the ab recurrence when a,b < 0.  We can then use a single
0616      // recurrence step to convert this to the ratio for the ab recursion and proceed largely as before.
0617      // The routine is quite insensitive to the size of z, but requires |a| < |5b| for accuracy.
0618      // Fortunately the accuracy outside this domain falls off steadily rather than suddenly switching
0619      // to a different behaviour.
0620      //
0621      template <class T, class Policy>
0622      T hypergeometric_1F1_from_function_ratio_negative_ab(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
0623      {
0624         BOOST_MATH_STD_USING
0625         //
0626         // Get the function ratio, M(a+1, b+1, z)/M(a,b,z):
0627         //
0628         std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
0629         boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> coef(a, b + 1, z);
0630         T ratio = boost::math::tools::function_ratio_from_backwards_recurrence(coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
0631         boost::math::policies::check_series_iterations<T>("boost::math::hypergeometric_1F1_from_function_ratio_negative_b_positive_a<%1%>(%1%,%1%,%1%)", max_iter, pol);
0632         //
0633         // We need to use A&S 13.4.3 to convert a ratio for M(a, b + 1, z) / M(a, b, z)
0634         // to M(a+1, b+1, z) / M(a, b, z)
0635         //
0636         // We have:        (a-b)M(a, b+1, z) - aM(a+1, b+1, z) + bM(a, b, z) = 0
0637         // and therefore:  M(a + 1, b + 1, z) / M(a, b, z) = ((a - b)M(a, b + 1, z) / M(a, b, z) + b) / a
0638         //
0639         ratio = ((a - b) * ratio + b) / a;
0640         //
0641         // Let M2 = M(1+a-b, 2-b, z)
0642         // This is going to be a mighty big number:
0643         //
0644         long long local_scaling = 0;
0645         T M2 = boost::math::detail::hypergeometric_1F1_imp(T(1 + a - b), T(2 - b), z, pol, local_scaling);
0646         log_scaling -= local_scaling; // all the M2 terms are in the denominator
0647         //
0648         // Let M3 = M(1+a-b + 1, 2-b + 1, z)
0649         // We don't use the ratio to get this as it's not clear that it's reliable:
0650         //
0651         long long local_scaling2 = 0;
0652         T M3 = boost::math::detail::hypergeometric_1F1_imp(T(2 + a - b), T(3 - b), z, pol, local_scaling2);
0653         //
0654         // M2 and M3 must be identically scaled:
0655         //
0656         if (local_scaling != local_scaling2)
0657         {
0658            M3 *= exp(T(local_scaling2 - local_scaling));
0659         }
0660         //
0661         // Get the RHS of the Wronksian:
0662         //
0663         long long fz = lltrunc(z);
0664         log_scaling += fz;
0665         T rhs = (1 - b) * exp(z - fz);
0666         //
0667         // We need to divide that by:
0668         // [(a-b+1)z/(2-b) M2(a+1,b+1,z) + (1-b) M2(a,b,z) - a/b z^b R(a,b,z) M2(a,b,z) ]
0669         // Note that at this stage, both M3 and M2 are scaled by exp(local_scaling).
0670         //
0671         T lhs = (a - b + 1) * z * M3 / (2 - b) + (1 - b) * M2 - a * z * ratio * M2 / b;
0672
0673         return rhs / lhs;
0674      }
0675
0676   } } } // namespaces
0677
0678 #endif // BOOST_HYPERGEOMETRIC_1F1_BY_RATIOS_HPP_