boost/math/special_functions/detail/hypergeometric_1F1_scaled_series.hpp

   1 ///////////////////////////////////////////////////////////////////////////////
   2 //  Copyright 2017 John Maddock
   3 //  Distributed under the Boost
   4 //  Software License, Version 1.0. (See accompanying file
   5 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   6 //
   7 #ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_SCALED_SERIES_HPP
   8 #define BOOST_MATH_HYPERGEOMETRIC_1F1_SCALED_SERIES_HPP
   9
  10 #include <boost/array.hpp>
  11
  12   namespace boost{ namespace math{ namespace detail{
  13
  14      template <class T, class Policy>
  15      T hypergeometric_1F1_scaled_series(const T& a, const T& b, T z, const Policy& pol, const char* function)
  16      {
  17         BOOST_MATH_STD_USING
  18         //
  19         // Result is returned scaled by e^-z.
  20         // Whenever the terms start becoming too large, we scale by some factor e^-n
  21         // and keep track of the integer scaling factor n.  At the end we can perform
  22         // an exact subtraction of n from z and scale the result:
  23         //
  24         T sum(0), term(1), upper_limit(sqrt(boost::math::tools::max_value<T>())), diff;
  25         unsigned n = 0;
  26         boost::intmax_t log_scaling_factor = 1 - itrunc(boost::math::tools::log_max_value<T>());
  27         T scaling_factor = exp(T(log_scaling_factor));
  28         boost::intmax_t current_scaling = 0;
  29
  30         do
  31         {
  32            sum += term;
  33            if (sum >= upper_limit)
  34            {
  35               sum *= scaling_factor;
  36               term *= scaling_factor;
  37               current_scaling += log_scaling_factor;
  38            }
  39            term *= (((a + n) / ((b + n) * (n + 1))) * z);
  40            if (n > boost::math::policies::get_max_series_iterations<Policy>())
  41               return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
  42            ++n;
  43            diff = fabs(term / sum);
  44         } while (diff > boost::math::policies::get_epsilon<T, Policy>());
  45
  46         z = -z - current_scaling;
  47         while (z < log_scaling_factor)
  48         {
  49            z -= log_scaling_factor;
  50            sum *= scaling_factor;
  51         }
  52         return sum * exp(z);
  53      }
  54
  55
  56
  57   } } } // namespaces
  58
  59 #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_SCALED_SERIES_HPP