boost/math/special_functions/jacobi.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_JACOBI_HPP
   7 #define BOOST_MATH_SPECIAL_JACOBI_HPP
   8
   9 #include <limits>
  10 #include <stdexcept>
  11
  12 namespace boost { namespace math {
  13
  14 template<typename Real>
  15 Real jacobi(unsigned n, Real alpha, Real beta, Real x)
  16 {
  17     static_assert(!std::is_integral<Real>::value, "Jacobi polynomials do not work with integer arguments.");
  18
  19     if (n == 0) {
  20         return Real(1);
  21     }
  22     Real y0 = 1;
  23     Real y1 = (alpha+1) + (alpha+beta+2)*(x-1)/Real(2);
  24
  25     Real yk = y1;
  26     Real k = 2;
  27     Real k_max = n*(1+std::numeric_limits<Real>::epsilon());
  28     while(k < k_max)
  29     {
  30         // Hoping for lots of common subexpression elimination by the compiler:
  31         Real denom = 2*k*(k+alpha+beta)*(2*k+alpha+beta-2);
  32         Real gamma1 = (2*k+alpha+beta-1)*( (2*k+alpha+beta)*(2*k+alpha+beta-2)*x + alpha*alpha -beta*beta);
  33         Real gamma0 = -2*(k+alpha-1)*(k+beta-1)*(2*k+alpha+beta);
  34         yk = (gamma1*y1 + gamma0*y0)/denom;
  35         y0 = y1;
  36         y1 = yk;
  37         k += 1;
  38     }
  39     return yk;
  40 }
  41
  42 template<typename Real>
  43 Real jacobi_derivative(unsigned n, Real alpha, Real beta, Real x, unsigned k)
  44 {
  45     if (k > n) {
  46         return Real(0);
  47     }
  48     Real scale = 1;
  49     for(unsigned j = 1; j <= k; ++j) {
  50         scale *= (alpha + beta + n + j)/2;
  51     }
  52
  53     return scale*jacobi<Real>(n-k, alpha + k, beta+k, x);
  54 }
  55
  56 template<typename Real>
  57 Real jacobi_prime(unsigned n, Real alpha, Real beta, Real x)
  58 {
  59     return jacobi_derivative<Real>(n, alpha, beta, x, 1);
  60 }
  61
  62 template<typename Real>
  63 Real jacobi_double_prime(unsigned n, Real alpha, Real beta, Real x)
  64 {
  65     return jacobi_derivative<Real>(n, alpha, beta, x, 2);
  66 }
  67
  68 }}
  69 #endif