boost/math/special_functions/gegenbauer.hpp

   1 //  (C) Copyright Nick Thompson 2019.
   2 //  Use, modification and distribution are subject to the
   3 //  Boost Software License, Version 1.0. (See accompanying file
   4 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   5
   6 #ifndef BOOST_MATH_SPECIAL_GEGENBAUER_HPP
   7 #define BOOST_MATH_SPECIAL_GEGENBAUER_HPP
   8
   9 #include <stdexcept>
  10 #include <type_traits>
  11
  12 namespace boost { namespace math {
  13
  14 template<typename Real>
  15 Real gegenbauer(unsigned n, Real lambda, Real x)
  16 {
  17     static_assert(!std::is_integral<Real>::value, "Gegenbauer polynomials required floating point arguments.");
  18     if (lambda <= -1/Real(2)) {
  19         throw std::domain_error("lambda > -1/2 is required.");
  20     }
  21     // 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.
  22     // I haven't observed this, but then again, I haven't managed to test an exhaustive number of parameters.
  23     // In any case, the routine is distinctly faster without this test:
  24     //if (x < 0) {
  25     //    if (n&1) {
  26     //        return -gegenbauer(n, lambda, -x);
  27     //    }
  28     //    return gegenbauer(n, lambda, -x);
  29     //}
  30
  31     if (n == 0) {
  32         return Real(1);
  33     }
  34     Real y0 = 1;
  35     Real y1 = 2*lambda*x;
  36
  37     Real yk = y1;
  38     Real k = 2;
  39     Real k_max = n*(1+std::numeric_limits<Real>::epsilon());
  40     Real gamma = 2*(lambda - 1);
  41     while(k < k_max)
  42     {
  43         yk = ( (2 + gamma/k)*x*y1 - (1+gamma/k)*y0);
  44         y0 = y1;
  45         y1 = yk;
  46         k += 1;
  47     }
  48     return yk;
  49 }
  50
  51
  52 template<typename Real>
  53 Real gegenbauer_derivative(unsigned n, Real lambda, Real x, unsigned k)
  54 {
  55     if (k > n) {
  56         return Real(0);
  57     }
  58     Real gegen = gegenbauer<Real>(n-k, lambda + k, x);
  59     Real scale = 1;
  60     for (unsigned j = 0; j < k; ++j) {
  61         scale *= 2*lambda;
  62         lambda += 1;
  63     }
  64     return scale*gegen;
  65 }
  66
  67 template<typename Real>
  68 Real gegenbauer_prime(unsigned n, Real lambda, Real x) {
  69     return gegenbauer_derivative<Real>(n, lambda, x, 1);
  70 }
  71
  72
  73 }}
  74 #endif