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