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

//  (C) Copyright Nick Thompson 2017.
//  Use, modification and distribution are subject to 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_SPECIAL_CHEBYSHEV_HPP
#define BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
#include <cmath>
#include <boost/math/special_functions/math_fwd.hpp>
#include <boost/math/policies/error_handling.hpp>
#include <boost/math/constants/constants.hpp>
#include <boost/math/tools/promotion.hpp>

#if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610))
#  define BOOST_MATH_CHEB_USE_STD_ACOSH
#endif

#ifndef BOOST_MATH_CHEB_USE_STD_ACOSH
#  include <boost/math/special_functions/acosh.hpp>
#endif

namespace boost { namespace math {

template<class T1, class T2, class T3>
inline typename tools::promote_args<T1, T2, T3>::type chebyshev_next(T1 const & x, T2 const & Tn, T3 const & Tn_1)
{
    return 2*x*Tn - Tn_1;
}

namespace detail {

template<class Real, bool second, class Policy>
inline Real chebyshev_imp(unsigned n, Real const & x, const Policy&)
{
#ifdef BOOST_MATH_CHEB_USE_STD_ACOSH
    using std::acosh;
#define BOOST_MATH_ACOSH_POLICY
#else
   using boost::math::acosh;
#define BOOST_MATH_ACOSH_POLICY , Policy()
#endif
    using std::cosh;
    using std::pow;
    using std::sqrt;
    Real T0 = 1;
    Real T1;
    if (second)
    {
        if (x > 1 || x < -1)
        {
            Real t = sqrt(x*x -1);
            return static_cast<Real>((pow(x+t, (int)(n+1)) - pow(x-t, (int)(n+1)))/(2*t));
        }
        T1 = 2*x;
    }
    else
    {
        if (x > 1)
        {
            return cosh(n*acosh(x BOOST_MATH_ACOSH_POLICY));
        }
        if (x < -1)
        {
            if (n & 1)
            {
                return -cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
            }
            else
            {
                return cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
            }
        }
        T1 = x;
    }

    if (n == 0)
    {
        return T0;
    }

    unsigned l = 1;
    while(l < n)
    {
       std::swap(T0, T1);
       T1 = boost::math::chebyshev_next(x, T0, T1);
       ++l;
    }
    return T1;
}
} // namespace detail

template <class Real, class Policy>
inline typename tools::promote_args<Real>::type
chebyshev_t(unsigned n, Real const & x, const Policy&)
{
   typedef typename tools::promote_args<Real>::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;

   return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, false>(n, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_t<%1%>(unsigned, %1%)");
}

template<class Real>
inline typename tools::promote_args<Real>::type chebyshev_t(unsigned n, Real const & x)
{
    return chebyshev_t(n, x, policies::policy<>());
}

template <class Real, class Policy>
inline typename tools::promote_args<Real>::type
chebyshev_u(unsigned n, Real const & x, const Policy&)
{
   typedef typename tools::promote_args<Real>::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;

   return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, true>(n, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_u<%1%>(unsigned, %1%)");
}

template<class Real>
inline typename tools::promote_args<Real>::type chebyshev_u(unsigned n, Real const & x)
{
    return chebyshev_u(n, x, policies::policy<>());
}

template <class Real, class Policy>
inline typename tools::promote_args<Real>::type
chebyshev_t_prime(unsigned n, Real const & x, const Policy&)
{
   typedef typename tools::promote_args<Real>::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;
   if (n == 0)
   {
      return result_type(0);
   }
   return policies::checked_narrowing_cast<result_type, Policy>(n * detail::chebyshev_imp<value_type, true>(n - 1, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_t_prime<%1%>(unsigned, %1%)");
}

template<class Real>
inline typename tools::promote_args<Real>::type chebyshev_t_prime(unsigned n, Real const & x)
{
   return chebyshev_t_prime(n, x, policies::policy<>());
}

/*
 * This is Algorithm 3.1 of
 * Gil, Amparo, Javier Segura, and Nico M. Temme.
 * Numerical methods for special functions.
 * Society for Industrial and Applied Mathematics, 2007.
 * https://www.siam.org/books/ot99/OT99SampleChapter.pdf
 * However, our definition of c0 differs by a factor of 1/2, as stated in the docs. . .
 */
template<class Real, class T2>
inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x)
{
    using boost::math::constants::half;
    if (length < 2)
    {
        if (length == 0)
        {
            return 0;
        }
        return c[0]/2;
    }
    Real b2 = 0;
    Real b1 = c[length -1];
    for(size_t j = length - 2; j >= 1; --j)
    {
        Real tmp = 2*x*b1 - b2 + c[j];
        b2 = b1;
        b1 = tmp;
    }
    return x*b1 - b2 + half<Real>()*c[0];
}


namespace detail {
template<class Real>
inline Real unchecked_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
{
    Real t;
    Real u;
    // This cutoff is not super well defined, but it's a good estimate.
    // See "An Error Analysis of the Modified Clenshaw Method for Evaluating Chebyshev and Fourier Series"
    // J. OLIVER, IMA Journal of Applied Mathematics, Volume 20, Issue 3, November 1977, Pages 379-391
    // https://doi.org/10.1093/imamat/20.3.379
    const Real cutoff = 0.6;
    if (x - a < b - x)
    {
        u = 2*(x-a)/(b-a);
        t = u - 1;
        if (t > -cutoff)
        {
            Real b2 = 0;
            Real b1 = c[length -1];
            for(size_t j = length - 2; j >= 1; --j)
            {
                Real tmp = 2*t*b1 - b2 + c[j];
                b2 = b1;
                b1 = tmp;
            }
            return t*b1 - b2 + c[0]/2;
        }
        else
        {
            Real b = c[length -1];
            Real d = b;
            Real b2 = 0;
            for (size_t r = length - 2; r >= 1; --r)
            {
                d = 2*u*b - d + c[r];
                b2 = b;
                b = d - b;
            }
            return t*b - b2 + c[0]/2;
        }
    }
    else
    {
        u = -2*(b-x)/(b-a);
        t = u + 1;
        if (t < cutoff)
        {
            Real b2 = 0;
            Real b1 = c[length -1];
            for(size_t j = length - 2; j >= 1; --j)
            {
                Real tmp = 2*t*b1 - b2 + c[j];
                b2 = b1;
                b1 = tmp;
            }
            return t*b1 - b2 + c[0]/2;
        }
        else
        {
            Real b = c[length -1];
            Real d = b;
            Real b2 = 0;
            for (size_t r = length - 2; r >= 1; --r)
            {
                d = 2*u*b + d + c[r];
                b2 = b;
                b = d + b;
            }
            return t*b - b2 + c[0]/2;
        }
    }
}

} // namespace detail

template<class Real>
inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
{
    if (x < a || x > b)
    {
        throw std::domain_error("x in [a, b] is required.");
    }
    if (length < 2)
    {
        if (length == 0)
        {
            return 0;
        }
        return c[0]/2;
    }
    return detail::unchecked_chebyshev_clenshaw_recurrence(c, length, a, b, x);
}


}}
#endif
Commit	Line	Data
b32b8144 FG	1	// (C) Copyright Nick Thompson 2017.
	2	// Use, modification and distribution are subject to the
	3	// Boost Software License, Version 1.0. (See accompanying file
	4	// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
	5
	6	#ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
	7	#define BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
	8	#include <cmath>
1e59de90	9	#include <boost/math/special_functions/math_fwd.hpp>
20effc67	10	#include <boost/math/policies/error_handling.hpp>
b32b8144	11	#include <boost/math/constants/constants.hpp>
1e59de90	12	#include <boost/math/tools/promotion.hpp>
b32b8144 FG	13
	14	#if (__cplusplus > 201103) \|\| (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610))
	15	# define BOOST_MATH_CHEB_USE_STD_ACOSH
	16	#endif
	17
	18	#ifndef BOOST_MATH_CHEB_USE_STD_ACOSH
	19	# include <boost/math/special_functions/acosh.hpp>
	20	#endif
	21
	22	namespace boost { namespace math {
	23
	24	template<class T1, class T2, class T3>
	25	inline typename tools::promote_args<T1, T2, T3>::type chebyshev_next(T1 const & x, T2 const & Tn, T3 const & Tn_1)
	26	{
	27	return 2xTn - Tn_1;
	28	}
	29
	30	namespace detail {
	31
20effc67 TL	32	template<class Real, bool second, class Policy>
20effc67 TL	33	inline Real chebyshev_imp(unsigned n, Real const & x, const Policy&)
b32b8144 FG	34	{
	35	#ifdef BOOST_MATH_CHEB_USE_STD_ACOSH
	36	using std::acosh;
20effc67	37	#define BOOST_MATH_ACOSH_POLICY
b32b8144 FG	38	#else
b32b8144 FG	39	using boost::math::acosh;
20effc67	40	#define BOOST_MATH_ACOSH_POLICY , Policy()
b32b8144 FG	41	#endif
	42	using std::cosh;
	43	using std::pow;
	44	using std::sqrt;
	45	Real T0 = 1;
	46	Real T1;
	47	if (second)
	48	{
	49	if (x > 1 \|\| x < -1)
	50	{
	51	Real t = sqrt(x*x -1);
	52	return static_cast<Real>((pow(x+t, (int)(n+1)) - pow(x-t, (int)(n+1)))/(2*t));
	53	}
	54	T1 = 2*x;
	55	}
	56	else
	57	{
	58	if (x > 1)
	59	{
20effc67	60	return cosh(n*acosh(x BOOST_MATH_ACOSH_POLICY));
b32b8144 FG	61	}
	62	if (x < -1)
	63	{
	64	if (n & 1)
	65	{
20effc67	66	return -cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
b32b8144 FG	67	}
	68	else
	69	{
20effc67	70	return cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
b32b8144 FG	71	}
	72	}
	73	T1 = x;
	74	}
	75
	76	if (n == 0)
	77	{
	78	return T0;
	79	}
	80
	81	unsigned l = 1;
	82	while(l < n)
	83	{
	84	std::swap(T0, T1);
	85	T1 = boost::math::chebyshev_next(x, T0, T1);
	86	++l;
	87	}
	88	return T1;
	89	}
	90	} // namespace detail
	91
	92	template <class Real, class Policy>
	93	inline typename tools::promote_args<Real>::type
	94	chebyshev_t(unsigned n, Real const & x, const Policy&)
	95	{
	96	typedef typename tools::promote_args<Real>::type result_type;
	97	typedef typename policies::evaluation<result_type, Policy>::type value_type;
20effc67 TL	98	typedef typename policies::normalise<
	99	Policy,
	100	policies::promote_float<false>,
	101	policies::promote_double<false>,
	102	policies::discrete_quantile<>,
	103	policies::assert_undefined<> >::type forwarding_policy;
	104
	105	return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, false>(n, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_t<%1%>(unsigned, %1%)");
b32b8144 FG	106	}
	107
	108	template<class Real>
	109	inline typename tools::promote_args<Real>::type chebyshev_t(unsigned n, Real const & x)
	110	{
	111	return chebyshev_t(n, x, policies::policy<>());
	112	}
	113
	114	template <class Real, class Policy>
	115	inline typename tools::promote_args<Real>::type
	116	chebyshev_u(unsigned n, Real const & x, const Policy&)
	117	{
	118	typedef typename tools::promote_args<Real>::type result_type;
	119	typedef typename policies::evaluation<result_type, Policy>::type value_type;
20effc67 TL	120	typedef typename policies::normalise<
	121	Policy,
	122	policies::promote_float<false>,
	123	policies::promote_double<false>,
	124	policies::discrete_quantile<>,
	125	policies::assert_undefined<> >::type forwarding_policy;
	126
	127	return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, true>(n, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_u<%1%>(unsigned, %1%)");
b32b8144 FG	128	}
	129
	130	template<class Real>
	131	inline typename tools::promote_args<Real>::type chebyshev_u(unsigned n, Real const & x)
	132	{
	133	return chebyshev_u(n, x, policies::policy<>());
	134	}
	135
	136	template <class Real, class Policy>
	137	inline typename tools::promote_args<Real>::type
	138	chebyshev_t_prime(unsigned n, Real const & x, const Policy&)
	139	{
	140	typedef typename tools::promote_args<Real>::type result_type;
	141	typedef typename policies::evaluation<result_type, Policy>::type value_type;
20effc67 TL	142	typedef typename policies::normalise<
	143	Policy,
	144	policies::promote_float<false>,
	145	policies::promote_double<false>,
	146	policies::discrete_quantile<>,
	147	policies::assert_undefined<> >::type forwarding_policy;
b32b8144 FG	148	if (n == 0)
	149	{
	150	return result_type(0);
	151	}
20effc67	152	return policies::checked_narrowing_cast<result_type, Policy>(n * detail::chebyshev_imp<value_type, true>(n - 1, static_cast<value_type>(x), forwarding_policy()), "boost::math::chebyshev_t_prime<%1%>(unsigned, %1%)");
b32b8144 FG	153	}
	154
	155	template<class Real>
	156	inline typename tools::promote_args<Real>::type chebyshev_t_prime(unsigned n, Real const & x)
	157	{
	158	return chebyshev_t_prime(n, x, policies::policy<>());
	159	}
	160
	161	/*
	162	* This is Algorithm 3.1 of
	163	* Gil, Amparo, Javier Segura, and Nico M. Temme.
	164	* Numerical methods for special functions.
	165	* Society for Industrial and Applied Mathematics, 2007.
	166	* https://www.siam.org/books/ot99/OT99SampleChapter.pdf
	167	* However, our definition of c0 differs by a factor of 1/2, as stated in the docs. . .
	168	*/
	169	template<class Real, class T2>
	170	inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x)
	171	{
	172	using boost::math::constants::half;
	173	if (length < 2)
	174	{
	175	if (length == 0)
	176	{
	177	return 0;
	178	}
	179	return c[0]/2;
	180	}
	181	Real b2 = 0;
	182	Real b1 = c[length -1];
	183	for(size_t j = length - 2; j >= 1; --j)
	184	{
	185	Real tmp = 2xb1 - b2 + c[j];
	186	b2 = b1;
	187	b1 = tmp;
	188	}
	189	return xb1 - b2 + half<Real>()c[0];
	190	}
	191
	192
20effc67 TL	193
	194	namespace detail {
	195	template<class Real>
	196	inline Real unchecked_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
	197	{
	198	Real t;
	199	Real u;
	200	// This cutoff is not super well defined, but it's a good estimate.
	201	// See "An Error Analysis of the Modified Clenshaw Method for Evaluating Chebyshev and Fourier Series"
1e59de90	202	// J. OLIVER, IMA Journal of Applied Mathematics, Volume 20, Issue 3, November 1977, Pages 379-391
20effc67 TL	203	// https://doi.org/10.1093/imamat/20.3.379
	204	const Real cutoff = 0.6;
	205	if (x - a < b - x)
	206	{
	207	u = 2*(x-a)/(b-a);
	208	t = u - 1;
	209	if (t > -cutoff)
	210	{
	211	Real b2 = 0;
	212	Real b1 = c[length -1];
	213	for(size_t j = length - 2; j >= 1; --j)
	214	{
	215	Real tmp = 2tb1 - b2 + c[j];
	216	b2 = b1;
	217	b1 = tmp;
	218	}
	219	return t*b1 - b2 + c[0]/2;
	220	}
	221	else
	222	{
	223	Real b = c[length -1];
	224	Real d = b;
	225	Real b2 = 0;
	226	for (size_t r = length - 2; r >= 1; --r)
	227	{
	228	d = 2ub - d + c[r];
	229	b2 = b;
	230	b = d - b;
	231	}
	232	return t*b - b2 + c[0]/2;
	233	}
	234	}
	235	else
	236	{
	237	u = -2*(b-x)/(b-a);
	238	t = u + 1;
	239	if (t < cutoff)
	240	{
	241	Real b2 = 0;
	242	Real b1 = c[length -1];
	243	for(size_t j = length - 2; j >= 1; --j)
	244	{
	245	Real tmp = 2tb1 - b2 + c[j];
	246	b2 = b1;
	247	b1 = tmp;
	248	}
	249	return t*b1 - b2 + c[0]/2;
	250	}
	251	else
	252	{
	253	Real b = c[length -1];
	254	Real d = b;
	255	Real b2 = 0;
	256	for (size_t r = length - 2; r >= 1; --r)
	257	{
	258	d = 2ub + d + c[r];
	259	b2 = b;
	260	b = d + b;
	261	}
	262	return t*b - b2 + c[0]/2;
	263	}
	264	}
	265	}
	266
267	} // namespace detail
268
269	template<class Real>
270	inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
271	{
272	if (x < a \|\| x > b)
273	{
274	throw std::domain_error("x in [a, b] is required.");
275	}
276	if (length < 2)
277	{
278	if (length == 0)
279	{
280	return 0;
281	}
282	return c[0]/2;
283	}
284	return detail::unchecked_chebyshev_clenshaw_recurrence(c, length, a, b, x);
285	}
286
287
b32b8144 FG	288	}}
b32b8144 FG	289	#endif