3 * GCD for polynomials in prime field. */
6 * GiNaC Copyright (C) 1999-2009 Johannes Gutenberg University Mainz, Germany
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.
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.
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
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"
34 primpart_content(ex& pp, ex& c, ex e, const exvector& vars, const long p);
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)
43 static const ex ex1(1);
56 // Checks for univariate polynomial
57 if (vars.size() == 1) {
58 ex ret = euclid_gcd(A, B, vars[0], p); // Univariate GCD
61 const ex& mainvar(vars.back());
63 // gcd of the contents
64 ex H = 0, Hprev = 0; // GCD candidate
65 ex newton_poly = 1; // for Newton Interpolation
67 // Contents and primparts of A and B
69 primpart_content(Aprim, contA, A, vars, p);
71 primpart_content(Bprim, contB, B, vars, p);
72 // gcd of univariate polynomials
73 const ex cont_gcd = euclid_gcd(contA, contB, mainvar, p);
75 exvector restvars = vars;
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);
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;
86 eval_point_finder find_eval_point(p);
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
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);
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);
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.
114 newton_poly = mainvar - numeric(b);
116 gcd_deg = img_gcd_deg;
119 if (img_gcd_deg > gcd_deg) {
120 // The degree of images GCD is too high, this
121 // evaluation point is bad. Skip it.
124 if (img_gcd_deg == 0)
127 // Image has the same degree as the previous one
128 // (or at least not higher than the limit)
130 H = newton_interp(Cb, b, H, newton_poly, mainvar, p);
131 newton_poly = newton_poly*(mainvar - b);
133 // try to reduce the number of division tests.
134 const ex H_lcoeff = lcoeff_wrt(H, restvars);
136 if (H_lcoeff.is_equal(lc_gcd)) {
137 if ((Hprev-H).expand().smod(pn).is_zero())
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);
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