special_functions/detail/hypergeometric_1F1_scaled_series.hpp

0001 ///////////////////////////////////////////////////////////////////////////////
0002 //  Copyright 2017 John Maddock
0003 //  Distributed under the Boost
0004 //  Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006 //
0007 #ifndef BOOST_MATH_HYPERGEOMETRIC_1F1_SCALED_SERIES_HPP
0008 #define BOOST_MATH_HYPERGEOMETRIC_1F1_SCALED_SERIES_HPP
0009
0010 #include <array>
0011 #include <cstdint>
0012
0013   namespace boost{ namespace math{ namespace detail{
0014
0015      template <class T, class Policy>
0016      T hypergeometric_1F1_scaled_series(const T& a, const T& b, T z, const Policy& pol, const char* function)
0017      {
0018         BOOST_MATH_STD_USING
0019         //
0020         // Result is returned scaled by e^-z.
0021         // Whenever the terms start becoming too large, we scale by some factor e^-n
0022         // and keep track of the integer scaling factor n.  At the end we can perform
0023         // an exact subtraction of n from z and scale the result:
0024         //
0025         T sum(0), term(1), upper_limit(sqrt(boost::math::tools::max_value<T>())), diff;
0026         unsigned n = 0;
0027         long long log_scaling_factor = 1 - lltrunc(boost::math::tools::log_max_value<T>());
0028         T scaling_factor = exp(T(log_scaling_factor));
0029         std::intmax_t current_scaling = 0;
0030
0031         do
0032         {
0033            sum += term;
0034            if (sum >= upper_limit)
0035            {
0036               sum *= scaling_factor;
0037               term *= scaling_factor;
0038               current_scaling += log_scaling_factor;
0039            }
0040            term *= (((a + n) / ((b + n) * (n + 1))) * z);
0041            if (n > boost::math::policies::get_max_series_iterations<Policy>())
0042               return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
0043            ++n;
0044            diff = fabs(term / sum);
0045         } while (diff > boost::math::policies::get_epsilon<T, Policy>());
0046
0047         z = -z - current_scaling;
0048         while (z < log_scaling_factor)
0049         {
0050            z -= log_scaling_factor;
0051            sum *= scaling_factor;
0052         }
0053         return sum * exp(z);
0054      }
0055
0056
0057
0058   } } } // namespaces
0059
0060 #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_SCALED_SERIES_HPP