ceph/src/boost/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 <array>
  11 #include <cstdint>
  12
  13   namespace boost{ namespace math{ namespace detail{
  14
  15      template <class T, class Policy>
  16      T hypergeometric_1F1_scaled_series(const T& a, const T& b, T z, const Policy& pol, const char* function)
  17      {
  18         BOOST_MATH_STD_USING
  19         //
  20         // Result is returned scaled by e^-z.
  21         // Whenever the terms start becoming too large, we scale by some factor e^-n
  22         // and keep track of the integer scaling factor n.  At the end we can perform
  23         // an exact subtraction of n from z and scale the result:
  24         //
  25         T sum(0), term(1), upper_limit(sqrt(boost::math::tools::max_value<T>())), diff;
  26         unsigned n = 0;
  27         long long log_scaling_factor = 1 - lltrunc(boost::math::tools::log_max_value<T>());
  28         T scaling_factor = exp(T(log_scaling_factor));
  29         std::intmax_t current_scaling = 0;
  30
  31         do
  32         {
  33            sum += term;
  34            if (sum >= upper_limit)
  35            {
  36               sum *= scaling_factor;
  37               term *= scaling_factor;
  38               current_scaling += log_scaling_factor;
  39            }
  40            term *= (((a + n) / ((b + n) * (n + 1))) * z);
  41            if (n > boost::math::policies::get_max_series_iterations<Policy>())
  42               return boost::math::policies::raise_evaluation_error(function, "Series did not converge, best value is %1%", sum, pol);
  43            ++n;
  44            diff = fabs(term / sum);
  45         } while (diff > boost::math::policies::get_epsilon<T, Policy>());
  46
  47         z = -z - current_scaling;
  48         while (z < log_scaling_factor)
  49         {
  50            z -= log_scaling_factor;
  51            sum *= scaling_factor;
  52         }
  53         return sum * exp(z);
  54      }
  55
  56
  57
  58   } } } // namespaces
  59
  60 #endif // BOOST_MATH_HYPERGEOMETRIC_1F1_SCALED_SERIES_HPP