math/special_functions/gegenbauer.hpp

0001 //  (C) Copyright Nick Thompson 2019.
0002 //  (C) Copyright Matt Borland 2024.
0003 //  Use, modification and distribution are subject to the
0004 //  Boost Software License, Version 1.0. (See accompanying file
0005 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0006
0007 #ifndef BOOST_MATH_SPECIAL_GEGENBAUER_HPP
0008 #define BOOST_MATH_SPECIAL_GEGENBAUER_HPP
0009
0010 #include <boost/math/tools/config.hpp>
0011 #include <boost/math/tools/type_traits.hpp>
0012 #include <boost/math/tools/numeric_limits.hpp>
0013
0014 #ifndef BOOST_MATH_NO_EXCEPTIONS
0015 #include <stdexcept>
0016 #endif
0017
0018 namespace boost { namespace math {
0019
0020 template<typename Real>
0021 BOOST_MATH_GPU_ENABLED Real gegenbauer(unsigned n, Real lambda, Real x)
0022 {
0023     static_assert(!boost::math::is_integral<Real>::value, "Gegenbauer polynomials required floating point arguments.");
0024     if (lambda <= -1/Real(2)) {
0025 #ifndef BOOST_MATH_NO_EXCEPTIONS
0026        throw std::domain_error("lambda > -1/2 is required.");
0027 #else
0028        return boost::math::numeric_limits<Real>::quiet_NaN();
0029 #endif
0030     }
0031     // 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.
0032     // I haven't observed this, but then again, I haven't managed to test an exhaustive number of parameters.
0033     // In any case, the routine is distinctly faster without this test:
0034     //if (x < 0) {
0035     //    if (n&1) {
0036     //        return -gegenbauer(n, lambda, -x);
0037     //    }
0038     //    return gegenbauer(n, lambda, -x);
0039     //}
0040
0041     if (n == 0) {
0042         return Real(1);
0043     }
0044     Real y0 = 1;
0045     Real y1 = 2*lambda*x;
0046
0047     Real yk = y1;
0048     Real k = 2;
0049     Real k_max = n*(1+boost::math::numeric_limits<Real>::epsilon());
0050     Real gamma = 2*(lambda - 1);
0051     while(k < k_max)
0052     {
0053         yk = ( (2 + gamma/k)*x*y1 - (1+gamma/k)*y0);
0054         y0 = y1;
0055         y1 = yk;
0056         k += 1;
0057     }
0058     return yk;
0059 }
0060
0061
0062 template<typename Real>
0063 BOOST_MATH_GPU_ENABLED Real gegenbauer_derivative(unsigned n, Real lambda, Real x, unsigned k)
0064 {
0065     if (k > n) {
0066         return Real(0);
0067     }
0068     Real gegen = gegenbauer<Real>(n-k, lambda + k, x);
0069     Real scale = 1;
0070     for (unsigned j = 0; j < k; ++j) {
0071         scale *= 2*lambda;
0072         lambda += 1;
0073     }
0074     return scale*gegen;
0075 }
0076
0077 template<typename Real>
0078 BOOST_MATH_GPU_ENABLED Real gegenbauer_prime(unsigned n, Real lambda, Real x) {
0079     return gegenbauer_derivative<Real>(n, lambda, x, 1);
0080 }
0081
0082
0083 }}
0084 #endif