315 lines · plain
1// (C) Copyright Nick Thompson 2017.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_HPP7#define BOOST_MATH_SPECIAL_CHEBYSHEV_HPP8#include <cmath>9#include <type_traits>10#include <boost/math/special_functions/math_fwd.hpp>11#include <boost/math/policies/error_handling.hpp>12#include <boost/math/constants/constants.hpp>13#include <boost/math/tools/promotion.hpp>14#include <boost/math/tools/throw_exception.hpp>15 16#if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610))17# define BOOST_MATH_CHEB_USE_STD_ACOSH18#endif19 20#ifndef BOOST_MATH_CHEB_USE_STD_ACOSH21# include <boost/math/special_functions/acosh.hpp>22#endif23 24namespace boost { namespace math {25 26template <class T1, class T2, class T3>27inline tools::promote_args_t<T1, T2, T3> chebyshev_next(T1 const & x, T2 const & Tn, T3 const & Tn_1)28{29 return 2*x*Tn - Tn_1;30}31 32namespace detail {33 34// https://stackoverflow.com/questions/5625431/efficient-way-to-compute-pq-exponentiation-where-q-is-an-integer35template <typename T, typename std::enable_if<std::is_arithmetic<T>::value, bool>::type = true>36T expt(T p, unsigned q)37{38 T r = 1;39 40 while (q != 0) {41 if (q % 2 == 1) { // q is odd42 r *= p;43 q--;44 }45 p *= p;46 q /= 2;47 }48 49 return r;50}51 52template <typename T, typename std::enable_if<!std::is_arithmetic<T>::value, bool>::type = true>53T expt(T p, unsigned q)54{55 using std::pow;56 return pow(p, static_cast<int>(q));57}58 59template<class Real, bool second, class Policy>60inline Real chebyshev_imp(unsigned n, Real const & x, const Policy&)61{62#ifdef BOOST_MATH_CHEB_USE_STD_ACOSH63 using std::acosh;64#define BOOST_MATH_ACOSH_POLICY65#else66 using boost::math::acosh;67#define BOOST_MATH_ACOSH_POLICY , Policy()68#endif69 using std::cosh;70 using std::pow;71 using std::sqrt;72 Real T0 = 1;73 Real T1;74 75 BOOST_MATH_IF_CONSTEXPR (second)76 {77 if (x > 1 || x < -1)78 {79 Real t = sqrt(x*x -1);80 return static_cast<Real>((expt(static_cast<Real>(x+t), n+1) - expt(static_cast<Real>(x-t), n+1))/(2*t));81 }82 T1 = 2*x;83 }84 else85 {86 if (x > 1)87 {88 return cosh(n*acosh(x BOOST_MATH_ACOSH_POLICY));89 }90 if (x < -1)91 {92 if (n & 1)93 {94 return -cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));95 }96 else97 {98 return cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));99 }100 }101 T1 = x;102 }103 104 if (n == 0)105 {106 return T0;107 }108 109 unsigned l = 1;110 while(l < n)111 {112 std::swap(T0, T1);113 T1 = static_cast<Real>(boost::math::chebyshev_next(x, T0, T1));114 ++l;115 }116 return T1;117}118} // namespace detail119 120template <class Real, class Policy>121inline tools::promote_args_t<Real> chebyshev_t(unsigned n, Real const & x, const Policy&)122{123 using result_type = tools::promote_args_t<Real>;124 using value_type = typename policies::evaluation<result_type, Policy>::type;125 using forwarding_policy = typename policies::normalise<126 Policy,127 policies::promote_float<false>,128 policies::promote_double<false>,129 policies::discrete_quantile<>,130 policies::assert_undefined<> >::type;131 132 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%)");133}134 135template <class Real>136inline tools::promote_args_t<Real> chebyshev_t(unsigned n, Real const & x)137{138 return chebyshev_t(n, x, policies::policy<>());139}140 141template <class Real, class Policy>142inline tools::promote_args_t<Real> chebyshev_u(unsigned n, Real const & x, const Policy&)143{144 using result_type = tools::promote_args_t<Real>;145 using value_type = typename policies::evaluation<result_type, Policy>::type;146 using forwarding_policy = typename policies::normalise<147 Policy,148 policies::promote_float<false>,149 policies::promote_double<false>,150 policies::discrete_quantile<>,151 policies::assert_undefined<> >::type;152 153 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%)");154}155 156template <class Real>157inline tools::promote_args_t<Real> chebyshev_u(unsigned n, Real const & x)158{159 return chebyshev_u(n, x, policies::policy<>());160}161 162template <class Real, class Policy>163inline tools::promote_args_t<Real> chebyshev_t_prime(unsigned n, Real const & x, const Policy&)164{165 using result_type = tools::promote_args_t<Real>;166 using value_type = typename policies::evaluation<result_type, Policy>::type;167 using forwarding_policy = typename policies::normalise<168 Policy,169 policies::promote_float<false>,170 policies::promote_double<false>,171 policies::discrete_quantile<>,172 policies::assert_undefined<> >::type;173 if (n == 0)174 {175 return result_type(0);176 }177 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%)");178}179 180template <class Real>181inline tools::promote_args_t<Real> chebyshev_t_prime(unsigned n, Real const & x)182{183 return chebyshev_t_prime(n, x, policies::policy<>());184}185 186/*187 * This is Algorithm 3.1 of188 * Gil, Amparo, Javier Segura, and Nico M. Temme.189 * Numerical methods for special functions.190 * Society for Industrial and Applied Mathematics, 2007.191 * https://www.siam.org/books/ot99/OT99SampleChapter.pdf192 * However, our definition of c0 differs by a factor of 1/2, as stated in the docs. . .193 */194template <class Real, class T2>195inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x)196{197 using boost::math::constants::half;198 if (length < 2)199 {200 if (length == 0)201 {202 return 0;203 }204 return c[0]/2;205 }206 Real b2 = 0;207 Real b1 = c[length -1];208 for(size_t j = length - 2; j >= 1; --j)209 {210 Real tmp = 2*x*b1 - b2 + c[j];211 b2 = b1;212 b1 = tmp;213 }214 return x*b1 - b2 + half<Real>()*c[0];215}216 217 218 219namespace detail {220template <class Real>221inline Real unchecked_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)222{223 Real t;224 Real u;225 // This cutoff is not super well defined, but it's a good estimate.226 // See "An Error Analysis of the Modified Clenshaw Method for Evaluating Chebyshev and Fourier Series"227 // J. OLIVER, IMA Journal of Applied Mathematics, Volume 20, Issue 3, November 1977, Pages 379-391228 // https://doi.org/10.1093/imamat/20.3.379229 const auto cutoff = static_cast<Real>(0.6L);230 if (x - a < b - x)231 {232 u = 2*(x-a)/(b-a);233 t = u - 1;234 if (t > -cutoff)235 {236 Real b2 = 0;237 Real b1 = c[length -1];238 for(size_t j = length - 2; j >= 1; --j)239 {240 Real tmp = 2*t*b1 - b2 + c[j];241 b2 = b1;242 b1 = tmp;243 }244 return t*b1 - b2 + c[0]/2;245 }246 else247 {248 Real b1 = c[length - 1];249 Real d = b1;250 Real b2 = 0;251 for (size_t r = length - 2; r >= 1; --r)252 {253 d = 2*u*b1 - d + c[r];254 b2 = b1;255 b1 = d - b1;256 }257 return t*b1 - b2 + c[0]/2;258 }259 }260 else261 {262 u = -2*(b-x)/(b-a);263 t = u + 1;264 if (t < cutoff)265 {266 Real b2 = 0;267 Real b1 = c[length -1];268 for(size_t j = length - 2; j >= 1; --j)269 {270 Real tmp = 2*t*b1 - b2 + c[j];271 b2 = b1;272 b1 = tmp;273 }274 return t*b1 - b2 + c[0]/2;275 }276 else277 {278 Real b1 = c[length - 1];279 Real d = b1;280 Real b2 = 0;281 for (size_t r = length - 2; r >= 1; --r)282 {283 d = 2*u*b1 + d + c[r];284 b2 = b1;285 b1 = d + b1;286 }287 return t*b1 - b2 + c[0]/2;288 }289 }290}291 292} // namespace detail293 294template <class Real>295inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)296{297 if (x < a || x > b)298 {299 BOOST_MATH_THROW_EXCEPTION(std::domain_error("x in [a, b] is required."));300 }301 if (length < 2)302 {303 if (length == 0)304 {305 return 0;306 }307 return c[0]/2;308 }309 return detail::unchecked_chebyshev_clenshaw_recurrence(c, length, a, b, x);310}311 312}} // Namespace boost::math313 314#endif // BOOST_MATH_SPECIAL_CHEBYSHEV_HPP315