1 // Compute the Legendre polynomials.
4 #include <cl_integer.h>
5 #include <cl_rational.h>
6 #include <cl_univpoly.h>
7 #include <cl_modinteger.h>
8 #include <cl_univpoly_rational.h>
9 #include <cl_univpoly_modint.h>
13 // Computes the n-th Legendre polynomial in R[x], using the formula
14 // P_n(x) = 1/(2^n n!) * (d/dx)^n (x^2-1)^n. (Assume n >= 0.)
16 const cl_UP_RA legendre (const cl_rational_ring& R, int n)
18 cl_univpoly_rational_ring PR = cl_find_univpoly_ring(R);
19 cl_UP_RA b = PR->create(2);
23 b.finalize(); // b is now x^2-1
24 cl_UP_RA p = (n==0 ? PR->one() : expt_pos(b,n));
25 for (int i = 0; i < n; i++)
27 cl_RA factor = recip(factorial(n)*ash(1,n));
28 for (int j = degree(p); j >= 0; j--)
29 p.set_coeff(j, coeff(p,j) * factor);
34 const cl_UP_MI legendre (const cl_modint_ring& R, int n)
36 cl_univpoly_modint_ring PR = cl_find_univpoly_ring(R);
37 cl_UP_MI b = PR->create(2);
38 b.set_coeff(2,R->canonhom(1));
39 b.set_coeff(1,R->canonhom(0));
40 b.set_coeff(0,R->canonhom(-1));
41 b.finalize(); // b is now x^2-1
42 cl_UP_MI p = (n==0 ? PR->one() : expt_pos(b,n));
43 for (int i = 0; i < n; i++)
45 cl_MI factor = recip(R->canonhom(factorial(n)*ash(1,n)));
46 for (int j = degree(p); j >= 0; j--)
47 p.set_coeff(j, coeff(p,j) * factor);
52 int main (int argc, char* argv[])
54 if (!(argc == 2 || argc == 3)) {
55 fprint(cl_stderr, "Usage: legendre n [m]\n");
58 int n = atoi(argv[1]);
60 fprint(cl_stderr, "Usage: legendre n [m] with n >= 0\n");
64 cl_UP p = legendre(cl_RA_ring,n);
68 cl_UP p = legendre(cl_find_modint_ring(m),n);
71 fprint(cl_stdout, "\n");