src/polynomial/misc/cl_UP_RA_legendre.cc

   1 // legendre().
   2
   3 // General includes.
   4 #include "cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/univpoly_rational.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/integer.h"
  13 #include "cln/rational.h"
  14
  15 namespace cln {
  16
  17 const cl_UP_RA legendre (sintL n)
  18 {
  19 // The Legendre polynomials P_n(x) are defined as
  20 //
  21 //                1   ( d  ) n          n
  22 //    P_n(x) = ------ (----)   (x^2 - 1)
  23 //             2^n n! ( dx )
  24 //
  25 // They satisfy the recurrence relation
  26 //
  27 //    P_0(x) = 1
  28 //    P_{n+1}(x) = 1/(n+1) * ((2n+1)x P_n(x) - n P_{n-1}(x)) for n >= 0.
  29 //
  30 // Theorem:
  31 //    P_n(x) satisfies the differential equation
  32 //    (1-x^2)*P_n''(x) - 2x*P_n'(x) + (n^2+n)*P_n(x) = 0.
  33 //
  34 // Proof: See elsewhere.
  35 //
  36 // Corollary:
  37 //    The coefficients c_{n,k} of P_n(x) = sum(k=0..n, c_{n,k} x^k)
  38 //    satisfy:
  39 //       c_{n,n} = (2n)!/(n!^2 2^n),
  40 //       c_{n,n-1} = 0,
  41 //       c_{n,k} = (k+1)(k+2)/(k-n)(k+n+1)*c_{n,k+2}
  42 //
  43 // It follows that for n>=0
  44 //
  45 // P_n(x) = sum(j=0..floor(n/2), (-1)^j (2n-2j)!/j!(n-2j)!(n-j)! 2^-n x^(n-2j))
  46 //
  47         var cl_univpoly_rational_ring R = find_univpoly_ring(cl_RA_ring);
  48         var cl_UP_RA p = R->create(n);
  49         var cl_I denom = ash(1,n);
  50         var sintL k = n;
  51         var cl_I c_k = binomial(2*n,n);
  52         for (;;) {
  53                 p.set_coeff(k,c_k/denom);
  54                 k = k-2;
  55                 if (k < 0)
  56                         break;
  57                 c_k = exquo((cl_I)(k+1) * (cl_I)(k+2) * c_k,
  58                             (cl_I)(k-n) * (cl_I)(k+n+1));
  59         }
  60         p.finalize();
  61         return p;
  62 }
  63
  64 }  // namespace cln