Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:40:15 

0001 
0002 ///////////////////////////////////////////////////////////////////////////////
0003 //  Copyright 2018 John Maddock
0004 //  Distributed under the Boost
0005 //  Software License, Version 1.0. (See accompanying file
0006 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0007 
0008 #ifndef BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP
0009 #define BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP
0010 
0011 #include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
0012 #include <boost/math/tools/throw_exception.hpp>
0013 #include <chrono>
0014 #include <initializer_list>
0015 
0016 namespace boost {
0017    namespace math {
0018 
0019       namespace detail {
0020 
0021          struct pFq_termination_exception : public std::runtime_error
0022          {
0023             pFq_termination_exception(const char* p) : std::runtime_error(p) {}
0024          };
0025 
0026          struct timed_iteration_terminator
0027          {
0028             timed_iteration_terminator(std::uintmax_t i, double t) : max_iter(i), max_time(t), start_time(std::chrono::system_clock::now()) {}
0029 
0030             bool operator()(std::uintmax_t iter)const
0031             {
0032                if (iter > max_iter)
0033                   BOOST_MATH_THROW_EXCEPTION(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted iterations."));
0034                if (std::chrono::duration<double>(std::chrono::system_clock::now() - start_time).count() > max_time)
0035                   BOOST_MATH_THROW_EXCEPTION(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted evaluation time."));
0036                return false;
0037             }
0038 
0039             std::uintmax_t max_iter;
0040             double max_time;
0041             std::chrono::system_clock::time_point start_time;
0042          };
0043 
0044       }
0045 
0046       template <class Seq, class Real, class Policy>
0047       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)
0048       {
0049          typedef typename tools::promote_args<Real, typename Seq::value_type>::type result_type;
0050          typedef typename policies::evaluation<result_type, Policy>::type value_type;
0051          typedef typename policies::normalise<
0052             Policy,
0053             policies::promote_float<false>,
0054             policies::promote_double<false>,
0055             policies::discrete_quantile<>,
0056             policies::assert_undefined<> >::type forwarding_policy;
0057 
0058          BOOST_MATH_STD_USING
0059 
0060          long long scale = 0;
0061          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);
0062          r.first *= exp(Real(scale));
0063          r.second *= exp(Real(scale));
0064          if (p_abs_error)
0065             *p_abs_error = static_cast<Real>(r.second) * boost::math::tools::epsilon<Real>();
0066          return policies::checked_narrowing_cast<result_type, Policy>(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)");
0067       }
0068 
0069       template <class Seq, class Real>
0070       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)
0071       {
0072          return hypergeometric_pFq(aj, bj, z, p_abs_error, boost::math::policies::policy<>());
0073       }
0074 
0075       template <class R, class Real, class Policy>
0076       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)
0077       {
0078          return hypergeometric_pFq<std::initializer_list<R>, Real, Policy>(aj, bj, z, p_abs_error, pol);
0079       }
0080 
0081       template <class R, class Real>
0082       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 = nullptr)
0083       {
0084          return hypergeometric_pFq<std::initializer_list<R>, Real>(aj, bj, z, p_abs_error);
0085       }
0086 
0087 #ifndef BOOST_NO_EXCEPTIONS
0088       template <class T>
0089       struct scoped_precision
0090       {
0091          scoped_precision(unsigned p)
0092          {
0093             old_p = T::default_precision();
0094             T::default_precision(p);
0095          }
0096          ~scoped_precision()
0097          {
0098             T::default_precision(old_p);
0099          }
0100          unsigned old_p;
0101       };
0102 
0103       template <class Seq, class Real, class Policy>
0104       Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, Real z, unsigned digits10, double timeout, const Policy& pol)
0105       {
0106          unsigned current_precision = digits10 + 5;
0107 
0108          for (auto ai = aj.begin(); ai != aj.end(); ++ai)
0109          {
0110             current_precision = (std::max)(current_precision, ai->precision());
0111          }
0112          for (auto bi = bj.begin(); bi != bj.end(); ++bi)
0113          {
0114             current_precision = (std::max)(current_precision, bi->precision());
0115          }
0116          current_precision = (std::max)(current_precision, z.precision());
0117 
0118          Real r, norm;
0119          std::vector<Real> aa(aj), bb(bj);
0120          do
0121          {
0122             scoped_precision<Real> p(current_precision);
0123             for (auto ai = aa.begin(); ai != aa.end(); ++ai)
0124                ai->precision(current_precision);
0125             for (auto bi = bb.begin(); bi != bb.end(); ++bi)
0126                bi->precision(current_precision);
0127             z.precision(current_precision);
0128             try
0129             {
0130                long long scale = 0;
0131                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);
0132                rp.first *= exp(Real(scale));
0133                rp.second *= exp(Real(scale));
0134 
0135                r = rp.first;
0136                norm = rp.second;
0137 
0138                unsigned cancellation;
0139                try {
0140                   cancellation = itrunc(log10(abs(norm / r)));
0141                }
0142                catch (const boost::math::rounding_error&)
0143                {
0144                   // Happens when r is near enough zero:
0145                   cancellation = UINT_MAX;
0146                }
0147                if (cancellation >= current_precision - 1)
0148                {
0149                   current_precision *= 2;
0150                   continue;
0151                }
0152                unsigned precision_obtained = current_precision - 1 - cancellation;
0153                if (precision_obtained < digits10)
0154                {
0155                   current_precision += digits10 - precision_obtained + 5;
0156                }
0157                else
0158                   break;
0159             }
0160             catch (const boost::math::evaluation_error&)
0161             {
0162                current_precision *= 2;
0163             }
0164             catch (const detail::pFq_termination_exception& e)
0165             {
0166                //
0167                // Either we have exhausted the number of series iterations, or the timeout.
0168                // Either way we quit now.
0169                throw boost::math::evaluation_error(e.what());
0170             }
0171          } while (true);
0172 
0173          return r;
0174       }
0175       template <class Seq, class Real>
0176       Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, const Real& z, unsigned digits10, double timeout = 0.5)
0177       {
0178          return hypergeometric_pFq_precision(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());
0179       }
0180 
0181       template <class Real, class Policy>
0182       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)
0183       {
0184          return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, pol);
0185       }
0186       template <class Real>
0187       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)
0188       {
0189          return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());
0190       }
0191 #endif
0192    }
0193 } // namespaces
0194 
0195 #endif // BOOST_MATH_BESSEL_ITERATORS_HPP
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!