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