]> www.ginac.de Git - cln.git/blob - examples/legendre.cc
Move GETVAL macro from acinclude.m4 to m4/getval.m4.
[cln.git] / 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
13 using namespace std;
14 using 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
19 const 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
37 const 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
55 int 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 }