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>
chebyshev_next(T1 const & x,T2 const & Tn,T3 const & Tn_1)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>
chebyshev_imp(unsigned n,Real const & x)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
chebyshev_t(unsigned n,Real const & x,const Policy &)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>
chebyshev_t(unsigned n,Real const & x)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
chebyshev_u(unsigned n,Real const & x,const Policy &)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>
chebyshev_u(unsigned n,Real const & x)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
chebyshev_t_prime(unsigned n,Real const & x,const Policy &)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>
chebyshev_t_prime(unsigned n,Real const & x)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>
chebyshev_clenshaw_recurrence(const Real * const c,size_t length,const T2 & x)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