85 lines · plain
1// (C) Copyright Nick Thompson 2019.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_SPECIAL_GEGENBAUER_HPP8#define BOOST_MATH_SPECIAL_GEGENBAUER_HPP9 10#include <boost/math/tools/config.hpp>11#include <boost/math/tools/type_traits.hpp>12#include <boost/math/tools/numeric_limits.hpp>13 14#ifndef BOOST_MATH_NO_EXCEPTIONS15#include <stdexcept>16#endif17 18namespace boost { namespace math {19 20template<typename Real>21BOOST_MATH_GPU_ENABLED Real gegenbauer(unsigned n, Real lambda, Real x)22{23 static_assert(!boost::math::is_integral<Real>::value, "Gegenbauer polynomials required floating point arguments.");24 if (lambda <= -1/Real(2)) {25#ifndef BOOST_MATH_NO_EXCEPTIONS26 throw std::domain_error("lambda > -1/2 is required.");27#else28 return boost::math::numeric_limits<Real>::quiet_NaN();29#endif30 }31 // 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.32 // I haven't observed this, but then again, I haven't managed to test an exhaustive number of parameters.33 // In any case, the routine is distinctly faster without this test:34 //if (x < 0) {35 // if (n&1) {36 // return -gegenbauer(n, lambda, -x);37 // }38 // return gegenbauer(n, lambda, -x);39 //}40 41 if (n == 0) {42 return Real(1);43 }44 Real y0 = 1;45 Real y1 = 2*lambda*x;46 47 Real yk = y1;48 Real k = 2;49 Real k_max = n*(1+boost::math::numeric_limits<Real>::epsilon());50 Real gamma = 2*(lambda - 1);51 while(k < k_max)52 {53 yk = ( (2 + gamma/k)*x*y1 - (1+gamma/k)*y0);54 y0 = y1;55 y1 = yk;56 k += 1;57 }58 return yk;59}60 61 62template<typename Real>63BOOST_MATH_GPU_ENABLED Real gegenbauer_derivative(unsigned n, Real lambda, Real x, unsigned k)64{65 if (k > n) {66 return Real(0);67 }68 Real gegen = gegenbauer<Real>(n-k, lambda + k, x);69 Real scale = 1;70 for (unsigned j = 0; j < k; ++j) {71 scale *= 2*lambda;72 lambda += 1;73 }74 return scale*gegen;75}76 77template<typename Real>78BOOST_MATH_GPU_ENABLED Real gegenbauer_prime(unsigned n, Real lambda, Real x) {79 return gegenbauer_derivative<Real>(n, lambda, x, 1);80}81 82 83}}84#endif85