boost/math/special_functions/hypergeometric_pFq.hpp

   1
   2 ///////////////////////////////////////////////////////////////////////////////
   3 //  Copyright 2018 John Maddock
   4 //  Distributed under the Boost
   5 //  Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP
   9 #define BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP
  10
  11 #include <boost/config.hpp>
  12
  13 #if defined(BOOST_NO_CXX11_AUTO_DECLARATIONS) || defined(BOOST_NO_CXX11_LAMBDAS) || defined(BOOST_NO_CXX11_UNIFIED_INITIALIZATION_SYNTAX) || defined(BOOST_NO_CXX11_HDR_INITIALIZER_LIST)
  14 # error "hypergeometric_pFq requires a C++11 compiler"
  15 #endif
  16
  17 #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
  18 #include <boost/chrono.hpp>
  19 #include <initializer_list>
  20
  21 namespace boost {
  22    namespace math {
  23
  24       namespace detail {
  25
  26          struct pFq_termination_exception : public std::runtime_error
  27          {
  28             pFq_termination_exception(const char* p) : std::runtime_error(p) {}
  29          };
  30
  31          struct timed_iteration_terminator
  32          {
  33             timed_iteration_terminator(boost::uintmax_t i, double t) : max_iter(i), max_time(t), start_time(boost::chrono::system_clock::now()) {}
  34
  35             bool operator()(boost::uintmax_t iter)const
  36             {
  37                if (iter > max_iter)
  38                   boost::throw_exception(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted iterations."));
  39                if (boost::chrono::duration<double>(boost::chrono::system_clock::now() - start_time).count() > max_time)
  40                   boost::throw_exception(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted evaluation time."));
  41                return false;
  42             }
  43
  44             boost::uintmax_t max_iter;
  45             double max_time;
  46             boost::chrono::system_clock::time_point start_time;
  47          };
  48
  49       }
  50
  51       template <class Seq, class Real, class Policy>
  52       inline typename tools::promote_args<Real, typename Seq::value_type>::type hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error, const Policy& pol)
  53       {
  54          typedef typename tools::promote_args<Real, typename Seq::value_type>::type result_type;
  55          typedef typename policies::evaluation<result_type, Policy>::type value_type;
  56          typedef typename policies::normalise<
  57             Policy,
  58             policies::promote_float<false>,
  59             policies::promote_double<false>,
  60             policies::discrete_quantile<>,
  61             policies::assert_undefined<> >::type forwarding_policy;
  62
  63          BOOST_MATH_STD_USING
  64
  65          int scale = 0;
  66          std::pair<value_type, value_type> r = boost::math::detail::hypergeometric_pFq_checked_series_impl(aj, bj, value_type(z), pol, boost::math::detail::iteration_terminator(boost::math::policies::get_max_series_iterations<forwarding_policy>()), scale);
  67          r.first *= exp(Real(scale));
  68          r.second *= exp(Real(scale));
  69          if (p_abs_error)
  70             *p_abs_error = static_cast<Real>(r.second) * boost::math::tools::epsilon<Real>();
  71          return policies::checked_narrowing_cast<result_type, Policy>(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)");
  72       }
  73
  74       template <class Seq, class Real>
  75       inline typename tools::promote_args<Real, typename Seq::value_type>::type hypergeometric_pFq(const Seq& aj, const Seq& bj, const Real& z, Real* p_abs_error = 0)
  76       {
  77          return hypergeometric_pFq(aj, bj, z, p_abs_error, boost::math::policies::policy<>());
  78       }
  79
  80       template <class R, class Real, class Policy>
  81       inline typename tools::promote_args<Real, R>::type hypergeometric_pFq(const std::initializer_list<R>& aj, const std::initializer_list<R>& bj, const Real& z, Real* p_abs_error, const Policy& pol)
  82       {
  83          return hypergeometric_pFq<std::initializer_list<R>, Real, Policy>(aj, bj, z, p_abs_error, pol);
  84       }
  85
  86       template <class R, class Real>
  87       inline typename tools::promote_args<Real, R>::type  hypergeometric_pFq(const std::initializer_list<R>& aj, const std::initializer_list<R>& bj, const Real& z, Real* p_abs_error = 0)
  88       {
  89          return hypergeometric_pFq<std::initializer_list<R>, Real>(aj, bj, z, p_abs_error);
  90       }
  91
  92       template <class T>
  93       struct scoped_precision
  94       {
  95          scoped_precision(unsigned p)
  96          {
  97             old_p = T::default_precision();
  98             T::default_precision(p);
  99          }
 100          ~scoped_precision()
 101          {
 102             T::default_precision(old_p);
 103          }
 104          unsigned old_p;
 105       };
 106
 107       template <class Seq, class Real, class Policy>
 108       Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, Real z, unsigned digits10, double timeout, const Policy& pol)
 109       {
 110          unsigned current_precision = digits10 + 5;
 111
 112          for (auto ai = aj.begin(); ai != aj.end(); ++ai)
 113          {
 114             current_precision = (std::max)(current_precision, ai->precision());
 115          }
 116          for (auto bi = bj.begin(); bi != bj.end(); ++bi)
 117          {
 118             current_precision = (std::max)(current_precision, bi->precision());
 119          }
 120          current_precision = (std::max)(current_precision, z.precision());
 121
 122          Real r, norm;
 123          std::vector<Real> aa(aj), bb(bj);
 124          do
 125          {
 126             scoped_precision<Real> p(current_precision);
 127             for (auto ai = aa.begin(); ai != aa.end(); ++ai)
 128                ai->precision(current_precision);
 129             for (auto bi = bb.begin(); bi != bb.end(); ++bi)
 130                bi->precision(current_precision);
 131             z.precision(current_precision);
 132             try
 133             {
 134                int scale = 0;
 135                std::pair<Real, Real> rp = boost::math::detail::hypergeometric_pFq_checked_series_impl(aa, bb, z, pol, boost::math::detail::timed_iteration_terminator(boost::math::policies::get_max_series_iterations<Policy>(), timeout), scale);
 136                rp.first *= exp(Real(scale));
 137                rp.second *= exp(Real(scale));
 138
 139                r = rp.first;
 140                norm = rp.second;
 141
 142                unsigned cancellation;
 143                try {
 144                   cancellation = itrunc(log10(abs(norm / r)));
 145                }
 146                catch (const boost::math::rounding_error&)
 147                {
 148                   // Happens when r is near enough zero:
 149                   cancellation = UINT_MAX;
 150                }
 151                if (cancellation >= current_precision - 1)
 152                {
 153                   current_precision *= 2;
 154                   continue;
 155                }
 156                unsigned precision_obtained = current_precision - 1 - cancellation;
 157                if (precision_obtained < digits10)
 158                {
 159                   current_precision += digits10 - precision_obtained + 5;
 160                }
 161                else
 162                   break;
 163             }
 164             catch (const boost::math::evaluation_error&)
 165             {
 166                current_precision *= 2;
 167             }
 168             catch (const detail::pFq_termination_exception& e)
 169             {
 170                //
 171                // Either we have exhausted the number of series iterations, or the timeout.
 172                // Either way we quit now.
 173                throw boost::math::evaluation_error(e.what());
 174             }
 175          } while (true);
 176
 177          return r;
 178       }
 179       template <class Seq, class Real>
 180       Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, const Real& z, unsigned digits10, double timeout = 0.5)
 181       {
 182          return hypergeometric_pFq_precision(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());
 183       }
 184
 185       template <class Real, class Policy>
 186       Real hypergeometric_pFq_precision(const std::initializer_list<Real>& aj, const std::initializer_list<Real>& bj, const Real& z, unsigned digits10, double timeout, const Policy& pol)
 187       {
 188          return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, pol);
 189       }
 190       template <class Real>
 191       Real hypergeometric_pFq_precision(const std::initializer_list<Real>& aj, const std::initializer_list<Real>& bj, const Real& z, unsigned digits10, double timeout = 0.5)
 192       {
 193          return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());
 194       }
 195
 196    }
 197 } // namespaces
 198
 199 #endif // BOOST_MATH_BESSEL_ITERATORS_HPP