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)
6 #ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
7 #define BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
9 #include <boost/math/special_functions/math_fwd.hpp>
10 #include <boost/math/policies/error_handling.hpp>
11 #include <boost/math/constants/constants.hpp>
12 #include <boost/math/tools/promotion.hpp>
14 #if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610))
15 # define BOOST_MATH_CHEB_USE_STD_ACOSH
18 #ifndef BOOST_MATH_CHEB_USE_STD_ACOSH
19 # include <boost/math/special_functions/acosh.hpp>
22 namespace boost { namespace math {
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)
32 template<class Real, bool second, class Policy>
33 inline Real chebyshev_imp(unsigned n, Real const & x, const Policy&)
35 #ifdef BOOST_MATH_CHEB_USE_STD_ACOSH
37 #define BOOST_MATH_ACOSH_POLICY
39 using boost::math::acosh;
40 #define BOOST_MATH_ACOSH_POLICY , Policy()
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));
60 return cosh(n*acosh(x BOOST_MATH_ACOSH_POLICY));
66 return -cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
70 return cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
85 T1 = boost::math::chebyshev_next(x, T0, T1);
92 template <class Real, class Policy>
93 inline typename tools::promote_args<Real>::type
94 chebyshev_t(unsigned n, Real const & x, const Policy&)
96 typedef typename tools::promote_args<Real>::type result_type;
97 typedef typename policies::evaluation<result_type, Policy>::type value_type;
98 typedef typename policies::normalise<
100 policies::promote_float<false>,
101 policies::promote_double<false>,
102 policies::discrete_quantile<>,
103 policies::assert_undefined<> >::type forwarding_policy;
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%)");
109 inline typename tools::promote_args<Real>::type chebyshev_t(unsigned n, Real const & x)
111 return chebyshev_t(n, x, policies::policy<>());
114 template <class Real, class Policy>
115 inline typename tools::promote_args<Real>::type
116 chebyshev_u(unsigned n, Real const & x, const Policy&)
118 typedef typename tools::promote_args<Real>::type result_type;
119 typedef typename policies::evaluation<result_type, Policy>::type value_type;
120 typedef typename policies::normalise<
122 policies::promote_float<false>,
123 policies::promote_double<false>,
124 policies::discrete_quantile<>,
125 policies::assert_undefined<> >::type forwarding_policy;
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%)");
131 inline typename tools::promote_args<Real>::type chebyshev_u(unsigned n, Real const & x)
133 return chebyshev_u(n, x, policies::policy<>());
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&)
140 typedef typename tools::promote_args<Real>::type result_type;
141 typedef typename policies::evaluation<result_type, Policy>::type value_type;
142 typedef typename policies::normalise<
144 policies::promote_float<false>,
145 policies::promote_double<false>,
146 policies::discrete_quantile<>,
147 policies::assert_undefined<> >::type forwarding_policy;
150 return result_type(0);
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%)");
156 inline typename tools::promote_args<Real>::type chebyshev_t_prime(unsigned n, Real const & x)
158 return chebyshev_t_prime(n, x, policies::policy<>());
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. . .
169 template<class Real, class T2>
170 inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x)
172 using boost::math::constants::half;
182 Real b1 = c[length -1];
183 for(size_t j = length - 2; j >= 1; --j)
185 Real tmp = 2*x*b1 - b2 + c[j];
189 return x*b1 - b2 + half<Real>()*c[0];
196 inline Real unchecked_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
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"
202 // J. OLIVER, IMA Journal of Applied Mathematics, Volume 20, Issue 3, November 1977, Pages 379-391
203 // https://doi.org/10.1093/imamat/20.3.379
204 const Real cutoff = 0.6;
212 Real b1 = c[length -1];
213 for(size_t j = length - 2; j >= 1; --j)
215 Real tmp = 2*t*b1 - b2 + c[j];
219 return t*b1 - b2 + c[0]/2;
223 Real b = c[length -1];
226 for (size_t r = length - 2; r >= 1; --r)
228 d = 2*u*b - d + c[r];
232 return t*b - b2 + c[0]/2;
242 Real b1 = c[length -1];
243 for(size_t j = length - 2; j >= 1; --j)
245 Real tmp = 2*t*b1 - b2 + c[j];
249 return t*b1 - b2 + c[0]/2;
253 Real b = c[length -1];
256 for (size_t r = length - 2; r >= 1; --r)
258 d = 2*u*b + d + c[r];
262 return t*b - b2 + c[0]/2;
267 } // namespace detail
270 inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
274 throw std::domain_error("x in [a, b] is required.");
284 return detail::unchecked_chebyshev_clenshaw_recurrence(c, length, a, b, x);