PageRenderTime 49ms CodeModel.GetById 19ms RepoModel.GetById 0ms app.codeStats 0ms

/cln-1.3.2/examples/legendre.cc

#
C++ | 75 lines | 66 code | 6 blank | 3 comment | 12 complexity | 1ac17740850628250d313b3d22adfeb4 MD5 | raw file
Possible License(s): GPL-2.0
  1. // Compute the Legendre polynomials.
  2. #include <cln/number.h>
  3. #include <cln/integer.h>
  4. #include <cln/rational.h>
  5. #include <cln/univpoly.h>
  6. #include <cln/modinteger.h>
  7. #include <cln/univpoly_rational.h>
  8. #include <cln/univpoly_modint.h>
  9. #include <cln/io.h>
  10. #include <cstdlib>
  11. using namespace std;
  12. using namespace cln;
  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.)
  15. const cl_UP_RA legendre (const cl_rational_ring& R, int n)
  16. {
  17. cl_univpoly_rational_ring PR = find_univpoly_ring(R);
  18. cl_UP_RA b = PR->create(2);
  19. b.set_coeff(2,1);
  20. b.set_coeff(1,0);
  21. b.set_coeff(0,-1);
  22. b.finalize(); // b is now x^2-1
  23. cl_UP_RA p = (n==0 ? PR->one() : expt_pos(b,n));
  24. for (int i = 0; i < n; i++)
  25. p = deriv(p);
  26. cl_RA factor = recip(factorial(n)*ash(1,n));
  27. for (int j = degree(p); j >= 0; j--)
  28. p.set_coeff(j, coeff(p,j) * factor);
  29. p.finalize();
  30. return p;
  31. }
  32. const cl_UP_MI legendre (const cl_modint_ring& R, int n)
  33. {
  34. cl_univpoly_modint_ring PR = find_univpoly_ring(R);
  35. cl_UP_MI b = PR->create(2);
  36. b.set_coeff(2,R->canonhom(1));
  37. b.set_coeff(1,R->canonhom(0));
  38. b.set_coeff(0,R->canonhom(-1));
  39. b.finalize(); // b is now x^2-1
  40. cl_UP_MI p = (n==0 ? PR->one() : expt_pos(b,n));
  41. for (int i = 0; i < n; i++)
  42. p = deriv(p);
  43. cl_MI factor = recip(R->canonhom(factorial(n)*ash(1,n)));
  44. for (int j = degree(p); j >= 0; j--)
  45. p.set_coeff(j, coeff(p,j) * factor);
  46. p.finalize();
  47. return p;
  48. }
  49. int main (int argc, char* argv[])
  50. {
  51. if (!(argc == 2 || argc == 3)) {
  52. cerr << "Usage: legendre n [m]" << endl;
  53. exit(1);
  54. }
  55. int n = atoi(argv[1]);
  56. if (!(n >= 0)) {
  57. cerr << "Usage: legendre n [m] with n >= 0" << endl;
  58. exit(1);
  59. }
  60. if (argc == 2) {
  61. cl_UP p = legendre(cl_RA_ring,n);
  62. cout << p << endl;
  63. } else {
  64. cl_I m = argv[2];
  65. cl_UP p = legendre(find_modint_ring(m),n);
  66. cout << p << endl;
  67. }
  68. return 0;
  69. }