examples/legendre.cc

   1 // Compute the Legendre polynomials.
   2
   3 #include <cln/number.h>
   4 #include <cln/integer.h>
   5 #include <cln/rational.h>
   6 #include <cln/univpoly.h>
   7 #include <cln/modinteger.h>
   8 #include <cln/univpoly_rational.h>
   9 #include <cln/univpoly_modint.h>
  10 #include <cln/io.h>
  11 #include <stdlib.h>
  12
  13 using namespace cln;
  14
  15 // Computes the n-th Legendre polynomial in R[x], using the formula
  16 // P_n(x) = 1/(2^n n!) * (d/dx)^n (x^2-1)^n. (Assume n >= 0.)
  17
  18 const cl_UP_RA legendre (const cl_rational_ring& R, int n)
  19 {
  20         cl_univpoly_rational_ring PR = find_univpoly_ring(R);
  21         cl_UP_RA b = PR->create(2);
  22         b.set_coeff(2,1);
  23         b.set_coeff(1,0);
  24         b.set_coeff(0,-1);
  25         b.finalize(); // b is now x^2-1
  26         cl_UP_RA p = (n==0 ? PR->one() : expt_pos(b,n));
  27         for (int i = 0; i < n; i++)
  28                 p = deriv(p);
  29         cl_RA factor = recip(factorial(n)*ash(1,n));
  30         for (int j = degree(p); j >= 0; j--)
  31                 p.set_coeff(j, coeff(p,j) * factor);
  32         p.finalize();
  33         return p;
  34 }
  35
  36 const cl_UP_MI legendre (const cl_modint_ring& R, int n)
  37 {
  38         cl_univpoly_modint_ring PR = find_univpoly_ring(R);
  39         cl_UP_MI b = PR->create(2);
  40         b.set_coeff(2,R->canonhom(1));
  41         b.set_coeff(1,R->canonhom(0));
  42         b.set_coeff(0,R->canonhom(-1));
  43         b.finalize(); // b is now x^2-1
  44         cl_UP_MI p = (n==0 ? PR->one() : expt_pos(b,n));
  45         for (int i = 0; i < n; i++)
  46                 p = deriv(p);
  47         cl_MI factor = recip(R->canonhom(factorial(n)*ash(1,n)));
  48         for (int j = degree(p); j >= 0; j--)
  49                 p.set_coeff(j, coeff(p,j) * factor);
  50         p.finalize();
  51         return p;
  52 }
  53
  54 int main (int argc, char* argv[])
  55 {
  56         if (!(argc == 2 || argc == 3)) {
  57                 fprint(stderr, "Usage: legendre n [m]\n");
  58                 exit(1);
  59         }
  60         int n = atoi(argv[1]);
  61         if (!(n >= 0)) {
  62                 fprint(stderr, "Usage: legendre n [m]  with n >= 0\n");
  63                 exit(1);
  64         }
  65         if (argc == 2) {
  66                 cl_UP p = legendre(cl_RA_ring,n);
  67                 fprint(stdout, p);
  68         } else {
  69                 cl_I m = argv[2];
  70                 cl_UP p = legendre(find_modint_ring(m),n);
  71                 fprint(stdout, p);
  72         }
  73         fprint(stdout, "\n");
  74 }