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/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>
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>
chebyshev_next(T1 const & x,T2 const & Tn,T3 const & Tn_1)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 2*x*Tn - Tn_1;
28 }
29 
30 namespace detail {
31 
32 template<class Real, bool second, class Policy>
chebyshev_imp(unsigned n,Real const & x,const Policy &)33 inline Real chebyshev_imp(unsigned n, Real const & x, const Policy&)
34 {
35 #ifdef BOOST_MATH_CHEB_USE_STD_ACOSH
36     using std::acosh;
37 #define BOOST_MATH_ACOSH_POLICY
38 #else
39    using boost::math::acosh;
40 #define BOOST_MATH_ACOSH_POLICY , Policy()
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         {
60             return cosh(n*acosh(x BOOST_MATH_ACOSH_POLICY));
61         }
62         if (x < -1)
63         {
64             if (n & 1)
65             {
66                 return -cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
67             }
68             else
69             {
70                 return cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
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
chebyshev_t(unsigned n,Real const & x,const Policy &)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;
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%)");
106 }
107 
108 template<class Real>
chebyshev_t(unsigned n,Real const & x)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
chebyshev_u(unsigned n,Real const & x,const Policy &)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;
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%)");
128 }
129 
130 template<class Real>
chebyshev_u(unsigned n,Real const & x)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
chebyshev_t_prime(unsigned n,Real const & x,const Policy &)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;
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;
148    if (n == 0)
149    {
150       return result_type(0);
151    }
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%)");
153 }
154 
155 template<class Real>
chebyshev_t_prime(unsigned n,Real const & x)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>
chebyshev_clenshaw_recurrence(const Real * const c,size_t length,const T2 & x)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 = 2*x*b1 - b2 + c[j];
186         b2 = b1;
187         b1 = tmp;
188     }
189     return x*b1 - b2 + half<Real>()*c[0];
190 }
191 
192 
193 
194 namespace detail {
195 template<class Real>
unchecked_chebyshev_clenshaw_recurrence(const Real * const c,size_t length,const Real & a,const Real & b,const Real & x)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"
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;
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 = 2*t*b1 - 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 = 2*u*b - 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 = 2*t*b1 - 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 = 2*u*b + 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>
chebyshev_clenshaw_recurrence(const Real * const c,size_t length,const Real & a,const Real & b,const Real & x)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 
288 }}
289 #endif
290