special_functions/detail/hypergeometric_asym.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_ASYM_HPP
0011 #define BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP
0012
0013 #include <boost/math/special_functions/gamma.hpp>
0014 #include <boost/math/special_functions/hypergeometric_2F0.hpp>
0015
0016 #ifdef _MSC_VER
0017 #pragma warning(push)
0018 #pragma warning(disable:4127)
0019 #endif
0020
0021   namespace boost { namespace math {
0022
0023   namespace detail {
0024
0025      //
0026      // Asymptotic series based on https://dlmf.nist.gov/13.7#E1
0027      //
0028      // Note that a and b must not be negative integers, in addition
0029      // we require z > 0 and so apply Kummer's relation for z < 0.
0030      //
0031      template <class T, class Policy>
0032      inline T hypergeometric_1F1_asym_large_z_series(T a, const T& b, T z, const Policy& pol, long long& log_scaling)
0033      {
0034         BOOST_MATH_STD_USING
0035         static const char* function = "boost::math::hypergeometric_1F1_asym_large_z_series<%1%>(%1%, %1%, %1%)";
0036         T prefix;
0037         long long e;
0038         int s;
0039         if (z < 0)
0040         {
0041            a = b - a;
0042            z = -z;
0043            prefix = 1;
0044         }
0045         else
0046         {
0047            e = z > static_cast<T>((std::numeric_limits<long long>::max)()) ? (std::numeric_limits<long long>::max)() : lltrunc(z, pol);
0048            log_scaling += e;
0049            prefix = exp(z - e);
0050         }
0051         if ((fabs(a) < 10) && (fabs(b) < 10))
0052         {
0053            prefix *= pow(z, a) * pow(z, -b) * boost::math::tgamma(b, pol) / boost::math::tgamma(a, pol);
0054         }
0055         else
0056         {
0057            T t = log(z) * (a - b);
0058            e = lltrunc(t, pol);
0059            log_scaling += e;
0060            prefix *= exp(t - e);
0061
0062            t = boost::math::lgamma(b, &s, pol);
0063            e = lltrunc(t, pol);
0064            log_scaling += e;
0065            prefix *= s * exp(t - e);
0066
0067            t = boost::math::lgamma(a, &s, pol);
0068            e = lltrunc(t, pol);
0069            log_scaling -= e;
0070            prefix /= s * exp(t - e);
0071         }
0072         //
0073         // Checked 2F0:
0074         //
0075         unsigned k = 0;
0076         T a1_poch(1 - a);
0077         T a2_poch(b - a);
0078         T z_mult(1 / z);
0079         T sum = 0;
0080         T abs_sum = 0;
0081         T term = 1;
0082         T last_term = 0;
0083         do
0084         {
0085            sum += term;
0086            last_term = term;
0087            abs_sum += fabs(sum);
0088            term *= a1_poch * a2_poch * z_mult;
0089            term /= ++k;
0090            a1_poch += 1;
0091            a2_poch += 1;
0092            if (fabs(sum) * boost::math::policies::get_epsilon<T, Policy>() > fabs(term))
0093               break;
0094            if(fabs(sum) / abs_sum < boost::math::policies::get_epsilon<T, Policy>())
0095               return boost::math::policies::raise_evaluation_error<T>(function, "Large-z asymptotic approximation to 1F1 has destroyed all the digits in the result due to cancellation.  Current best guess is %1%",
0096                  prefix * sum, Policy());
0097            if(k > boost::math::policies::get_max_series_iterations<Policy>())
0098               return boost::math::policies::raise_evaluation_error<T>(function, "1F1: Unable to locate solution in a reasonable time:"
0099                  " large-z asymptotic approximation.  Current best guess is %1%", prefix * sum, Policy());
0100            if((k > 10) && (fabs(term) > fabs(last_term)))
0101               return boost::math::policies::raise_evaluation_error<T>(function, "Large-z asymptotic approximation to 1F1 is divergent.  Current best guess is %1%", prefix * sum, Policy());
0102         } while (true);
0103
0104         return prefix * sum;
0105      }
0106
0107
0108   // experimental range
0109   template <class T, class Policy>
0110   inline bool hypergeometric_1F1_asym_region(const T& a, const T& b, const T& z, const Policy&)
0111   {
0112     BOOST_MATH_STD_USING
0113     int half_digits = policies::digits<T, Policy>() / 2;
0114     bool in_region = false;
0115
0116     if (fabs(a) < 0.001f)
0117        return false; // Haven't been able to make this work, why not?  TODO!
0118
0119     //
0120     // We use the following heuristic, if after we have had half_digits terms
0121     // of the 2F0 series, we require terms to be decreasing in size by a factor
0122     // of at least 0.7.  Assuming the earlier terms were converging much faster
0123     // than this, then this should be enough to achieve convergence before the
0124     // series shoots off to infinity.
0125     //
0126     if (z > 0)
0127     {
0128        T one_minus_a = 1 - a;
0129        T b_minus_a = b - a;
0130        if (fabs((one_minus_a + half_digits) * (b_minus_a + half_digits) / (half_digits * z)) < 0.7)
0131        {
0132           in_region = true;
0133           //
0134           // double check that we are not divergent at the start if a,b < 0:
0135           //
0136           if ((one_minus_a < 0) || (b_minus_a < 0))
0137           {
0138              if (fabs(one_minus_a * b_minus_a / z) > 0.5)
0139                 in_region = false;
0140           }
0141        }
0142     }
0143     else if (fabs((1 - (b - a) + half_digits) * (a + half_digits) / (half_digits * z)) < 0.7)
0144     {
0145        if ((floor(b - a) == (b - a)) && (b - a < 0))
0146           return false;  // Can't have a negative integer b-a.
0147        in_region = true;
0148        //
0149        // double check that we are not divergent at the start if a,b < 0:
0150        //
0151        T a1 = 1 - (b - a);
0152        if ((a1 < 0) || (a < 0))
0153        {
0154           if (fabs(a1 * a / z) > 0.5)
0155              in_region = false;
0156        }
0157     }
0158     //
0159     // Check for a and b negative integers as these aren't supported by the approximation:
0160     //
0161     if (in_region)
0162     {
0163        if ((a < 0) && (floor(a) == a))
0164           in_region = false;
0165        if ((b < 0) && (floor(b) == b))
0166           in_region = false;
0167        if (fabs(z) < 40)
0168           in_region = false;
0169     }
0170     return in_region;
0171   }
0172
0173   } } } // namespaces
0174
0175 #ifdef _MSC_VER
0176 #pragma warning(pop)
0177 #endif
0178
0179 #endif // BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP