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