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/constants/constants.hpp>
10 
11 #if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610))
12 #  define BOOST_MATH_CHEB_USE_STD_ACOSH
13 #endif
14 
15 #ifndef BOOST_MATH_CHEB_USE_STD_ACOSH
16 #  include <boost/math/special_functions/acosh.hpp>
17 #endif
18 
19 namespace boost { namespace math {
20 
21 template<class T1, class T2, class T3>
22 inline typename tools::promote_args<T1, T2, T3>::type chebyshev_next(T1 const & x, T2 const & Tn, T3 const & Tn_1)
23 {
24     return 2*x*Tn - Tn_1;
25 }
26 
27 namespace detail {
28 
29 template<class Real, bool second>
30 inline Real chebyshev_imp(unsigned n, Real const & x)
31 {
32 #ifdef BOOST_MATH_CHEB_USE_STD_ACOSH
33     using std::acosh;
34 #else
35    using boost::math::acosh;
36 #endif
37     using std::cosh;
38     using std::pow;
39     using std::sqrt;
40     Real T0 = 1;
41     Real T1;
42     if (second)
43     {
44         if (x > 1 || x < -1)
45         {
46             Real t = sqrt(x*x -1);
47             return static_cast<Real>((pow(x+t, (int)(n+1)) - pow(x-t, (int)(n+1)))/(2*t));
48         }
49         T1 = 2*x;
50     }
51     else
52     {
53         if (x > 1)
54         {
55             return cosh(n*acosh(x));
56         }
57         if (x < -1)
58         {
59             if (n & 1)
60             {
61                 return -cosh(n*acosh(-x));
62             }
63             else
64             {
65                 return cosh(n*acosh(-x));
66             }
67         }
68         T1 = x;
69     }
70 
71     if (n == 0)
72     {
73         return T0;
74     }
75 
76     unsigned l = 1;
77     while(l < n)
78     {
79        std::swap(T0, T1);
80        T1 = boost::math::chebyshev_next(x, T0, T1);
81        ++l;
82     }
83     return T1;
84 }
85 } // namespace detail
86 
87 template <class Real, class Policy>
88 inline typename tools::promote_args<Real>::type
89 chebyshev_t(unsigned n, Real const & x, const Policy&)
90 {
91    typedef typename tools::promote_args<Real>::type result_type;
92    typedef typename policies::evaluation<result_type, Policy>::type value_type;
93    return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, false>(n, static_cast<value_type>(x)), "boost::math::chebyshev_t<%1%>(unsigned, %1%)");
94 }
95 
96 template<class Real>
97 inline typename tools::promote_args<Real>::type chebyshev_t(unsigned n, Real const & x)
98 {
99     return chebyshev_t(n, x, policies::policy<>());
100 }
101 
102 template <class Real, class Policy>
103 inline typename tools::promote_args<Real>::type
104 chebyshev_u(unsigned n, Real const & x, const Policy&)
105 {
106    typedef typename tools::promote_args<Real>::type result_type;
107    typedef typename policies::evaluation<result_type, Policy>::type value_type;
108    return policies::checked_narrowing_cast<result_type, Policy>(detail::chebyshev_imp<value_type, true>(n, static_cast<value_type>(x)), "boost::math::chebyshev_u<%1%>(unsigned, %1%)");
109 }
110 
111 template<class Real>
112 inline typename tools::promote_args<Real>::type chebyshev_u(unsigned n, Real const & x)
113 {
114     return chebyshev_u(n, x, policies::policy<>());
115 }
116 
117 template <class Real, class Policy>
118 inline typename tools::promote_args<Real>::type
119 chebyshev_t_prime(unsigned n, Real const & x, const Policy&)
120 {
121    typedef typename tools::promote_args<Real>::type result_type;
122    typedef typename policies::evaluation<result_type, Policy>::type value_type;
123    if (n == 0)
124    {
125       return result_type(0);
126    }
127    return policies::checked_narrowing_cast<result_type, Policy>(n * detail::chebyshev_imp<value_type, true>(n - 1, static_cast<value_type>(x)), "boost::math::chebyshev_t_prime<%1%>(unsigned, %1%)");
128 }
129 
130 template<class Real>
131 inline typename tools::promote_args<Real>::type chebyshev_t_prime(unsigned n, Real const & x)
132 {
133    return chebyshev_t_prime(n, x, policies::policy<>());
134 }
135 
136 /*
137  * This is Algorithm 3.1 of
138  * Gil, Amparo, Javier Segura, and Nico M. Temme.
139  * Numerical methods for special functions.
140  * Society for Industrial and Applied Mathematics, 2007.
141  * https://www.siam.org/books/ot99/OT99SampleChapter.pdf
142  * However, our definition of c0 differs by a factor of 1/2, as stated in the docs. . .
143  */
144 template<class Real, class T2>
145 inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x)
146 {
147     using boost::math::constants::half;
148     if (length < 2)
149     {
150         if (length == 0)
151         {
152             return 0;
153         }
154         return c[0]/2;
155     }
156     Real b2 = 0;
157     Real b1 = c[length -1];
158     for(size_t j = length - 2; j >= 1; --j)
159     {
160         Real tmp = 2*x*b1 - b2 + c[j];
161         b2 = b1;
162         b1 = tmp;
163     }
164     return x*b1 - b2 + half<Real>()*c[0];
165 }
166 
167 
168 }}
169 #endif
170