3 * Implementation of modular GCD. */
6 * GiNaC Copyright (C) 1999-2018 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 "gcd_euclid.h"
25 #include "cra_garner.h"
28 #include <cln/numtheory.h>
29 #include <cln/random.h>
36 * @brief Remove the integer content from univariate polynomials A and B.
38 * As a byproduct compute the GCD of contents.
40 static void remove_content(upoly& A, upoly& B, upoly::value_type& c)
43 upoly::value_type acont, bcont;
44 normalize_in_ring(A, &acont);
45 normalize_in_ring(B, &bcont);
46 c = gcd(acont, bcont);
49 /// Check if @a candidate divides both @a A and @a B
51 do_division_check(const upoly& A, const upoly& B, const upoly& candidate)
54 remainder_in_ring(r1, A, candidate);
59 remainder_in_ring(r2, B, candidate);
67 * Given two GCD candidates H \in Z/q[x], and C \in Z/p[x] (where p is a prime)
68 * compute the candidate in Z/(q*p) with CRA (chinise remainder algorithm)
70 * @param H \in Z/q[x] GCD candidate, will be updated by this function
71 * @param q modulus of H, will NOT be updated by this function
72 * @param C \in Z/p[x] GCD candidate, will be updated by this function
73 * @param p modulus of C
76 update_the_candidate(upoly& H, const upoly::value_type& q,
78 const upoly::value_type& p,
79 const cln::cl_modint_ring& R)
81 typedef upoly::value_type ring_t;
82 std::vector<ring_t> moduli(2);
85 if (H.size() < C.size())
88 for (std::size_t i = C.size(); i-- != 0; ) {
89 std::vector<ring_t> coeffs(2);
91 coeffs[1] = R->retract(C[i]);
92 H[i] = integer_cra(coeffs, moduli);
96 /// Convert Z/p[x] -> Z[x]
97 static void retract(upoly& p, const umodpoly& cp,
98 const cln::cl_modint_ring& Rp)
101 for (std::size_t i = cp.size(); i-- != 0; )
102 p[i] = Rp->retract(cp[i]);
106 /// Find the prime which is > p, and does NOT divide g
107 static void find_next_prime(cln::cl_I& p, const cln::cl_I& g)
111 p = nextprobprime(p);
112 } while (zerop(mod(g, p)));
115 /// Compute the GCD of univariate polynomials A, B \in Z[x]
116 void mod_gcd(upoly& result, upoly A, upoly B)
118 typedef upoly::value_type ring_t;
120 // remove the integer content
121 remove_content(A, B, content_gcd);
123 // compute the coefficient bound for GCD(a, b)
124 ring_t g = gcd(lcoeff(A), lcoeff(B));
125 std::size_t max_gcd_degree = std::min(degree(A), degree(B));
126 ring_t limit = (ring_t(1) << max_gcd_degree)*g*
127 std::min(max_coeff(A), max_coeff(B));
132 const ring_t p_threshold = ring_t(1) << (8*sizeof(void *));
133 ring_t p = isqrt(std::min(max_coeff(A), max_coeff(B)));
143 find_next_prime(p, g);
145 // Map the polynomials onto Z/p[x]
146 cln::cl_modint_ring Rp = cln::find_modint_ring(p);
147 cln::cl_MI gp = Rp->canonhom(g);
148 umodpoly ap(A.size()), bp(B.size());
149 make_umodpoly(ap, A, Rp);
150 make_umodpoly(bp, B, Rp);
152 // Compute the GCD in Z/p[x]
154 gcd_euclid(cp, ap, bp);
155 bug_on(cp.size() == 0, "gcd(ap, bp) = 0, with ap = " <<
156 ap << ", and bp = " << bp);
159 // Normalize the candidate so that its leading coefficient
161 umodpoly::value_type norm_factor = gp*recip(lcoeff(cp));
162 bug_on(zerop(norm_factor), "division in a field give 0");
165 for (std::size_t k = cp.size() - 1; k-- != 0; )
166 cp[k] = cp[k]*norm_factor;
169 // check for unlucky homomorphisms
170 if (degree(cp) < max_gcd_degree) {
172 max_gcd_degree = degree(cp);
175 update_the_candidate(H, q, cp, p, Rp);
181 normalize_in_ring(result);
182 // if H divides both A and B it's a GCD
183 if (do_division_check(A, B, result)) {
184 result *= content_gcd;
187 // H does not divide A and/or B, look for
189 } else if (degree(cp) == 0) {
190 // Polynomials are relatively prime
192 result[0] = content_gcd;