ginac/polynomial/pgcd.cpp

   1 #include "pgcd.h"
   2 #include "collect_vargs.h"
   3 #include "smod_helpers.h"
   4 #include "euclid_gcd_wrap.h"
   5 #include "eval_point_finder.h"
   6 #include "newton_interpolate.h"
   7 #include "divide_in_z_p.h"
   8
   9 namespace GiNaC
  10 {
  11
  12 extern void
  13 primpart_content(ex& pp, ex& c, ex e, const exvector& vars, const long p);
  14
  15 // Computes the GCD of two polynomials over a prime field.
  16 // Based on Algorithm 7.2 from "Algorithms for Computer Algebra"
  17 // A and B are considered as Z_p[x_n][x_0, \ldots, x_{n-1}], that is,
  18 // as a polynomials in variables x_0, \ldots x_{n-1} having coefficients
  19 // from the ring Z_p[x_n]
  20 ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
  21 {
  22         static const ex ex1(1);
  23         if (A.is_zero())
  24                 return B;
  25
  26         if (B.is_zero())
  27                 return A;
  28
  29         if (is_a<numeric>(A))
  30                 return ex1;
  31
  32         if (is_a<numeric>(B))
  33                 return ex1;
  34
  35         // Checks for univariate polynomial
  36         if (vars.size() == 1) {
  37                 ex ret = euclid_gcd(A, B, vars[0], p); // Univariate GCD
  38                 return ret;
  39         }
  40         const ex& mainvar(vars.back());
  41
  42         // gcd of the contents
  43         ex H = 0, Hprev = 0; // GCD candidate
  44         ex newton_poly = 1;  // for Newton Interpolation
  45
  46         // Contents and primparts of A and B
  47         ex Aprim, contA;
  48         primpart_content(Aprim, contA, A, vars, p);
  49         ex Bprim, contB;
  50         primpart_content(Bprim, contB, B, vars, p);
  51         // gcd of univariate polynomials
  52         const ex cont_gcd = euclid_gcd(contA, contB, mainvar, p);
  53
  54         exvector restvars = vars;
  55         restvars.pop_back();
  56         const ex AL = lcoeff_wrt(Aprim, restvars);
  57         const ex BL = lcoeff_wrt(Bprim, restvars);
  58         // gcd of univariate polynomials
  59         const ex lc_gcd = euclid_gcd(AL, BL, mainvar, p);
  60
  61         // The estimate of degree of the gcd of Ab and Bb
  62         int gcd_deg = std::min(degree(Aprim, mainvar), degree(Bprim, mainvar));
  63         eval_point_finder::value_type b;
  64
  65         eval_point_finder find_eval_point(p);
  66         const numeric pn(p);
  67         do {
  68                 // Find a `good' evaluation point b.
  69                 bool has_more_pts = find_eval_point(b, lc_gcd, mainvar);
  70                 // If there are no more possible evaluation points, bail out
  71                 if (!has_more_pts)
  72                         throw pgcd_failed();
  73
  74                 const numeric bn(b);
  75                 // Evaluate the polynomials in b
  76                 ex Ab = Aprim.subs(mainvar == bn).smod(pn);
  77                 ex Bb = Bprim.subs(mainvar == bn).smod(pn);
  78                 ex Cb = pgcd(Ab, Bb, restvars, p);
  79
  80                 // Set the correct the leading coefficient
  81                 const cln::cl_I lcb_gcd =
  82                         smod(to_cl_I(lc_gcd.subs(mainvar == bn)), p);
  83                 const cln::cl_I Cblc = integer_lcoeff(Cb, restvars);
  84                 const cln::cl_I correct_lc = smod(lcb_gcd*recip(Cblc, p), p);
  85                 Cb = (Cb*numeric(correct_lc)).smod(pn);
  86
  87                 // Test for unlucky homomorphisms
  88                 const int img_gcd_deg = Cb.degree(restvars.back());
  89                 if (img_gcd_deg < gcd_deg) {
  90                         // The degree decreased, previous homomorphisms were
  91                         // bad, so we have to start it all over.
  92                         H = Cb;
  93                         newton_poly = mainvar - numeric(b);
  94                         Hprev = 0;
  95                         gcd_deg  = img_gcd_deg;
  96                         continue;
  97                 }
  98                 if (img_gcd_deg > gcd_deg) {
  99                         // The degree of images GCD is too high, this
 100                         // evaluation point is bad. Skip it.
 101                         continue;
 102                 }
 103
 104                 // Image has the same degree as the previous one
 105                 // (or at least not higher than the limit)
 106                 Hprev = H;
 107                 H = newton_interp(Cb, b, H, newton_poly, mainvar, p);
 108                 newton_poly = newton_poly*(mainvar - b);
 109
 110                 // try to reduce the number of division tests.
 111                 const ex H_lcoeff = lcoeff_wrt(H, restvars);
 112
 113                 if (H_lcoeff.is_equal(lc_gcd)) {
 114                         if ((Hprev-H).expand().smod(pn).is_zero())
 115                                 continue;
 116                         ex C /* primitive part of H */, contH /* dummy */;
 117                         primpart_content(C, contH, H, vars, p);
 118                         // Normalize GCD so that leading coefficient is 1
 119                         const cln::cl_I Clc = recip(integer_lcoeff(C, vars), p);
 120                         C = (C*numeric(Clc)).expand().smod(pn);
 121
 122                         ex dummy1, dummy2;
 123
 124                         if (divide_in_z_p(Aprim, C, dummy1, vars, p) &&
 125                                         divide_in_z_p(Bprim, C, dummy2, vars, p))
 126                                 return (cont_gcd*C).expand().smod(pn);
 127                         else if (img_gcd_deg == 0)
 128                                 return cont_gcd;
 129                         // else continue building the candidate
 130                 }
 131         } while(true);
 132         throw pgcd_failed();
 133 }
 134
 135 } // namespace GiNaC
 136