[ceph.git] / ceph / src / boost / boost / math / special_functions / hypergeometric_pFq.hpp


///////////////////////////////////////////////////////////////////////////////
//  Copyright 2018 John Maddock
//  Distributed under the Boost
//  Software License, Version 1.0. (See accompanying file
//  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)

#ifndef BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP
#define BOOST_MATH_HYPERGEOMETRIC_PFQ_HPP

#include <boost/config.hpp>

#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) || defined(BOOST_NO_CXX11_HDR_CHRONO)
# error "hypergeometric_pFq requires a C++11 compiler"
#endif

#include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
#include <chrono>
#include <initializer_list>

namespace boost {
   namespace math {

      namespace detail {

         struct pFq_termination_exception : public std::runtime_error
         {
            pFq_termination_exception(const char* p) : std::runtime_error(p) {}
         };

         struct timed_iteration_terminator
         {
            timed_iteration_terminator(boost::uintmax_t i, double t) : max_iter(i), max_time(t), start_time(std::chrono::system_clock::now()) {}

            bool operator()(boost::uintmax_t iter)const
            {
               if (iter > max_iter)
                  boost::throw_exception(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted iterations."));
               if (std::chrono::duration<double>(std::chrono::system_clock::now() - start_time).count() > max_time)
                  boost::throw_exception(boost::math::detail::pFq_termination_exception("pFq exceeded maximum permitted evaluation time."));
               return false;
            }

            boost::uintmax_t max_iter;
            double max_time;
            std::chrono::system_clock::time_point start_time;
         };

      }

      template <class Seq, class Real, class Policy>
      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)
      {
         typedef typename tools::promote_args<Real, typename Seq::value_type>::type result_type;
         typedef typename policies::evaluation<result_type, Policy>::type value_type;
         typedef typename policies::normalise<
            Policy,
            policies::promote_float<false>,
            policies::promote_double<false>,
            policies::discrete_quantile<>,
            policies::assert_undefined<> >::type forwarding_policy;

         BOOST_MATH_STD_USING

         int scale = 0;
         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);
         r.first *= exp(Real(scale));
         r.second *= exp(Real(scale));
         if (p_abs_error)
            *p_abs_error = static_cast<Real>(r.second) * boost::math::tools::epsilon<Real>();
         return policies::checked_narrowing_cast<result_type, Policy>(r.first, "boost::math::hypergeometric_pFq<%1%>(%1%,%1%,%1%)");
      }

      template <class Seq, class Real>
      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)
      {
         return hypergeometric_pFq(aj, bj, z, p_abs_error, boost::math::policies::policy<>());
      }

      template <class R, class Real, class Policy>
      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)
      {
         return hypergeometric_pFq<std::initializer_list<R>, Real, Policy>(aj, bj, z, p_abs_error, pol);
      }
      
      template <class R, class Real>
      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)
      {
         return hypergeometric_pFq<std::initializer_list<R>, Real>(aj, bj, z, p_abs_error);
      }

      template <class T>
      struct scoped_precision
      {
         scoped_precision(unsigned p)
         {
            old_p = T::default_precision();
            T::default_precision(p);
         }
         ~scoped_precision()
         {
            T::default_precision(old_p);
         }
         unsigned old_p;
      };

      template <class Seq, class Real, class Policy>
      Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, Real z, unsigned digits10, double timeout, const Policy& pol)
      {
         unsigned current_precision = digits10 + 5;

         for (auto ai = aj.begin(); ai != aj.end(); ++ai)
         {
            current_precision = (std::max)(current_precision, ai->precision());
         }
         for (auto bi = bj.begin(); bi != bj.end(); ++bi)
         {
            current_precision = (std::max)(current_precision, bi->precision());
         }
         current_precision = (std::max)(current_precision, z.precision());

         Real r, norm;
         std::vector<Real> aa(aj), bb(bj);
         do
         {
            scoped_precision<Real> p(current_precision);
            for (auto ai = aa.begin(); ai != aa.end(); ++ai)
               ai->precision(current_precision);
            for (auto bi = bb.begin(); bi != bb.end(); ++bi)
               bi->precision(current_precision);
            z.precision(current_precision);
            try
            {
               int scale = 0;
               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);
               rp.first *= exp(Real(scale));
               rp.second *= exp(Real(scale));

               r = rp.first;
               norm = rp.second;

               unsigned cancellation;
               try {
                  cancellation = itrunc(log10(abs(norm / r)));
               }
               catch (const boost::math::rounding_error&)
               {
                  // Happens when r is near enough zero:
                  cancellation = UINT_MAX;
               }
               if (cancellation >= current_precision - 1)
               {
                  current_precision *= 2;
                  continue;
               }
               unsigned precision_obtained = current_precision - 1 - cancellation;
               if (precision_obtained < digits10)
               {
                  current_precision += digits10 - precision_obtained + 5;
               }
               else
                  break;
            }
            catch (const boost::math::evaluation_error&)
            {
               current_precision *= 2;
            }
            catch (const detail::pFq_termination_exception& e)
            {
               //
               // Either we have exhausted the number of series iterations, or the timeout.
               // Either way we quit now.
               throw boost::math::evaluation_error(e.what());
            }
         } while (true);

         return r;
      }
      template <class Seq, class Real>
      Real hypergeometric_pFq_precision(const Seq& aj, const Seq& bj, const Real& z, unsigned digits10, double timeout = 0.5)
      {
         return hypergeometric_pFq_precision(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());
      }

      template <class Real, class Policy>
      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)
      {
         return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, pol);
      }
      template <class Real>
      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)
      {
         return hypergeometric_pFq_precision< std::initializer_list<Real>, Real>(aj, bj, z, digits10, timeout, boost::math::policies::policy<>());
      }

   }
} // namespaces

#endif // BOOST_MATH_BESSEL_ITERATORS_HPP
Commit	Line	Data
92f5a8d4 TL	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
f67539c2	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) \|\| defined(BOOST_NO_CXX11_HDR_CHRONO)
92f5a8d4 TL	14	# error "hypergeometric_pFq requires a C++11 compiler"
	15	#endif
	16
	17	#include <boost/math/special_functions/detail/hypergeometric_pFq_checked_series.hpp>
f67539c2	18	#include <chrono>
92f5a8d4 TL	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	{
f67539c2	33	timed_iteration_terminator(boost::uintmax_t i, double t) : max_iter(i), max_time(t), start_time(std::chrono::system_clock::now()) {}
92f5a8d4 TL	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."));
f67539c2	39	if (std::chrono::duration<double>(std::chrono::system_clock::now() - start_time).count() > max_time)
92f5a8d4 TL	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;
f67539c2	46	std::chrono::system_clock::time_point start_time;
92f5a8d4 TL	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