86834001fcaa231de6d24306e3a746735f656abd
[ginac.git] / ginac / polynomial / mod_gcd.cpp
1 /** @file mod_gcd.cpp
2  *
3  *  Implementation of modular GCD. */
4
5 /*
6  *  GiNaC Copyright (C) 1999-2015 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 "upoly.h"
24 #include "gcd_euclid.h"
25 #include "cra_garner.h"
26 #include "debug.h"
27
28 #include <cln/numtheory.h>
29 #include <cln/random.h>
30 #include <stdexcept>
31
32 namespace GiNaC {
33
34 /**
35  * @brief Remove the integer content from univariate polynomials A and B.
36  *
37  * As a byproduct compute the GCD of contents.
38  */
39 static void remove_content(upoly& A, upoly& B, upoly::value_type& c)
40 {
41         // remove the integer
42         upoly::value_type acont, bcont;
43         normalize_in_ring(A, &acont);
44         normalize_in_ring(B, &bcont);
45         c = gcd(acont, bcont);
46 };
47
48 /// Check if @a candidate divides both @a A and @a B
49 static bool
50 do_division_check(const upoly& A, const upoly& B, const upoly& candidate)
51 {
52         upoly r1;
53         remainder_in_ring(r1, A, candidate);
54         if (r1.size() != 0)
55                 return false;
56
57         upoly r2;
58         remainder_in_ring(r2, B, candidate);
59         if (r2.size() != 0)
60                 return false;
61
62         return true;
63 }
64
65 /**
66  * Given two GCD candidates H \in Z/q[x], and C \in Z/p[x] (where p is a prime)
67  * compute the candidate in Z/(q*p) with CRA (chinise remainder algorithm)
68  *
69  * @param H \in Z/q[x] GCD candidate, will be updated by this function
70  * @param q modulus of H, will NOT be updated by this function
71  * @param C \in Z/p[x] GCD candidate, will be updated by this function
72  * @param p modulus of C
73  */
74 static void
75 update_the_candidate(upoly& H, const upoly::value_type& q,
76                      const umodpoly& C,
77                      const upoly::value_type& p,
78                      const cln::cl_modint_ring& R)
79 {
80         typedef upoly::value_type ring_t;
81         std::vector<ring_t> moduli(2);
82         moduli[0] = q;
83         moduli[1] = p;
84         if (H.size() < C.size())
85                 H.resize(C.size());
86
87         for (std::size_t  i = C.size(); i-- != 0; ) {
88                 std::vector<ring_t> coeffs(2);
89                 coeffs[0] = H[i];
90                 coeffs[1] = R->retract(C[i]); 
91                 H[i] = integer_cra(coeffs, moduli);
92         }
93 }
94
95 /// Convert Z/p[x] -> Z[x]
96 static void retract(upoly& p, const umodpoly& cp,
97                     const cln::cl_modint_ring& Rp)
98 {
99         p.resize(cp.size());
100         for (std::size_t i = cp.size(); i-- != 0; )
101                 p[i] = Rp->retract(cp[i]);
102 }
103
104
105 /// Find the prime which is > p, and does NOT divide g
106 static void find_next_prime(cln::cl_I& p, const cln::cl_I& g)
107 {
108         do {
109                 ++p;
110                 p = nextprobprime(p);
111         } while (zerop(mod(g, p)));
112 }
113
114 /// Compute the GCD of univariate polynomials A, B \in Z[x]
115 void mod_gcd(upoly& result, upoly A, upoly B)
116 {
117         typedef upoly::value_type ring_t;
118         ring_t content_gcd;
119         // remove the integer content
120         remove_content(A, B, content_gcd);
121
122         // compute the coefficient bound for GCD(a, b)
123         ring_t g = gcd(lcoeff(A), lcoeff(B));
124         std::size_t max_gcd_degree = std::min(degree(A), degree(B));
125         ring_t limit = (ring_t(1) << max_gcd_degree)*g*
126                        std::min(max_coeff(A), max_coeff(B));
127         ring_t q(0);
128         upoly H;
129
130         int count = 0;
131         const ring_t p_threshold = ring_t(1) << (8*sizeof(void *));
132         ring_t p = isqrt(std::min(max_coeff(A), max_coeff(B)));
133         while (true) {
134                 if (count >= 8) {
135                         count = 0;
136                         if (p < p_threshold)
137                                 p <<= 1;
138                         else
139                                 p = p + (p >> 4);
140                 } else 
141                         ++count;
142                 find_next_prime(p, g);
143
144                 // Map the polynomials onto Z/p[x]
145                 cln::cl_modint_ring Rp = cln::find_modint_ring(p);
146                 cln::cl_MI gp = Rp->canonhom(g);
147                 umodpoly ap(A.size()), bp(B.size());
148                 make_umodpoly(ap, A, Rp);
149                 make_umodpoly(bp, B, Rp);
150
151                 // Compute the GCD in Z/p[x]
152                 umodpoly cp;
153                 gcd_euclid(cp, ap, bp);
154                 bug_on(cp.size() == 0, "gcd(ap, bp) = 0, with ap = " <<
155                                         ap << ", and bp = " << bp);
156
157
158                 // Normalize the candidate so that its leading coefficient
159                 // is g mod p
160                 umodpoly::value_type norm_factor = gp*recip(lcoeff(cp));
161                 bug_on(zerop(norm_factor), "division in a field give 0");
162
163                 lcoeff(cp) = gp;
164                 for (std::size_t k = cp.size() - 1; k-- != 0; )
165                         cp[k] = cp[k]*norm_factor;
166
167
168                 // check for unlucky homomorphisms
169                 if (degree(cp) < max_gcd_degree) {
170                         q = p;
171                         max_gcd_degree = degree(cp);
172                         retract(H, cp, Rp);
173                 } else {
174                         update_the_candidate(H, q, cp, p, Rp);
175                         q = q*p;
176                 }
177
178                 if (q > limit) {
179                         result = H;
180                         normalize_in_ring(result);
181                         // if H divides both A and B it's a GCD
182                         if (do_division_check(A, B, result)) {
183                                 result *= content_gcd;
184                                 break; 
185                         }
186                         // H does not divide A and/or B, look for
187                         // another one
188                 } else if (degree(cp) == 0) {
189                         // Polynomials are relatively prime
190                         result.resize(1);
191                         result[0] = content_gcd;
192                         break;
193                 }
194         }
195 }
196
197 } // namespace GiNaC