math/special_functions/chebyshev.hpp

0001 //  (C) Copyright Nick Thompson 2017.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
0007 #define BOOST_MATH_SPECIAL_CHEBYSHEV_HPP
0008 #include <cmath>
0009 #include <type_traits>
0010 #include <boost/math/special_functions/math_fwd.hpp>
0011 #include <boost/math/policies/error_handling.hpp>
0012 #include <boost/math/constants/constants.hpp>
0013 #include <boost/math/tools/promotion.hpp>
0014 #include <boost/math/tools/throw_exception.hpp>
0015
0016 #if (__cplusplus > 201103) || (defined(_CPPLIB_VER) && (_CPPLIB_VER >= 610))
0017 #  define BOOST_MATH_CHEB_USE_STD_ACOSH
0018 #endif
0019
0020 #ifndef BOOST_MATH_CHEB_USE_STD_ACOSH
0021 #  include <boost/math/special_functions/acosh.hpp>
0022 #endif
0023
0024 namespace boost { namespace math {
0025
0026 template <class T1, class T2, class T3>
0027 inline tools::promote_args_t<T1, T2, T3> chebyshev_next(T1 const & x, T2 const & Tn, T3 const & Tn_1)
0028 {
0029     return 2*x*Tn - Tn_1;
0030 }
0031
0032 namespace detail {
0033
0034 // https://stackoverflow.com/questions/5625431/efficient-way-to-compute-pq-exponentiation-where-q-is-an-integer
0035 template <typename T, typename std::enable_if<std::is_arithmetic<T>::value, bool>::type = true>
0036 T expt(T p, unsigned q)
0037 {
0038     T r = 1;
0039
0040     while (q != 0) {
0041         if (q % 2 == 1) {    // q is odd
0042             r *= p;
0043             q--;
0044         }
0045         p *= p;
0046         q /= 2;
0047     }
0048
0049     return r;
0050 }
0051
0052 template <typename T, typename std::enable_if<!std::is_arithmetic<T>::value, bool>::type = true>
0053 T expt(T p, unsigned q)
0054 {
0055     using std::pow;
0056     return pow(p, static_cast<int>(q));
0057 }
0058
0059 template<class Real, bool second, class Policy>
0060 inline Real chebyshev_imp(unsigned n, Real const & x, const Policy&)
0061 {
0062 #ifdef BOOST_MATH_CHEB_USE_STD_ACOSH
0063     using std::acosh;
0064 #define BOOST_MATH_ACOSH_POLICY
0065 #else
0066    using boost::math::acosh;
0067 #define BOOST_MATH_ACOSH_POLICY , Policy()
0068 #endif
0069     using std::cosh;
0070     using std::pow;
0071     using std::sqrt;
0072     Real T0 = 1;
0073     Real T1;
0074
0075     BOOST_IF_CONSTEXPR (second)
0076     {
0077         if (x > 1 || x < -1)
0078         {
0079             Real t = sqrt(x*x -1);
0080             return static_cast<Real>((expt(static_cast<Real>(x+t), n+1) - expt(static_cast<Real>(x-t), n+1))/(2*t));
0081         }
0082         T1 = 2*x;
0083     }
0084     else
0085     {
0086         if (x > 1)
0087         {
0088             return cosh(n*acosh(x BOOST_MATH_ACOSH_POLICY));
0089         }
0090         if (x < -1)
0091         {
0092             if (n & 1)
0093             {
0094                 return -cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
0095             }
0096             else
0097             {
0098                 return cosh(n*acosh(-x BOOST_MATH_ACOSH_POLICY));
0099             }
0100         }
0101         T1 = x;
0102     }
0103
0104     if (n == 0)
0105     {
0106         return T0;
0107     }
0108
0109     unsigned l = 1;
0110     while(l < n)
0111     {
0112        std::swap(T0, T1);
0113        T1 = static_cast<Real>(boost::math::chebyshev_next(x, T0, T1));
0114        ++l;
0115     }
0116     return T1;
0117 }
0118 } // namespace detail
0119
0120 template <class Real, class Policy>
0121 inline tools::promote_args_t<Real> chebyshev_t(unsigned n, Real const & x, const Policy&)
0122 {
0123    using result_type = tools::promote_args_t<Real>;
0124    using value_type = typename policies::evaluation<result_type, Policy>::type;
0125    using forwarding_policy = typename policies::normalise<
0126                                                             Policy,
0127                                                             policies::promote_float<false>,
0128                                                             policies::promote_double<false>,
0129                                                             policies::discrete_quantile<>,
0130                                                             policies::assert_undefined<> >::type;
0131
0132    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%)");
0133 }
0134
0135 template <class Real>
0136 inline tools::promote_args_t<Real> chebyshev_t(unsigned n, Real const & x)
0137 {
0138     return chebyshev_t(n, x, policies::policy<>());
0139 }
0140
0141 template <class Real, class Policy>
0142 inline tools::promote_args_t<Real> chebyshev_u(unsigned n, Real const & x, const Policy&)
0143 {
0144    using result_type = tools::promote_args_t<Real>;
0145    using value_type = typename policies::evaluation<result_type, Policy>::type;
0146    using forwarding_policy =  typename policies::normalise<
0147                                                             Policy,
0148                                                             policies::promote_float<false>,
0149                                                             policies::promote_double<false>,
0150                                                             policies::discrete_quantile<>,
0151                                                             policies::assert_undefined<> >::type;
0152
0153    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%)");
0154 }
0155
0156 template <class Real>
0157 inline tools::promote_args_t<Real> chebyshev_u(unsigned n, Real const & x)
0158 {
0159     return chebyshev_u(n, x, policies::policy<>());
0160 }
0161
0162 template <class Real, class Policy>
0163 inline tools::promote_args_t<Real> chebyshev_t_prime(unsigned n, Real const & x, const Policy&)
0164 {
0165    using result_type = tools::promote_args_t<Real>;
0166    using value_type = typename policies::evaluation<result_type, Policy>::type;
0167    using forwarding_policy = typename policies::normalise<
0168                                                             Policy,
0169                                                             policies::promote_float<false>,
0170                                                             policies::promote_double<false>,
0171                                                             policies::discrete_quantile<>,
0172                                                             policies::assert_undefined<> >::type;
0173    if (n == 0)
0174    {
0175       return result_type(0);
0176    }
0177    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%)");
0178 }
0179
0180 template <class Real>
0181 inline tools::promote_args_t<Real> chebyshev_t_prime(unsigned n, Real const & x)
0182 {
0183    return chebyshev_t_prime(n, x, policies::policy<>());
0184 }
0185
0186 /*
0187  * This is Algorithm 3.1 of
0188  * Gil, Amparo, Javier Segura, and Nico M. Temme.
0189  * Numerical methods for special functions.
0190  * Society for Industrial and Applied Mathematics, 2007.
0191  * https://www.siam.org/books/ot99/OT99SampleChapter.pdf
0192  * However, our definition of c0 differs by a factor of 1/2, as stated in the docs. . .
0193  */
0194 template <class Real, class T2>
0195 inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const T2& x)
0196 {
0197     using boost::math::constants::half;
0198     if (length < 2)
0199     {
0200         if (length == 0)
0201         {
0202             return 0;
0203         }
0204         return c[0]/2;
0205     }
0206     Real b2 = 0;
0207     Real b1 = c[length -1];
0208     for(size_t j = length - 2; j >= 1; --j)
0209     {
0210         Real tmp = 2*x*b1 - b2 + c[j];
0211         b2 = b1;
0212         b1 = tmp;
0213     }
0214     return x*b1 - b2 + half<Real>()*c[0];
0215 }
0216
0217
0218
0219 namespace detail {
0220 template <class Real>
0221 inline Real unchecked_chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
0222 {
0223     Real t;
0224     Real u;
0225     // This cutoff is not super well defined, but it's a good estimate.
0226     // See "An Error Analysis of the Modified Clenshaw Method for Evaluating Chebyshev and Fourier Series"
0227     // J. OLIVER, IMA Journal of Applied Mathematics, Volume 20, Issue 3, November 1977, Pages 379-391
0228     // https://doi.org/10.1093/imamat/20.3.379
0229     const auto cutoff = static_cast<Real>(0.6L);
0230     if (x - a < b - x)
0231     {
0232         u = 2*(x-a)/(b-a);
0233         t = u - 1;
0234         if (t > -cutoff)
0235         {
0236             Real b2 = 0;
0237             Real b1 = c[length -1];
0238             for(size_t j = length - 2; j >= 1; --j)
0239             {
0240                 Real tmp = 2*t*b1 - b2 + c[j];
0241                 b2 = b1;
0242                 b1 = tmp;
0243             }
0244             return t*b1 - b2 + c[0]/2;
0245         }
0246         else
0247         {
0248             Real b1 = c[length - 1];
0249             Real d = b1;
0250             Real b2 = 0;
0251             for (size_t r = length - 2; r >= 1; --r)
0252             {
0253                 d = 2*u*b1 - d + c[r];
0254                 b2 = b1;
0255                 b1 = d - b1;
0256             }
0257             return t*b1 - b2 + c[0]/2;
0258         }
0259     }
0260     else
0261     {
0262         u = -2*(b-x)/(b-a);
0263         t = u + 1;
0264         if (t < cutoff)
0265         {
0266             Real b2 = 0;
0267             Real b1 = c[length -1];
0268             for(size_t j = length - 2; j >= 1; --j)
0269             {
0270                 Real tmp = 2*t*b1 - b2 + c[j];
0271                 b2 = b1;
0272                 b1 = tmp;
0273             }
0274             return t*b1 - b2 + c[0]/2;
0275         }
0276         else
0277         {
0278             Real b1 = c[length - 1];
0279             Real d = b1;
0280             Real b2 = 0;
0281             for (size_t r = length - 2; r >= 1; --r)
0282             {
0283                 d = 2*u*b1 + d + c[r];
0284                 b2 = b1;
0285                 b1 = d + b1;
0286             }
0287             return t*b1 - b2 + c[0]/2;
0288         }
0289     }
0290 }
0291
0292 } // namespace detail
0293
0294 template <class Real>
0295 inline Real chebyshev_clenshaw_recurrence(const Real* const c, size_t length, const Real & a, const Real & b, const Real& x)
0296 {
0297     if (x < a || x > b)
0298     {
0299        BOOST_MATH_THROW_EXCEPTION(std::domain_error("x in [a, b] is required."));
0300     }
0301     if (length < 2)
0302     {
0303         if (length == 0)
0304         {
0305             return 0;
0306         }
0307         return c[0]/2;
0308     }
0309     return detail::unchecked_chebyshev_clenshaw_recurrence(c, length, a, b, x);
0310 }
0311
0312 }} // Namespace boost::math
0313
0314 #endif // BOOST_MATH_SPECIAL_CHEBYSHEV_HPP