Fixed a bug in squarefree(). Some polynomials were erroneously classified as
[ginac.git] / ginac / factor.cpp
1 /** @file factor.cpp
2  *
3  *  Polynomial factorization code (implementation).
4  *
5  *  Algorithms used can be found in
6  *    [W1]  An Improved Multivariate Polynomial Factoring Algorithm,
7  *          P.S.Wang, Mathematics of Computation, Vol. 32, No. 144 (1978) 1215--1231.
8  *    [GCL] Algorithms for Computer Algebra,
9  *          K.O.Geddes, S.R.Czapor, G.Labahn, Springer Verlag, 1992.
10  */
11
12 /*
13  *  GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
14  *
15  *  This program is free software; you can redistribute it and/or modify
16  *  it under the terms of the GNU General Public License as published by
17  *  the Free Software Foundation; either version 2 of the License, or
18  *  (at your option) any later version.
19  *
20  *  This program is distributed in the hope that it will be useful,
21  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
22  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
23  *  GNU General Public License for more details.
24  *
25  *  You should have received a copy of the GNU General Public License
26  *  along with this program; if not, write to the Free Software
27  *  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
28  */
29
30 //#define DEBUGFACTOR
31
32 #include "factor.h"
33
34 #include "ex.h"
35 #include "numeric.h"
36 #include "operators.h"
37 #include "inifcns.h"
38 #include "symbol.h"
39 #include "relational.h"
40 #include "power.h"
41 #include "mul.h"
42 #include "normal.h"
43 #include "add.h"
44
45 #include <algorithm>
46 #include <cmath>
47 #include <limits>
48 #include <list>
49 #include <vector>
50 #ifdef DEBUGFACTOR
51 #include <ostream>
52 #endif
53 using namespace std;
54
55 #include <cln/cln.h>
56 using namespace cln;
57
58 namespace GiNaC {
59
60 #ifdef DEBUGFACTOR
61 #define DCOUT(str) cout << #str << endl
62 #define DCOUTVAR(var) cout << #var << ": " << var << endl
63 #define DCOUT2(str,var) cout << #str << ": " << var << endl
64 #else
65 #define DCOUT(str)
66 #define DCOUTVAR(var)
67 #define DCOUT2(str,var)
68 #endif
69
70 // anonymous namespace to hide all utility functions
71 namespace {
72
73 typedef vector<cl_MI> mvec;
74 #ifdef DEBUGFACTOR
75 ostream& operator<<(ostream& o, const vector<int>& v)
76 {
77         vector<int>::const_iterator i = v.begin(), end = v.end();
78         while ( i != end ) {
79                 o << *i++ << " ";
80         }
81         return o;
82 }
83 ostream& operator<<(ostream& o, const vector<cl_I>& v)
84 {
85         vector<cl_I>::const_iterator i = v.begin(), end = v.end();
86         while ( i != end ) {
87                 o << *i << "[" << i-v.begin() << "]" << " ";
88                 ++i;
89         }
90         return o;
91 }
92 ostream& operator<<(ostream& o, const vector<cl_MI>& v)
93 {
94         vector<cl_MI>::const_iterator i = v.begin(), end = v.end();
95         while ( i != end ) {
96                 o << *i << "[" << i-v.begin() << "]" << " ";
97                 ++i;
98         }
99         return o;
100 }
101 ostream& operator<<(ostream& o, const vector< vector<cl_MI> >& v)
102 {
103         vector< vector<cl_MI> >::const_iterator i = v.begin(), end = v.end();
104         while ( i != end ) {
105                 o << i-v.begin() << ": " << *i << endl;
106                 ++i;
107         }
108         return o;
109 }
110 #endif
111
112 ////////////////////////////////////////////////////////////////////////////////
113 // modular univariate polynomial code
114
115 typedef std::vector<cln::cl_MI> umodpoly;
116 typedef std::vector<cln::cl_I> upoly;
117 typedef vector<umodpoly> upvec;
118
119 // COPY FROM UPOLY.HPP
120
121 // CHANGED size_t -> int !!!
122 template<typename T> static int degree(const T& p)
123 {
124         return p.size() - 1;
125 }
126
127 template<typename T> static typename T::value_type lcoeff(const T& p)
128 {
129         return p[p.size() - 1];
130 }
131
132 static bool normalize_in_field(umodpoly& a)
133 {
134         if (a.size() == 0)
135                 return true;
136         if ( lcoeff(a) == a[0].ring()->one() ) {
137                 return true;
138         }
139
140         const cln::cl_MI lc_1 = recip(lcoeff(a));
141         for (std::size_t k = a.size(); k-- != 0; )
142                 a[k] = a[k]*lc_1;
143         return false;
144 }
145
146 template<typename T> static void
147 canonicalize(T& p, const typename T::size_type hint = std::numeric_limits<typename T::size_type>::max())
148 {
149         if (p.empty())
150                 return;
151
152         std::size_t i = p.size() - 1;
153         // Be fast if the polynomial is already canonicalized
154         if (!zerop(p[i]))
155                 return;
156
157         if (hint < p.size())
158                 i = hint;
159
160         bool is_zero = false;
161         do {
162                 if (!zerop(p[i])) {
163                         ++i;
164                         break;
165                 }
166                 if (i == 0) {
167                         is_zero = true;
168                         break;
169                 }
170                 --i;
171         } while (true);
172
173         if (is_zero) {
174                 p.clear();
175                 return;
176         }
177
178         p.erase(p.begin() + i, p.end());
179 }
180
181 // END COPY FROM UPOLY.HPP
182
183 static void expt_pos(umodpoly& a, unsigned int q)
184 {
185         if ( a.empty() ) return;
186         cl_MI zero = a[0].ring()->zero(); 
187         int deg = degree(a);
188         a.resize(degree(a)*q+1, zero);
189         for ( int i=deg; i>0; --i ) {
190                 a[i*q] = a[i];
191                 a[i] = zero;
192         }
193 }
194
195 template<typename T>
196 static T operator+(const T& a, const T& b)
197 {
198         int sa = a.size();
199         int sb = b.size();
200         if ( sa >= sb ) {
201                 T r(sa);
202                 int i = 0;
203                 for ( ; i<sb; ++i ) {
204                         r[i] = a[i] + b[i];
205                 }
206                 for ( ; i<sa; ++i ) {
207                         r[i] = a[i];
208                 }
209                 canonicalize(r);
210                 return r;
211         }
212         else {
213                 T r(sb);
214                 int i = 0;
215                 for ( ; i<sa; ++i ) {
216                         r[i] = a[i] + b[i];
217                 }
218                 for ( ; i<sb; ++i ) {
219                         r[i] = b[i];
220                 }
221                 canonicalize(r);
222                 return r;
223         }
224 }
225
226 template<typename T>
227 static T operator-(const T& a, const T& b)
228 {
229         int sa = a.size();
230         int sb = b.size();
231         if ( sa >= sb ) {
232                 T r(sa);
233                 int i = 0;
234                 for ( ; i<sb; ++i ) {
235                         r[i] = a[i] - b[i];
236                 }
237                 for ( ; i<sa; ++i ) {
238                         r[i] = a[i];
239                 }
240                 canonicalize(r);
241                 return r;
242         }
243         else {
244                 T r(sb);
245                 int i = 0;
246                 for ( ; i<sa; ++i ) {
247                         r[i] = a[i] - b[i];
248                 }
249                 for ( ; i<sb; ++i ) {
250                         r[i] = -b[i];
251                 }
252                 canonicalize(r);
253                 return r;
254         }
255 }
256
257 static upoly operator*(const upoly& a, const upoly& b)
258 {
259         upoly c;
260         if ( a.empty() || b.empty() ) return c;
261
262         int n = degree(a) + degree(b);
263         c.resize(n+1, 0);
264         for ( int i=0 ; i<=n; ++i ) {
265                 for ( int j=0 ; j<=i; ++j ) {
266                         if ( j > degree(a) || (i-j) > degree(b) ) continue;
267                         c[i] = c[i] + a[j] * b[i-j];
268                 }
269         }
270         canonicalize(c);
271         return c;
272 }
273
274 static umodpoly operator*(const umodpoly& a, const umodpoly& b)
275 {
276         umodpoly c;
277         if ( a.empty() || b.empty() ) return c;
278
279         int n = degree(a) + degree(b);
280         c.resize(n+1, a[0].ring()->zero());
281         for ( int i=0 ; i<=n; ++i ) {
282                 for ( int j=0 ; j<=i; ++j ) {
283                         if ( j > degree(a) || (i-j) > degree(b) ) continue;
284                         c[i] = c[i] + a[j] * b[i-j];
285                 }
286         }
287         canonicalize(c);
288         return c;
289 }
290
291 static upoly operator*(const upoly& a, const cl_I& x)
292 {
293         if ( zerop(x) ) {
294                 upoly r;
295                 return r;
296         }
297         upoly r(a.size());
298         for ( size_t i=0; i<a.size(); ++i ) {
299                 r[i] = a[i] * x;
300         }
301         return r;
302 }
303
304 static upoly operator/(const upoly& a, const cl_I& x)
305 {
306         if ( zerop(x) ) {
307                 upoly r;
308                 return r;
309         }
310         upoly r(a.size());
311         for ( size_t i=0; i<a.size(); ++i ) {
312                 r[i] = exquo(a[i],x);
313         }
314         return r;
315 }
316
317 static umodpoly operator*(const umodpoly& a, const cl_MI& x)
318 {
319         umodpoly r(a.size());
320         for ( size_t i=0; i<a.size(); ++i ) {
321                 r[i] = a[i] * x;
322         }
323         canonicalize(r);
324         return r;
325 }
326
327 static void upoly_from_ex(upoly& up, const ex& e, const ex& x)
328 {
329         // assert: e is in Z[x]
330         int deg = e.degree(x);
331         up.resize(deg+1);
332         int ldeg = e.ldegree(x);
333         for ( ; deg>=ldeg; --deg ) {
334                 up[deg] = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
335         }
336         for ( ; deg>=0; --deg ) {
337                 up[deg] = 0;
338         }
339         canonicalize(up);
340 }
341
342 static void umodpoly_from_upoly(umodpoly& ump, const upoly& e, const cl_modint_ring& R)
343 {
344         int deg = degree(e);
345         ump.resize(deg+1);
346         for ( ; deg>=0; --deg ) {
347                 ump[deg] = R->canonhom(e[deg]);
348         }
349         canonicalize(ump);
350 }
351
352 static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_modint_ring& R)
353 {
354         // assert: e is in Z[x]
355         int deg = e.degree(x);
356         ump.resize(deg+1);
357         int ldeg = e.ldegree(x);
358         for ( ; deg>=ldeg; --deg ) {
359                 cl_I coeff = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
360                 ump[deg] = R->canonhom(coeff);
361         }
362         for ( ; deg>=0; --deg ) {
363                 ump[deg] = R->zero();
364         }
365         canonicalize(ump);
366 }
367
368 static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_I& modulus)
369 {
370         umodpoly_from_ex(ump, e, x, find_modint_ring(modulus));
371 }
372
373 static ex upoly_to_ex(const upoly& a, const ex& x)
374 {
375         if ( a.empty() ) return 0;
376         ex e;
377         for ( int i=degree(a); i>=0; --i ) {
378                 e += numeric(a[i]) * pow(x, i);
379         }
380         return e;
381 }
382
383 static ex umodpoly_to_ex(const umodpoly& a, const ex& x)
384 {
385         if ( a.empty() ) return 0;
386         cl_modint_ring R = a[0].ring();
387         cl_I mod = R->modulus;
388         cl_I halfmod = (mod-1) >> 1;
389         ex e;
390         for ( int i=degree(a); i>=0; --i ) {
391                 cl_I n = R->retract(a[i]);
392                 if ( n > halfmod ) {
393                         e += numeric(n-mod) * pow(x, i);
394                 } else {
395                         e += numeric(n) * pow(x, i);
396                 }
397         }
398         return e;
399 }
400
401 static upoly umodpoly_to_upoly(const umodpoly& a)
402 {
403         upoly e(a.size());
404         if ( a.empty() ) return e;
405         cl_modint_ring R = a[0].ring();
406         cl_I mod = R->modulus;
407         cl_I halfmod = (mod-1) >> 1;
408         for ( int i=degree(a); i>=0; --i ) {
409                 cl_I n = R->retract(a[i]);
410                 if ( n > halfmod ) {
411                         e[i] = n-mod;
412                 } else {
413                         e[i] = n;
414                 }
415         }
416         return e;
417 }
418
419 static umodpoly umodpoly_to_umodpoly(const umodpoly& a, const cl_modint_ring& R, unsigned int m)
420 {
421         umodpoly e;
422         if ( a.empty() ) return e;
423         cl_modint_ring oldR = a[0].ring();
424         size_t sa = a.size();
425         e.resize(sa+m, R->zero());
426         for ( size_t i=0; i<sa; ++i ) {
427                 e[i+m] = R->canonhom(oldR->retract(a[i]));
428         }
429         canonicalize(e);
430         return e;
431 }
432
433 /** Divides all coefficients of the polynomial a by the integer x.
434  *  All coefficients are supposed to be divisible by x. If they are not, the
435  *  the<cl_I> cast will raise an exception.
436  *
437  *  @param[in,out] a  polynomial of which the coefficients will be reduced by x
438  *  @param[in]     x  integer that divides the coefficients
439  */
440 static void reduce_coeff(umodpoly& a, const cl_I& x)
441 {
442         if ( a.empty() ) return;
443
444         cl_modint_ring R = a[0].ring();
445         umodpoly::iterator i = a.begin(), end = a.end();
446         for ( ; i!=end; ++i ) {
447                 // cln cannot perform this division in the modular field
448                 cl_I c = R->retract(*i);
449                 *i = cl_MI(R, the<cl_I>(c / x));
450         }
451 }
452
453 /** Calculates remainder of a/b.
454  *  Assertion: a and b not empty.
455  *
456  *  @param[in]  a  polynomial dividend
457  *  @param[in]  b  polynomial divisor
458  *  @param[out] r  polynomial remainder
459  */
460 static void rem(const umodpoly& a, const umodpoly& b, umodpoly& r)
461 {
462         int k, n;
463         n = degree(b);
464         k = degree(a) - n;
465         r = a;
466         if ( k < 0 ) return;
467
468         do {
469                 cl_MI qk = div(r[n+k], b[n]);
470                 if ( !zerop(qk) ) {
471                         for ( int i=0; i<n; ++i ) {
472                                 unsigned int j = n + k - 1 - i;
473                                 r[j] = r[j] - qk * b[j-k];
474                         }
475                 }
476         } while ( k-- );
477
478         fill(r.begin()+n, r.end(), a[0].ring()->zero());
479         canonicalize(r);
480 }
481
482 /** Calculates quotient of a/b.
483  *  Assertion: a and b not empty.
484  *
485  *  @param[in]  a  polynomial dividend
486  *  @param[in]  b  polynomial divisor
487  *  @param[out] q  polynomial quotient
488  */
489 static void div(const umodpoly& a, const umodpoly& b, umodpoly& q)
490 {
491         int k, n;
492         n = degree(b);
493         k = degree(a) - n;
494         q.clear();
495         if ( k < 0 ) return;
496
497         umodpoly r = a;
498         q.resize(k+1, a[0].ring()->zero());
499         do {
500                 cl_MI qk = div(r[n+k], b[n]);
501                 if ( !zerop(qk) ) {
502                         q[k] = qk;
503                         for ( int i=0; i<n; ++i ) {
504                                 unsigned int j = n + k - 1 - i;
505                                 r[j] = r[j] - qk * b[j-k];
506                         }
507                 }
508         } while ( k-- );
509
510         canonicalize(q);
511 }
512
513 /** Calculates quotient and remainder of a/b.
514  *  Assertion: a and b not empty.
515  *
516  *  @param[in]  a  polynomial dividend
517  *  @param[in]  b  polynomial divisor
518  *  @param[out] r  polynomial remainder
519  *  @param[out] q  polynomial quotient
520  */
521 static void remdiv(const umodpoly& a, const umodpoly& b, umodpoly& r, umodpoly& q)
522 {
523         int k, n;
524         n = degree(b);
525         k = degree(a) - n;
526         q.clear();
527         r = a;
528         if ( k < 0 ) return;
529
530         q.resize(k+1, a[0].ring()->zero());
531         do {
532                 cl_MI qk = div(r[n+k], b[n]);
533                 if ( !zerop(qk) ) {
534                         q[k] = qk;
535                         for ( int i=0; i<n; ++i ) {
536                                 unsigned int j = n + k - 1 - i;
537                                 r[j] = r[j] - qk * b[j-k];
538                         }
539                 }
540         } while ( k-- );
541
542         fill(r.begin()+n, r.end(), a[0].ring()->zero());
543         canonicalize(r);
544         canonicalize(q);
545 }
546
547 /** Calculates the GCD of polynomial a and b.
548  *
549  *  @param[in]  a  polynomial
550  *  @param[in]  b  polynomial
551  *  @param[out] c  GCD
552  */
553 static void gcd(const umodpoly& a, const umodpoly& b, umodpoly& c)
554 {
555         if ( degree(a) < degree(b) ) return gcd(b, a, c);
556
557         c = a;
558         normalize_in_field(c);
559         umodpoly d = b;
560         normalize_in_field(d);
561         umodpoly r;
562         while ( !d.empty() ) {
563                 rem(c, d, r);
564                 c = d;
565                 d = r;
566         }
567         normalize_in_field(c);
568 }
569
570 /** Calculates the derivative of the polynomial a.
571  *  
572  *  @param[in]  a  polynomial of which to take the derivative
573  *  @param[out] d  result/derivative
574  */
575 static void deriv(const umodpoly& a, umodpoly& d)
576 {
577         d.clear();
578         if ( a.size() <= 1 ) return;
579
580         d.insert(d.begin(), a.begin()+1, a.end());
581         int max = d.size();
582         for ( int i=1; i<max; ++i ) {
583                 d[i] = d[i] * (i+1);
584         }
585         canonicalize(d);
586 }
587
588 static bool unequal_one(const umodpoly& a)
589 {
590         if ( a.empty() ) return true;
591         return ( a.size() != 1 || a[0] != a[0].ring()->one() );
592 }
593
594 static bool equal_one(const umodpoly& a)
595 {
596         return ( a.size() == 1 && a[0] == a[0].ring()->one() );
597 }
598
599 /** Returns true if polynomial a is square free.
600  *
601  *  @param[in] a  polynomial to check
602  *  @return       true if polynomial is square free, false otherwise
603  */
604 static bool squarefree(const umodpoly& a)
605 {
606         umodpoly b;
607         deriv(a, b);
608         if ( b.empty() ) {
609                 return false;
610         }
611         umodpoly c;
612         gcd(a, b, c);
613         return equal_one(c);
614 }
615
616 // END modular univariate polynomial code
617 ////////////////////////////////////////////////////////////////////////////////
618
619 ////////////////////////////////////////////////////////////////////////////////
620 // modular matrix
621
622 class modular_matrix
623 {
624         friend ostream& operator<<(ostream& o, const modular_matrix& m);
625 public:
626         modular_matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
627         {
628                 m.resize(c*r, init);
629         }
630         size_t rowsize() const { return r; }
631         size_t colsize() const { return c; }
632         cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
633         cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
634         void mul_col(size_t col, const cl_MI x)
635         {
636                 mvec::iterator i = m.begin() + col;
637                 for ( size_t rc=0; rc<r; ++rc ) {
638                         *i = *i * x;
639                         i += c;
640                 }
641         }
642         void sub_col(size_t col1, size_t col2, const cl_MI fac)
643         {
644                 mvec::iterator i1 = m.begin() + col1;
645                 mvec::iterator i2 = m.begin() + col2;
646                 for ( size_t rc=0; rc<r; ++rc ) {
647                         *i1 = *i1 - *i2 * fac;
648                         i1 += c;
649                         i2 += c;
650                 }
651         }
652         void switch_col(size_t col1, size_t col2)
653         {
654                 cl_MI buf;
655                 mvec::iterator i1 = m.begin() + col1;
656                 mvec::iterator i2 = m.begin() + col2;
657                 for ( size_t rc=0; rc<r; ++rc ) {
658                         buf = *i1; *i1 = *i2; *i2 = buf;
659                         i1 += c;
660                         i2 += c;
661                 }
662         }
663         void mul_row(size_t row, const cl_MI x)
664         {
665                 vector<cl_MI>::iterator i = m.begin() + row*c;
666                 for ( size_t cc=0; cc<c; ++cc ) {
667                         *i = *i * x;
668                         ++i;
669                 }
670         }
671         void sub_row(size_t row1, size_t row2, const cl_MI fac)
672         {
673                 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
674                 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
675                 for ( size_t cc=0; cc<c; ++cc ) {
676                         *i1 = *i1 - *i2 * fac;
677                         ++i1;
678                         ++i2;
679                 }
680         }
681         void switch_row(size_t row1, size_t row2)
682         {
683                 cl_MI buf;
684                 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
685                 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
686                 for ( size_t cc=0; cc<c; ++cc ) {
687                         buf = *i1; *i1 = *i2; *i2 = buf;
688                         ++i1;
689                         ++i2;
690                 }
691         }
692         bool is_col_zero(size_t col) const
693         {
694                 mvec::const_iterator i = m.begin() + col;
695                 for ( size_t rr=0; rr<r; ++rr ) {
696                         if ( !zerop(*i) ) {
697                                 return false;
698                         }
699                         i += c;
700                 }
701                 return true;
702         }
703         bool is_row_zero(size_t row) const
704         {
705                 mvec::const_iterator i = m.begin() + row*c;
706                 for ( size_t cc=0; cc<c; ++cc ) {
707                         if ( !zerop(*i) ) {
708                                 return false;
709                         }
710                         ++i;
711                 }
712                 return true;
713         }
714         void set_row(size_t row, const vector<cl_MI>& newrow)
715         {
716                 mvec::iterator i1 = m.begin() + row*c;
717                 mvec::const_iterator i2 = newrow.begin(), end = newrow.end();
718                 for ( ; i2 != end; ++i1, ++i2 ) {
719                         *i1 = *i2;
720                 }
721         }
722         mvec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
723         mvec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
724 private:
725         size_t r, c;
726         mvec m;
727 };
728
729 #ifdef DEBUGFACTOR
730 modular_matrix operator*(const modular_matrix& m1, const modular_matrix& m2)
731 {
732         const unsigned int r = m1.rowsize();
733         const unsigned int c = m2.colsize();
734         modular_matrix o(r,c,m1(0,0));
735
736         for ( size_t i=0; i<r; ++i ) {
737                 for ( size_t j=0; j<c; ++j ) {
738                         cl_MI buf;
739                         buf = m1(i,0) * m2(0,j);
740                         for ( size_t k=1; k<c; ++k ) {
741                                 buf = buf + m1(i,k)*m2(k,j);
742                         }
743                         o(i,j) = buf;
744                 }
745         }
746         return o;
747 }
748
749 ostream& operator<<(ostream& o, const modular_matrix& m)
750 {
751         cl_modint_ring R = m(0,0).ring();
752         o << "{";
753         for ( size_t i=0; i<m.rowsize(); ++i ) {
754                 o << "{";
755                 for ( size_t j=0; j<m.colsize()-1; ++j ) {
756                         o << R->retract(m(i,j)) << ",";
757                 }
758                 o << R->retract(m(i,m.colsize()-1)) << "}";
759                 if ( i != m.rowsize()-1 ) {
760                         o << ",";
761                 }
762         }
763         o << "}";
764         return o;
765 }
766 #endif // def DEBUGFACTOR
767
768 // END modular matrix
769 ////////////////////////////////////////////////////////////////////////////////
770
771 static void q_matrix(const umodpoly& a_, modular_matrix& Q)
772 {
773         umodpoly a = a_;
774         normalize_in_field(a);
775
776         int n = degree(a);
777         unsigned int q = cl_I_to_uint(a[0].ring()->modulus);
778         umodpoly r(n, a[0].ring()->zero());
779         r[0] = a[0].ring()->one();
780         Q.set_row(0, r);
781         unsigned int max = (n-1) * q;
782         for ( size_t m=1; m<=max; ++m ) {
783                 cl_MI rn_1 = r.back();
784                 for ( size_t i=n-1; i>0; --i ) {
785                         r[i] = r[i-1] - (rn_1 * a[i]);
786                 }
787                 r[0] = -rn_1 * a[0];
788                 if ( (m % q) == 0 ) {
789                         Q.set_row(m/q, r);
790                 }
791         }
792 }
793
794 static void nullspace(modular_matrix& M, vector<mvec>& basis)
795 {
796         const size_t n = M.rowsize();
797         const cl_MI one = M(0,0).ring()->one();
798         for ( size_t i=0; i<n; ++i ) {
799                 M(i,i) = M(i,i) - one;
800         }
801         for ( size_t r=0; r<n; ++r ) {
802                 size_t cc = 0;
803                 for ( ; cc<n; ++cc ) {
804                         if ( !zerop(M(r,cc)) ) {
805                                 if ( cc < r ) {
806                                         if ( !zerop(M(cc,cc)) ) {
807                                                 continue;
808                                         }
809                                         M.switch_col(cc, r);
810                                 }
811                                 else if ( cc > r ) {
812                                         M.switch_col(cc, r);
813                                 }
814                                 break;
815                         }
816                 }
817                 if ( cc < n ) {
818                         M.mul_col(r, recip(M(r,r)));
819                         for ( cc=0; cc<n; ++cc ) {
820                                 if ( cc != r ) {
821                                         M.sub_col(cc, r, M(r,cc));
822                                 }
823                         }
824                 }
825         }
826
827         for ( size_t i=0; i<n; ++i ) {
828                 M(i,i) = M(i,i) - one;
829         }
830         for ( size_t i=0; i<n; ++i ) {
831                 if ( !M.is_row_zero(i) ) {
832                         mvec nu(M.row_begin(i), M.row_end(i));
833                         basis.push_back(nu);
834                 }
835         }
836 }
837
838 static void berlekamp(const umodpoly& a, upvec& upv)
839 {
840         cl_modint_ring R = a[0].ring();
841         umodpoly one(1, R->one());
842
843         modular_matrix Q(degree(a), degree(a), R->zero());
844         q_matrix(a, Q);
845         vector<mvec> nu;
846         nullspace(Q, nu);
847
848         const unsigned int k = nu.size();
849         if ( k == 1 ) {
850                 return;
851         }
852
853         list<umodpoly> factors;
854         factors.push_back(a);
855         unsigned int size = 1;
856         unsigned int r = 1;
857         unsigned int q = cl_I_to_uint(R->modulus);
858
859         list<umodpoly>::iterator u = factors.begin();
860
861         while ( true ) {
862                 for ( unsigned int s=0; s<q; ++s ) {
863                         umodpoly nur = nu[r];
864                         nur[0] = nur[0] - cl_MI(R, s);
865                         canonicalize(nur);
866                         umodpoly g;
867                         gcd(nur, *u, g);
868                         if ( unequal_one(g) && g != *u ) {
869                                 umodpoly uo;
870                                 div(*u, g, uo);
871                                 if ( equal_one(uo) ) {
872                                         throw logic_error("berlekamp: unexpected divisor.");
873                                 }
874                                 else {
875                                         *u = uo;
876                                 }
877                                 factors.push_back(g);
878                                 size = 0;
879                                 list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
880                                 while ( i != end ) {
881                                         if ( degree(*i) ) ++size; 
882                                         ++i;
883                                 }
884                                 if ( size == k ) {
885                                         list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
886                                         while ( i != end ) {
887                                                 upv.push_back(*i++);
888                                         }
889                                         return;
890                                 }
891                         }
892                 }
893                 if ( ++r == k ) {
894                         r = 1;
895                         ++u;
896                 }
897         }
898 }
899
900 static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
901 {
902         size_t newdeg = degree(a)/prime;
903         ap.resize(newdeg+1);
904         ap[0] = a[0];
905         for ( size_t i=1; i<=newdeg; ++i ) {
906                 ap[i] = a[i*prime];
907         }
908 }
909
910 static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
911 {
912         const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus);
913         int i = 1;
914         umodpoly b;
915         deriv(a, b);
916         if ( b.size() ) {
917                 umodpoly c;
918                 gcd(a, b, c);
919                 umodpoly w;
920                 div(a, c, w);
921                 while ( unequal_one(w) ) {
922                         umodpoly y;
923                         gcd(w, c, y);
924                         umodpoly z;
925                         div(w, y, z);
926                         factors.push_back(z);
927                         mult.push_back(i);
928                         ++i;
929                         w = y;
930                         umodpoly buf;
931                         div(c, y, buf);
932                         c = buf;
933                 }
934                 if ( unequal_one(c) ) {
935                         umodpoly cp;
936                         expt_1_over_p(c, prime, cp);
937                         size_t previ = mult.size();
938                         modsqrfree(cp, factors, mult);
939                         for ( size_t i=previ; i<mult.size(); ++i ) {
940                                 mult[i] *= prime;
941                         }
942                 }
943         }
944         else {
945                 umodpoly ap;
946                 expt_1_over_p(a, prime, ap);
947                 size_t previ = mult.size();
948                 modsqrfree(ap, factors, mult);
949                 for ( size_t i=previ; i<mult.size(); ++i ) {
950                         mult[i] *= prime;
951                 }
952         }
953 }
954
955 static void distinct_degree_factor(const umodpoly& a_, vector<int>& degrees, upvec& ddfactors)
956 {
957         umodpoly a = a_;
958
959         cl_modint_ring R = a[0].ring();
960         int q = cl_I_to_int(R->modulus);
961         int nhalf = degree(a)/2;
962
963         int i = 1;
964         umodpoly w(2);
965         w[0] = R->zero();
966         w[1] = R->one();
967         umodpoly x = w;
968
969         while ( i <= nhalf ) {
970                 expt_pos(w, q);
971                 umodpoly buf;
972                 rem(w, a, buf);
973                 w = buf;
974                 umodpoly wx = w - x;
975                 gcd(a, wx, buf);
976                 if ( unequal_one(buf) ) {
977                         degrees.push_back(i);
978                         ddfactors.push_back(buf);
979                 }
980                 if ( unequal_one(buf) ) {
981                         umodpoly buf2;
982                         div(a, buf, buf2);
983                         a = buf2;
984                         nhalf = degree(a)/2;
985                         rem(w, a, buf);
986                         w = buf;
987                 }
988                 ++i;
989         }
990         if ( unequal_one(a) ) {
991                 degrees.push_back(degree(a));
992                 ddfactors.push_back(a);
993         }
994 }
995
996 static void same_degree_factor(const umodpoly& a, upvec& upv)
997 {
998         cl_modint_ring R = a[0].ring();
999
1000         vector<int> degrees;
1001         upvec ddfactors;
1002         distinct_degree_factor(a, degrees, ddfactors);
1003
1004         for ( size_t i=0; i<degrees.size(); ++i ) {
1005                 if ( degrees[i] == degree(ddfactors[i]) ) {
1006                         upv.push_back(ddfactors[i]);
1007                 }
1008                 else {
1009                         berlekamp(ddfactors[i], upv);
1010                 }
1011         }
1012 }
1013
1014 #define USE_SAME_DEGREE_FACTOR
1015
1016 static void factor_modular(const umodpoly& p, upvec& upv)
1017 {
1018 #ifdef USE_SAME_DEGREE_FACTOR
1019         same_degree_factor(p, upv);
1020 #else
1021         berlekamp(p, upv);
1022 #endif
1023 }
1024
1025 /** Calculates polynomials s and t such that a*s+b*t==1.
1026  *  Assertion: a and b are relatively prime and not zero.
1027  *
1028  *  @param[in]  a  polynomial
1029  *  @param[in]  b  polynomial
1030  *  @param[out] s  polynomial
1031  *  @param[out] t  polynomial
1032  */
1033 static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpoly& t)
1034 {
1035         if ( degree(a) < degree(b) ) {
1036                 exteuclid(b, a, t, s);
1037                 return;
1038         }
1039
1040         umodpoly one(1, a[0].ring()->one());
1041         umodpoly c = a; normalize_in_field(c);
1042         umodpoly d = b; normalize_in_field(d);
1043         s = one;
1044         t.clear();
1045         umodpoly d1;
1046         umodpoly d2 = one;
1047         umodpoly q;
1048         while ( true ) {
1049                 div(c, d, q);
1050                 umodpoly r = c - q * d;
1051                 umodpoly r1 = s - q * d1;
1052                 umodpoly r2 = t - q * d2;
1053                 c = d;
1054                 s = d1;
1055                 t = d2;
1056                 if ( r.empty() ) break;
1057                 d = r;
1058                 d1 = r1;
1059                 d2 = r2;
1060         }
1061         cl_MI fac = recip(lcoeff(a) * lcoeff(c));
1062         umodpoly::iterator i = s.begin(), end = s.end();
1063         for ( ; i!=end; ++i ) {
1064                 *i = *i * fac;
1065         }
1066         canonicalize(s);
1067         fac = recip(lcoeff(b) * lcoeff(c));
1068         i = t.begin(), end = t.end();
1069         for ( ; i!=end; ++i ) {
1070                 *i = *i * fac;
1071         }
1072         canonicalize(t);
1073 }
1074
1075 static upoly replace_lc(const upoly& poly, const cl_I& lc)
1076 {
1077         if ( poly.empty() ) return poly;
1078         upoly r = poly;
1079         r.back() = lc;
1080         return r;
1081 }
1082
1083 static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
1084 {
1085         cl_I maxcoeff = 0;
1086         cl_R coeff = 0;
1087         for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
1088                 cl_I aa = abs(the<cl_I>(ex_to<numeric>(a.coeff(x, i)).to_cl_N()));
1089                 if ( aa > maxcoeff ) maxcoeff = aa;
1090                 coeff = coeff + square(aa);
1091         }
1092         cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
1093         cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1094         return ( B > maxcoeff ) ? B : maxcoeff;
1095 }
1096
1097 static inline cl_I calc_bound(const upoly& a, int maxdeg)
1098 {
1099         cl_I maxcoeff = 0;
1100         cl_R coeff = 0;
1101         for ( int i=degree(a); i>=0; --i ) {
1102                 cl_I aa = abs(a[i]);
1103                 if ( aa > maxcoeff ) maxcoeff = aa;
1104                 coeff = coeff + square(aa);
1105         }
1106         cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
1107         cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1108         return ( B > maxcoeff ) ? B : maxcoeff;
1109 }
1110
1111 static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w)
1112 {
1113         upoly a = a_;
1114         const cl_modint_ring& R = u1_[0].ring();
1115
1116         // calc bound B
1117         int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
1118         cl_I maxmodulus = 2*calc_bound(a, maxdeg);
1119
1120         // step 1
1121         cl_I alpha = lcoeff(a);
1122         a = a * alpha;
1123         umodpoly nu1 = u1_;
1124         normalize_in_field(nu1);
1125         umodpoly nw1 = w1_;
1126         normalize_in_field(nw1);
1127         upoly phi;
1128         phi = umodpoly_to_upoly(nu1) * alpha;
1129         umodpoly u1;
1130         umodpoly_from_upoly(u1, phi, R);
1131         phi = umodpoly_to_upoly(nw1) * alpha;
1132         umodpoly w1;
1133         umodpoly_from_upoly(w1, phi, R);
1134
1135         // step 2
1136         umodpoly s;
1137         umodpoly t;
1138         exteuclid(u1, w1, s, t);
1139
1140         // step 3
1141         u = replace_lc(umodpoly_to_upoly(u1), alpha);
1142         w = replace_lc(umodpoly_to_upoly(w1), alpha);
1143         upoly e = a - u * w;
1144         cl_I modulus = p;
1145
1146         // step 4
1147         while ( !e.empty() && modulus < maxmodulus ) {
1148                 upoly c = e / modulus;
1149                 phi = umodpoly_to_upoly(s) * c;
1150                 umodpoly sigmatilde;
1151                 umodpoly_from_upoly(sigmatilde, phi, R);
1152                 phi = umodpoly_to_upoly(t) * c;
1153                 umodpoly tautilde;
1154                 umodpoly_from_upoly(tautilde, phi, R);
1155                 umodpoly r, q;
1156                 remdiv(sigmatilde, w1, r, q);
1157                 umodpoly sigma = r;
1158                 phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1);
1159                 umodpoly tau;
1160                 umodpoly_from_upoly(tau, phi, R);
1161                 u = u + umodpoly_to_upoly(tau) * modulus;
1162                 w = w + umodpoly_to_upoly(sigma) * modulus;
1163                 e = a - u * w;
1164                 modulus = modulus * p;
1165         }
1166
1167         // step 5
1168         if ( e.empty() ) {
1169                 cl_I g = u[0];
1170                 for ( size_t i=1; i<u.size(); ++i ) {
1171                         g = gcd(g, u[i]);
1172                         if ( g == 1 ) break;
1173                 }
1174                 if ( g != 1 ) {
1175                         u = u / g;
1176                         w = w * g;
1177                 }
1178                 if ( alpha != 1 ) {
1179                         w = w / alpha;
1180                 }
1181         }
1182         else {
1183                 u.clear();
1184         }
1185 }
1186
1187 static unsigned int next_prime(unsigned int p)
1188 {
1189         static vector<unsigned int> primes;
1190         if ( primes.size() == 0 ) {
1191                 primes.push_back(3); primes.push_back(5); primes.push_back(7);
1192         }
1193         vector<unsigned int>::const_iterator it = primes.begin();
1194         if ( p >= primes.back() ) {
1195                 unsigned int candidate = primes.back() + 2;
1196                 while ( true ) {
1197                         size_t n = primes.size()/2;
1198                         for ( size_t i=0; i<n; ++i ) {
1199                                 if ( candidate % primes[i] ) continue;
1200                                 candidate += 2;
1201                                 i=-1;
1202                         }
1203                         primes.push_back(candidate);
1204                         if ( candidate > p ) break;
1205                 }
1206                 return candidate;
1207         }
1208         vector<unsigned int>::const_iterator end = primes.end();
1209         for ( ; it!=end; ++it ) {
1210                 if ( *it > p ) {
1211                         return *it;
1212                 }
1213         }
1214         throw logic_error("next_prime: should not reach this point!");
1215 }
1216
1217 class factor_partition
1218 {
1219 public:
1220         factor_partition(const upvec& factors_) : factors(factors_)
1221         {
1222                 n = factors.size();
1223                 k.resize(n, 0);
1224                 k[0] = 1;
1225                 cache.resize(n-1);
1226                 one.resize(1, factors.front()[0].ring()->one());
1227                 len = 1;
1228                 last = 0;
1229                 split();
1230         }
1231         int operator[](size_t i) const { return k[i]; }
1232         size_t size() const { return n; }
1233         size_t size_left() const { return n-len; }
1234         size_t size_right() const { return len; }
1235 #ifdef DEBUGFACTOR
1236         void get() const { DCOUTVAR(k); }
1237 #endif
1238         bool next()
1239         {
1240                 if ( last == n-1 ) {
1241                         int rem = len - 1;
1242                         int p = last - 1;
1243                         while ( rem ) {
1244                                 if ( k[p] ) {
1245                                         --rem;
1246                                         --p;
1247                                         continue;
1248                                 }
1249                                 last = p - 1;
1250                                 while ( k[last] == 0 ) { --last; }
1251                                 if ( last == 0 && n == 2*len ) return false;
1252                                 k[last++] = 0;
1253                                 for ( size_t i=0; i<=len-rem; ++i ) {
1254                                         k[last] = 1;
1255                                         ++last;
1256                                 }
1257                                 fill(k.begin()+last, k.end(), 0);
1258                                 --last;
1259                                 split();
1260                                 return true;
1261                         }
1262                         last = len;
1263                         ++len;
1264                         if ( len > n/2 ) return false;
1265                         fill(k.begin(), k.begin()+len, 1);
1266                         fill(k.begin()+len+1, k.end(), 0);
1267                 }
1268                 else {
1269                         k[last++] = 0;
1270                         k[last] = 1;
1271                 }
1272                 split();
1273                 return true;
1274         }
1275         umodpoly& left() { return lr[0]; }
1276         umodpoly& right() { return lr[1]; }
1277 private:
1278         void split_cached()
1279         {
1280                 size_t i = 0;
1281                 do {
1282                         size_t pos = i;
1283                         int group = k[i++];
1284                         size_t d = 0;
1285                         while ( i < n && k[i] == group ) { ++d; ++i; }
1286                         if ( d ) {
1287                                 if ( cache[pos].size() >= d ) {
1288                                         lr[group] = lr[group] * cache[pos][d-1];
1289                                 }
1290                                 else {
1291                                         if ( cache[pos].size() == 0 ) {
1292                                                 cache[pos].push_back(factors[pos] * factors[pos+1]);
1293                                         }
1294                                         size_t j = pos + cache[pos].size() + 1;
1295                                         d -= cache[pos].size();
1296                                         while ( d ) {
1297                                                 umodpoly buf = cache[pos].back() * factors[j];
1298                                                 cache[pos].push_back(buf);
1299                                                 --d;
1300                                                 ++j;
1301                                         }
1302                                         lr[group] = lr[group] * cache[pos].back();
1303                                 }
1304                         }
1305                         else {
1306                                 lr[group] = lr[group] * factors[pos];
1307                         }
1308                 } while ( i < n );
1309         }
1310         void split()
1311         {
1312                 lr[0] = one;
1313                 lr[1] = one;
1314                 if ( n > 6 ) {
1315                         split_cached();
1316                 }
1317                 else {
1318                         for ( size_t i=0; i<n; ++i ) {
1319                                 lr[k[i]] = lr[k[i]] * factors[i];
1320                         }
1321                 }
1322         }
1323 private:
1324         umodpoly lr[2];
1325         vector< vector<umodpoly> > cache;
1326         upvec factors;
1327         umodpoly one;
1328         size_t n;
1329         size_t len;
1330         size_t last;
1331         vector<int> k;
1332 };
1333
1334 struct ModFactors
1335 {
1336         upoly poly;
1337         upvec factors;
1338 };
1339
1340 static ex factor_univariate(const ex& poly, const ex& x)
1341 {
1342         ex unit, cont, prim_ex;
1343         poly.unitcontprim(x, unit, cont, prim_ex);
1344         upoly prim;
1345         upoly_from_ex(prim, prim_ex, x);
1346
1347         // determine proper prime and minimize number of modular factors
1348         unsigned int p = 3, lastp = 3;
1349         cl_modint_ring R;
1350         unsigned int trials = 0;
1351         unsigned int minfactors = 0;
1352         cl_I lc = lcoeff(prim);
1353         upvec factors;
1354         while ( trials < 2 ) {
1355                 umodpoly modpoly;
1356                 while ( true ) {
1357                         p = next_prime(p);
1358                         if ( !zerop(rem(lc, p)) ) {
1359                                 R = find_modint_ring(p);
1360                                 umodpoly_from_upoly(modpoly, prim, R);
1361                                 if ( squarefree(modpoly) ) break;
1362                         }
1363                 }
1364
1365                 // do modular factorization
1366                 upvec trialfactors;
1367                 factor_modular(modpoly, trialfactors);
1368                 if ( trialfactors.size() <= 1 ) {
1369                         // irreducible for sure
1370                         return poly;
1371                 }
1372
1373                 if ( minfactors == 0 || trialfactors.size() < minfactors ) {
1374                         factors = trialfactors;
1375                         minfactors = trialfactors.size();
1376                         lastp = p;
1377                         trials = 1;
1378                 }
1379                 else {
1380                         ++trials;
1381                 }
1382         }
1383         p = lastp;
1384         R = find_modint_ring(p);
1385
1386         // lift all factor combinations
1387         stack<ModFactors> tocheck;
1388         ModFactors mf;
1389         mf.poly = prim;
1390         mf.factors = factors;
1391         tocheck.push(mf);
1392         upoly f1, f2;
1393         ex result = 1;
1394         while ( tocheck.size() ) {
1395                 const size_t n = tocheck.top().factors.size();
1396                 factor_partition part(tocheck.top().factors);
1397                 while ( true ) {
1398                         hensel_univar(tocheck.top().poly, p, part.left(), part.right(), f1, f2);
1399                         if ( !f1.empty() ) {
1400                                 if ( part.size_left() == 1 ) {
1401                                         if ( part.size_right() == 1 ) {
1402                                                 result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1403                                                 tocheck.pop();
1404                                                 break;
1405                                         }
1406                                         result *= upoly_to_ex(f1, x);
1407                                         tocheck.top().poly = f2;
1408                                         for ( size_t i=0; i<n; ++i ) {
1409                                                 if ( part[i] == 0 ) {
1410                                                         tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1411                                                         break;
1412                                                 }
1413                                         }
1414                                         break;
1415                                 }
1416                                 else if ( part.size_right() == 1 ) {
1417                                         if ( part.size_left() == 1 ) {
1418                                                 result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1419                                                 tocheck.pop();
1420                                                 break;
1421                                         }
1422                                         result *= upoly_to_ex(f2, x);
1423                                         tocheck.top().poly = f1;
1424                                         for ( size_t i=0; i<n; ++i ) {
1425                                                 if ( part[i] == 1 ) {
1426                                                         tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1427                                                         break;
1428                                                 }
1429                                         }
1430                                         break;
1431                                 }
1432                                 else {
1433                                         upvec newfactors1(part.size_left()), newfactors2(part.size_right());
1434                                         upvec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1435                                         for ( size_t i=0; i<n; ++i ) {
1436                                                 if ( part[i] ) {
1437                                                         *i2++ = tocheck.top().factors[i];
1438                                                 }
1439                                                 else {
1440                                                         *i1++ = tocheck.top().factors[i];
1441                                                 }
1442                                         }
1443                                         tocheck.top().factors = newfactors1;
1444                                         tocheck.top().poly = f1;
1445                                         ModFactors mf;
1446                                         mf.factors = newfactors2;
1447                                         mf.poly = f2;
1448                                         tocheck.push(mf);
1449                                         break;
1450                                 }
1451                         }
1452                         else {
1453                                 if ( !part.next() ) {
1454                                         result *= upoly_to_ex(tocheck.top().poly, x);
1455                                         tocheck.pop();
1456                                         break;
1457                                 }
1458                         }
1459                 }
1460         }
1461
1462         return unit * cont * result;
1463 }
1464
1465 struct EvalPoint
1466 {
1467         ex x;
1468         int evalpoint;
1469 };
1470
1471 // forward declaration
1472 vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex& c, const vector<EvalPoint>& I, unsigned int d, unsigned int p, unsigned int k);
1473
1474 upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k)
1475 {
1476         const size_t r = a.size();
1477         cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1478         upvec q(r-1);
1479         q[r-2] = a[r-1];
1480         for ( size_t j=r-2; j>=1; --j ) {
1481                 q[j-1] = a[j] * q[j];
1482         }
1483         umodpoly beta(1, R->one());
1484         upvec s;
1485         for ( size_t j=1; j<r; ++j ) {
1486                 vector<ex> mdarg(2);
1487                 mdarg[0] = umodpoly_to_ex(q[j-1], x);
1488                 mdarg[1] = umodpoly_to_ex(a[j-1], x);
1489                 vector<EvalPoint> empty;
1490                 vector<ex> exsigma = multivar_diophant(mdarg, x, umodpoly_to_ex(beta, x), empty, 0, p, k);
1491                 umodpoly sigma1;
1492                 umodpoly_from_ex(sigma1, exsigma[0], x, R);
1493                 umodpoly sigma2;
1494                 umodpoly_from_ex(sigma2, exsigma[1], x, R);
1495                 beta = sigma1;
1496                 s.push_back(sigma2);
1497         }
1498         s.push_back(beta);
1499         return s;
1500 }
1501
1502 /**
1503  *  Assert: a not empty.
1504  */
1505 void change_modulus(const cl_modint_ring& R, umodpoly& a)
1506 {
1507         if ( a.empty() ) return;
1508         cl_modint_ring oldR = a[0].ring();
1509         umodpoly::iterator i = a.begin(), end = a.end();
1510         for ( ; i!=end; ++i ) {
1511                 *i = R->canonhom(oldR->retract(*i));
1512         }
1513         canonicalize(a);
1514 }
1515
1516 void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_)
1517 {
1518         cl_modint_ring R = find_modint_ring(p);
1519         umodpoly amod = a;
1520         change_modulus(R, amod);
1521         umodpoly bmod = b;
1522         change_modulus(R, bmod);
1523
1524         umodpoly smod;
1525         umodpoly tmod;
1526         exteuclid(amod, bmod, smod, tmod);
1527
1528         cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1529         umodpoly s = smod;
1530         change_modulus(Rpk, s);
1531         umodpoly t = tmod;
1532         change_modulus(Rpk, t);
1533
1534         cl_I modulus(p);
1535         umodpoly one(1, Rpk->one());
1536         for ( size_t j=1; j<k; ++j ) {
1537                 umodpoly e = one - a * s - b * t;
1538                 reduce_coeff(e, modulus);
1539                 umodpoly c = e;
1540                 change_modulus(R, c);
1541                 umodpoly sigmabar = smod * c;
1542                 umodpoly taubar = tmod * c;
1543                 umodpoly sigma, q;
1544                 remdiv(sigmabar, bmod, sigma, q);
1545                 umodpoly tau = taubar + q * amod;
1546                 umodpoly sadd = sigma;
1547                 change_modulus(Rpk, sadd);
1548                 cl_MI modmodulus(Rpk, modulus);
1549                 s = s + sadd * modmodulus;
1550                 umodpoly tadd = tau;
1551                 change_modulus(Rpk, tadd);
1552                 t = t + tadd * modmodulus;
1553                 modulus = modulus * p;
1554         }
1555
1556         s_ = s; t_ = t;
1557 }
1558
1559 upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1560 {
1561         cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1562
1563         const size_t r = a.size();
1564         upvec result;
1565         if ( r > 2 ) {
1566                 upvec s = multiterm_eea_lift(a, x, p, k);
1567                 for ( size_t j=0; j<r; ++j ) {
1568                         umodpoly bmod = umodpoly_to_umodpoly(s[j], R, m);
1569                         umodpoly buf;
1570                         rem(bmod, a[j], buf);
1571                         result.push_back(buf);
1572                 }
1573         }
1574         else {
1575                 umodpoly s, t;
1576                 eea_lift(a[1], a[0], x, p, k, s, t);
1577                 umodpoly bmod = umodpoly_to_umodpoly(s, R, m);
1578                 umodpoly buf, q;
1579                 remdiv(bmod, a[0], buf, q);
1580                 result.push_back(buf);
1581                 umodpoly t1mod = umodpoly_to_umodpoly(t, R, m);
1582                 buf = t1mod + q * a[1];
1583                 result.push_back(buf);
1584         }
1585
1586         return result;
1587 }
1588
1589 struct make_modular_map : public map_function {
1590         cl_modint_ring R;
1591         make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1592         ex operator()(const ex& e)
1593         {
1594                 if ( is_a<add>(e) || is_a<mul>(e) ) {
1595                         return e.map(*this);
1596                 }
1597                 else if ( is_a<numeric>(e) ) {
1598                         numeric mod(R->modulus);
1599                         numeric halfmod = (mod-1)/2;
1600                         cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1601                         numeric n(R->retract(emod));
1602                         if ( n > halfmod ) {
1603                                 return n-mod;
1604                         }
1605                         else {
1606                                 return n;
1607                         }
1608                 }
1609                 return e;
1610         }
1611 };
1612
1613 static ex make_modular(const ex& e, const cl_modint_ring& R)
1614 {
1615         make_modular_map map(R);
1616         return map(e.expand());
1617 }
1618
1619 vector<ex> multivar_diophant(const vector<ex>& a_, const ex& x, const ex& c, const vector<EvalPoint>& I, unsigned int d, unsigned int p, unsigned int k)
1620 {
1621         vector<ex> a = a_;
1622
1623         const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1624         const size_t r = a.size();
1625         const size_t nu = I.size() + 1;
1626
1627         vector<ex> sigma;
1628         if ( nu > 1 ) {
1629                 ex xnu = I.back().x;
1630                 int alphanu = I.back().evalpoint;
1631
1632                 ex A = 1;
1633                 for ( size_t i=0; i<r; ++i ) {
1634                         A *= a[i];
1635                 }
1636                 vector<ex> b(r);
1637                 for ( size_t i=0; i<r; ++i ) {
1638                         b[i] = normal(A / a[i]);
1639                 }
1640
1641                 vector<ex> anew = a;
1642                 for ( size_t i=0; i<r; ++i ) {
1643                         anew[i] = anew[i].subs(xnu == alphanu);
1644                 }
1645                 ex cnew = c.subs(xnu == alphanu);
1646                 vector<EvalPoint> Inew = I;
1647                 Inew.pop_back();
1648                 sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1649
1650                 ex buf = c;
1651                 for ( size_t i=0; i<r; ++i ) {
1652                         buf -= sigma[i] * b[i];
1653                 }
1654                 ex e = make_modular(buf, R);
1655
1656                 ex monomial = 1;
1657                 for ( size_t m=1; !e.is_zero() && e.has(xnu) && m<=d; ++m ) {
1658                         monomial *= (xnu - alphanu);
1659                         monomial = expand(monomial);
1660                         ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1661                         cm = make_modular(cm, R);
1662                         if ( !cm.is_zero() ) {
1663                                 vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1664                                 ex buf = e;
1665                                 for ( size_t j=0; j<delta_s.size(); ++j ) {
1666                                         delta_s[j] *= monomial;
1667                                         sigma[j] += delta_s[j];
1668                                         buf -= delta_s[j] * b[j];
1669                                 }
1670                                 e = make_modular(buf, R);
1671                         }
1672                 }
1673         }
1674         else {
1675                 upvec amod;
1676                 for ( size_t i=0; i<a.size(); ++i ) {
1677                         umodpoly up;
1678                         umodpoly_from_ex(up, a[i], x, R);
1679                         amod.push_back(up);
1680                 }
1681
1682                 sigma.insert(sigma.begin(), r, 0);
1683                 size_t nterms;
1684                 ex z;
1685                 if ( is_a<add>(c) ) {
1686                         nterms = c.nops();
1687                         z = c.op(0);
1688                 }
1689                 else {
1690                         nterms = 1;
1691                         z = c;
1692                 }
1693                 for ( size_t i=0; i<nterms; ++i ) {
1694                         int m = z.degree(x);
1695                         cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1696                         upvec delta_s = univar_diophant(amod, x, m, p, k);
1697                         cl_MI modcm;
1698                         cl_I poscm = cm;
1699                         while ( poscm < 0 ) {
1700                                 poscm = poscm + expt_pos(cl_I(p),k);
1701                         }
1702                         modcm = cl_MI(R, poscm);
1703                         for ( size_t j=0; j<delta_s.size(); ++j ) {
1704                                 delta_s[j] = delta_s[j] * modcm;
1705                                 sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j], x);
1706                         }
1707                         if ( nterms > 1 ) {
1708                                 z = c.op(i+1);
1709                         }
1710                 }
1711         }
1712
1713         for ( size_t i=0; i<sigma.size(); ++i ) {
1714                 sigma[i] = make_modular(sigma[i], R);
1715         }
1716
1717         return sigma;
1718 }
1719
1720 #ifdef DEBUGFACTOR
1721 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1722 {
1723         for ( size_t i=0; i<v.size(); ++i ) {
1724                 o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1725         }
1726         return o;
1727 }
1728 #endif // def DEBUGFACTOR
1729
1730 ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I, unsigned int p, const cl_I& l, const upvec& u, const vector<ex>& lcU)
1731 {
1732         const size_t nu = I.size() + 1;
1733         const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1734
1735         vector<ex> A(nu);
1736         A[nu-1] = a;
1737
1738         for ( size_t j=nu; j>=2; --j ) {
1739                 ex x = I[j-2].x;
1740                 int alpha = I[j-2].evalpoint;
1741                 A[j-2] = A[j-1].subs(x==alpha);
1742                 A[j-2] = make_modular(A[j-2], R);
1743         }
1744
1745         int maxdeg = a.degree(I.front().x);
1746         for ( size_t i=1; i<I.size(); ++i ) {
1747                 int maxdeg2 = a.degree(I[i].x);
1748                 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
1749         }
1750
1751         const size_t n = u.size();
1752         vector<ex> U(n);
1753         for ( size_t i=0; i<n; ++i ) {
1754                 U[i] = umodpoly_to_ex(u[i], x);
1755         }
1756
1757         for ( size_t j=2; j<=nu; ++j ) {
1758                 vector<ex> U1 = U;
1759                 ex monomial = 1;
1760                 for ( size_t m=0; m<n; ++m) {
1761                         if ( lcU[m] != 1 ) {
1762                                 ex coef = lcU[m];
1763                                 for ( size_t i=j-1; i<nu-1; ++i ) {
1764                                         coef = coef.subs(I[i].x == I[i].evalpoint);
1765                                 }
1766                                 coef = make_modular(coef, R);
1767                                 int deg = U[m].degree(x);
1768                                 U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
1769                         }
1770                 }
1771                 ex Uprod = 1;
1772                 for ( size_t i=0; i<n; ++i ) {
1773                         Uprod *= U[i];
1774                 }
1775                 ex e = expand(A[j-1] - Uprod);
1776
1777                 vector<EvalPoint> newI;
1778                 for ( size_t i=1; i<=j-2; ++i ) {
1779                         newI.push_back(I[i-1]);
1780                 }
1781
1782                 ex xj = I[j-2].x;
1783                 int alphaj = I[j-2].evalpoint;
1784                 size_t deg = A[j-1].degree(xj);
1785                 for ( size_t k=1; k<=deg; ++k ) {
1786                         if ( !e.is_zero() ) {
1787                                 monomial *= (xj - alphaj);
1788                                 monomial = expand(monomial);
1789                                 ex dif = e.diff(ex_to<symbol>(xj), k);
1790                                 ex c = dif.subs(xj==alphaj) / factorial(k);
1791                                 if ( !c.is_zero() ) {
1792                                         vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
1793                                         for ( size_t i=0; i<n; ++i ) {
1794                                                 deltaU[i] *= monomial;
1795                                                 U[i] += deltaU[i];
1796                                                 U[i] = make_modular(U[i], R);
1797                                         }
1798                                         ex Uprod = 1;
1799                                         for ( size_t i=0; i<n; ++i ) {
1800                                                 Uprod *= U[i];
1801                                         }
1802                                         e = A[j-1] - Uprod;
1803                                         e = make_modular(e, R);
1804                                 }
1805                         }
1806                 }
1807         }
1808
1809         ex acand = 1;
1810         for ( size_t i=0; i<U.size(); ++i ) {
1811                 acand *= U[i];
1812         }
1813         if ( expand(a-acand).is_zero() ) {
1814                 lst res;
1815                 for ( size_t i=0; i<U.size(); ++i ) {
1816                         res.append(U[i]);
1817                 }
1818                 return res;
1819         }
1820         else {
1821                 lst res;
1822                 return lst();
1823         }
1824 }
1825
1826 static ex put_factors_into_lst(const ex& e)
1827 {
1828         lst result;
1829
1830         if ( is_a<numeric>(e) ) {
1831                 result.append(e);
1832                 return result;
1833         }
1834         if ( is_a<power>(e) ) {
1835                 result.append(1);
1836                 result.append(e.op(0));
1837                 result.append(e.op(1));
1838                 return result;
1839         }
1840         if ( is_a<symbol>(e) || is_a<add>(e) ) {
1841                 result.append(1);
1842                 result.append(e);
1843                 result.append(1);
1844                 return result;
1845         }
1846         if ( is_a<mul>(e) ) {
1847                 ex nfac = 1;
1848                 for ( size_t i=0; i<e.nops(); ++i ) {
1849                         ex op = e.op(i);
1850                         if ( is_a<numeric>(op) ) {
1851                                 nfac = op;
1852                         }
1853                         if ( is_a<power>(op) ) {
1854                                 result.append(op.op(0));
1855                                 result.append(op.op(1));
1856                         }
1857                         if ( is_a<symbol>(op) || is_a<add>(op) ) {
1858                                 result.append(op);
1859                                 result.append(1);
1860                         }
1861                 }
1862                 result.prepend(nfac);
1863                 return result;
1864         }
1865         throw runtime_error("put_factors_into_lst: bad term.");
1866 }
1867
1868 #ifdef DEBUGFACTOR
1869 ostream& operator<<(ostream& o, const vector<numeric>& v)
1870 {
1871         for ( size_t i=0; i<v.size(); ++i ) {
1872                 o << v[i] << " ";
1873         }
1874         return o;
1875 }
1876 #endif // def DEBUGFACTOR
1877
1878 static bool checkdivisors(const lst& f, vector<numeric>& d)
1879 {
1880         const int k = f.nops()-2;
1881         numeric q, r;
1882         d[0] = ex_to<numeric>(f.op(0) * f.op(f.nops()-1));
1883         if ( d[0] == 1 && k == 1 && abs(f.op(1)) != 1 ) {
1884                 return false;
1885         }
1886         for ( int i=1; i<=k; ++i ) {
1887                 q = ex_to<numeric>(abs(f.op(i)));
1888                 for ( int j=i-1; j>=0; --j ) {
1889                         r = d[j];
1890                         do {
1891                                 r = gcd(r, q);
1892                                 q = q/r;
1893                         } while ( r != 1 );
1894                         if ( q == 1 ) {
1895                                 return true;
1896                         }
1897                 }
1898                 d[i] = q;
1899         }
1900         return false;
1901 }
1902
1903 static bool generate_set(const ex& u, const ex& vn, const exset& syms, const ex& f, const numeric& modulus, vector<numeric>& a, vector<numeric>& d)
1904 {
1905         // computation of d is actually not necessary
1906         const ex& x = *syms.begin();
1907         bool trying = true;
1908         do {
1909                 ex u0 = u;
1910                 ex vna = vn;
1911                 ex vnatry;
1912                 exset::const_iterator s = syms.begin();
1913                 ++s;
1914                 for ( size_t i=0; i<a.size(); ++i ) {
1915                         do {
1916                                 a[i] = mod(numeric(rand()), 2*modulus) - modulus;
1917                                 vnatry = vna.subs(*s == a[i]);
1918                         } while ( vnatry == 0 );
1919                         vna = vnatry;
1920                         u0 = u0.subs(*s == a[i]);
1921                         ++s;
1922                 }
1923                 if ( gcd(u0,u0.diff(ex_to<symbol>(x))) != 1 ) {
1924                         continue;
1925                 }
1926                 if ( is_a<numeric>(vn) ) {
1927                         trying = false;
1928                 }
1929                 else {
1930                         lst fnum;
1931                         lst::const_iterator i = ex_to<lst>(f).begin();
1932                         fnum.append(*i++);
1933                         bool problem = false;
1934                         while ( i!=ex_to<lst>(f).end() ) {
1935                                 ex fs = *i;
1936                                 if ( !is_a<numeric>(fs) ) {
1937                                         s = syms.begin();
1938                                         ++s;
1939                                         for ( size_t j=0; j<a.size(); ++j ) {
1940                                                 fs = fs.subs(*s == a[j]);
1941                                                 ++s;
1942                                         }
1943                                         if ( abs(fs) == 1 ) {
1944                                                 problem = true;
1945                                                 break;
1946                                         }
1947                                 }
1948                                 fnum.append(fs);
1949                                 ++i; ++i;
1950                         }
1951                         if ( problem ) {
1952                                 return true;
1953                         }
1954                         ex con = u0.content(x);
1955                         fnum.append(con);
1956                         trying = checkdivisors(fnum, d);
1957                 }
1958         } while ( trying );
1959         return false;
1960 }
1961
1962 static ex factor_multivariate(const ex& poly, const exset& syms)
1963 {
1964         exset::const_iterator s;
1965         const ex& x = *syms.begin();
1966
1967         /* make polynomial primitive */
1968         ex p = poly.expand().collect(x);
1969         ex cont = p.lcoeff(x);
1970         for ( numeric i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
1971                 cont = gcd(cont, p.coeff(x,ex_to<numeric>(i).to_int()));
1972                 if ( cont == 1 ) break;
1973         }
1974         ex pp = expand(normal(p / cont));
1975         if ( !is_a<numeric>(cont) ) {
1976                 return factor(cont) * factor(pp);
1977         }
1978
1979         /* factor leading coefficient */
1980         pp = pp.collect(x);
1981         ex vn = pp.lcoeff(x);
1982         pp = pp.expand();
1983         ex vnlst;
1984         if ( is_a<numeric>(vn) ) {
1985                 vnlst = lst(vn);
1986         }
1987         else {
1988                 ex vnfactors = factor(vn);
1989                 vnlst = put_factors_into_lst(vnfactors);
1990         }
1991
1992         const numeric maxtrials = 3;
1993         numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3;
1994         numeric minimalr = -1;
1995         vector<numeric> a(syms.size()-1, 0);
1996         vector<numeric> d((vnlst.nops()-1)/2+1, 0);
1997
1998         while ( true ) {
1999                 numeric trialcount = 0;
2000                 ex u, delta;
2001                 unsigned int prime = 3;
2002                 size_t factor_count = 0;
2003                 ex ufac;
2004                 ex ufaclst;
2005                 while ( trialcount < maxtrials ) {
2006                         bool problem = generate_set(pp, vn, syms, vnlst, modulus, a, d);
2007                         if ( problem ) {
2008                                 ++modulus;
2009                                 continue;
2010                         }
2011                         u = pp;
2012                         s = syms.begin();
2013                         ++s;
2014                         for ( size_t i=0; i<a.size(); ++i ) {
2015                                 u = u.subs(*s == a[i]);
2016                                 ++s;
2017                         }
2018                         delta = u.content(x);
2019
2020                         // determine proper prime
2021                         prime = 3;
2022                         cl_modint_ring R = find_modint_ring(prime);
2023                         while ( true ) {
2024                                 if ( irem(ex_to<numeric>(u.lcoeff(x)), prime) != 0 ) {
2025                                         umodpoly modpoly;
2026                                         umodpoly_from_ex(modpoly, u, x, R);
2027                                         if ( squarefree(modpoly) ) break;
2028                                 }
2029                                 prime = next_prime(prime);
2030                                 R = find_modint_ring(prime);
2031                         }
2032
2033                         ufac = factor(u);
2034                         ufaclst = put_factors_into_lst(ufac);
2035                         factor_count = (ufaclst.nops()-1)/2;
2036
2037                         // veto factorization for which gcd(u_i, u_j) != 1 for all i,j
2038                         upvec tryu;
2039                         for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) {
2040                                 umodpoly newu;
2041                                 umodpoly_from_ex(newu, ufaclst.op(i*2+1), x, R);
2042                                 tryu.push_back(newu);
2043                         }
2044                         bool veto = false;
2045                         for ( size_t i=0; i<tryu.size()-1; ++i ) {
2046                                 for ( size_t j=i+1; j<tryu.size(); ++j ) {
2047                                         umodpoly tryg;
2048                                         gcd(tryu[i], tryu[j], tryg);
2049                                         if ( unequal_one(tryg) ) {
2050                                                 veto = true;
2051                                                 goto escape_quickly;
2052                                         }
2053                                 }
2054                         }
2055                         escape_quickly: ;
2056                         if ( veto ) {
2057                                 continue;
2058                         }
2059
2060                         if ( factor_count <= 1 ) {
2061                                 return poly;
2062                         }
2063
2064                         if ( minimalr < 0 ) {
2065                                 minimalr = factor_count;
2066                         }
2067                         else if ( minimalr == factor_count ) {
2068                                 ++trialcount;
2069                                 ++modulus;
2070                         }
2071                         else if ( minimalr > factor_count ) {
2072                                 minimalr = factor_count;
2073                                 trialcount = 0;
2074                         }
2075                         if ( minimalr <= 1 ) {
2076                                 return poly;
2077                         }
2078                 }
2079
2080                 vector<numeric> ftilde((vnlst.nops()-1)/2+1);
2081                 ftilde[0] = ex_to<numeric>(vnlst.op(0));
2082                 for ( size_t i=1; i<ftilde.size(); ++i ) {
2083                         ex ft = vnlst.op((i-1)*2+1);
2084                         s = syms.begin();
2085                         ++s;
2086                         for ( size_t j=0; j<a.size(); ++j ) {
2087                                 ft = ft.subs(*s == a[j]);
2088                                 ++s;
2089                         }
2090                         ftilde[i] = ex_to<numeric>(ft);
2091                 }
2092
2093                 vector<bool> used_flag((vnlst.nops()-1)/2+1, false);
2094                 vector<ex> D(factor_count, 1);
2095                 for ( size_t i=0; i<=factor_count; ++i ) {
2096                         numeric prefac;
2097                         if ( i == 0 ) {
2098                                 prefac = ex_to<numeric>(ufaclst.op(0));
2099                                 ftilde[0] = ftilde[0] / prefac;
2100                                 vnlst.let_op(0) = vnlst.op(0) / prefac;
2101                                 continue;
2102                         }
2103                         else {
2104                                 prefac = ex_to<numeric>(ufaclst.op(2*(i-1)+1).lcoeff(x));
2105                         }
2106                         for ( size_t j=(vnlst.nops()-1)/2+1; j>0; --j ) {
2107                                 if ( abs(ftilde[j-1]) == 1 ) {
2108                                         used_flag[j-1] = true;
2109                                         continue;
2110                                 }
2111                                 numeric g = gcd(prefac, ftilde[j-1]);
2112                                 if ( g != 1 ) {
2113                                         prefac = prefac / g;
2114                                         numeric count = abs(iquo(g, ftilde[j-1]));
2115                                         used_flag[j-1] = true;
2116                                         if ( i > 0 ) {
2117                                                 if ( j == 1 ) {
2118                                                         D[i-1] = D[i-1] * pow(vnlst.op(0), count);
2119                                                 }
2120                                                 else {
2121                                                         D[i-1] = D[i-1] * pow(vnlst.op(2*(j-2)+1), count);
2122                                                 }
2123                                         }
2124                                         else {
2125                                                 ftilde[j-1] = ftilde[j-1] / prefac;
2126                                                 break;
2127                                         }
2128                                         ++j;
2129                                 }
2130                         }
2131                 }
2132
2133                 bool some_factor_unused = false;
2134                 for ( size_t i=0; i<used_flag.size(); ++i ) {
2135                         if ( !used_flag[i] ) {
2136                                 some_factor_unused = true;
2137                                 break;
2138                         }
2139                 }
2140                 if ( some_factor_unused ) {
2141                         continue;
2142                 }
2143
2144                 vector<ex> C(factor_count);
2145                 if ( delta == 1 ) {
2146                         for ( size_t i=0; i<D.size(); ++i ) {
2147                                 ex Dtilde = D[i];
2148                                 s = syms.begin();
2149                                 ++s;
2150                                 for ( size_t j=0; j<a.size(); ++j ) {
2151                                         Dtilde = Dtilde.subs(*s == a[j]);
2152                                         ++s;
2153                                 }
2154                                 C[i] = D[i] * (ufaclst.op(2*i+1).lcoeff(x) / Dtilde);
2155                         }
2156                 }
2157                 else {
2158                         for ( size_t i=0; i<D.size(); ++i ) {
2159                                 ex Dtilde = D[i];
2160                                 s = syms.begin();
2161                                 ++s;
2162                                 for ( size_t j=0; j<a.size(); ++j ) {
2163                                         Dtilde = Dtilde.subs(*s == a[j]);
2164                                         ++s;
2165                                 }
2166                                 ex ui;
2167                                 if ( i == 0 ) {
2168                                         ui = ufaclst.op(0);
2169                                 }
2170                                 else {
2171                                         ui = ufaclst.op(2*(i-1)+1);
2172                                 }
2173                                 while ( true ) {
2174                                         ex d = gcd(ui.lcoeff(x), Dtilde);
2175                                         C[i] = D[i] * ( ui.lcoeff(x) / d );
2176                                         ui = ui * ( Dtilde[i] / d );
2177                                         delta = delta / ( Dtilde[i] / d );
2178                                         if ( delta == 1 ) break;
2179                                         ui = delta * ui;
2180                                         C[i] = delta * C[i];
2181                                         pp = pp * pow(delta, D.size()-1);
2182                                 }
2183                         }
2184                 }
2185
2186                 EvalPoint ep;
2187                 vector<EvalPoint> epv;
2188                 s = syms.begin();
2189                 ++s;
2190                 for ( size_t i=0; i<a.size(); ++i ) {
2191                         ep.x = *s++;
2192                         ep.evalpoint = a[i].to_int();
2193                         epv.push_back(ep);
2194                 }
2195
2196                 // calc bound p^l
2197                 int maxdeg = 0;
2198                 for ( size_t i=0; i<factor_count; ++i ) {
2199                         if ( ufaclst[2*i+1].degree(x) > maxdeg ) {
2200                                 maxdeg = ufaclst[2*i+1].degree(x);
2201                         }
2202                 }
2203                 cl_I B = 2*calc_bound(u, x, maxdeg);
2204                 cl_I l = 1;
2205                 cl_I pl = prime;
2206                 while ( pl < B ) {
2207                         l = l + 1;
2208                         pl = pl * prime;
2209                 }
2210
2211                 upvec uvec;
2212                 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
2213                 for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) {
2214                         umodpoly newu;
2215                         umodpoly_from_ex(newu, ufaclst.op(i*2+1), x, R);
2216                         uvec.push_back(newu);
2217                 }
2218
2219                 ex res = hensel_multivar(ufaclst.op(0)*pp, x, epv, prime, l, uvec, C);
2220                 if ( res != lst() ) {
2221                         ex result = cont * ufaclst.op(0);
2222                         for ( size_t i=0; i<res.nops(); ++i ) {
2223                                 result *= res.op(i).content(x) * res.op(i).unit(x);
2224                                 result *= res.op(i).primpart(x);
2225                         }
2226                         return result;
2227                 }
2228         }
2229 }
2230
2231 struct find_symbols_map : public map_function {
2232         exset syms;
2233         ex operator()(const ex& e)
2234         {
2235                 if ( is_a<symbol>(e) ) {
2236                         syms.insert(e);
2237                         return e;
2238                 }
2239                 return e.map(*this);
2240         }
2241 };
2242
2243 static ex factor_sqrfree(const ex& poly)
2244 {
2245         // determine all symbols in poly
2246         find_symbols_map findsymbols;
2247         findsymbols(poly);
2248         if ( findsymbols.syms.size() == 0 ) {
2249                 return poly;
2250         }
2251
2252         if ( findsymbols.syms.size() == 1 ) {
2253                 // univariate case
2254                 const ex& x = *(findsymbols.syms.begin());
2255                 if ( poly.ldegree(x) > 0 ) {
2256                         int ld = poly.ldegree(x);
2257                         ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2258                         return res * pow(x,ld);
2259                 }
2260                 else {
2261                         ex res = factor_univariate(poly, x);
2262                         return res;
2263                 }
2264         }
2265
2266         // multivariate case
2267         ex res = factor_multivariate(poly, findsymbols.syms);
2268         return res;
2269 }
2270
2271 struct apply_factor_map : public map_function {
2272         unsigned options;
2273         apply_factor_map(unsigned options_) : options(options_) { }
2274         ex operator()(const ex& e)
2275         {
2276                 if ( e.info(info_flags::polynomial) ) {
2277                         return factor(e, options);
2278                 }
2279                 if ( is_a<add>(e) ) {
2280                         ex s1, s2;
2281                         for ( size_t i=0; i<e.nops(); ++i ) {
2282                                 if ( e.op(i).info(info_flags::polynomial) ) {
2283                                         s1 += e.op(i);
2284                                 }
2285                                 else {
2286                                         s2 += e.op(i);
2287                                 }
2288                         }
2289                         s1 = s1.eval();
2290                         s2 = s2.eval();
2291                         return factor(s1, options) + s2.map(*this);
2292                 }
2293                 return e.map(*this);
2294         }
2295 };
2296
2297 } // anonymous namespace
2298
2299 ex factor(const ex& poly, unsigned options)
2300 {
2301         // check arguments
2302         if ( !poly.info(info_flags::polynomial) ) {
2303                 if ( options & factor_options::all ) {
2304                         options &= ~factor_options::all;
2305                         apply_factor_map factor_map(options);
2306                         return factor_map(poly);
2307                 }
2308                 return poly;
2309         }
2310
2311         // determine all symbols in poly
2312         find_symbols_map findsymbols;
2313         findsymbols(poly);
2314         if ( findsymbols.syms.size() == 0 ) {
2315                 return poly;
2316         }
2317         lst syms;
2318         exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
2319         for ( ; i!=end; ++i ) {
2320                 syms.append(*i);
2321         }
2322
2323         // make poly square free
2324         ex sfpoly = sqrfree(poly, syms);
2325
2326         // factorize the square free components
2327         if ( is_a<power>(sfpoly) ) {
2328                 // case: (polynomial)^exponent
2329                 const ex& base = sfpoly.op(0);
2330                 if ( !is_a<add>(base) ) {
2331                         // simple case: (monomial)^exponent
2332                         return sfpoly;
2333                 }
2334                 ex f = factor_sqrfree(base);
2335                 return pow(f, sfpoly.op(1));
2336         }
2337         if ( is_a<mul>(sfpoly) ) {
2338                 // case: multiple factors
2339                 ex res = 1;
2340                 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
2341                         const ex& t = sfpoly.op(i);
2342                         if ( is_a<power>(t) ) {
2343                                 const ex& base = t.op(0);
2344                                 if ( !is_a<add>(base) ) {
2345                                         res *= t;
2346                                 }
2347                                 else {
2348                                         ex f = factor_sqrfree(base);
2349                                         res *= pow(f, t.op(1));
2350                                 }
2351                         }
2352                         else if ( is_a<add>(t) ) {
2353                                 ex f = factor_sqrfree(t);
2354                                 res *= f;
2355                         }
2356                         else {
2357                                 res *= t;
2358                         }
2359                 }
2360                 return res;
2361         }
2362         if ( is_a<symbol>(sfpoly) ) {
2363                 return poly;
2364         }
2365         // case: (polynomial)
2366         ex f = factor_sqrfree(sfpoly);
2367         return f;
2368 }
2369
2370 } // namespace GiNaC
2371
2372 #ifdef DEBUGFACTOR
2373 #include "test.h"
2374 #endif