ginac/polynomial/pgcd.cpp

   1 /** @file pgcd.cpp
   2  *
   3  *  GCD for polynomials in prime field. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2010 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  21  */
  22
  23 #include "pgcd.h"
  24 #include "collect_vargs.h"
  25 #include "smod_helpers.h"
  26 #include "euclid_gcd_wrap.h"
  27 #include "eval_point_finder.h"
  28 #include "newton_interpolate.h"
  29 #include "divide_in_z_p.h"
  30
  31 namespace GiNaC {
  32
  33 extern void
  34 primpart_content(ex& pp, ex& c, ex e, const exvector& vars, const long p);
  35
  36 // Computes the GCD of two polynomials over a prime field.
  37 // Based on Algorithm 7.2 from "Algorithms for Computer Algebra"
  38 // A and B are considered as Z_p[x_n][x_0, \ldots, x_{n-1}], that is,
  39 // as a polynomials in variables x_0, \ldots x_{n-1} having coefficients
  40 // from the ring Z_p[x_n]
  41 ex pgcd(const ex& A, const ex& B, const exvector& vars, const long p)
  42 {
  43         static const ex ex1(1);
  44         if (A.is_zero())
  45                 return B;
  46
  47         if (B.is_zero())
  48                 return A;
  49
  50         if (is_a<numeric>(A))
  51                 return ex1;
  52
  53         if (is_a<numeric>(B))
  54                 return ex1;
  55
  56         // Checks for univariate polynomial
  57         if (vars.size() == 1) {
  58                 ex ret = euclid_gcd(A, B, vars[0], p); // Univariate GCD
  59                 return ret;
  60         }
  61         const ex& mainvar(vars.back());
  62
  63         // gcd of the contents
  64         ex H = 0, Hprev = 0; // GCD candidate
  65         ex newton_poly = 1;  // for Newton Interpolation
  66
  67         // Contents and primparts of A and B
  68         ex Aprim, contA;
  69         primpart_content(Aprim, contA, A, vars, p);
  70         ex Bprim, contB;
  71         primpart_content(Bprim, contB, B, vars, p);
  72         // gcd of univariate polynomials
  73         const ex cont_gcd = euclid_gcd(contA, contB, mainvar, p);
  74
  75         exvector restvars = vars;
  76         restvars.pop_back();
  77         const ex AL = lcoeff_wrt(Aprim, restvars);
  78         const ex BL = lcoeff_wrt(Bprim, restvars);
  79         // gcd of univariate polynomials
  80         const ex lc_gcd = euclid_gcd(AL, BL, mainvar, p);
  81
  82         // The estimate of degree of the gcd of Ab and Bb
  83         int gcd_deg = std::min(degree(Aprim, mainvar), degree(Bprim, mainvar));
  84         eval_point_finder::value_type b;
  85
  86         eval_point_finder find_eval_point(p);
  87         const numeric pn(p);
  88         do {
  89                 // Find a `good' evaluation point b.
  90                 bool has_more_pts = find_eval_point(b, lc_gcd, mainvar);
  91                 // If there are no more possible evaluation points, bail out
  92                 if (!has_more_pts)
  93                         throw pgcd_failed();
  94
  95                 const numeric bn(b);
  96                 // Evaluate the polynomials in b
  97                 ex Ab = Aprim.subs(mainvar == bn).smod(pn);
  98                 ex Bb = Bprim.subs(mainvar == bn).smod(pn);
  99                 ex Cb = pgcd(Ab, Bb, restvars, p);
 100
 101                 // Set the correct the leading coefficient
 102                 const cln::cl_I lcb_gcd =
 103                         smod(to_cl_I(lc_gcd.subs(mainvar == bn)), p);
 104                 const cln::cl_I Cblc = integer_lcoeff(Cb, restvars);
 105                 const cln::cl_I correct_lc = smod(lcb_gcd*recip(Cblc, p), p);
 106                 Cb = (Cb*numeric(correct_lc)).smod(pn);
 107
 108                 // Test for unlucky homomorphisms
 109                 const int img_gcd_deg = Cb.degree(restvars.back());
 110                 if (img_gcd_deg < gcd_deg) {
 111                         // The degree decreased, previous homomorphisms were
 112                         // bad, so we have to start it all over.
 113                         H = Cb;
 114                         newton_poly = mainvar - numeric(b);
 115                         Hprev = 0;
 116                         gcd_deg  = img_gcd_deg;
 117                         continue;
 118                 }
 119                 if (img_gcd_deg > gcd_deg) {
 120                         // The degree of images GCD is too high, this
 121                         // evaluation point is bad. Skip it.
 122                         continue;
 123                 }
 124                 if (img_gcd_deg == 0)
 125                         return cont_gcd;
 126
 127                 // Image has the same degree as the previous one
 128                 // (or at least not higher than the limit)
 129                 Hprev = H;
 130                 H = newton_interp(Cb, b, H, newton_poly, mainvar, p);
 131                 newton_poly = newton_poly*(mainvar - b);
 132
 133                 // try to reduce the number of division tests.
 134                 const ex H_lcoeff = lcoeff_wrt(H, restvars);
 135
 136                 if (H_lcoeff.is_equal(lc_gcd)) {
 137                         if ((Hprev-H).expand().smod(pn).is_zero())
 138                                 continue;
 139                         ex C /* primitive part of H */, contH /* dummy */;
 140                         primpart_content(C, contH, H, vars, p);
 141                         // Normalize GCD so that leading coefficient is 1
 142                         const cln::cl_I Clc = recip(integer_lcoeff(C, vars), p);
 143                         C = (C*numeric(Clc)).expand().smod(pn);
 144
 145                         ex dummy1, dummy2;
 146
 147                         if (divide_in_z_p(Aprim, C, dummy1, vars, p) &&
 148                                         divide_in_z_p(Bprim, C, dummy2, vars, p))
 149                                 return (cont_gcd*C).expand().smod(pn);
 150                         // else continue building the candidate
 151                 }
 152         } while(true);
 153         throw pgcd_failed();
 154 }
 155
 156 } // namespace GiNaC