src/polynomial/misc/cl_UP_I_tchebychev.cc

   1 // tschebychev().
   2
   3 // General includes.
   4 #include "base/cl_sysdep.h"
   5
   6 // Specification.
   7 #include "cln/univpoly_integer.h"
   8
   9
  10 // Implementation.
  11
  12 #include "cln/integer.h"
  13
  14 namespace cln {
  15
  16 const cl_UP_I tschebychev (sintL n)
  17 {
  18 // The Tschebychev polynomials (of the 1st kind) T_n(x) are defined
  19 // through the recurrence relation
  20 //
  21 //    T_0(x) = 1
  22 //    T_1(x) = x
  23 //    T_{n+2}(x) = 2x T_{n+1}(x) - T_n(x) for n >= 0.
  24 //
  25 // Theorem:
  26 //    T_n(x) satisfies the differential equation
  27 //    (x^2-1)*T_n''(x) + x*T_n'(x) - n^2*T_n(x) = 0.
  28 //
  29 // Proof: See elsewhere.
  30 //
  31 // Corollary:
  32 //    The coefficients c_{n,k} of T_n(x) = sum(k=0..n, c_{n,k} x^k)
  33 //    satisfy:
  34 //       c_{n,n} = 2^(n-1) for n>=1, 1 for n=0,
  35 //       c_{n,n-1} = 0,
  36 //       c_{n,k} = (k+1)(k+2)/(k^2-n^2)*c_{n,k+2}
  37 //
  38 // It follows that for n>0
  39 //
  40 // T_n(x) = sum(j=0..floor(n/2), (-1)^j (n-j-1)!n/j!(n-2j)! 2^(n-2j-1) x^(n-2j))
  41 //
  42         var cl_univpoly_integer_ring R = find_univpoly_ring(cl_I_ring);
  43         if (n == 0)
  44                 return R->one();
  45         var cl_UP_I t = R->create(n);
  46         var sintL k = n;
  47         var cl_I c_k = ash(1,n-1);
  48         for (;;) {
  49                 t.set_coeff(k,c_k);
  50                 k = k-2;
  51                 if (k < 0)
  52                         break;
  53                 c_k = exquo((cl_I)(k+1) * (cl_I)(k+2) * c_k,
  54                             (cl_I)(k-n) * (cl_I)(k+n));
  55         }
  56         t.finalize();
  57         return t;
  58 }
  59
  60 }  // namespace cln