boost/math/special_functions/detail/hypergeometric_asym.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2014 Anton Bikineev
   3 //  Copyright 2014 Christopher Kormanyos
   4 //  Copyright 2014 John Maddock
   5 //  Copyright 2014 Paul Bristow
   6 //  Distributed under the Boost
   7 //  Software License, Version 1.0. (See accompanying file
   8 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   9 //
  10 #ifndef BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP
  11 #define BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP
  12
  13 #include <boost/math/special_functions/gamma.hpp>
  14 #include <boost/math/special_functions/hypergeometric_2F0.hpp>
  15
  16 #ifdef BOOST_MSVC
  17 #pragma warning(push)
  18 #pragma warning(disable:4127)
  19 #endif
  20
  21   namespace boost { namespace math {
  22
  23   namespace detail {
  24
  25      //
  26      // Asymptotic series based on https://dlmf.nist.gov/13.7#E1
  27      //
  28      // Note that a and b must not be negative integers, in addition
  29      // we require z > 0 and so apply Kummer's relation for z < 0.
  30      //
  31      template <class T, class Policy>
  32      inline T hypergeometric_1F1_asym_large_z_series(T a, const T& b, T z, const Policy& pol, int& log_scaling)
  33      {
  34         BOOST_MATH_STD_USING
  35         static const char* function = "boost::math::hypergeometric_1F1_asym_large_z_series<%1%>(%1%, %1%, %1%)";
  36         T prefix;
  37         int e, s;
  38         if (z < 0)
  39         {
  40            a = b - a;
  41            z = -z;
  42            prefix = 1;
  43         }
  44         else
  45         {
  46            e = z > INT_MAX ? INT_MAX : itrunc(z, pol);
  47            log_scaling += e;
  48            prefix = exp(z - e);
  49         }
  50         if ((fabs(a) < 10) && (fabs(b) < 10))
  51         {
  52            prefix *= pow(z, a) * pow(z, -b) * boost::math::tgamma(b, pol) / boost::math::tgamma(a, pol);
  53         }
  54         else
  55         {
  56            T t = log(z) * (a - b);
  57            e = itrunc(t, pol);
  58            log_scaling += e;
  59            prefix *= exp(t - e);
  60
  61            t = boost::math::lgamma(b, &s, pol);
  62            e = itrunc(t, pol);
  63            log_scaling += e;
  64            prefix *= s * exp(t - e);
  65
  66            t = boost::math::lgamma(a, &s, pol);
  67            e = itrunc(t, pol);
  68            log_scaling -= e;
  69            prefix /= s * exp(t - e);
  70         }
  71         //
  72         // Checked 2F0:
  73         //
  74         unsigned k = 0;
  75         T a1_poch(1 - a);
  76         T a2_poch(b - a);
  77         T z_mult(1 / z);
  78         T sum = 0;
  79         T abs_sum = 0;
  80         T term = 1;
  81         T last_term = 0;
  82         do
  83         {
  84            sum += term;
  85            last_term = term;
  86            abs_sum += fabs(sum);
  87            term *= a1_poch * a2_poch * z_mult;
  88            term /= ++k;
  89            a1_poch += 1;
  90            a2_poch += 1;
  91            if (fabs(sum) * boost::math::policies::get_epsilon<T, Policy>() > fabs(term))
  92               break;
  93            if(fabs(sum) / abs_sum < boost::math::policies::get_epsilon<T, Policy>())
  94               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%",
  95                  prefix * sum, Policy());
  96            if(k > boost::math::policies::get_max_series_iterations<Policy>())
  97               return boost::math::policies::raise_evaluation_error<T>(function, "1F1: Unable to locate solution in a reasonable time:"
  98                  " large-z asymptotic approximation.  Current best guess is %1%", prefix * sum, Policy());
  99            if((k > 10) && (fabs(term) > fabs(last_term)))
 100               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());
 101         } while (true);
 102
 103         return prefix * sum;
 104      }
 105
 106
 107   // experimental range
 108   template <class T, class Policy>
 109   inline bool hypergeometric_1F1_asym_region(const T& a, const T& b, const T& z, const Policy&)
 110   {
 111     BOOST_MATH_STD_USING
 112     int half_digits = policies::digits<T, Policy>() / 2;
 113     bool in_region = false;
 114
 115     if (fabs(a) < 0.001f)
 116        return false; // Haven't been able to make this work, why not?  TODO!
 117
 118     //
 119     // We use the following heuristic, if after we have had half_digits terms
 120     // of the 2F0 series, we require terms to be decreasing in size by a factor
 121     // of at least 0.7.  Assuming the earlier terms were converging much faster
 122     // than this, then this should be enough to achieve convergence before the
 123     // series shoots off to infinity.
 124     //
 125     if (z > 0)
 126     {
 127        T one_minus_a = 1 - a;
 128        T b_minus_a = b - a;
 129        if (fabs((one_minus_a + half_digits) * (b_minus_a + half_digits) / (half_digits * z)) < 0.7)
 130        {
 131           in_region = true;
 132           //
 133           // double check that we are not divergent at the start if a,b < 0:
 134           //
 135           if ((one_minus_a < 0) || (b_minus_a < 0))
 136           {
 137              if (fabs(one_minus_a * b_minus_a / z) > 0.5)
 138                 in_region = false;
 139           }
 140        }
 141     }
 142     else if (fabs((1 - (b - a) + half_digits) * (a + half_digits) / (half_digits * z)) < 0.7)
 143     {
 144        if ((floor(b - a) == (b - a)) && (b - a < 0))
 145           return false;  // Can't have a negative integer b-a.
 146        in_region = true;
 147        //
 148        // double check that we are not divergent at the start if a,b < 0:
 149        //
 150        T a1 = 1 - (b - a);
 151        if ((a1 < 0) || (a < 0))
 152        {
 153           if (fabs(a1 * a / z) > 0.5)
 154              in_region = false;
 155        }
 156     }
 157     //
 158     // Check for a and b negative integers as these aren't supported by the approximation:
 159     //
 160     if (in_region)
 161     {
 162        if ((a < 0) && (floor(a) == a))
 163           in_region = false;
 164        if ((b < 0) && (floor(b) == b))
 165           in_region = false;
 166        if (fabs(z) < 40)
 167           in_region = false;
 168     }
 169     return in_region;
 170   }
 171
 172   } } } // namespaces
 173
 174 #ifdef BOOST_MSVC
 175 #pragma warning(pop)
 176 #endif
 177
 178 #endif // BOOST_MATH_HYPERGEOMETRIC_ASYM_HPP