CoolVLViewer/include/boost/math/special_functions/detail/hypergeometric_1F1_by_ratios.hpp

///////////////////////////////////////////////////////////////////////////////
//  Copyright 2018 John Maddock
//  Distributed under the Boost
//  Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_HYPERGEOMETRIC_1F1_BY_RATIOS_HPP_
#define BOOST_HYPERGEOMETRIC_1F1_BY_RATIOS_HPP_

#include <boost/math/tools/recurrence.hpp>
#include <boost/math/policies/error_handling.hpp>

  namespace boost { namespace math { namespace detail {

     template <class T, class Policy>
     T hypergeometric_1F1_imp(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling);

     /*
      Evaluation by method of ratios for domain b < 0 < a,z

      We first convert the recurrence relation into a ratio
      of M(a+1, b+1, z) / M(a, b, z) using Shintan's equivalence
      between solving a recurrence relation using Miller's method
      and continued fractions.  The continued fraction is VERY rapid
      to converge (typically < 10 terms), but may converge to entirely
      the wrong value if we're in a bad part of the domain.  Strangely
      it seems to matter not whether we use recurrence on a, b or a and b
      they all work or not work about the same, so we might as well make
      life easy for ourselves and use the a and b recurrence to avoid
      having to apply one extra recurrence to convert from an a or b
      recurrence to an a+b one.

      See: H. Shintan, Note on Miller's recurrence algorithm, J. Sci. Hiroshima Univ. Ser. A-I Math., 29 (1965), pp. 121-133.
      Also: Computational Aspects of Three Term Recurrence Relations, SIAM Review, January 1967.

      The following table lists by experiment, how large z can be in order to
      ensure the continued fraction converges to the correct value:

          a         b    max  z
         13,      -130,      22
         13,     -1300,     335
         13,    -13000,    3585
        130,      -130,       8
        130,     -1300,     263
        130,   - 13000,    3420
       1300,      -130,       1
       1300,     -1300,      90
       1300,    -13000,    2650
      13000,       -13,       -
      13000,      -130,       -
      13000,     -1300,      13
      13000,    -13000,     750

      So try z_limit = -b / (4 - 5 * sqrt(log(a)) * a / b);
      or     z_limit = -b / (4 - 5 * (log(a)) * a / b)  for a < 100
      
      This still isn't quite right for both a and b small, but we'll be using a Bessel approximation
      in that region anyway.

      Normalization using wronksian {1,2} from A&S 13.1.20 (also 13.1.12, 13.1.13):

      W{ M(a,b,z), z^(1-b)M(1+a-b, 2-b, z) } = (1-b)z^-b e^z

       = M(a,b,z) M2'(a,b,z) - M'(a,b,z) M2(a,b,z)
       = 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)
       = 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)
       = 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) ]
       so:
       (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) ]

      */

     template <class T>
     inline bool is_in_hypergeometric_1F1_from_function_ratio_negative_b_region(const T& a, const T& b, const T& z)
     {
        BOOST_MATH_STD_USING
        if (a < 100)
           return z < -b / (4 - 5 * (log(a)) * a / b);
        else
           return z < -b / (4 - 5 * sqrt(log(a)) * a / b);
     }

     template <class T, class Policy>
     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)
     {
        BOOST_MATH_STD_USING
        //
        // Let M2 = M(1+a-b, 2-b, z)
        // This is going to be a mighty big number:
        //
        long long local_scaling = 0;
        T M2 = boost::math::detail::hypergeometric_1F1_imp(T(1 + a - b), T(2 - b), z, pol, local_scaling);
        log_scaling -= local_scaling; // all the M2 terms are in the denominator
        //
        // Since a, b and z are all likely to be large we need the Wronksian
        // calculation below to not overflow, so scale everything right down:
        //
        if (fabs(M2) > 1)
        {
           long long s = lltrunc(log(fabs(M2)));
           log_scaling -= s;  // M2 will be in the denominator, so subtract the scaling!
           local_scaling += s;
           M2 *= exp(T(-s));
        }
        //
        // Let M3 = M(1+a-b + 1, 2-b + 1, z)
        // we can get to this from the ratio which is cheaper to calculate:
        //
        std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
        boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<T> coef2(1 + a - b + 1, 2 - b + 1, z);
        T M3 = boost::math::tools::function_ratio_from_backwards_recurrence(coef2, boost::math::policies::get_epsilon<T, Policy>(), max_iter) * M2;
        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);
        //
        // Get the RHS of the Wronksian:
        //
        long long fz = lltrunc(z);
        log_scaling += fz;
        T rhs = (1 - b) * exp(z - fz);
        //
        // We need to divide that by:
        // [(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) ]
        // Note that at this stage, both M3 and M2 are scaled by exp(local_scaling).
        //
        T lhs = (a - b + 1) * z * M3 / (2 - b) + (1 - b) * M2 - a * z * ratio * M2 / b;

        return rhs / lhs;
     }

     template <class T, class Policy>
     T hypergeometric_1F1_from_function_ratio_negative_b(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
     {
        BOOST_MATH_STD_USING
        //
        // Get the function ratio, M(a+1, b+1, z)/M(a,b,z):
        //
        std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
        boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<T> coef(a + 1, b + 1, z);
        T ratio = boost::math::tools::function_ratio_from_backwards_recurrence(coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
        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);
        return hypergeometric_1F1_from_function_ratio_negative_b(a, b, z, pol, log_scaling, ratio);
     }
     //
     // And over again, this time via forwards recurrence when z is large enough:
     //
     template <class T>
     bool hypergeometric_1F1_is_in_forwards_recurence_for_negative_b_region(const T& a, const T& b, const T& z)
     {
        //
        // There's no easy relation between a, b and z that tells us whether we're in the region
        // where forwards recursion is stable, so use a lookup table, note that the minimum
        // permissible z-value is decreasing with a, and increasing with |b|:
        //
        static const float data[][3] = {
           {7.500e+00f, -7.500e+00f, 8.409e+00f },
           {7.500e+00f, -1.125e+01f, 8.409e+00f },
           {7.500e+00f, -1.688e+01f, 9.250e+00f },
           {7.500e+00f, -2.531e+01f, 1.119e+01f },
           {7.500e+00f, -3.797e+01f, 1.354e+01f },
           {7.500e+00f, -5.695e+01f, 1.983e+01f },
           {7.500e+00f, -8.543e+01f, 2.639e+01f },
           {7.500e+00f, -1.281e+02f, 3.864e+01f },
           {7.500e+00f, -1.922e+02f, 5.657e+01f },
           {7.500e+00f, -2.883e+02f, 8.283e+01f },
           {7.500e+00f, -4.325e+02f, 1.213e+02f },
           {7.500e+00f, -6.487e+02f, 1.953e+02f },
           {7.500e+00f, -9.731e+02f, 2.860e+02f },
           {7.500e+00f, -1.460e+03f, 4.187e+02f },
           {7.500e+00f, -2.189e+03f, 6.130e+02f },
           {7.500e+00f, -3.284e+03f, 9.872e+02f },
           {7.500e+00f, -4.926e+03f, 1.445e+03f },
           {7.500e+00f, -7.389e+03f, 2.116e+03f },
           {7.500e+00f, -1.108e+04f, 3.098e+03f },
           {7.500e+00f, -1.663e+04f, 4.990e+03f },
           {1.125e+01f, -7.500e+00f, 6.318e+00f },
           {1.125e+01f, -1.125e+01f, 6.950e+00f },
           {1.125e+01f, -1.688e+01f, 7.645e+00f },
           {1.125e+01f, -2.531e+01f, 9.250e+00f },
           {1.125e+01f, -3.797e+01f, 1.231e+01f },
           {1.125e+01f, -5.695e+01f, 1.639e+01f },
           {1.125e+01f, -8.543e+01f, 2.399e+01f },
           {1.125e+01f, -1.281e+02f, 3.513e+01f },
           {1.125e+01f, -1.922e+02f, 5.657e+01f },
           {1.125e+01f, -2.883e+02f, 8.283e+01f },
           {1.125e+01f, -4.325e+02f, 1.213e+02f },
           {1.125e+01f, -6.487e+02f, 1.776e+02f },
           {1.125e+01f, -9.731e+02f, 2.860e+02f },
           {1.125e+01f, -1.460e+03f, 4.187e+02f },
           {1.125e+01f, -2.189e+03f, 6.130e+02f },
           {1.125e+01f, -3.284e+03f, 9.872e+02f },
           {1.125e+01f, -4.926e+03f, 1.445e+03f },
           {1.125e+01f, -7.389e+03f, 2.116e+03f },
           {1.125e+01f, -1.108e+04f, 3.098e+03f },
           {1.125e+01f, -1.663e+04f, 4.990e+03f },
           {1.688e+01f, -7.500e+00f, 4.747e+00f },
           {1.688e+01f, -1.125e+01f, 5.222e+00f },
           {1.688e+01f, -1.688e+01f, 5.744e+00f },
           {1.688e+01f, -2.531e+01f, 7.645e+00f },
           {1.688e+01f, -3.797e+01f, 1.018e+01f },
           {1.688e+01f, -5.695e+01f, 1.490e+01f },
           {1.688e+01f, -8.543e+01f, 2.181e+01f },
           {1.688e+01f, -1.281e+02f, 3.193e+01f },
           {1.688e+01f, -1.922e+02f, 5.143e+01f },
           {1.688e+01f, -2.883e+02f, 7.530e+01f },
           {1.688e+01f, -4.325e+02f, 1.213e+02f },
           {1.688e+01f, -6.487e+02f, 1.776e+02f },
           {1.688e+01f, -9.731e+02f, 2.600e+02f },
           {1.688e+01f, -1.460e+03f, 4.187e+02f },
           {1.688e+01f, -2.189e+03f, 6.130e+02f },
           {1.688e+01f, -3.284e+03f, 9.872e+02f },
           {1.688e+01f, -4.926e+03f, 1.445e+03f },
           {1.688e+01f, -7.389e+03f, 2.116e+03f },
           {1.688e+01f, -1.108e+04f, 3.098e+03f },
           {1.688e+01f, -1.663e+04f, 4.990e+03f },
           {2.531e+01f, -7.500e+00f, 3.242e+00f },
           {2.531e+01f, -1.125e+01f, 3.566e+00f },
           {2.531e+01f, -1.688e+01f, 4.315e+00f },
           {2.531e+01f, -2.531e+01f, 5.744e+00f },
           {2.531e+01f, -3.797e+01f, 7.645e+00f },
           {2.531e+01f, -5.695e+01f, 1.231e+01f },
           {2.531e+01f, -8.543e+01f, 1.803e+01f },
           {2.531e+01f, -1.281e+02f, 2.903e+01f },
           {2.531e+01f, -1.922e+02f, 4.676e+01f },
           {2.531e+01f, -2.883e+02f, 6.845e+01f },
           {2.531e+01f, -4.325e+02f, 1.102e+02f },
           {2.531e+01f, -6.487e+02f, 1.776e+02f },
           {2.531e+01f, -9.731e+02f, 2.600e+02f },
           {2.531e+01f, -1.460e+03f, 4.187e+02f },
           {2.531e+01f, -2.189e+03f, 6.130e+02f },
           {2.531e+01f, -3.284e+03f, 8.974e+02f },
           {2.531e+01f, -4.926e+03f, 1.445e+03f },
           {2.531e+01f, -7.389e+03f, 2.116e+03f },
           {2.531e+01f, -1.108e+04f, 3.098e+03f },
           {2.531e+01f, -1.663e+04f, 4.990e+03f },
           {3.797e+01f, -7.500e+00f, 2.214e+00f },
           {3.797e+01f, -1.125e+01f, 2.679e+00f },
           {3.797e+01f, -1.688e+01f, 3.242e+00f },
           {3.797e+01f, -2.531e+01f, 4.315e+00f },
           {3.797e+01f, -3.797e+01f, 6.318e+00f },
           {3.797e+01f, -5.695e+01f, 9.250e+00f },
           {3.797e+01f, -8.543e+01f, 1.490e+01f },
           {3.797e+01f, -1.281e+02f, 2.399e+01f },
           {3.797e+01f, -1.922e+02f, 3.864e+01f },
           {3.797e+01f, -2.883e+02f, 6.223e+01f },
           {3.797e+01f, -4.325e+02f, 1.002e+02f },
           {3.797e+01f, -6.487e+02f, 1.614e+02f },
           {3.797e+01f, -9.731e+02f, 2.600e+02f },
           {3.797e+01f, -1.460e+03f, 3.806e+02f },
           {3.797e+01f, -2.189e+03f, 6.130e+02f },
           {3.797e+01f, -3.284e+03f, 8.974e+02f },
           {3.797e+01f, -4.926e+03f, 1.445e+03f },
           {3.797e+01f, -7.389e+03f, 2.116e+03f },
           {3.797e+01f, -1.108e+04f, 3.098e+03f },
           {3.797e+01f, -1.663e+04f, 4.990e+03f },
           {5.695e+01f, -7.500e+00f, 1.513e+00f },
           {5.695e+01f, -1.125e+01f, 1.830e+00f },
           {5.695e+01f, -1.688e+01f, 2.214e+00f },
           {5.695e+01f, -2.531e+01f, 3.242e+00f },
           {5.695e+01f, -3.797e+01f, 4.315e+00f },
           {5.695e+01f, -5.695e+01f, 7.645e+00f },
           {5.695e+01f, -8.543e+01f, 1.231e+01f },
           {5.695e+01f, -1.281e+02f, 1.983e+01f },
           {5.695e+01f, -1.922e+02f, 3.513e+01f },
           {5.695e+01f, -2.883e+02f, 5.657e+01f },
           {5.695e+01f, -4.325e+02f, 9.111e+01f },
           {5.695e+01f, -6.487e+02f, 1.467e+02f },
           {5.695e+01f, -9.731e+02f, 2.363e+02f },
           {5.695e+01f, -1.460e+03f, 3.806e+02f },
           {5.695e+01f, -2.189e+03f, 5.572e+02f },
           {5.695e+01f, -3.284e+03f, 8.974e+02f },
           {5.695e+01f, -4.926e+03f, 1.314e+03f },
           {5.695e+01f, -7.389e+03f, 2.116e+03f },
           {5.695e+01f, -1.108e+04f, 3.098e+03f },
           {5.695e+01f, -1.663e+04f, 4.990e+03f },
           {8.543e+01f, -7.500e+00f, 1.250e+00f },
           {8.543e+01f, -1.125e+01f, 1.250e+00f },
           {8.543e+01f, -1.688e+01f, 1.513e+00f },
           {8.543e+01f, -2.531e+01f, 2.214e+00f },
           {8.543e+01f, -3.797e+01f, 3.242e+00f },
           {8.543e+01f, -5.695e+01f, 5.222e+00f },
           {8.543e+01f, -8.543e+01f, 9.250e+00f },
           {8.543e+01f, -1.281e+02f, 1.639e+01f },
           {8.543e+01f, -1.922e+02f, 2.903e+01f },
           {8.543e+01f, -2.883e+02f, 5.143e+01f },
           {8.543e+01f, -4.325e+02f, 8.283e+01f },
           {8.543e+01f, -6.487e+02f, 1.334e+02f },
           {8.543e+01f, -9.731e+02f, 2.148e+02f },
           {8.543e+01f, -1.460e+03f, 3.460e+02f },
           {8.543e+01f, -2.189e+03f, 5.572e+02f },
           {8.543e+01f, -3.284e+03f, 8.974e+02f },
           {8.543e+01f, -4.926e+03f, 1.314e+03f },
           {8.543e+01f, -7.389e+03f, 2.116e+03f },
           {8.543e+01f, -1.108e+04f, 3.098e+03f },
           {8.543e+01f, -1.663e+04f, 4.536e+03f },
           {1.281e+02f, -7.500e+00f, 1.250e+00f },
           {1.281e+02f, -1.125e+01f, 1.250e+00f },
           {1.281e+02f, -1.688e+01f, 1.250e+00f },
           {1.281e+02f, -2.531e+01f, 1.513e+00f },
           {1.281e+02f, -3.797e+01f, 2.214e+00f },
           {1.281e+02f, -5.695e+01f, 3.923e+00f },
           {1.281e+02f, -8.543e+01f, 6.950e+00f },
           {1.281e+02f, -1.281e+02f, 1.231e+01f },
           {1.281e+02f, -1.922e+02f, 2.181e+01f },
           {1.281e+02f, -2.883e+02f, 4.250e+01f },
           {1.281e+02f, -4.325e+02f, 6.845e+01f },
           {1.281e+02f, -6.487e+02f, 1.213e+02f },
           {1.281e+02f, -9.731e+02f, 1.953e+02f },
           {1.281e+02f, -1.460e+03f, 3.460e+02f },
           {1.281e+02f, -2.189e+03f, 5.572e+02f },
           {1.281e+02f, -3.284e+03f, 8.159e+02f },
           {1.281e+02f, -4.926e+03f, 1.314e+03f },
           {1.281e+02f, -7.389e+03f, 1.924e+03f },
           {1.281e+02f, -1.108e+04f, 3.098e+03f },
           {1.281e+02f, -1.663e+04f, 4.536e+03f },
           {1.922e+02f, -7.500e+00f, 1.250e+00f },
           {1.922e+02f, -1.125e+01f, 1.250e+00f },
           {1.922e+02f, -1.688e+01f, 1.250e+00f },
           {1.922e+02f, -2.531e+01f, 1.250e+00f },
           {1.922e+02f, -3.797e+01f, 1.664e+00f },
           {1.922e+02f, -5.695e+01f, 2.679e+00f },
           {1.922e+02f, -8.543e+01f, 5.222e+00f },
           {1.922e+02f, -1.281e+02f, 9.250e+00f },
           {1.922e+02f, -1.922e+02f, 1.803e+01f },
           {1.922e+02f, -2.883e+02f, 3.193e+01f },
           {1.922e+02f, -4.325e+02f, 5.657e+01f },
           {1.922e+02f, -6.487e+02f, 1.002e+02f },
           {1.922e+02f, -9.731e+02f, 1.776e+02f },
           {1.922e+02f, -1.460e+03f, 3.145e+02f },
           {1.922e+02f, -2.189e+03f, 5.066e+02f },
           {1.922e+02f, -3.284e+03f, 8.159e+02f },
           {1.922e+02f, -4.926e+03f, 1.194e+03f },
           {1.922e+02f, -7.389e+03f, 1.924e+03f },
           {1.922e+02f, -1.108e+04f, 3.098e+03f },
           {1.922e+02f, -1.663e+04f, 4.536e+03f },
           {2.883e+02f, -7.500e+00f, 1.250e+00f },
           {2.883e+02f, -1.125e+01f, 1.250e+00f },
           {2.883e+02f, -1.688e+01f, 1.250e+00f },
           {2.883e+02f, -2.531e+01f, 1.250e+00f },
           {2.883e+02f, -3.797e+01f, 1.250e+00f },
           {2.883e+02f, -5.695e+01f, 2.013e+00f },
           {2.883e+02f, -8.543e+01f, 3.566e+00f },
           {2.883e+02f, -1.281e+02f, 6.950e+00f },
           {2.883e+02f, -1.922e+02f, 1.354e+01f },
           {2.883e+02f, -2.883e+02f, 2.399e+01f },
           {2.883e+02f, -4.325e+02f, 4.676e+01f },
           {2.883e+02f, -6.487e+02f, 8.283e+01f },
           {2.883e+02f, -9.731e+02f, 1.614e+02f },
           {2.883e+02f, -1.460e+03f, 2.600e+02f },
           {2.883e+02f, -2.189e+03f, 4.605e+02f },
           {2.883e+02f, -3.284e+03f, 7.417e+02f },
           {2.883e+02f, -4.926e+03f, 1.194e+03f },
           {2.883e+02f, -7.389e+03f, 1.924e+03f },
           {2.883e+02f, -1.108e+04f, 2.817e+03f },
           {2.883e+02f, -1.663e+04f, 4.536e+03f },
           {4.325e+02f, -7.500e+00f, 1.250e+00f },
           {4.325e+02f, -1.125e+01f, 1.250e+00f },
           {4.325e+02f, -1.688e+01f, 1.250e+00f },
           {4.325e+02f, -2.531e+01f, 1.250e+00f },
           {4.325e+02f, -3.797e+01f, 1.250e+00f },
           {4.325e+02f, -5.695e+01f, 1.375e+00f },
           {4.325e+02f, -8.543e+01f, 2.436e+00f },
           {4.325e+02f, -1.281e+02f, 4.747e+00f },
           {4.325e+02f, -1.922e+02f, 9.250e+00f },
           {4.325e+02f, -2.883e+02f, 1.803e+01f },
           {4.325e+02f, -4.325e+02f, 3.513e+01f },
           {4.325e+02f, -6.487e+02f, 6.845e+01f },
           {4.325e+02f, -9.731e+02f, 1.334e+02f },
           {4.325e+02f, -1.460e+03f, 2.363e+02f },
           {4.325e+02f, -2.189e+03f, 3.806e+02f },
           {4.325e+02f, -3.284e+03f, 6.743e+02f },
           {4.325e+02f, -4.926e+03f, 1.086e+03f },
           {4.325e+02f, -7.389e+03f, 1.749e+03f },
           {4.325e+02f, -1.108e+04f, 2.817e+03f },
           {4.325e+02f, -1.663e+04f, 4.536e+03f },
           {6.487e+02f, -7.500e+00f, 1.250e+00f },
           {6.487e+02f, -1.125e+01f, 1.250e+00f },
           {6.487e+02f, -1.688e+01f, 1.250e+00f },
           {6.487e+02f, -2.531e+01f, 1.250e+00f },
           {6.487e+02f, -3.797e+01f, 1.250e+00f },
           {6.487e+02f, -5.695e+01f, 1.250e+00f },
           {6.487e+02f, -8.543e+01f, 1.664e+00f },
           {6.487e+02f, -1.281e+02f, 3.242e+00f },
           {6.487e+02f, -1.922e+02f, 6.950e+00f },
           {6.487e+02f, -2.883e+02f, 1.354e+01f },
           {6.487e+02f, -4.325e+02f, 2.639e+01f },
           {6.487e+02f, -6.487e+02f, 5.143e+01f },
           {6.487e+02f, -9.731e+02f, 1.002e+02f },
           {6.487e+02f, -1.460e+03f, 1.953e+02f },
           {6.487e+02f, -2.189e+03f, 3.460e+02f },
           {6.487e+02f, -3.284e+03f, 6.130e+02f },
           {6.487e+02f, -4.926e+03f, 9.872e+02f },
           {6.487e+02f, -7.389e+03f, 1.590e+03f },
           {6.487e+02f, -1.108e+04f, 2.561e+03f },
           {6.487e+02f, -1.663e+04f, 4.124e+03f },
           {9.731e+02f, -7.500e+00f, 1.250e+00f },
           {9.731e+02f, -1.125e+01f, 1.250e+00f },
           {9.731e+02f, -1.688e+01f, 1.250e+00f },
           {9.731e+02f, -2.531e+01f, 1.250e+00f },
           {9.731e+02f, -3.797e+01f, 1.250e+00f },
           {9.731e+02f, -5.695e+01f, 1.250e+00f },
           {9.731e+02f, -8.543e+01f, 1.250e+00f },
           {9.731e+02f, -1.281e+02f, 2.214e+00f },
           {9.731e+02f, -1.922e+02f, 4.747e+00f },
           {9.731e+02f, -2.883e+02f, 9.250e+00f },
           {9.731e+02f, -4.325e+02f, 1.983e+01f },
           {9.731e+02f, -6.487e+02f, 3.864e+01f },
           {9.731e+02f, -9.731e+02f, 7.530e+01f },
           {9.731e+02f, -1.460e+03f, 1.467e+02f },
           {9.731e+02f, -2.189e+03f, 2.860e+02f },
           {9.731e+02f, -3.284e+03f, 5.066e+02f },
           {9.731e+02f, -4.926e+03f, 8.974e+02f },
           {9.731e+02f, -7.389e+03f, 1.445e+03f },
           {9.731e+02f, -1.108e+04f, 2.561e+03f },
           {9.731e+02f, -1.663e+04f, 4.124e+03f },
           {1.460e+03f, -7.500e+00f, 1.250e+00f },
           {1.460e+03f, -1.125e+01f, 1.250e+00f },
           {1.460e+03f, -1.688e+01f, 1.250e+00f },
           {1.460e+03f, -2.531e+01f, 1.250e+00f },
           {1.460e+03f, -3.797e+01f, 1.250e+00f },
           {1.460e+03f, -5.695e+01f, 1.250e+00f },
           {1.460e+03f, -8.543e+01f, 1.250e+00f },
           {1.460e+03f, -1.281e+02f, 1.513e+00f },
           {1.460e+03f, -1.922e+02f, 3.242e+00f },
           {1.460e+03f, -2.883e+02f, 6.950e+00f },
           {1.460e+03f, -4.325e+02f, 1.354e+01f },
           {1.460e+03f, -6.487e+02f, 2.903e+01f },
           {1.460e+03f, -9.731e+02f, 5.657e+01f },
           {1.460e+03f, -1.460e+03f, 1.213e+02f },
           {1.460e+03f, -2.189e+03f, 2.148e+02f },
           {1.460e+03f, -3.284e+03f, 4.187e+02f },
           {1.460e+03f, -4.926e+03f, 7.417e+02f },
           {1.460e+03f, -7.389e+03f, 1.314e+03f },
           {1.460e+03f, -1.108e+04f, 2.328e+03f },
           {1.460e+03f, -1.663e+04f, 3.749e+03f },
           {2.189e+03f, -7.500e+00f, 1.250e+00f },
           {2.189e+03f, -1.125e+01f, 1.250e+00f },
           {2.189e+03f, -1.688e+01f, 1.250e+00f },
           {2.189e+03f, -2.531e+01f, 1.250e+00f },
           {2.189e+03f, -3.797e+01f, 1.250e+00f },
           {2.189e+03f, -5.695e+01f, 1.250e+00f },
           {2.189e+03f, -8.543e+01f, 1.250e+00f },
           {2.189e+03f, -1.281e+02f, 1.250e+00f },
           {2.189e+03f, -1.922e+02f, 2.214e+00f },
           {2.189e+03f, -2.883e+02f, 4.747e+00f },
           {2.189e+03f, -4.325e+02f, 9.250e+00f },
           {2.189e+03f, -6.487e+02f, 1.983e+01f },
           {2.189e+03f, -9.731e+02f, 4.250e+01f },
           {2.189e+03f, -1.460e+03f, 8.283e+01f },
           {2.189e+03f, -2.189e+03f, 1.776e+02f },
           {2.189e+03f, -3.284e+03f, 3.460e+02f },
           {2.189e+03f, -4.926e+03f, 6.130e+02f },
           {2.189e+03f, -7.389e+03f, 1.086e+03f },
           {2.189e+03f, -1.108e+04f, 1.924e+03f },
           {2.189e+03f, -1.663e+04f, 3.408e+03f },
           {3.284e+03f, -7.500e+00f, 1.250e+00f },
           {3.284e+03f, -1.125e+01f, 1.250e+00f },
           {3.284e+03f, -1.688e+01f, 1.250e+00f },
           {3.284e+03f, -2.531e+01f, 1.250e+00f },
           {3.284e+03f, -3.797e+01f, 1.250e+00f },
           {3.284e+03f, -5.695e+01f, 1.250e+00f },
           {3.284e+03f, -8.543e+01f, 1.250e+00f },
           {3.284e+03f, -1.281e+02f, 1.250e+00f },
           {3.284e+03f, -1.922e+02f, 1.513e+00f },
           {3.284e+03f, -2.883e+02f, 3.242e+00f },
           {3.284e+03f, -4.325e+02f, 6.318e+00f },
           {3.284e+03f, -6.487e+02f, 1.354e+01f },
           {3.284e+03f, -9.731e+02f, 2.903e+01f },
           {3.284e+03f, -1.460e+03f, 6.223e+01f },
           {3.284e+03f, -2.189e+03f, 1.334e+02f },
           {3.284e+03f, -3.284e+03f, 2.600e+02f },
           {3.284e+03f, -4.926e+03f, 5.066e+02f },
           {3.284e+03f, -7.389e+03f, 9.872e+02f },
           {3.284e+03f, -1.108e+04f, 1.749e+03f },
           {3.284e+03f, -1.663e+04f, 3.098e+03f },
           {4.926e+03f, -7.500e+00f, 1.250e+00f },
           {4.926e+03f, -1.125e+01f, 1.250e+00f },
           {4.926e+03f, -1.688e+01f, 1.250e+00f },
           {4.926e+03f, -2.531e+01f, 1.250e+00f },
           {4.926e+03f, -3.797e+01f, 1.250e+00f },
           {4.926e+03f, -5.695e+01f, 1.250e+00f },
           {4.926e+03f, -8.543e+01f, 1.250e+00f },
           {4.926e+03f, -1.281e+02f, 1.250e+00f },
           {4.926e+03f, -1.922e+02f, 1.250e+00f },
           {4.926e+03f, -2.883e+02f, 2.013e+00f },
           {4.926e+03f, -4.325e+02f, 4.315e+00f },
           {4.926e+03f, -6.487e+02f, 9.250e+00f },
           {4.926e+03f, -9.731e+02f, 2.181e+01f },
           {4.926e+03f, -1.460e+03f, 4.250e+01f },
           {4.926e+03f, -2.189e+03f, 9.111e+01f },
           {4.926e+03f, -3.284e+03f, 1.953e+02f },
           {4.926e+03f, -4.926e+03f, 3.806e+02f },
           {4.926e+03f, -7.389e+03f, 7.417e+02f },
           {4.926e+03f, -1.108e+04f, 1.445e+03f },
           {4.926e+03f, -1.663e+04f, 2.561e+03f },
           {7.389e+03f, -7.500e+00f, 1.250e+00f },
           {7.389e+03f, -1.125e+01f, 1.250e+00f },
           {7.389e+03f, -1.688e+01f, 1.250e+00f },
           {7.389e+03f, -2.531e+01f, 1.250e+00f },
           {7.389e+03f, -3.797e+01f, 1.250e+00f },
           {7.389e+03f, -5.695e+01f, 1.250e+00f },
           {7.389e+03f, -8.543e+01f, 1.250e+00f },
           {7.389e+03f, -1.281e+02f, 1.250e+00f },
           {7.389e+03f, -1.922e+02f, 1.250e+00f },
           {7.389e+03f, -2.883e+02f, 1.375e+00f },
           {7.389e+03f, -4.325e+02f, 2.947e+00f },
           {7.389e+03f, -6.487e+02f, 6.318e+00f },
           {7.389e+03f, -9.731e+02f, 1.490e+01f },
           {7.389e+03f, -1.460e+03f, 3.193e+01f },
           {7.389e+03f, -2.189e+03f, 6.845e+01f },
           {7.389e+03f, -3.284e+03f, 1.334e+02f },
           {7.389e+03f, -4.926e+03f, 2.860e+02f },
           {7.389e+03f, -7.389e+03f, 5.572e+02f },
           {7.389e+03f, -1.108e+04f, 1.086e+03f },
           {7.389e+03f, -1.663e+04f, 2.116e+03f },
           {1.108e+04f, -7.500e+00f, 1.250e+00f },
           {1.108e+04f, -1.125e+01f, 1.250e+00f },
           {1.108e+04f, -1.688e+01f, 1.250e+00f },
           {1.108e+04f, -2.531e+01f, 1.250e+00f },
           {1.108e+04f, -3.797e+01f, 1.250e+00f },
           {1.108e+04f, -5.695e+01f, 1.250e+00f },
           {1.108e+04f, -8.543e+01f, 1.250e+00f },
           {1.108e+04f, -1.281e+02f, 1.250e+00f },
           {1.108e+04f, -1.922e+02f, 1.250e+00f },
           {1.108e+04f, -2.883e+02f, 1.250e+00f },
           {1.108e+04f, -4.325e+02f, 2.013e+00f },
           {1.108e+04f, -6.487e+02f, 4.315e+00f },
           {1.108e+04f, -9.731e+02f, 1.018e+01f },
           {1.108e+04f, -1.460e+03f, 2.181e+01f },
           {1.108e+04f, -2.189e+03f, 4.676e+01f },
           {1.108e+04f, -3.284e+03f, 1.002e+02f },
           {1.108e+04f, -4.926e+03f, 2.148e+02f },
           {1.108e+04f, -7.389e+03f, 4.187e+02f },
           {1.108e+04f, -1.108e+04f, 8.974e+02f },
           {1.108e+04f, -1.663e+04f, 1.749e+03f },
           {1.663e+04f, -7.500e+00f, 1.250e+00f },
           {1.663e+04f, -1.125e+01f, 1.250e+00f },
           {1.663e+04f, -1.688e+01f, 1.250e+00f },
           {1.663e+04f, -2.531e+01f, 1.250e+00f },
           {1.663e+04f, -3.797e+01f, 1.250e+00f },
           {1.663e+04f, -5.695e+01f, 1.250e+00f },
           {1.663e+04f, -8.543e+01f, 1.250e+00f },
           {1.663e+04f, -1.281e+02f, 1.250e+00f },
           {1.663e+04f, -1.922e+02f, 1.250e+00f },
           {1.663e+04f, -2.883e+02f, 1.250e+00f },
           {1.663e+04f, -4.325e+02f, 1.375e+00f },
           {1.663e+04f, -6.487e+02f, 2.947e+00f },
           {1.663e+04f, -9.731e+02f, 6.318e+00f },
           {1.663e+04f, -1.460e+03f, 1.490e+01f },
           {1.663e+04f, -2.189e+03f, 3.193e+01f },
           {1.663e+04f, -3.284e+03f, 6.845e+01f },
           {1.663e+04f, -4.926e+03f, 1.467e+02f },
           {1.663e+04f, -7.389e+03f, 3.145e+02f },
           {1.663e+04f, -1.108e+04f, 6.743e+02f },
           {1.663e+04f, -1.663e+04f, 1.314e+03f },
        };
        if ((a > 1.663e+04) || (-b > 1.663e+04))
           return z > -b;  // Way overly conservative?
        if (a < data[0][0])
           return false;
        int index = 0;
        while (data[index][0] < a)
           ++index;
        if(a != data[index][0])
           --index;
        while ((data[index][1] < b) && (data[index][2] > 1.25))
           --index;
        ++index;
        BOOST_MATH_ASSERT(a > data[index][0]);
        BOOST_MATH_ASSERT(-b < -data[index][1]);
        return z > data[index][2];
     }
     template <class T, class Policy>
     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)
     {
        BOOST_MATH_STD_USING
        //
        // Get the function ratio, M(a+1, b+1, z)/M(a,b,z):
        //
        std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
        boost::math::detail::hypergeometric_1F1_recurrence_a_and_b_coefficients<T> coef(a, b, z);
        T ratio = 1 / boost::math::tools::function_ratio_from_forwards_recurrence(coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
        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);
        //
        // We can't normalise via the Wronksian as the subtraction in the Wronksian will suffer an exquisite amount of cancellation - 
        // potentially many hundreds of digits in this region.  However, if forwards iteration is stable at this point
        // 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
        // use a reference value just over the origin to normalise:
        //
        long long scale = 0;
        int steps = itrunc(ceil(-b));
        T reference_value = hypergeometric_1F1_imp(T(a + steps), T(b + steps), z, pol, log_scaling);
        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);
        log_scaling -= scale;
        if ((fabs(reference_value) < 1) && (fabs(reference_value) < tools::min_value<T>() * fabs(found)))
        {
           // Possible underflow, rescale
           long long s = lltrunc(tools::log_max_value<T>());
           log_scaling -= s;
           reference_value *= exp(T(s));
        }
        else if ((fabs(found) < 1) && (fabs(reference_value) > tools::max_value<T>() * fabs(found)))
        {
           // Overflow, rescale:
           long long s = lltrunc(tools::log_max_value<T>());
           log_scaling += s;
           reference_value /= exp(T(s));
        }
        return reference_value / found;
     }



     //
     // This next version is largely the same as above, but calculates the ratio for the b recurrence relation
     // which has a larger area of stability than the ab recurrence when a,b < 0.  We can then use a single
     // recurrence step to convert this to the ratio for the ab recursion and proceed largely as before.
     // The routine is quite insensitive to the size of z, but requires |a| < |5b| for accuracy.
     // Fortunately the accuracy outside this domain falls off steadily rather than suddenly switching
     // to a different behaviour.
     //
     template <class T, class Policy>
     T hypergeometric_1F1_from_function_ratio_negative_ab(const T& a, const T& b, const T& z, const Policy& pol, long long& log_scaling)
     {
        BOOST_MATH_STD_USING
        //
        // Get the function ratio, M(a+1, b+1, z)/M(a,b,z):
        //
        std::uintmax_t max_iter = boost::math::policies::get_max_series_iterations<Policy>();
        boost::math::detail::hypergeometric_1F1_recurrence_b_coefficients<T> coef(a, b + 1, z);
        T ratio = boost::math::tools::function_ratio_from_backwards_recurrence(coef, boost::math::policies::get_epsilon<T, Policy>(), max_iter);
        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);
        //
        // We need to use A&S 13.4.3 to convert a ratio for M(a, b + 1, z) / M(a, b, z)
        // to M(a+1, b+1, z) / M(a, b, z)
        //
        // We have:        (a-b)M(a, b+1, z) - aM(a+1, b+1, z) + bM(a, b, z) = 0
        // 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
        //
        ratio = ((a - b) * ratio + b) / a;
        //
        // Let M2 = M(1+a-b, 2-b, z)
        // This is going to be a mighty big number:
        //
        long long local_scaling = 0;
        T M2 = boost::math::detail::hypergeometric_1F1_imp(T(1 + a - b), T(2 - b), z, pol, local_scaling);
        log_scaling -= local_scaling; // all the M2 terms are in the denominator
        //
        // Let M3 = M(1+a-b + 1, 2-b + 1, z)
        // We don't use the ratio to get this as it's not clear that it's reliable:
        //
        long long local_scaling2 = 0;
        T M3 = boost::math::detail::hypergeometric_1F1_imp(T(2 + a - b), T(3 - b), z, pol, local_scaling2);
        //
        // M2 and M3 must be identically scaled:
        //
        if (local_scaling != local_scaling2)
        {
           M3 *= exp(T(local_scaling2 - local_scaling));
        }
        //
        // Get the RHS of the Wronksian:
        //
        long long fz = lltrunc(z);
        log_scaling += fz;
        T rhs = (1 - b) * exp(z - fz);
        //
        // We need to divide that by:
        // [(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) ]
        // Note that at this stage, both M3 and M2 are scaled by exp(local_scaling).
        //
        T lhs = (a - b + 1) * z * M3 / (2 - b) + (1 - b) * M2 - a * z * ratio * M2 / b;

        return rhs / lhs;
     }

  } } } // namespaces

#endif // BOOST_HYPERGEOMETRIC_1F1_BY_RATIOS_HPP_