math/special_functions/gegenbauer.hpp

0001 //  (C) Copyright Nick Thompson 2019.
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_GEGENBAUER_HPP
0007 #define BOOST_MATH_SPECIAL_GEGENBAUER_HPP
0008
0009 #include <limits>
0010 #include <stdexcept>
0011 #include <type_traits>
0012
0013 namespace boost { namespace math {
0014
0015 template<typename Real>
0016 Real gegenbauer(unsigned n, Real lambda, Real x)
0017 {
0018     static_assert(!std::is_integral<Real>::value, "Gegenbauer polynomials required floating point arguments.");
0019     if (lambda <= -1/Real(2)) {
0020 #ifndef BOOST_NO_EXCEPTIONS
0021        throw std::domain_error("lambda > -1/2 is required.");
0022 #else
0023        return std::numeric_limits<Real>::quiet_NaN();
0024 #endif
0025     }
0026     // The only reason to do this is because of some instability that could be present for x < 0 that is not present for x > 0.
0027     // I haven't observed this, but then again, I haven't managed to test an exhaustive number of parameters.
0028     // In any case, the routine is distinctly faster without this test:
0029     //if (x < 0) {
0030     //    if (n&1) {
0031     //        return -gegenbauer(n, lambda, -x);
0032     //    }
0033     //    return gegenbauer(n, lambda, -x);
0034     //}
0035
0036     if (n == 0) {
0037         return Real(1);
0038     }
0039     Real y0 = 1;
0040     Real y1 = 2*lambda*x;
0041
0042     Real yk = y1;
0043     Real k = 2;
0044     Real k_max = n*(1+std::numeric_limits<Real>::epsilon());
0045     Real gamma = 2*(lambda - 1);
0046     while(k < k_max)
0047     {
0048         yk = ( (2 + gamma/k)*x*y1 - (1+gamma/k)*y0);
0049         y0 = y1;
0050         y1 = yk;
0051         k += 1;
0052     }
0053     return yk;
0054 }
0055
0056
0057 template<typename Real>
0058 Real gegenbauer_derivative(unsigned n, Real lambda, Real x, unsigned k)
0059 {
0060     if (k > n) {
0061         return Real(0);
0062     }
0063     Real gegen = gegenbauer<Real>(n-k, lambda + k, x);
0064     Real scale = 1;
0065     for (unsigned j = 0; j < k; ++j) {
0066         scale *= 2*lambda;
0067         lambda += 1;
0068     }
0069     return scale*gegen;
0070 }
0071
0072 template<typename Real>
0073 Real gegenbauer_prime(unsigned n, Real lambda, Real x) {
0074     return gegenbauer_derivative<Real>(n, lambda, x, 1);
0075 }
0076
0077
0078 }}
0079 #endif