1 /*++
2 Copyright (c) 2011 Microsoft Corporation
3
4 Module Name:
5
6 polynomial_factorization.cpp
7
8 Abstract:
9
10 Implementation of polynomial factorization.
11
12 Author:
13
14 Dejan (t-dejanj) 2011-11-15
15
16 Notes:
17
18 [1] Elwyn Ralph Berlekamp. Factoring Polynomials over Finite Fields. Bell System Technical Journal,
19 46(8-10):1853-1859, 1967.
20 [2] Donald Ervin Knuth. The Art of Computer Programming, volume 2: Seminumerical Algorithms. Addison Wesley, third
21 edition, 1997.
22 [3] Henri Cohen. A Course in Computational Algebraic Number Theory. Springer Verlag, 1993.
23
24 --*/
25 #include "util/trace.h"
26 #include "util/util.h"
27 #include "math/polynomial/upolynomial_factorization_int.h"
28 #include "util/prime_generator.h"
29
30 using namespace std;
31
32 namespace upolynomial {
33
34 // get the prime as unsigned while throwing exceptions if anything goes bad`
get_p_from_manager(zp_numeral_manager const & zp_nm)35 unsigned get_p_from_manager(zp_numeral_manager const & zp_nm) {
36 z_numeral_manager & nm = zp_nm.m();
37 numeral const & p = zp_nm.p();
38 if (!nm.is_uint64(p)) {
39 throw upolynomial_exception("The prime number attempted in factorization is too big!");
40 }
41 uint64_t p_uint64 = nm.get_uint64(p);
42 unsigned p_uint = static_cast<unsigned>(p_uint64);
43 if (((uint64_t)p_uint) != p_uint64) {
44 throw upolynomial_exception("The prime number attempted in factorization is too big!");
45 }
46 return p_uint;
47 }
48
49 /**
50 \brief The Q-I matrix from Berelkamp's algorithm [1,2].
51
52 Given a polynomial f = \sum f_k x^k of degree n, with f_i in Z_p[x], the i-th row of Q is a representation of
53 x^(p*i) modulo f, i.e.
54
55 x^(p*i) modulo f = \sum_j q[i][j] x^j
56
57 If f is of degree n, the matrix is square nxn. When the this matrix is constructed we obtain Q - I, because
58 this is what we need in the algorithm. After construction, the null space vectors can be extracted one-by one using
59 the next_null_space_vector method.
60 */
61 class berlekamp_matrix {
62
63 zp_manager & m_upm;
64 mpzzp_manager & m_zpm;
65
66 svector<mpz> m_matrix;
67 unsigned m_size;
68
69 unsigned m_null_row; // 0, ..., m_size - 1, state for null vectors
70 svector<int> m_column_pivot; // position of pivots in the columns
71 svector<int> m_row_pivot; // position of pivots in the rows
72
get(unsigned i,unsigned j)73 mpz & get(unsigned i, unsigned j) {
74 SASSERT(i < m_size && j < m_size);
75 return m_matrix[i*m_size + j];
76 }
77
get(unsigned i,unsigned j) const78 mpz const & get(unsigned i, unsigned j) const {
79 SASSERT(i < m_size && j < m_size);
80 return m_matrix[i*m_size + j];
81 }
82
83 public:
84
85 /**
86 \brief Construct the matrix as explained above, f should be in the Z_p.
87 */
berlekamp_matrix(zp_manager & upm,numeral_vector const & f)88 berlekamp_matrix(zp_manager & upm, numeral_vector const & f)
89 : m_upm(upm),
90 m_zpm(upm.m()),
91 m_size(m_upm.degree(f)),
92 m_null_row(1),
93 m_column_pivot(m_size, -1),
94 m_row_pivot(m_size, -1) {
95 unsigned p = get_p_from_manager(m_zpm);
96
97 TRACE("polynomial::factorization::bughunt", tout << "polynomial::berlekamp_matrix("; m_upm.display(tout, f); tout << ", " << p << ")" << endl;);
98
99 // the first row is always the vector [1, 0, ..., 0], since x^0 = 0 (modulo f)
100 m_matrix.push_back(1);
101 for (unsigned j = 0; j < m_size; ++ j) {
102 m_matrix.push_back(0);
103 }
104
105 // the rest of the rows, we can get as follows, given f = x^n + f_{n-1} x^{n-1} + ... + f_1 x + f_0
106 // if x^k = \sum a_{k,j} x^j (modulo p), hence 0 <= j <= n-1 then
107 // x^{k+1} = a_{k,n-1}(-f_{n-1} x^{n-1} + ... + f_0) + a_{k,n-2} x^{n-1} + ... + a_{k, 0} x
108 // so we can compute a_{k+1,j} = a_{k, j-1} - a_{k,n-1}*f_j
109 // every p-th row we add to the matrix
110 scoped_numeral tmp(m_zpm);
111 unsigned row = 0, previous_row = 0;
112 for (unsigned k = 1; true; previous_row = row, ++ k) {
113
114 // add another row if we need it
115 if (k % p == 1) {
116 if (++ row >= m_size) {
117 break;
118 }
119 for (unsigned j = 0; j < m_size; ++ j) {
120 m_matrix.push_back(0);
121 }
122 }
123
124 // the multiplier
125 m_zpm.set(tmp, get(previous_row, m_size - 1));
126
127 // go down the row and shift it
128 for (unsigned j = m_size - 1; j > 0; -- j) {
129 m_zpm.submul(get(previous_row, j-1), tmp, f[j], get(row, j));
130 }
131
132 // add the 0 element (same formula with a_{k,-1} = 0)
133 m_zpm.mul(f[0], tmp, get(row, 0));
134 m_zpm.neg(get(row, 0));
135 }
136
137 // do Q - I
138 for (unsigned i = 0; i < m_size; ++ i) {
139 m_zpm.dec(get(i, i));
140 }
141
142 TRACE("polynomial::factorization::bughunt", tout << "polynomial::berlekamp_matrix():" << endl; display(tout); tout << endl;);
143 }
144
145 /**
146 \brief Destructor, just removing the numbers
147 */
~berlekamp_matrix()148 ~berlekamp_matrix() {
149 for (unsigned k = 0; k < m_matrix.size(); ++ k) {
150 m_zpm.del(m_matrix[k]);
151 }
152 }
153
154 /**
155 \brief 'Diagonalizes' the matrix using only column operations. The resulting matrix will have -1 at pivot
156 elements. Returns the rank of the null space.
157 */
diagonalize()158 unsigned diagonalize() {
159
160 scoped_numeral multiplier(m_zpm);
161
162 unsigned null_rank = 0;
163 for (unsigned i = 0; i < m_size; ++ i) {
164
165 // get the first non-zero entry in the m_null_row row
166 bool column_found = false;
167 for (unsigned j = 0; j < m_size; ++ j) {
168 if (m_column_pivot[j] < 0 && !m_zpm.is_zero(get(i, j))) {
169 column_found = true;
170 m_column_pivot[j] = i;
171 m_row_pivot[i] = j;
172
173 // found a pivot, to make it -1 we compute the multiplier -p^-1
174 m_zpm.set(multiplier, get(i, j));
175 m_zpm.inv(multiplier);
176 m_zpm.neg(multiplier);
177
178 // multiply the pivot column with the multiplier
179 for (unsigned k = m_null_row; k < m_size; ++ k) {
180 m_zpm.mul(get(k, j), multiplier, get(k, j));
181 }
182 // pivot is -1 so we can add it to the rest of the columns to eliminate the row
183 for (unsigned other_j = 0; other_j < m_size; ++ other_j) {
184 if (j != other_j) {
185 m_zpm.set(multiplier, get(i, other_j));
186 for (unsigned k = m_null_row; k < m_size; ++ k) {
187 m_zpm.addmul(get(k, other_j), multiplier, get(k, j), get(k, other_j));
188 }
189 }
190 }
191 }
192 }
193 if (!column_found) {
194 null_rank ++;
195 }
196 }
197
198 TRACE("polynomial::factorization::bughunt", tout << "polynomial::diagonalize():" << endl; display(tout); tout << endl;);
199
200 return null_rank;
201 }
202
203 /**
204 If rank of the matrix is n - r, we are interested in linearly independent vectors v_1, ..., v_r (the basis of
205 the null space), such that v_k A = 0. This method will give one at a time. The method returns true if vector has
206 been computed properly. The first vector [1, 0, ..., 0] is ignored (m_null_row starts from 1).
207 */
next_null_space_vector(numeral_vector & v)208 bool next_null_space_vector(numeral_vector & v) {
209 SASSERT(v.size() <= m_size);
210 v.resize(m_size);
211 for (; m_null_row < m_size; ++ m_null_row) {
212 if (m_row_pivot[m_null_row] < 0) {
213 // output the vector
214 for (unsigned j = 0; j < m_size; ++ j) {
215 if (m_row_pivot[j] >= 0) {
216 m_zpm.set(v[j], get(m_null_row, m_row_pivot[j]));
217 }
218 else {
219 if (j == m_null_row) {
220 m_zpm.set(v[j], 1);
221 }
222 else {
223 m_zpm.set(v[j], 0);
224 }
225 }
226 }
227 ++ m_null_row;
228 return true;
229 }
230 }
231 // didn't find the vector
232 return false;
233 }
234
235 /**
236 \brief Display the matrix on the output stream.
237 */
display(std::ostream & out) const238 void display(std::ostream & out) const {
239 for (unsigned i = 0; i < m_matrix.size() / m_size; ++ i) {
240 for (unsigned j = 0; j < m_size; ++ j) {
241 out << m_zpm.to_string(get(i, j)) << "\t";
242 }
243 out << endl;
244 }
245 }
246 };
247
248 /**
249 See [3] p. 125.
250 */
zp_square_free_factor(zp_manager & upm,numeral_vector const & f,zp_factors & sq_free_factors)251 void zp_square_free_factor(zp_manager & upm, numeral_vector const & f, zp_factors & sq_free_factors) {
252
253 zp_numeral_manager & nm = upm.m();
254 unsigned p = get_p_from_manager(upm.m());
255
256 TRACE("polynomial::factorization", tout << "polynomial::square_free_factor("; upm.display(tout, f); tout << ") over Z_" << p << endl;);
257
258 scoped_numeral_vector div_tmp(nm);
259
260 // [initialize] T_0 = f, e = 1
261 // trim and get the make it monic if not already
262 SASSERT(f.size() > 1);
263 scoped_numeral_vector T_0(nm);
264 upm.set(f.size(), f.c_ptr(), T_0);
265 scoped_numeral constant(nm);
266 upm.mk_monic(T_0.size(), T_0.c_ptr(), constant);
267 sq_free_factors.set_constant(constant);
268 TRACE("polynomial::factorization::bughunt",
269 tout << "Initial factors: " << sq_free_factors << endl;
270 tout << "R.<x> = GF(" << p << ")['x']" << endl;
271 tout << "T_0 = "; upm.display(tout, T_0); tout << endl;
272 );
273 unsigned e = 1;
274
275 // we repeat until we get a constant
276 scoped_numeral_vector T_0_d(nm);
277 scoped_numeral_vector T(nm);
278 scoped_numeral_vector V(nm);
279 scoped_numeral_vector W(nm);
280 scoped_numeral_vector A_ek(nm);
281 while (T_0.size() > 1)
282 {
283 // [initialize e-loop] T = gcd(T_0, T_0'), V / T_0/T, k = 0
284 unsigned k = 0;
285 TRACE("polynomial::factorization::bughunt", tout << "k = 0" << endl;);
286
287 // T_0_d = T_0'
288 upm.derivative(T_0.size(), T_0.c_ptr(), T_0_d);
289 TRACE("polynomial::factorization::bughunt",
290 tout << "T_0_d = T_0.derivative(x)" << endl;
291 tout << "T_0_d == "; upm.display(tout, T_0_d); tout << endl;
292 );
293
294 // T = gcd(T_0, T_0')
295 upm.gcd(T_0.size(), T_0.c_ptr(), T_0_d.size(), T_0_d.c_ptr(), T);
296 TRACE("polynomial::factorization::bughunt",
297 tout << "T = T_0.gcd(T_0_d)" << endl;
298 tout << "T == "; upm.display(tout, T); tout << endl;
299 );
300
301 // V = T_0 / T
302 upm.div(T_0.size(), T_0.c_ptr(), T.size(), T.c_ptr(), V);
303 TRACE("polynomial::factorization::bughunt",
304 tout << "V = T_0.quo_rem(T)[0]" << endl;
305 tout << "V == "; upm.display(tout, V); tout << endl;
306 );
307
308 while (V.size() > 1) {
309 // [special case]
310 if ((++k) % p == 0) {
311 ++ k;
312 // T = T/V
313 upm.div(T.size(), T.c_ptr(), V.size(), V.c_ptr(), T);
314 TRACE("polynomial::factorization::bughunt",
315 tout << "T = T.quo_rem(V)[0]" << endl;
316 tout << "T == "; upm.display(tout, T); tout << endl;
317 );
318 }
319
320 // [compute A_ek]
321
322 // W = gcd(T, V)
323 upm.gcd(T.size(), T.c_ptr(), V.size(), V.c_ptr(), W);
324 TRACE("polynomial::factorization::bughunt",
325 tout << "W = T.gcd(V)" << endl;
326 upm.display(tout, W); tout << endl;
327 );
328
329 // A_ek = V/W
330 upm.div(V.size(), V.c_ptr(), W.size(), W.c_ptr(), A_ek);
331 TRACE("polynomial::factorization::bughunt",
332 tout << "A_ek = V.quo_rem(W)[0]" << endl;
333 tout << "A_ek == "; upm.display(tout, A_ek); tout << endl;
334 );
335
336 // V = W
337 V.swap(W);
338 TRACE("polynomial::factorization::bughunt",
339 tout << "V = W" << endl;
340 tout << "V == "; upm.display(tout, V); tout << endl;
341 );
342
343 // T = T/V
344 upm.div(T.size(), T.c_ptr(), V.size(), V.c_ptr(), T);
345 TRACE("polynomial::factorization::bughunt",
346 tout << "T = T.quo_rem(V)[0]" << endl;
347 tout << "T == "; upm.display(tout, T); tout << endl;
348 );
349
350 // if not constant, we output it
351 if (A_ek.size() > 1) {
352 TRACE("polynomial::factorization::bughunt", tout << "Factor: ("; upm.display(tout, A_ek); tout << ")^" << e*k << endl;);
353 sq_free_factors.push_back(A_ek, e*k);
354 }
355 }
356
357 // [finished e-loop] T_0 = \sum_{p div j} t_j x^j, set T_0 \sum_{p div j} t_j x^{j/p}, e = pe
358 e *= p;
359 T_0.reset();
360 for (unsigned deg_p = 0; deg_p < T.size(); deg_p += p) {
361 T_0.push_back(numeral());
362 nm.set(T_0.back(), T[deg_p]);
363 }
364 TRACE("polynomial::factorization::bughunt",
365 tout << "T_0 = "; upm.display(tout, T_0); tout << endl;
366 );
367 }
368
369 TRACE("polynomial::factorization", tout << "polynomial::square_free_factor("; upm.display(tout, f); tout << ") => " << sq_free_factors << endl;);
370 }
371
zp_factor(zp_manager & upm,numeral_vector const & f,zp_factors & factors)372 bool zp_factor(zp_manager & upm, numeral_vector const & f, zp_factors & factors) {
373
374
375 TRACE("polynomial::factorization",
376 unsigned p = get_p_from_manager(upm.m());
377 tout << "polynomial::factor("; upm.display(tout, f); tout << ") over Z_" << p << endl;);
378
379 // get the sq-free parts (all of them will be monic)
380 zp_factors sq_free_factors(upm);
381 zp_square_free_factor(upm, f, sq_free_factors);
382
383 // factor the sq-free parts individually
384 for (unsigned i = 0; i < sq_free_factors.distinct_factors(); ++ i) {
385 unsigned j = factors.distinct_factors();
386 if (upm.degree(sq_free_factors[i]) > 1) {
387 zp_factor_square_free(upm, sq_free_factors[i], factors); // monic from aq-free decomposition
388 for (; j < factors.distinct_factors(); ++ j) {
389 factors.set_degree(j, sq_free_factors.get_degree(i)*factors.get_degree(j));
390 }
391 }
392 else {
393 factors.push_back(sq_free_factors[i], sq_free_factors.get_degree(i));
394 }
395 }
396 // add the constant
397 factors.set_constant(sq_free_factors.get_constant());
398
399 TRACE("polynomial::factorization", tout << "polynomial::factor("; upm.display(tout, f); tout << ") => " << factors << endl;);
400
401 return factors.total_factors() > 1;
402 }
403
zp_factor_square_free(zp_manager & upm,numeral_vector const & f,zp_factors & factors)404 bool zp_factor_square_free(zp_manager & upm, numeral_vector const & f, zp_factors & factors) {
405 return zp_factor_square_free_berlekamp(upm, f, factors, false);
406 }
407
zp_factor_square_free_berlekamp(zp_manager & upm,numeral_vector const & f,zp_factors & factors,bool randomized)408 bool zp_factor_square_free_berlekamp(zp_manager & upm, numeral_vector const & f, zp_factors & factors, bool randomized) {
409 SASSERT(upm.degree(f) > 1);
410
411 mpzzp_manager & zpm = upm.m();
412 unsigned p = get_p_from_manager(zpm);
413
414 TRACE("polynomial::factorization", tout << "upolynomial::factor_square_free_berlekamp("; upm.display(tout, f); tout << ", " << p << ")" << endl;);
415 SASSERT(zpm.is_one(f.back()));
416
417 // construct the berlekamp Q matrix to get the null space
418 berlekamp_matrix Q_I(upm, f);
419
420 // copy the initial polynomial to factors
421 unsigned first_factor = factors.distinct_factors();
422 factors.push_back(f, 1);
423
424 // rank of the null-space (and the number of factors)
425 unsigned r = Q_I.diagonalize();
426 if (r == 1) {
427 // since r == 1 == number of factors, then f is irreducible
428 TRACE("polynomial::factorization", tout << "upolynomial::factor_square_free_berlekamp("; upm.display(tout, f); tout << ", " << p << ") => " << factors << endl;);
429 return false;
430 }
431
432 TRACE("polynomial::factorization::bughunt", tout << "upolynomial::factor_square_free_berlekamp(): computing factors, expecting " << r << endl;);
433
434 scoped_numeral_vector gcd(zpm);
435 scoped_numeral_vector div(zpm);
436
437 // process the null space vectors (skip first one, it's [1, 0, ..., 0]) while generating the factors
438 unsigned d = upm.degree(f);
439 (void)d;
440 scoped_numeral_vector v_k(zpm);
441 while (Q_I.next_null_space_vector(v_k)) {
442
443 TRACE("polynomial::factorization::bughunt",
444 tout << "null vector: ";
445 for(unsigned j = 0; j < d; ++ j) {
446 tout << zpm.to_string(v_k[j]) << " ";
447 }
448 tout << endl;
449 );
450
451 upm.trim(v_k);
452 // TRACE("polynomial::factorization", tout << "v_k = "; upm.display(tout, v_k); tout << endl;);
453
454 unsigned current_factor_end = factors.distinct_factors();
455 for (unsigned current_factor_i = first_factor; current_factor_i < current_factor_end; ++ current_factor_i) {
456
457 // we have v such that vQ = v, viewing v as a polynomial, we get that v^n - v = 0 (mod f)
458 // since v^n -v = v*(v-1)*...*(v - p-1) we compute the gcd(v - s, f) to extract the
459 // factors. it also holds that g = \prod gcd(v - s, f), so we just accumulate them
460
461 // if it's of degree 1, we're done (have to index the array as we are adding to it), as
462 if (factors[current_factor_i].size() == 2) {
463 continue;
464 }
465
466 for (unsigned s = 0; s < p; ++ s) {
467
468 numeral_vector const & current_factor = factors[current_factor_i];
469
470 // we just take one off v_k each time to get all of them
471 zpm.dec(v_k[0]);
472
473 // get the gcd
474 upm.gcd(v_k.size(), v_k.c_ptr(), current_factor.size(), current_factor.c_ptr(), gcd);
475
476 // if the gcd is 1, or the gcd is f, we just ignore it
477 if (gcd.size() != 1 && gcd.size() != current_factor.size()) {
478
479 // get the divisor also (no need to normalize the div, both are monic)
480 upm.div(current_factor.size(), current_factor.c_ptr(), gcd.size(), gcd.c_ptr(), div);
481
482 TRACE("polynomial::factorization::bughunt",
483 tout << "current_factor = "; upm.display(tout, current_factor); tout << endl;
484 tout << "gcd_norm = "; upm.display(tout, gcd); tout << endl;
485 tout << "div = "; upm.display(tout, div); tout << endl;
486 );
487
488 // continue with the rest
489 factors.swap_factor(current_factor_i, div);
490
491 // add the new factor(s)
492 factors.push_back(gcd, 1);
493
494 }
495
496 // at the point where we have all the factors, we are done
497 if (factors.distinct_factors() - first_factor == r) {
498 TRACE("polynomial::factorization", tout << "polynomial::factor("; upm.display(tout, f); tout << ", " << p << ") => " << factors << " of degree " << factors.get_degree() << endl;);
499 return true;
500 }
501 }
502 }
503 }
504
505 // should never get here
506 SASSERT(false);
507 return true;
508 }
509
510 /**
511 Check if the hensel lifting was correct, i.e. that C = A*B (mod br).
512 */
check_hansel_lift(z_manager & upm,numeral_vector const & C,numeral const & a,numeral const & b,numeral const & r,numeral_vector const & A,numeral_vector const & B,numeral_vector const & A_lifted,numeral_vector const & B_lifted)513 bool check_hansel_lift(z_manager & upm, numeral_vector const & C,
514 numeral const & a, numeral const & b, numeral const & r,
515 numeral_vector const & A, numeral_vector const & B,
516 numeral_vector const & A_lifted, numeral_vector const & B_lifted)
517 {
518 z_numeral_manager & nm = upm.zm();
519
520 scoped_mpz br(nm);
521 nm.mul(b, r, br);
522
523 zp_manager br_upm(upm.lim(), upm.zm());
524 br_upm.set_zp(br);
525
526 if (A_lifted.size() != A.size()) return false;
527 if (B_lifted.size() != B.size()) return false;
528 if (!nm.eq(A.back(), A_lifted.back())) return false;
529
530 // test1: check that C = A_lifted * B_lifted (mod b*r)
531 scoped_mpz_vector test1(nm);
532 upm.mul(A_lifted.size(), A_lifted.c_ptr(), B_lifted.size(), B_lifted.c_ptr(), test1);
533 upm.sub(C.size(), C.c_ptr(), test1.size(), test1.c_ptr(), test1);
534 to_zp_manager(br_upm, test1);
535 if (!test1.empty()) {
536 TRACE("polynomial::factorization::bughunt",
537 tout << "sage: R.<x> = ZZ['x']" << endl;
538 tout << "sage: A = "; upm.display(tout, A); tout << endl;
539 tout << "sage: B = "; upm.display(tout, B); tout << endl;
540 tout << "sage: C = "; upm.display(tout, C); tout << endl;
541 tout << "sage: test1 = C - AB" << endl;
542 tout << "sage: print test1.change_ring(GF(" << nm.to_string(br) << "))" << endl;
543 tout << "sage: print 'expected 0'" << endl;
544 );
545 return false;
546 }
547
548 zp_manager b_upm(upm.lim(), nm);
549 b_upm.set_zp(b);
550
551 // test2: A_lifted = A (mod b)
552 scoped_mpz_vector test2a(nm), test2b(nm);
553 to_zp_manager(b_upm, A, test2a);
554 to_zp_manager(b_upm, A_lifted, test2b);
555 if (!upm.eq(test2a, test2b)) {
556 return false;
557 }
558
559 // test3: B_lifted = B (mod b)
560 scoped_mpz_vector test3a(nm), test3b(nm);
561 to_zp_manager(b_upm, B, test3a);
562 to_zp_manager(b_upm, B_lifted, test3b);
563 if (!upm.eq(test3a, test3b)) {
564 return false;
565 }
566
567 return true;
568 }
569
570 /**
571 Performs a Hensel lift of A and B in Z_a to Z_b, where p is prime and a = p^{a_k}, b = p^{b_k},
572 r = (a, b), with the following assumptions:
573
574 (1) UA + VB = 1 (mod a)
575 (2) C = A*B (mod b)
576 (3) (l(A), r) = 1 (important in order to divide by A, i.e. to invert l(A))
577 (4) deg(A) + deg(B) = deg(C)
578
579 The output of is two polynomials A1, B1 such that A1 = A (mod b), B1 = B (mod b),
580 l(A1) = l(A), deg(A1) = deg(A), deg(B1) = deg(B) and C = A1 B1 (mod b*r). Such A1, B1 are unique if
581 r is prime. See [3] p. 138.
582 */
hensel_lift(z_manager & upm,numeral const & a,numeral const & b,numeral const & r,numeral_vector const & U,numeral_vector const & A,numeral_vector const & V,numeral_vector const & B,numeral_vector const & C,numeral_vector & A_lifted,numeral_vector & B_lifted)583 void hensel_lift(z_manager & upm, numeral const & a, numeral const & b, numeral const & r,
584 numeral_vector const & U, numeral_vector const & A, numeral_vector const & V, numeral_vector const & B,
585 numeral_vector const & C, numeral_vector & A_lifted, numeral_vector & B_lifted) {
586
587 z_numeral_manager & nm = upm.zm();
588
589 TRACE("polynomial::factorization::bughunt",
590 tout << "polynomial::hensel_lift(";
591 tout << "a = " << nm.to_string(a) << ", ";
592 tout << "b = " << nm.to_string(b) << ", ";
593 tout << "r = " << nm.to_string(r) << ", ";
594 tout << "U = "; upm.display(tout, U); tout << ", ";
595 tout << "A = "; upm.display(tout, A); tout << ", ";
596 tout << "V = "; upm.display(tout, V); tout << ", ";
597 tout << "B = "; upm.display(tout, B); tout << ", ";
598 tout << "C = "; upm.display(tout, C); tout << ")" << endl;
599 );
600
601 zp_manager r_upm(upm.lim(), nm);
602 r_upm.set_zp(r);
603
604 SASSERT(upm.degree(C) == upm.degree(A) + upm.degree(B));
605 SASSERT(upm.degree(U) < upm.degree(B) && upm.degree(V) < upm.degree(A));
606
607 // by (2) C = AB (mod b), hence (C - AB) is divisible by b
608 // define thus let f = (C - AB)/b in Z_r
609 scoped_numeral_vector f(upm.m());
610 upm.mul(A.size(), A.c_ptr(), B.size(), B.c_ptr(), f);
611 upm.sub(C.size(), C.c_ptr(), f.size(), f.c_ptr(), f);
612 upm.div(f, b);
613 to_zp_manager(r_upm, f);
614 TRACE("polynomial::factorization",
615 tout << "f = "; upm.display(tout, f); tout << endl;
616 );
617
618 // we need to get A1 = A (mod b), B1 = B (mode b) so we know that we need
619 // A1 = A + b*S, B1 = B + b*T in Z[x] for some S and T with deg(S) <= deg(A), deg(T) <= deg(B)
620 // we also need (mod b*r) C = A1*B1 = (A + b*S)*(B + b*T) = AB + b(AT + BS) + b^2*ST
621 // if we find S and T, we will have found our A1 and B1
622 // since r divides b, then b^2 contains b*r and we know that it must be that C = AB + b(AT + BS) (mod b*r)
623 // which is equivalent to
624 // (5) f = (C - AB)/b = AT + BS (mod r)
625 // having (1) AU + BV = 1 (mod r) and (5) AT + BS = f (mod r), we know that
626 // A*(fU) + B*(fV) = f (mod r), i.e. T = fU, S = fV is a solution
627 // but we also know that we need an S with deg(S) <= deg(A) so we can do the following
628 // we know that l(A) is invertible so we can find the exact remainder of fV with A, i.e. find the quotient
629 // t in the division and set
630 // A*(fU + tB) + B*(fV - tA) = f
631 // T = fU + tB, S = fU - tA
632 // since l(A) is invertible in Z_r, we can (in Z_r) use exact division to get Vf = At + R with deg(R) < A
633 // we now know that deg(A+bS) = deg(A), but we also know (4) which will guarantee that deg(B+bT) = deg(B)
634
635 // compute the S, T (compute in Z_r[x])
636 scoped_numeral_vector Vf(r_upm.m()), t(r_upm.m()), S(r_upm.m());
637 TRACE("polynomial::factorization::bughunt",
638 tout << "V == "; upm.display(tout, V); tout << endl;
639 );
640 r_upm.mul(V.size(), V.c_ptr(), f.size(), f.c_ptr(), Vf);
641 TRACE("polynomial::factorization::bughunt",
642 tout << "Vf = V*f" << endl;
643 tout << "Vf == "; upm.display(tout, Vf); tout << endl;
644 );
645 r_upm.div_rem(Vf.size(), Vf.c_ptr(), A.size(), A.c_ptr(), t, S);
646 TRACE("polynomial::factorization::bughunt",
647 tout << "[t, S] = Vf.quo_rem(A)" << endl;
648 tout << "t == "; upm.display(tout, t); tout << endl;
649 tout << "S == "; upm.display(tout, S); tout << endl;
650 );
651 scoped_numeral_vector T(r_upm.m()), tmp(r_upm.m());
652 r_upm.mul(U.size(), U.c_ptr(), f.size(), f.c_ptr(), T); // T = fU
653 TRACE("polynomial::factorization::bughunt",
654 tout << "T == U*f" << endl;
655 tout << "T == "; upm.display(tout, T); tout << endl;
656 );
657 r_upm.mul(B.size(), B.c_ptr(), t.size(), t.c_ptr(), tmp); // tmp = Bt
658 TRACE("polynomial::factorization::bughunt",
659 tout << "tmp = B*t" << endl;
660 tout << "tmp == "; upm.display(tout, tmp); tout << endl;
661 );
662 r_upm.add(T.size(), T.c_ptr(), tmp.size(), tmp.c_ptr(), T); // T = Uf + Bt
663 TRACE("polynomial::factorization::bughunt",
664 tout << "T = B*tmp" << endl;
665 tout << "T == "; upm.display(tout, T); tout << endl;
666 );
667
668 // set the result, A1 = A + b*S, B1 = B + b*T (now we compute in Z[x])
669 upm.mul(S, b);
670 upm.mul(T, b);
671 upm.add(A.size(), A.c_ptr(), S.size(), S.c_ptr(), A_lifted);
672 upm.add(B.size(), B.c_ptr(), T.size(), T.c_ptr(), B_lifted);
673
674 CASSERT("polynomial::factorizatio::bughunt", check_hansel_lift(upm, C, a, b, r, A, B, A_lifted, B_lifted));
675 }
676
check_quadratic_hensel(zp_manager & zpe_upm,numeral_vector const & U,numeral_vector const & A,numeral_vector const & V,numeral_vector const & B)677 bool check_quadratic_hensel(zp_manager & zpe_upm, numeral_vector const & U, numeral_vector const & A, numeral_vector const & V, numeral_vector const & B) {
678 z_numeral_manager & nm = zpe_upm.zm();
679
680 // compute UA+BV expecting to get 1 (in Z_pe[x])
681 scoped_mpz_vector tmp1(nm);
682 scoped_mpz_vector tmp2(nm);
683 zpe_upm.mul(U.size(), U.c_ptr(), A.size(), A.c_ptr(), tmp1);
684 zpe_upm.mul(V.size(), V.c_ptr(), B.size(), B.c_ptr(), tmp2);
685 scoped_mpz_vector one(nm);
686 zpe_upm.add(tmp1.size(), tmp1.c_ptr(), tmp2.size(), tmp2.c_ptr(), one);
687 if (one.size() != 1 || !nm.is_one(one[0])) {
688 TRACE("polynomial::factorization::bughunt",
689 tout << "sage: R.<x> = Zmod(" << nm.to_string(zpe_upm.m().p()) << ")['x']" << endl;
690 tout << "sage: A = "; zpe_upm.display(tout, A); tout << endl;
691 tout << "sage: B = "; zpe_upm.display(tout, B); tout << endl;
692 tout << "sage: U = "; zpe_upm.display(tout, U); tout << endl;
693 tout << "sage: V = "; zpe_upm.display(tout, V); tout << endl;
694 tout << "sage: print (1 - UA - VB)" << endl;
695 tout << "sage: print 'expected 0'" << endl;
696 );
697 return false;
698 }
699
700 return true;
701 }
702
703 /**
704 Lift C = A*B from Z_p[x] to Z_{p^e}[x] such that:
705 * A = A_lift (mod p)
706 * B = B_lift (mod p)
707 * C = A*B (mod p^e)
708 */
hensel_lift_quadratic(z_manager & upm,numeral_vector const & C,zp_manager & zpe_upm,numeral_vector & A,numeral_vector & B,unsigned e)709 void hensel_lift_quadratic(z_manager& upm, numeral_vector const & C,
710 zp_manager & zpe_upm, numeral_vector & A, numeral_vector & B, unsigned e) {
711 z_numeral_manager & nm = upm.zm();
712
713 TRACE("polynomial::factorization::bughunt",
714 tout << "polynomial::hansel_lift_quadratic(";
715 tout << "A = "; upm.display(tout, A); tout << ", ";
716 tout << "B = "; upm.display(tout, B); tout << ", ";
717 tout << "C = "; upm.display(tout, C); tout << ", ";
718 tout << "p = " << nm.to_string(zpe_upm.m().p()) << ", e = " << e << ")" << endl;
719 );
720
721 // we create a new Z_p manager, since we'll be changing the input one
722 zp_manager zp_upm(upm.lim(), nm);
723 zp_upm.set_zp(zpe_upm.m().p());
724
725 // get the U, V, such that A*U + B*V = 1 (mod p)
726 scoped_mpz_vector U(nm), V(nm), D(nm);
727 zp_upm.ext_gcd(A.size(), A.c_ptr(), B.size(), B.c_ptr(), U, V, D);
728 SASSERT(D.size() == 1 && zp_upm.m().is_one(D[0]));
729
730 // we start lifting from (a = p, b = p, r = p)
731 scoped_mpz_vector A_lifted(nm), B_lifted(nm);
732 for (unsigned k = 1; k < e; k *= 2) {
733 upm.checkpoint();
734 // INVARIANT(a = p^e, b = p^e, r = gcd(a, b) = p^e):
735 // C = AB (mod b), UA + VB = 1 (mod a)
736 // deg(U) < deg(B), dev(V) < deg(A), deg(C) = deg(A) + deg(V)
737 // gcd(l(A), r) = 1
738
739 // regular hensel lifting from a to b*r, here from pe -> pk*pk = p^{k*k}
740 numeral const & pe = zpe_upm.m().p();
741
742 hensel_lift(upm, pe, pe, pe, U, A, V, B, C, A_lifted, B_lifted);
743 // now we have C = AB (mod b*r)
744 TRACE("polynomial::factorization::bughunt",
745 tout << "A_lifted = "; upm.display(tout, A_lifted); tout << endl;
746 tout << "B_lifted = "; upm.display(tout, B_lifted); tout << endl;
747 tout << "C = "; upm.display(tout, C); tout << endl;
748 );
749
750 // we employ similar reasoning as in the regular hansel lemma now
751 // we need to lift UA + VB = 1 (mod a)
752 // since after the lift, we still have UA + VB = 1 (mod a) we know that (1 - UA - VB)/a is in Z[x]
753 // so we can compute g = (1 - UA - VB)/a
754 // we need U1 and V1 such that U1A + V1B = 1 (mod a^2), with U1 = U mod a, V1 = V mod a
755 // hence U1 = U + aS, V1 = V + aT and we need
756 // (U + aS)A + (V + aT)B = 1 (mod a^2) same as
757 // UA + VB + a(SA + TB) = 1 (mod a^2) same as
758 // SA + TB = g (mod a) hence
759 // (gU + tB)A + (gV - tA)B = g (mod a) will be a solution and we pick t such that deg(gV - tA) < deg(A)
760
761 // compute g
762 scoped_mpz_vector tmp1(nm), g(nm);
763 g.push_back(numeral());
764 nm.set(g.back(), 1); // g = 1
765 upm.mul(A_lifted.size(), A_lifted.c_ptr(), U.size(), U.c_ptr(), tmp1); // tmp1 = AU
766 upm.sub(g.size(), g.c_ptr(), tmp1.size(), tmp1.c_ptr(), g); // g = 1 - UA
767 upm.mul(B_lifted.size(), B_lifted.c_ptr(), V.size(), V.c_ptr(), tmp1); // tmp1 = BV
768 upm.sub(g.size(), g.c_ptr(), tmp1.size(), tmp1.c_ptr(), g); // g = 1 - UA - VB
769 upm.div(g, pe);
770 to_zp_manager(zpe_upm, g);
771 TRACE("polynomial::factorization::bughunt",
772 tout << "g = (1 - A_lifted*U - B_lifted*V)/" << nm.to_string(pe) << endl;
773 tout << "g == "; upm.display(tout, g); tout << endl;
774 );
775
776 // compute the S, T
777 scoped_mpz_vector S(nm), T(nm), t(nm), tmp2(nm);
778 zpe_upm.mul(g.size(), g.c_ptr(), V.size(), V.c_ptr(), tmp1); // tmp1 = gV
779 zpe_upm.div_rem(tmp1.size(), tmp1.c_ptr(), A.size(), A.c_ptr(), t, T); // T = gV - tA, deg(T) < deg(A)
780 zpe_upm.mul(g.size(), g.c_ptr(), U.size(), U.c_ptr(), tmp1); // tmp1 = gU
781 zpe_upm.mul(t.size(), t.c_ptr(), B.size(), B.c_ptr(), tmp2); // tmp2 = tB
782 zpe_upm.add(tmp1.size(), tmp1.c_ptr(), tmp2.size(), tmp2.c_ptr(), S);
783
784 // now update U = U + a*S and V = V + a*T
785 upm.mul(S.size(), S.c_ptr(), pe);
786 upm.mul(T.size(), T.c_ptr(), pe);
787 upm.add(U.size(), U.c_ptr(), S.size(), S.c_ptr(), U);
788 upm.add(V.size(), V.c_ptr(), T.size(), T.c_ptr(), V); // deg(V) < deg(A), deg(T) < deg(A) => deg(V') < deg(A)
789
790 // we go quadratic
791 zpe_upm.m().set_p_sq();
792 to_zp_manager(zpe_upm, U);
793 to_zp_manager(zpe_upm, V);
794 to_zp_manager(zpe_upm, A_lifted);
795 to_zp_manager(zpe_upm, B_lifted);
796
797 // at this point we have INVARIANT(a = (p^e)^2, b = (p^e)^2, r = (p^e)^2)
798 A.swap(A_lifted);
799 B.swap(B_lifted);
800
801 CASSERT("polynomial::factorizatio::bughunt", check_quadratic_hensel(zpe_upm, U, A, V, B));
802 }
803 }
804
check_hensel_lift(z_manager & upm,numeral_vector const & f,zp_factors const & zp_fs,zp_factors const & zpe_fs,unsigned e)805 bool check_hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors const & zp_fs, zp_factors const & zpe_fs, unsigned e) {
806 numeral_manager & nm(upm.m());
807
808 zp_manager & zp_upm = zp_fs.upm();
809 zp_manager & zpe_upm = zpe_fs.upm();
810
811 numeral const & p = zp_fs.nm().p();
812 numeral const & pe = zpe_fs.nm().p();
813
814 scoped_numeral power(nm);
815 nm.power(p, e, power);
816 if (!nm.ge(pe, power)) {
817 return false;
818 }
819
820 // check f = lc(f) * zp_fs (mod p)
821 scoped_numeral_vector mult_zp(nm), f_zp(nm);
822 zp_fs.multiply(mult_zp);
823 to_zp_manager(zp_upm, f, f_zp);
824 zp_upm.mul(mult_zp, f_zp.back());
825 if (!upm.eq(mult_zp, f_zp)) {
826 TRACE("polynomial::factorization::bughunt",
827 tout << "f = "; upm.display(tout, f); tout << endl;
828 tout << "zp_fs = " << zp_fs << endl;
829 tout << "sage: R.<x> = Zmod(" << nm.to_string(p) << ")['x']" << endl;
830 tout << "sage: mult_zp = "; upm.display(tout, mult_zp); tout << endl;
831 tout << "sage: f_zp = "; upm.display(tout, f_zp); tout << endl;
832 tout << "sage: mult_zp == f_zp" << endl;
833 );
834 return false;
835 }
836
837 // check individual factors
838 if (zpe_fs.distinct_factors() != zp_fs.distinct_factors()) {
839 return false;
840 }
841
842 // check f = lc(f) * zpe_fs (mod p^e)
843 scoped_numeral_vector mult_zpe(nm), f_zpe(nm);
844 zpe_fs.multiply(mult_zpe);
845 to_zp_manager(zpe_upm, f, f_zpe);
846 zpe_upm.mul(mult_zpe, f_zpe.back());
847 if (!upm.eq(mult_zpe, f_zpe)) {
848 TRACE("polynomial::factorization::bughunt",
849 tout << "f = "; upm.display(tout, f); tout << endl;
850 tout << "zpe_fs = " << zpe_fs << endl;
851 tout << "sage: R.<x> = Zmod(" << nm.to_string(pe) << ")['x']" << endl;
852 tout << "sage: mult_zpe = "; upm.display(tout, mult_zpe); tout << endl;
853 tout << "sage: f_zpe = "; upm.display(tout, f_zpe); tout << endl;
854 tout << "sage: mult_zpe == f_zpe" << endl;
855 );
856 return false;
857 }
858
859 return true;
860 }
861
check_individual_lift(zp_manager & zp_upm,numeral_vector const & A_p,zp_manager & zpe_upm,numeral_vector const & A_pe)862 bool check_individual_lift(zp_manager & zp_upm, numeral_vector const & A_p, zp_manager & zpe_upm, numeral_vector const & A_pe) {
863 scoped_numeral_vector A_pe_p(zp_upm.m());
864 to_zp_manager(zp_upm, A_pe, A_pe_p);
865
866 if (!zp_upm.eq(A_p, A_pe_p)) {
867 return false;
868 }
869
870 return true;
871 }
872
hensel_lift(z_manager & upm,numeral_vector const & f,zp_factors const & zp_fs,unsigned e,zp_factors & zpe_fs)873 void hensel_lift(z_manager & upm, numeral_vector const & f, zp_factors const & zp_fs, unsigned e, zp_factors & zpe_fs) {
874 SASSERT(zp_fs.total_factors() > 0);
875
876 zp_numeral_manager & zp_nm = zp_fs.nm();
877 zp_manager & zp_upm = zp_fs.upm();
878 z_numeral_manager & nm = zp_nm.m();
879
880 SASSERT(nm.is_one(zp_fs.get_constant()));
881
882 zp_numeral_manager & zpe_nm = zpe_fs.nm();
883 zp_manager & zpe_upm = zpe_fs.upm();
884 zpe_nm.set_zp(zp_nm.p());
885
886 TRACE("polynomial::factorization::bughunt",
887 tout << "polynomial::hansel_lift("; upm.display(tout, f); tout << ", " << zp_fs << ")" << endl;
888 );
889
890 // lift the factors one by one
891 scoped_mpz_vector A(nm), B(nm), C(nm), f_parts(nm); // these will all be in Z_p
892
893 // copy of f, that we'll be cutting parts of
894 upm.set(f.size(), f.c_ptr(), f_parts);
895
896 // F_k are factors mod Z_p, A_k the factors mod p^e
897 // the invariant we keep is that:
898 // (1) f_parts = C = lc(f) * \prod_{k = i}^n F_k, C in Z_p[x]
899 // (2) A_k = F_k (mod p), for k < i
900 // (3) f = (\prod_{k < i} A_k) * f_parts (mod p^e)
901 for (int i = 0, i_end = zp_fs.distinct_factors()-1; i < i_end; ++ i) {
902 SASSERT(zp_fs.get_degree(i) == 1); // p was chosen so that f is square-free
903
904 // F_i = A (mod Z_p)
905 zp_upm.set(zp_fs[i].size(), zp_fs[i].c_ptr(), A);
906 TRACE("polynomial::factorization::bughunt",
907 tout << "A = "; upm.display(tout, A); tout << endl;
908 );
909
910 // C = f_parts (mod Z_p)
911 if (i > 0) {
912 to_zp_manager(zp_upm, f_parts, C);
913 }
914 else {
915 // first time around, we don't have p^e yet, so first time we just compute C
916 zp_fs.multiply(C);
917 scoped_mpz lc(nm);
918 zp_nm.set(lc, f.back());
919 zp_upm.mul(C, lc);
920 }
921 TRACE("polynomial::factorization::bughunt",
922 tout << "C = "; upm.display(tout, C); tout << endl;
923 );
924
925 // we take B to be what's left from C and A
926 zp_upm.div(C.size(), C.c_ptr(), A.size(), A.c_ptr(), B);
927 TRACE("polynomial::factorization::bughunt",
928 tout << "B = "; upm.display(tout, B); tout << endl;
929 );
930
931 // lift A and B to p^e (this will change p^e and put it zp_upm)
932 zpe_nm.set_zp(zp_nm.p());
933 hensel_lift_quadratic(upm, f_parts, zpe_upm, A, B, e);
934 CASSERT("polynomial::factorizatio::bughunt", check_individual_lift(zp_upm, zp_fs[i], zpe_upm, A));
935 TRACE("polynomial::factorization",
936 tout << "lifted to " << nm.to_string(zpe_upm.m().p()) << endl;
937 tout << "A = "; upm.display(tout, A); tout << endl;
938 tout << "B = "; upm.display(tout, B); tout << endl;
939 );
940
941 // if this is the first time round, we also construct f_parts (now we have correct p^e)
942 if (i == 0) {
943 to_zp_manager(zpe_upm, f, f_parts);
944 }
945
946 // take the lifted A out of f_parts
947 zpe_upm.div(f_parts.size(), f_parts.c_ptr(), A.size(), A.c_ptr(), f_parts);
948
949 // add the lifted factor (kills A)
950 zpe_fs.push_back_swap(A, 1);
951 }
952
953 // we have one last factor, but it also contains lc(f), so take it out
954 scoped_mpz lc_inv(nm);
955 zpe_nm.set(lc_inv, f.back());
956 zpe_nm.inv(lc_inv);
957 zpe_upm.mul(B, lc_inv);
958 zpe_fs.push_back_swap(B, 1);
959
960 TRACE("polynomial::factorization::bughunt",
961 tout << "polynomial::hansel_lift("; upm.display(tout, f); tout << ", " << zp_fs << ") => " << zpe_fs << endl;
962 );
963
964 CASSERT("polynomial::factorizatio::bughunt", check_hensel_lift(upm, f, zp_fs, zpe_fs, e));
965 }
966
967 // get a bound on B for the factors of f with degree less or equal to deg(f)/2
968 // and then choose e to be smallest such that p^e > 2*lc(f)*B, we use the mignotte
969 //
970 // from [3, pg 134]
971 // |p| = sqrt(\sum |p_i|^2). If a = \sum a_i x^i, b = \sum b_j, deg b = n, and b divides a, then for all j
972 //
973 // |b_j| <= (n-1 over j)|a| + (n-1 over j-1)|lc(a)|
974 //
975 // when factoring a polynomial, we find a bound B for a factor of f of degree <= deg(f)/2
976 // allowing both positive and negative as coefficients of b, we now want a power p^e such that all coefficients
977 // can be represented, i.e. [-B, B] \subset [-\ceil(p^e/2), \ceil(p^e)/2], so we peek e, such that p^e >= 2B
978
mignotte_bound(z_manager & upm,numeral_vector const & f,numeral const & p)979 static unsigned mignotte_bound(z_manager & upm, numeral_vector const & f, numeral const & p) {
980 numeral_manager & nm = upm.m();
981
982 SASSERT(upm.degree(f) >= 2);
983 unsigned n = upm.degree(f)/2; // >= 1
984
985 // get the approximation for the norm of a
986 scoped_numeral f_norm(nm);
987 for (unsigned i = 0; i < f.size(); ++ i) {
988 if (!nm.is_zero(f[i])) {
989 nm.addmul(f_norm, f[i], f[i], f_norm);
990 }
991 }
992 nm.root(f_norm, 2);
993
994 // by above we can pick (n-1 over (n-1/2))|a| + (n-1 over (n-1)/2)lc(a)
995 // we approximate both binomial-coefficients with 2^(n-1), so to get 2B we use 2^n(f_norm + lc(f))
996 scoped_numeral bound(nm);
997 nm.set(bound, 1);
998 nm.mul2k(bound, n, bound);
999 scoped_numeral tmp(nm);
1000 nm.set(tmp, f.back());
1001 nm.abs(tmp);
1002 nm.add(f_norm, tmp, f_norm);
1003 nm.mul(bound, f_norm, bound);
1004
1005 // we need e such that p^e >= B
1006 nm.set(tmp, p);
1007 unsigned e;
1008 for (e = 1; nm.le(tmp, bound); e *= 2) {
1009 nm.mul(tmp, tmp, tmp);
1010 }
1011
1012 return e;
1013 }
1014
1015 /**
1016 \brief Given f from Z[x] that is square free, it factors it.
1017 This method also assumes f is primitive.
1018 */
factor_square_free(z_manager & upm,numeral_vector const & f,factors & fs,unsigned k,factor_params const & params)1019 bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs, unsigned k, factor_params const & params) {
1020 TRACE("polynomial::factorization::bughunt",
1021 tout << "sage: f = "; upm.display(tout, f); tout << endl;
1022 tout << "sage: if (not f.is_squarefree()): print 'Error, f is not square-free'" << endl;
1023 tout << "sage: print 'Factoring :', f" << endl;
1024 tout << "sage: print 'Expected factors: ', f.factor()" << endl;
1025 );
1026
1027 numeral_manager & nm = upm.m();
1028
1029 // This method assumes f is primitive. Thus, the content of f must be one
1030 DEBUG_CODE({
1031 scoped_numeral f_cont(nm);
1032 nm.gcd(f.size(), f.c_ptr(), f_cont);
1033 SASSERT(f.size() == 0 || nm.is_one(f_cont));
1034 });
1035
1036 scoped_numeral_vector f_pp(nm);
1037 upm.set(f.size(), f.c_ptr(), f_pp);
1038
1039 // make sure the leading coefficient is positive
1040 if (!f_pp.empty() && nm.is_neg(f_pp[f_pp.size() - 1])) {
1041 for (unsigned i = 0; i < f_pp.size(); i++)
1042 nm.neg(f_pp[i]);
1043 // flip sign constant if k is odd
1044 if (k % 2 == 1) {
1045 scoped_numeral c(nm);
1046 nm.set(c, fs.get_constant());
1047 nm.neg(c);
1048 fs.set_constant(c);
1049 }
1050 }
1051
1052 TRACE("polynomial::factorization::bughunt",
1053 tout << "sage: f_pp = "; upm.display(tout, f_pp); tout << endl;
1054 tout << "sage: if (not (f_pp * 1 == f): print 'Error, content computation wrong'" << endl;
1055 );
1056
1057 // the variables we'll be using and updating in Z_p
1058 scoped_numeral p(nm);
1059 nm.set(p, 2);
1060 zp_manager zp_upm(upm.lim(), nm.m());
1061 zp_upm.set_zp(p);
1062 zp_factors zp_fs(zp_upm);
1063 scoped_numeral zp_fs_p(nm); nm.set(zp_fs_p, 2);
1064
1065 // we keep all the possible sets of degrees of factors in this set
1066 factorization_degree_set degree_set(zp_upm);
1067
1068 // we try get some number of factorizations in Z_p, for some primes
1069 // get the prime p such that
1070 // (1) (f_prim mod p) stays square-free
1071 // (2) l(f_prim) mod p doesn't vanish, i.e. we don't get a polynomial of smaller degree
1072 prime_iterator prime_it;
1073 scoped_numeral gcd_tmp(nm);
1074 unsigned trials = 0;
1075 TRACE("polynomial::factorization::bughunt", tout << "trials: " << params.m_p_trials << "\n";);
1076 while (trials <= params.m_p_trials) {
1077 upm.checkpoint();
1078 // construct prime to check
1079 uint64_t next_prime = prime_it.next();
1080 if (next_prime > params.m_max_p) {
1081 fs.push_back(f_pp, k);
1082 return false;
1083 }
1084 nm.set(p, next_prime);
1085 zp_upm.set_zp(p);
1086
1087 // we need gcd(lc(f_pp), p) = 1
1088 nm.gcd(p, f_pp.back(), gcd_tmp);
1089 TRACE("polynomial::factorization::bughunt",
1090 tout << "sage: if (not (gcd(" << nm.to_string(p) << ", " << nm.to_string(f_pp.back()) << ")) == " <<
1091 nm.to_string(gcd_tmp) << "): print 'Error, wrong gcd'" << endl;
1092 );
1093 if (!nm.is_one(gcd_tmp)) {
1094 continue;
1095 }
1096
1097 // if it's not square free, we also try something else
1098 scoped_numeral_vector f_pp_zp(nm);
1099 to_zp_manager(zp_upm, f_pp, f_pp_zp);
1100
1101 TRACE("polynomial::factorization::bughunt",
1102 tout << "sage: Rp.<x_p> = GF(" << nm.to_string(p) << ")['x_p']"; tout << endl;
1103 tout << "sage: f_pp_zp = "; zp_upm.display(tout, f_pp_zp, "x_p"); tout << endl;
1104 );
1105
1106 if (!zp_upm.is_square_free(f_pp_zp.size(), f_pp_zp.c_ptr()))
1107 continue;
1108
1109 // we make it monic
1110 zp_upm.mk_monic(f_pp_zp.size(), f_pp_zp.c_ptr());
1111
1112 // found a candidate, factorize in Z_p and add back the constant
1113 zp_factors current_fs(zp_upm);
1114 bool factored = zp_factor_square_free(zp_upm, f_pp_zp, current_fs);
1115 if (!factored) {
1116 fs.push_back(f_pp, k);
1117 return true;
1118 }
1119
1120 // get the factors degrees set
1121 factorization_degree_set current_degree_set(current_fs);
1122 if (degree_set.max_degree() == 0) {
1123 // first time, initialize
1124 degree_set.swap(current_degree_set);
1125 }
1126 else {
1127 degree_set.intersect(current_degree_set);
1128 }
1129
1130 // if the degree set is trivial, we are done
1131 if (degree_set.is_trivial()) {
1132 fs.push_back(f_pp, k);
1133 return true;
1134 }
1135
1136 // we found a candidate, lets keep it if it has less factors than the current best
1137 trials ++;
1138 if (zp_fs.distinct_factors() == 0 || zp_fs.total_factors() > current_fs.total_factors()) {
1139 zp_fs.swap(current_fs);
1140 nm.set(zp_fs_p, p);
1141
1142 TRACE("polynomial::factorization::bughunt",
1143 tout << "best zp factorization (Z_" << nm.to_string(zp_fs_p) << "): ";
1144 tout << zp_fs << endl;
1145 tout << "best degree set: "; degree_set.display(tout); tout << endl;
1146 );
1147 }
1148 }
1149 #ifndef _EXTERNAL_RELEASE
1150 IF_VERBOSE(FACTOR_VERBOSE_LVL, verbose_stream() << "(polynomial-factorization :at GF_" << nm.to_string(zp_upm.p()) << ")" << std::endl;);
1151 #endif
1152
1153 // make sure to set the zp_manager back to modulo zp_fs_p
1154 zp_upm.set_zp(zp_fs_p);
1155
1156 TRACE("polynomial::factorization::bughunt",
1157 tout << "best zp factorization (Z_" << nm.to_string(zp_fs_p) << "): " << zp_fs << endl;
1158 tout << "best degree set: "; degree_set.display(tout); tout << endl;
1159 );
1160
1161 // get a bound on B for the factors of f_pp with degree less or equal to deg(f)/2
1162 // and then choose e to be smallest such that p^e > 2*lc(f)*B, we use the mignotte
1163 unsigned e = mignotte_bound(upm, f_pp, zp_fs_p);
1164 TRACE("polynomial::factorization::bughunt",
1165 tout << "out p = " << nm.to_string(zp_fs_p) << ", and we'll work p^e for e = " << e << endl;
1166 );
1167
1168 // we got a prime factoring, so we do the lifting now
1169 zp_manager zpe_upm(upm.lim(), nm.m());
1170 zpe_upm.set_zp(zp_fs_p);
1171 zp_numeral_manager & zpe_nm = zpe_upm.m();
1172
1173 zp_factors zpe_fs(zpe_upm);
1174 // this might give something bigger than p^e, but the lifting procedure will update the zpe_nm
1175 // zp factors are monic, so will be the zpe factors, i.e. f_pp = zpe_fs * lc(f_pp) (mod p^e)
1176 hensel_lift(upm, f_pp, zp_fs, e, zpe_fs);
1177
1178 #ifndef _EXTERNAL_RELEASE
1179 IF_VERBOSE(FACTOR_VERBOSE_LVL, verbose_stream() << "(polynomial-factorization :num-candidate-factors " << zpe_fs.distinct_factors() << ")" << std::endl;);
1180 #endif
1181
1182 // the leading coefficient of f_pp mod p^e
1183 scoped_numeral f_pp_lc(nm);
1184 zpe_nm.set(f_pp_lc, f_pp.back());
1185
1186 // we always keep in f_pp the actual primitive part f_pp*lc(f_pp)
1187 upm.mul(f_pp, f_pp_lc);
1188
1189 // now we go through the combinations of factors to check construct the factorization
1190 ufactorization_combination_iterator it(zpe_fs, degree_set);
1191 scoped_numeral_vector trial_factor(nm), trial_factor_quo(nm);
1192 scoped_numeral trial_factor_cont(nm);
1193 TRACE("polynomial::factorization::bughunt",
1194 tout << "STARTING TRIAL DIVISION" << endl;
1195 tout << "zpe_fs" << zpe_fs << endl;
1196 tout << "degree_set = "; degree_set.display(tout); tout << endl;
1197 );
1198 bool result = true;
1199 bool remove = false;
1200 unsigned counter = 0;
1201 while (it.next(remove)) {
1202 upm.checkpoint();
1203 counter++;
1204 if (counter > params.m_max_search_size) {
1205 // stop search
1206 result = false;
1207 break;
1208 }
1209 //
1210 // our bound ensures we can extract the right factors of degree at most 1/2 of the original
1211 // so, if out trial factor has degree bigger than 1/2, we need to take the rest of the factors
1212 // but, if we take the rest and it works, it doesn't mean that the rest is factorized, so we still take out
1213 // the original factor
1214 bool using_left = it.current_degree() <= zp_fs.get_degree()/2;
1215 if (using_left) {
1216 // do a quick check first
1217 scoped_numeral tmp(nm);
1218 it.get_left_tail_coeff(f_pp_lc, tmp);
1219 if (!nm.divides(tmp, f_pp[0])) {
1220 // don't remove this combination
1221 remove = false;
1222 continue;
1223 }
1224 it.left(trial_factor);
1225 }
1226 else {
1227 // do a quick check first
1228 scoped_numeral tmp(nm);
1229 it.get_right_tail_coeff(f_pp_lc, tmp);
1230 if (!nm.divides(tmp, f_pp[0])) {
1231 // don't remove this combination
1232 remove = false;
1233 continue;
1234 }
1235 it.right(trial_factor);
1236 }
1237
1238 // add the lc(f_pp) to the trial divisor
1239 zpe_upm.mul(trial_factor, f_pp_lc);
1240 TRACE("polynomial::factorization::bughunt",
1241 tout << "f_pp*lc(f_pp) = "; upm.display(tout, f_pp); tout << endl;
1242 tout << "trial_factor = "; upm.display(tout, trial_factor); tout << endl;
1243 );
1244
1245 bool true_factor = upm.exact_div(f_pp, trial_factor, trial_factor_quo);
1246
1247 TRACE("polynomial::factorization::bughunt",
1248 tout << "trial_factor = "; upm.display(tout, trial_factor); tout << endl;
1249 tout << "trial_factor_quo = "; upm.display(tout, trial_factor_quo); tout << endl;
1250 tout << "result = " << (true_factor ? "true" : "false") << endl;
1251 );
1252
1253 // if division is precise we have a factor
1254 if (true_factor) {
1255 if (!using_left) {
1256 // as noted above, we still use the original factor
1257 trial_factor.swap(trial_factor_quo);
1258 }
1259 // We need to get the content out of the factor
1260 upm.get_primitive_and_content(trial_factor, trial_factor, trial_factor_cont);
1261 // add the factor
1262 fs.push_back(trial_factor, k);
1263 // we continue with the quotient (with the content added back)
1264 // but we also have to keep lc(f_pp)*f_pp
1265 upm.get_primitive_and_content(trial_factor_quo, f_pp, trial_factor_cont);
1266 nm.set(f_pp_lc, f_pp.back());
1267 upm.mul(f_pp, f_pp_lc);
1268 // but we also remove it from the iterator
1269 remove = true;
1270 }
1271 else {
1272 // don't remove this combination
1273 remove = false;
1274 }
1275 TRACE("polynomial::factorization::bughunt",
1276 tout << "factors = " << fs << endl;
1277 tout << "f_pp*lc(f_pp) = "; upm.display(tout, f_pp); tout << endl;
1278 tout << "lc(f_pp) = " << f_pp_lc << endl;
1279 );
1280 }
1281 #ifndef _EXTERNAL_RELEASE
1282 IF_VERBOSE(FACTOR_VERBOSE_LVL, verbose_stream() << "(polynomial-factorization :search-size " << counter << ")" << std::endl;);
1283 #endif
1284
1285 // add the what's left to the factors (if not a constant)
1286 if (f_pp.size() > 1) {
1287 upm.div(f_pp, f_pp_lc);
1288 fs.push_back(f_pp, k);
1289 }
1290 else {
1291 // if a constant it must be 1 (it was primitive)
1292 SASSERT(f_pp.size() == 1 && nm.is_one(f_pp.back()));
1293 }
1294
1295 return result;
1296 }
1297
factor_square_free(z_manager & upm,numeral_vector const & f,factors & fs,factor_params const & params)1298 bool factor_square_free(z_manager & upm, numeral_vector const & f, factors & fs, factor_params const & params) {
1299 return factor_square_free(upm, f, fs, 1, params);
1300 }
1301
1302 }; // end upolynomial namespace
1303