Added modular square free factorization.
[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 ( int 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 true;
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         bool nontrivial = false;
970         while ( i <= nhalf ) {
971                 expt_pos(w, q);
972                 umodpoly buf;
973                 rem(w, a, buf);
974                 w = buf;
975                 umodpoly wx = w - x;
976                 gcd(a, wx, buf);
977                 if ( unequal_one(buf) ) {
978                         degrees.push_back(i);
979                         ddfactors.push_back(buf);
980                 }
981                 if ( unequal_one(buf) ) {
982                         umodpoly buf2;
983                         div(a, buf, buf2);
984                         a = buf2;
985                         nhalf = degree(a)/2;
986                         rem(w, a, buf);
987                         w = buf;
988                 }
989                 ++i;
990         }
991         if ( unequal_one(a) ) {
992                 degrees.push_back(degree(a));
993                 ddfactors.push_back(a);
994         }
995 }
996
997 static void same_degree_factor(const umodpoly& a, upvec& upv)
998 {
999         cl_modint_ring R = a[0].ring();
1000         int deg = degree(a);
1001
1002         vector<int> degrees;
1003         upvec ddfactors;
1004         distinct_degree_factor(a, degrees, ddfactors);
1005
1006         for ( size_t i=0; i<degrees.size(); ++i ) {
1007                 if ( degrees[i] == degree(ddfactors[i]) ) {
1008                         upv.push_back(ddfactors[i]);
1009                 }
1010                 else {
1011                         berlekamp(ddfactors[i], upv);
1012                 }
1013         }
1014 }
1015
1016 static void factor_modular(const umodpoly& p, upvec& upv)
1017 {
1018         upvec factors;
1019         vector<int> mult;
1020         modsqrfree(p, factors, mult);
1021
1022 #define USE_SAME_DEGREE_FACTOR
1023 #ifdef USE_SAME_DEGREE_FACTOR
1024         for ( size_t i=0; i<factors.size(); ++i ) {
1025                 upvec upvbuf;
1026                 same_degree_factor(factors[i], upvbuf);
1027                 for ( int j=mult[i]; j>0; --j ) {
1028                         upv.insert(upv.end(), upvbuf.begin(), upvbuf.end());
1029                 }
1030         }
1031 #else
1032         for ( size_t i=0; i<factors.size(); ++i ) {
1033                 upvec upvbuf;
1034                 berlekamp(factors[i], upvbuf);
1035                 if ( upvbuf.size() ) {
1036                         for ( size_t j=0; j<upvbuf.size(); ++j ) {
1037                                 upv.insert(upv.end(), mult[i], upvbuf[j]);
1038                         }
1039                 }
1040                 else {
1041                         for ( int j=mult[i]; j>0; --j ) {
1042                                 upv.push_back(factors[i]);
1043                         }
1044                 }
1045         }
1046 #endif
1047 }
1048
1049 /** Calculates polynomials s and t such that a*s+b*t==1.
1050  *  Assertion: a and b are relatively prime and not zero.
1051  *
1052  *  @param[in]  a  polynomial
1053  *  @param[in]  b  polynomial
1054  *  @param[out] s  polynomial
1055  *  @param[out] t  polynomial
1056  */
1057 static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpoly& t)
1058 {
1059         if ( degree(a) < degree(b) ) {
1060                 exteuclid(b, a, t, s);
1061                 return;
1062         }
1063
1064         umodpoly one(1, a[0].ring()->one());
1065         umodpoly c = a; normalize_in_field(c);
1066         umodpoly d = b; normalize_in_field(d);
1067         s = one;
1068         t.clear();
1069         umodpoly d1;
1070         umodpoly d2 = one;
1071         umodpoly q;
1072         while ( true ) {
1073                 div(c, d, q);
1074                 umodpoly r = c - q * d;
1075                 umodpoly r1 = s - q * d1;
1076                 umodpoly r2 = t - q * d2;
1077                 c = d;
1078                 s = d1;
1079                 t = d2;
1080                 if ( r.empty() ) break;
1081                 d = r;
1082                 d1 = r1;
1083                 d2 = r2;
1084         }
1085         cl_MI fac = recip(lcoeff(a) * lcoeff(c));
1086         umodpoly::iterator i = s.begin(), end = s.end();
1087         for ( ; i!=end; ++i ) {
1088                 *i = *i * fac;
1089         }
1090         canonicalize(s);
1091         fac = recip(lcoeff(b) * lcoeff(c));
1092         i = t.begin(), end = t.end();
1093         for ( ; i!=end; ++i ) {
1094                 *i = *i * fac;
1095         }
1096         canonicalize(t);
1097 }
1098
1099 static upoly replace_lc(const upoly& poly, const cl_I& lc)
1100 {
1101         if ( poly.empty() ) return poly;
1102         upoly r = poly;
1103         r.back() = lc;
1104         return r;
1105 }
1106
1107 static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w)
1108 {
1109         upoly a = a_;
1110         const cl_modint_ring& R = u1_[0].ring();
1111
1112         // calc bound B
1113         cl_R maxcoeff = 0;
1114         for ( int i=degree(a); i>=0; --i ) {
1115                 maxcoeff = maxcoeff + square(abs(a[i]));
1116         }
1117         cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(maxcoeff)));
1118         cl_I maxdegree = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
1119         cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
1120
1121         // step 1
1122         cl_I alpha = lcoeff(a);
1123         a = a * alpha;
1124         umodpoly nu1 = u1_;
1125         normalize_in_field(nu1);
1126         umodpoly nw1 = w1_;
1127         normalize_in_field(nw1);
1128         upoly phi;
1129         phi = umodpoly_to_upoly(nu1) * alpha;
1130         umodpoly u1;
1131         umodpoly_from_upoly(u1, phi, R);
1132         phi = umodpoly_to_upoly(nw1) * alpha;
1133         umodpoly w1;
1134         umodpoly_from_upoly(w1, phi, R);
1135
1136         // step 2
1137         umodpoly s;
1138         umodpoly t;
1139         exteuclid(u1, w1, s, t);
1140
1141         // step 3
1142         u = replace_lc(umodpoly_to_upoly(u1), alpha);
1143         w = replace_lc(umodpoly_to_upoly(w1), alpha);
1144         upoly e = a - u * w;
1145         cl_I modulus = p;
1146         const cl_I maxmodulus = 2*B*abs(alpha);
1147
1148         // step 4
1149         while ( !e.empty() && modulus < maxmodulus ) {
1150                 // ad-hoc divisablity check
1151                 for ( size_t k=0; k<e.size(); ++k ) {
1152                         if ( !zerop(mod(e[k], modulus)) ) {
1153                                 goto quickexit;
1154                         }
1155                 }
1156                 upoly c = e / modulus;
1157                 phi = umodpoly_to_upoly(s) * c;
1158                 umodpoly sigmatilde;
1159                 umodpoly_from_upoly(sigmatilde, phi, R);
1160                 phi = umodpoly_to_upoly(t) * c;
1161                 umodpoly tautilde;
1162                 umodpoly_from_upoly(tautilde, phi, R);
1163                 umodpoly r, q;
1164                 remdiv(sigmatilde, w1, r, q);
1165                 umodpoly sigma = r;
1166                 phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1);
1167                 umodpoly tau;
1168                 umodpoly_from_upoly(tau, phi, R);
1169                 u = u + umodpoly_to_upoly(tau) * modulus;
1170                 w = w + umodpoly_to_upoly(sigma) * modulus;
1171                 e = a - u * w;
1172                 modulus = modulus * p;
1173         }
1174 quickexit: ;
1175
1176         // step 5
1177         if ( e.empty() ) {
1178                 cl_I g = u[0];
1179                 for ( size_t i=1; i<u.size(); ++i ) {
1180                         g = gcd(g, u[i]);
1181                         if ( g == 1 ) break;
1182                 }
1183                 if ( g != 1 ) {
1184                         u = u / g;
1185                         w = w * g;
1186                 }
1187                 if ( alpha != 1 ) {
1188                         w = w / alpha;
1189                 }
1190         }
1191         else {
1192                 u.clear();
1193         }
1194 }
1195
1196 static unsigned int next_prime(unsigned int p)
1197 {
1198         static vector<unsigned int> primes;
1199         if ( primes.size() == 0 ) {
1200                 primes.push_back(3); primes.push_back(5); primes.push_back(7);
1201         }
1202         vector<unsigned int>::const_iterator it = primes.begin();
1203         if ( p >= primes.back() ) {
1204                 unsigned int candidate = primes.back() + 2;
1205                 while ( true ) {
1206                         size_t n = primes.size()/2;
1207                         for ( size_t i=0; i<n; ++i ) {
1208                                 if ( candidate % primes[i] ) continue;
1209                                 candidate += 2;
1210                                 i=-1;
1211                         }
1212                         primes.push_back(candidate);
1213                         if ( candidate > p ) break;
1214                 }
1215                 return candidate;
1216         }
1217         vector<unsigned int>::const_iterator end = primes.end();
1218         for ( ; it!=end; ++it ) {
1219                 if ( *it > p ) {
1220                         return *it;
1221                 }
1222         }
1223         throw logic_error("next_prime: should not reach this point!");
1224 }
1225
1226 class factor_partition
1227 {
1228 public:
1229         factor_partition(const upvec& factors_) : factors(factors_)
1230         {
1231                 n = factors.size();
1232                 k.resize(n, 1);
1233                 k[0] = 0;
1234                 sum = n-1;
1235                 one.resize(1, factors.front()[0].ring()->one());
1236                 split();
1237         }
1238         int operator[](size_t i) const { return k[i]; }
1239         size_t size() const { return n; }
1240         size_t size_first() const { return n-sum; }
1241         size_t size_second() const { return sum; }
1242 #ifdef DEBUGFACTOR
1243         void get() const { DCOUTVAR(k); }
1244 #endif
1245         bool next()
1246         {
1247                 for ( size_t i=n-1; i>=1; --i ) {
1248                         if ( k[i] ) {
1249                                 --k[i];
1250                                 --sum;
1251                                 if ( sum > 0 ) {
1252                                         split();
1253                                         return true;
1254                                 }
1255                                 else {
1256                                         return false;
1257                                 }
1258                         }
1259                         ++k[i];
1260                         ++sum;
1261                 }
1262                 return false;
1263         }
1264         void split()
1265         {
1266                 left = one;
1267                 right = one;
1268                 for ( size_t i=0; i<n; ++i ) {
1269                         if ( k[i] ) {
1270                                 right = right * factors[i];
1271                         }
1272                         else {
1273                                 left = left * factors[i];
1274                         }
1275                 }
1276         }
1277 public:
1278         umodpoly left, right;
1279 private:
1280         upvec factors;
1281         umodpoly one;
1282         size_t n, sum;
1283         vector<int> k;
1284 };
1285
1286 struct ModFactors
1287 {
1288         upoly poly;
1289         upvec factors;
1290 };
1291
1292 static ex factor_univariate(const ex& poly, const ex& x)
1293 {
1294         ex unit, cont, prim_ex;
1295         poly.unitcontprim(x, unit, cont, prim_ex);
1296         upoly prim;
1297         upoly_from_ex(prim, prim_ex, x);
1298
1299         // determine proper prime and minimize number of modular factors
1300         unsigned int p = 3, lastp = 3;
1301         cl_modint_ring R;
1302         unsigned int trials = 0;
1303         unsigned int minfactors = 0;
1304         cl_I lc = lcoeff(prim);
1305         upvec factors;
1306         while ( trials < 2 ) {
1307                 umodpoly modpoly;
1308                 while ( true ) {
1309                         p = next_prime(p);
1310                         if ( !zerop(rem(lc, p)) ) {
1311                                 R = find_modint_ring(p);
1312                                 umodpoly_from_upoly(modpoly, prim, R);
1313                                 if ( squarefree(modpoly) ) break;
1314                         }
1315                 }
1316
1317                 // do modular factorization
1318                 upvec trialfactors;
1319                 factor_modular(modpoly, trialfactors);
1320                 if ( trialfactors.size() <= 1 ) {
1321                         // irreducible for sure
1322                         return poly;
1323                 }
1324
1325                 if ( minfactors == 0 || trialfactors.size() < minfactors ) {
1326                         factors = trialfactors;
1327                         minfactors = factors.size();
1328                         lastp = p;
1329                         trials = 1;
1330                 }
1331                 else {
1332                         ++trials;
1333                 }
1334         }
1335         p = lastp;
1336         R = find_modint_ring(p);
1337
1338         // lift all factor combinations
1339         stack<ModFactors> tocheck;
1340         ModFactors mf;
1341         mf.poly = prim;
1342         mf.factors = factors;
1343         tocheck.push(mf);
1344         upoly f1, f2;
1345         ex result = 1;
1346         while ( tocheck.size() ) {
1347                 const size_t n = tocheck.top().factors.size();
1348                 factor_partition part(tocheck.top().factors);
1349                 while ( true ) {
1350                         hensel_univar(tocheck.top().poly, p, part.left, part.right, f1, f2);
1351                         if ( !f1.empty() ) {
1352                                 if ( part.size_first() == 1 ) {
1353                                         if ( part.size_second() == 1 ) {
1354                                                 result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1355                                                 tocheck.pop();
1356                                                 break;
1357                                         }
1358                                         result *= upoly_to_ex(f1, x);
1359                                         tocheck.top().poly = f2;
1360                                         for ( size_t i=0; i<n; ++i ) {
1361                                                 if ( part[i] == 0 ) {
1362                                                         tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1363                                                         break;
1364                                                 }
1365                                         }
1366                                         break;
1367                                 }
1368                                 else if ( part.size_second() == 1 ) {
1369                                         if ( part.size_first() == 1 ) {
1370                                                 result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
1371                                                 tocheck.pop();
1372                                                 break;
1373                                         }
1374                                         result *= upoly_to_ex(f2, x);
1375                                         tocheck.top().poly = f1;
1376                                         for ( size_t i=0; i<n; ++i ) {
1377                                                 if ( part[i] == 1 ) {
1378                                                         tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1379                                                         break;
1380                                                 }
1381                                         }
1382                                         break;
1383                                 }
1384                                 else {
1385                                         upvec newfactors1(part.size_first()), newfactors2(part.size_second());
1386                                         upvec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1387                                         for ( size_t i=0; i<n; ++i ) {
1388                                                 if ( part[i] ) {
1389                                                         *i2++ = tocheck.top().factors[i];
1390                                                 }
1391                                                 else {
1392                                                         *i1++ = tocheck.top().factors[i];
1393                                                 }
1394                                         }
1395                                         tocheck.top().factors = newfactors1;
1396                                         tocheck.top().poly = f1;
1397                                         ModFactors mf;
1398                                         mf.factors = newfactors2;
1399                                         mf.poly = f2;
1400                                         tocheck.push(mf);
1401                                         break;
1402                                 }
1403                         }
1404                         else {
1405                                 if ( !part.next() ) {
1406                                         result *= upoly_to_ex(tocheck.top().poly, x);
1407                                         tocheck.pop();
1408                                         break;
1409                                 }
1410                         }
1411                 }
1412         }
1413
1414         return unit * cont * result;
1415 }
1416
1417 struct EvalPoint
1418 {
1419         ex x;
1420         int evalpoint;
1421 };
1422
1423 // forward declaration
1424 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);
1425
1426 upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k)
1427 {
1428         const size_t r = a.size();
1429         cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1430         upvec q(r-1);
1431         q[r-2] = a[r-1];
1432         for ( size_t j=r-2; j>=1; --j ) {
1433                 q[j-1] = a[j] * q[j];
1434         }
1435         umodpoly beta(1, R->one());
1436         upvec s;
1437         for ( size_t j=1; j<r; ++j ) {
1438                 vector<ex> mdarg(2);
1439                 mdarg[0] = umodpoly_to_ex(q[j-1], x);
1440                 mdarg[1] = umodpoly_to_ex(a[j-1], x);
1441                 vector<EvalPoint> empty;
1442                 vector<ex> exsigma = multivar_diophant(mdarg, x, umodpoly_to_ex(beta, x), empty, 0, p, k);
1443                 umodpoly sigma1;
1444                 umodpoly_from_ex(sigma1, exsigma[0], x, R);
1445                 umodpoly sigma2;
1446                 umodpoly_from_ex(sigma2, exsigma[1], x, R);
1447                 beta = sigma1;
1448                 s.push_back(sigma2);
1449         }
1450         s.push_back(beta);
1451         return s;
1452 }
1453
1454 /**
1455  *  Assert: a not empty.
1456  */
1457 void change_modulus(const cl_modint_ring& R, umodpoly& a)
1458 {
1459         if ( a.empty() ) return;
1460         cl_modint_ring oldR = a[0].ring();
1461         umodpoly::iterator i = a.begin(), end = a.end();
1462         for ( ; i!=end; ++i ) {
1463                 *i = R->canonhom(oldR->retract(*i));
1464         }
1465         canonicalize(a);
1466 }
1467
1468 void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_)
1469 {
1470         cl_modint_ring R = find_modint_ring(p);
1471         umodpoly amod = a;
1472         change_modulus(R, amod);
1473         umodpoly bmod = b;
1474         change_modulus(R, bmod);
1475
1476         umodpoly smod;
1477         umodpoly tmod;
1478         exteuclid(amod, bmod, smod, tmod);
1479
1480         cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1481         umodpoly s = smod;
1482         change_modulus(Rpk, s);
1483         umodpoly t = tmod;
1484         change_modulus(Rpk, t);
1485
1486         cl_I modulus(p);
1487         umodpoly one(1, Rpk->one());
1488         for ( size_t j=1; j<k; ++j ) {
1489                 umodpoly e = one - a * s - b * t;
1490                 reduce_coeff(e, modulus);
1491                 umodpoly c = e;
1492                 change_modulus(R, c);
1493                 umodpoly sigmabar = smod * c;
1494                 umodpoly taubar = tmod * c;
1495                 umodpoly sigma, q;
1496                 remdiv(sigmabar, bmod, sigma, q);
1497                 umodpoly tau = taubar + q * amod;
1498                 umodpoly sadd = sigma;
1499                 change_modulus(Rpk, sadd);
1500                 cl_MI modmodulus(Rpk, modulus);
1501                 s = s + sadd * modmodulus;
1502                 umodpoly tadd = tau;
1503                 change_modulus(Rpk, tadd);
1504                 t = t + tadd * modmodulus;
1505                 modulus = modulus * p;
1506         }
1507
1508         s_ = s; t_ = t;
1509 }
1510
1511 upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1512 {
1513         cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1514
1515         const size_t r = a.size();
1516         upvec result;
1517         if ( r > 2 ) {
1518                 upvec s = multiterm_eea_lift(a, x, p, k);
1519                 for ( size_t j=0; j<r; ++j ) {
1520                         umodpoly bmod = umodpoly_to_umodpoly(s[j], R, m);
1521                         umodpoly buf;
1522                         rem(bmod, a[j], buf);
1523                         result.push_back(buf);
1524                 }
1525         }
1526         else {
1527                 umodpoly s, t;
1528                 eea_lift(a[1], a[0], x, p, k, s, t);
1529                 umodpoly bmod = umodpoly_to_umodpoly(s, R, m);
1530                 umodpoly buf, q;
1531                 remdiv(bmod, a[0], buf, q);
1532                 result.push_back(buf);
1533                 umodpoly t1mod = umodpoly_to_umodpoly(t, R, m);
1534                 buf = t1mod + q * a[1];
1535                 result.push_back(buf);
1536         }
1537
1538         return result;
1539 }
1540
1541 struct make_modular_map : public map_function {
1542         cl_modint_ring R;
1543         make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1544         ex operator()(const ex& e)
1545         {
1546                 if ( is_a<add>(e) || is_a<mul>(e) ) {
1547                         return e.map(*this);
1548                 }
1549                 else if ( is_a<numeric>(e) ) {
1550                         numeric mod(R->modulus);
1551                         numeric halfmod = (mod-1)/2;
1552                         cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1553                         numeric n(R->retract(emod));
1554                         if ( n > halfmod ) {
1555                                 return n-mod;
1556                         }
1557                         else {
1558                                 return n;
1559                         }
1560                 }
1561                 return e;
1562         }
1563 };
1564
1565 static ex make_modular(const ex& e, const cl_modint_ring& R)
1566 {
1567         make_modular_map map(R);
1568         return map(e.expand());
1569 }
1570
1571 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)
1572 {
1573         vector<ex> a = a_;
1574
1575         const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1576         const size_t r = a.size();
1577         const size_t nu = I.size() + 1;
1578
1579         vector<ex> sigma;
1580         if ( nu > 1 ) {
1581                 ex xnu = I.back().x;
1582                 int alphanu = I.back().evalpoint;
1583
1584                 ex A = 1;
1585                 for ( size_t i=0; i<r; ++i ) {
1586                         A *= a[i];
1587                 }
1588                 vector<ex> b(r);
1589                 for ( size_t i=0; i<r; ++i ) {
1590                         b[i] = normal(A / a[i]);
1591                 }
1592
1593                 vector<ex> anew = a;
1594                 for ( size_t i=0; i<r; ++i ) {
1595                         anew[i] = anew[i].subs(xnu == alphanu);
1596                 }
1597                 ex cnew = c.subs(xnu == alphanu);
1598                 vector<EvalPoint> Inew = I;
1599                 Inew.pop_back();
1600                 sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1601
1602                 ex buf = c;
1603                 for ( size_t i=0; i<r; ++i ) {
1604                         buf -= sigma[i] * b[i];
1605                 }
1606                 ex e = make_modular(buf, R);
1607
1608                 ex monomial = 1;
1609                 for ( size_t m=1; m<=d; ++m ) {
1610                         while ( !e.is_zero() && e.has(xnu) ) {
1611                                 monomial *= (xnu - alphanu);
1612                                 monomial = expand(monomial);
1613                                 ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1614                                 cm = make_modular(cm, R);
1615                                 if ( !cm.is_zero() ) {
1616                                         vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1617                                         ex buf = e;
1618                                         for ( size_t j=0; j<delta_s.size(); ++j ) {
1619                                                 delta_s[j] *= monomial;
1620                                                 sigma[j] += delta_s[j];
1621                                                 buf -= delta_s[j] * b[j];
1622                                         }
1623                                         e = make_modular(buf, R);
1624                                 }
1625                         }
1626                 }
1627         }
1628         else {
1629                 upvec amod;
1630                 for ( size_t i=0; i<a.size(); ++i ) {
1631                         umodpoly up;
1632                         umodpoly_from_ex(up, a[i], x, R);
1633                         amod.push_back(up);
1634                 }
1635
1636                 sigma.insert(sigma.begin(), r, 0);
1637                 size_t nterms;
1638                 ex z;
1639                 if ( is_a<add>(c) ) {
1640                         nterms = c.nops();
1641                         z = c.op(0);
1642                 }
1643                 else {
1644                         nterms = 1;
1645                         z = c;
1646                 }
1647                 for ( size_t i=0; i<nterms; ++i ) {
1648                         int m = z.degree(x);
1649                         cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1650                         upvec delta_s = univar_diophant(amod, x, m, p, k);
1651                         cl_MI modcm;
1652                         cl_I poscm = cm;
1653                         while ( poscm < 0 ) {
1654                                 poscm = poscm + expt_pos(cl_I(p),k);
1655                         }
1656                         modcm = cl_MI(R, poscm);
1657                         for ( size_t j=0; j<delta_s.size(); ++j ) {
1658                                 delta_s[j] = delta_s[j] * modcm;
1659                                 sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j], x);
1660                         }
1661                         if ( nterms > 1 ) {
1662                                 z = c.op(i+1);
1663                         }
1664                 }
1665         }
1666
1667         for ( size_t i=0; i<sigma.size(); ++i ) {
1668                 sigma[i] = make_modular(sigma[i], R);
1669         }
1670
1671         return sigma;
1672 }
1673
1674 #ifdef DEBUGFACTOR
1675 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1676 {
1677         for ( size_t i=0; i<v.size(); ++i ) {
1678                 o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1679         }
1680         return o;
1681 }
1682 #endif // def DEBUGFACTOR
1683
1684 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)
1685 {
1686         const size_t nu = I.size() + 1;
1687         const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1688
1689         vector<ex> A(nu);
1690         A[nu-1] = a;
1691
1692         for ( size_t j=nu; j>=2; --j ) {
1693                 ex x = I[j-2].x;
1694                 int alpha = I[j-2].evalpoint;
1695                 A[j-2] = A[j-1].subs(x==alpha);
1696                 A[j-2] = make_modular(A[j-2], R);
1697         }
1698
1699         int maxdeg = a.degree(I.front().x);
1700         for ( size_t i=1; i<I.size(); ++i ) {
1701                 int maxdeg2 = a.degree(I[i].x);
1702                 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
1703         }
1704
1705         const size_t n = u.size();
1706         vector<ex> U(n);
1707         for ( size_t i=0; i<n; ++i ) {
1708                 U[i] = umodpoly_to_ex(u[i], x);
1709         }
1710
1711         for ( size_t j=2; j<=nu; ++j ) {
1712                 vector<ex> U1 = U;
1713                 ex monomial = 1;
1714                 for ( size_t m=0; m<n; ++m) {
1715                         if ( lcU[m] != 1 ) {
1716                                 ex coef = lcU[m];
1717                                 for ( size_t i=j-1; i<nu-1; ++i ) {
1718                                         coef = coef.subs(I[i].x == I[i].evalpoint);
1719                                 }
1720                                 coef = make_modular(coef, R);
1721                                 int deg = U[m].degree(x);
1722                                 U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
1723                         }
1724                 }
1725                 ex Uprod = 1;
1726                 for ( size_t i=0; i<n; ++i ) {
1727                         Uprod *= U[i];
1728                 }
1729                 ex e = expand(A[j-1] - Uprod);
1730
1731                 vector<EvalPoint> newI;
1732                 for ( size_t i=1; i<=j-2; ++i ) {
1733                         newI.push_back(I[i-1]);
1734                 }
1735
1736                 ex xj = I[j-2].x;
1737                 int alphaj = I[j-2].evalpoint;
1738                 size_t deg = A[j-1].degree(xj);
1739                 for ( size_t k=1; k<=deg; ++k ) {
1740                         if ( !e.is_zero() ) {
1741                                 monomial *= (xj - alphaj);
1742                                 monomial = expand(monomial);
1743                                 ex dif = e.diff(ex_to<symbol>(xj), k);
1744                                 ex c = dif.subs(xj==alphaj) / factorial(k);
1745                                 if ( !c.is_zero() ) {
1746                                         vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
1747                                         for ( size_t i=0; i<n; ++i ) {
1748                                                 deltaU[i] *= monomial;
1749                                                 U[i] += deltaU[i];
1750                                                 U[i] = make_modular(U[i], R);
1751                                         }
1752                                         ex Uprod = 1;
1753                                         for ( size_t i=0; i<n; ++i ) {
1754                                                 Uprod *= U[i];
1755                                         }
1756                                         e = A[j-1] - Uprod;
1757                                         e = make_modular(e, R);
1758                                 }
1759                         }
1760                 }
1761         }
1762
1763         ex acand = 1;
1764         for ( size_t i=0; i<U.size(); ++i ) {
1765                 acand *= U[i];
1766         }
1767         if ( expand(a-acand).is_zero() ) {
1768                 lst res;
1769                 for ( size_t i=0; i<U.size(); ++i ) {
1770                         res.append(U[i]);
1771                 }
1772                 return res;
1773         }
1774         else {
1775                 lst res;
1776                 return lst();
1777         }
1778 }
1779
1780 static ex put_factors_into_lst(const ex& e)
1781 {
1782         lst result;
1783
1784         if ( is_a<numeric>(e) ) {
1785                 result.append(e);
1786                 return result;
1787         }
1788         if ( is_a<power>(e) ) {
1789                 result.append(1);
1790                 result.append(e.op(0));
1791                 result.append(e.op(1));
1792                 return result;
1793         }
1794         if ( is_a<symbol>(e) || is_a<add>(e) ) {
1795                 result.append(1);
1796                 result.append(e);
1797                 result.append(1);
1798                 return result;
1799         }
1800         if ( is_a<mul>(e) ) {
1801                 ex nfac = 1;
1802                 for ( size_t i=0; i<e.nops(); ++i ) {
1803                         ex op = e.op(i);
1804                         if ( is_a<numeric>(op) ) {
1805                                 nfac = op;
1806                         }
1807                         if ( is_a<power>(op) ) {
1808                                 result.append(op.op(0));
1809                                 result.append(op.op(1));
1810                         }
1811                         if ( is_a<symbol>(op) || is_a<add>(op) ) {
1812                                 result.append(op);
1813                                 result.append(1);
1814                         }
1815                 }
1816                 result.prepend(nfac);
1817                 return result;
1818         }
1819         throw runtime_error("put_factors_into_lst: bad term.");
1820 }
1821
1822 #ifdef DEBUGFACTOR
1823 ostream& operator<<(ostream& o, const vector<numeric>& v)
1824 {
1825         for ( size_t i=0; i<v.size(); ++i ) {
1826                 o << v[i] << " ";
1827         }
1828         return o;
1829 }
1830 #endif // def DEBUGFACTOR
1831
1832 static bool checkdivisors(const lst& f, vector<numeric>& d)
1833 {
1834         const int k = f.nops()-2;
1835         numeric q, r;
1836         d[0] = ex_to<numeric>(f.op(0) * f.op(f.nops()-1));
1837         if ( d[0] == 1 && k == 1 && abs(f.op(1)) != 1 ) {
1838                 return false;
1839         }
1840         for ( int i=1; i<=k; ++i ) {
1841                 q = ex_to<numeric>(abs(f.op(i)));
1842                 for ( int j=i-1; j>=0; --j ) {
1843                         r = d[j];
1844                         do {
1845                                 r = gcd(r, q);
1846                                 q = q/r;
1847                         } while ( r != 1 );
1848                         if ( q == 1 ) {
1849                                 return true;
1850                         }
1851                 }
1852                 d[i] = q;
1853         }
1854         return false;
1855 }
1856
1857 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)
1858 {
1859         // computation of d is actually not necessary
1860         const ex& x = *syms.begin();
1861         bool trying = true;
1862         do {
1863                 ex u0 = u;
1864                 ex vna = vn;
1865                 ex vnatry;
1866                 exset::const_iterator s = syms.begin();
1867                 ++s;
1868                 for ( size_t i=0; i<a.size(); ++i ) {
1869                         do {
1870                                 a[i] = mod(numeric(rand()), 2*modulus) - modulus;
1871                                 vnatry = vna.subs(*s == a[i]);
1872                         } while ( vnatry == 0 );
1873                         vna = vnatry;
1874                         u0 = u0.subs(*s == a[i]);
1875                         ++s;
1876                 }
1877                 if ( gcd(u0,u0.diff(ex_to<symbol>(x))) != 1 ) {
1878                         continue;
1879                 }
1880                 if ( is_a<numeric>(vn) ) {
1881                         trying = false;
1882                 }
1883                 else {
1884                         lst fnum;
1885                         lst::const_iterator i = ex_to<lst>(f).begin();
1886                         fnum.append(*i++);
1887                         bool problem = false;
1888                         while ( i!=ex_to<lst>(f).end() ) {
1889                                 ex fs = *i;
1890                                 if ( !is_a<numeric>(fs) ) {
1891                                         s = syms.begin();
1892                                         ++s;
1893                                         for ( size_t j=0; j<a.size(); ++j ) {
1894                                                 fs = fs.subs(*s == a[j]);
1895                                                 ++s;
1896                                         }
1897                                         if ( abs(fs) == 1 ) {
1898                                                 problem = true;
1899                                                 break;
1900                                         }
1901                                 }
1902                                 fnum.append(fs);
1903                                 ++i; ++i;
1904                         }
1905                         if ( problem ) {
1906                                 return true;
1907                         }
1908                         ex con = u0.content(x);
1909                         fnum.append(con);
1910                         trying = checkdivisors(fnum, d);
1911                 }
1912         } while ( trying );
1913         return false;
1914 }
1915
1916 static ex factor_multivariate(const ex& poly, const exset& syms)
1917 {
1918         exset::const_iterator s;
1919         const ex& x = *syms.begin();
1920
1921         /* make polynomial primitive */
1922         ex p = poly.expand().collect(x);
1923         ex cont = p.lcoeff(x);
1924         for ( numeric i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
1925                 cont = gcd(cont, p.coeff(x,ex_to<numeric>(i).to_int()));
1926                 if ( cont == 1 ) break;
1927         }
1928         ex pp = expand(normal(p / cont));
1929         if ( !is_a<numeric>(cont) ) {
1930                 return factor(cont) * factor(pp);
1931         }
1932
1933         /* factor leading coefficient */
1934         pp = pp.collect(x);
1935         ex vn = pp.lcoeff(x);
1936         pp = pp.expand();
1937         ex vnlst;
1938         if ( is_a<numeric>(vn) ) {
1939                 vnlst = lst(vn);
1940         }
1941         else {
1942                 ex vnfactors = factor(vn);
1943                 vnlst = put_factors_into_lst(vnfactors);
1944         }
1945
1946         const numeric maxtrials = 3;
1947         numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3;
1948         numeric minimalr = -1;
1949         vector<numeric> a(syms.size()-1, 0);
1950         vector<numeric> d((vnlst.nops()-1)/2+1, 0);
1951
1952         while ( true ) {
1953                 numeric trialcount = 0;
1954                 ex u, delta;
1955                 unsigned int prime = 3;
1956                 size_t factor_count = 0;
1957                 ex ufac;
1958                 ex ufaclst;
1959                 while ( trialcount < maxtrials ) {
1960                         bool problem = generate_set(pp, vn, syms, vnlst, modulus, a, d);
1961                         if ( problem ) {
1962                                 ++modulus;
1963                                 continue;
1964                         }
1965                         u = pp;
1966                         s = syms.begin();
1967                         ++s;
1968                         for ( size_t i=0; i<a.size(); ++i ) {
1969                                 u = u.subs(*s == a[i]);
1970                                 ++s;
1971                         }
1972                         delta = u.content(x);
1973
1974                         // determine proper prime
1975                         prime = 3;
1976                         cl_modint_ring R = find_modint_ring(prime);
1977                         while ( true ) {
1978                                 if ( irem(ex_to<numeric>(u.lcoeff(x)), prime) != 0 ) {
1979                                         umodpoly modpoly;
1980                                         umodpoly_from_ex(modpoly, u, x, R);
1981                                         if ( squarefree(modpoly) ) break;
1982                                 }
1983                                 prime = next_prime(prime);
1984                                 R = find_modint_ring(prime);
1985                         }
1986
1987                         ufac = factor(u);
1988                         ufaclst = put_factors_into_lst(ufac);
1989                         factor_count = (ufaclst.nops()-1)/2;
1990
1991                         // veto factorization for which gcd(u_i, u_j) != 1 for all i,j
1992                         upvec tryu;
1993                         for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) {
1994                                 umodpoly newu;
1995                                 umodpoly_from_ex(newu, ufaclst.op(i*2+1), x, R);
1996                                 tryu.push_back(newu);
1997                         }
1998                         bool veto = false;
1999                         for ( size_t i=0; i<tryu.size()-1; ++i ) {
2000                                 for ( size_t j=i+1; j<tryu.size(); ++j ) {
2001                                         umodpoly tryg;
2002                                         gcd(tryu[i], tryu[j], tryg);
2003                                         if ( unequal_one(tryg) ) {
2004                                                 veto = true;
2005                                                 goto escape_quickly;
2006                                         }
2007                                 }
2008                         }
2009                         escape_quickly: ;
2010                         if ( veto ) {
2011                                 continue;
2012                         }
2013
2014                         if ( factor_count <= 1 ) {
2015                                 return poly;
2016                         }
2017
2018                         if ( minimalr < 0 ) {
2019                                 minimalr = factor_count;
2020                         }
2021                         else if ( minimalr == factor_count ) {
2022                                 ++trialcount;
2023                                 ++modulus;
2024                         }
2025                         else if ( minimalr > factor_count ) {
2026                                 minimalr = factor_count;
2027                                 trialcount = 0;
2028                         }
2029                         if ( minimalr <= 1 ) {
2030                                 return poly;
2031                         }
2032                 }
2033
2034                 vector<numeric> ftilde((vnlst.nops()-1)/2+1);
2035                 ftilde[0] = ex_to<numeric>(vnlst.op(0));
2036                 for ( size_t i=1; i<ftilde.size(); ++i ) {
2037                         ex ft = vnlst.op((i-1)*2+1);
2038                         s = syms.begin();
2039                         ++s;
2040                         for ( size_t j=0; j<a.size(); ++j ) {
2041                                 ft = ft.subs(*s == a[j]);
2042                                 ++s;
2043                         }
2044                         ftilde[i] = ex_to<numeric>(ft);
2045                 }
2046
2047                 vector<bool> used_flag((vnlst.nops()-1)/2+1, false);
2048                 vector<ex> D(factor_count, 1);
2049                 for ( size_t i=0; i<=factor_count; ++i ) {
2050                         numeric prefac;
2051                         if ( i == 0 ) {
2052                                 prefac = ex_to<numeric>(ufaclst.op(0));
2053                                 ftilde[0] = ftilde[0] / prefac;
2054                                 vnlst.let_op(0) = vnlst.op(0) / prefac;
2055                                 continue;
2056                         }
2057                         else {
2058                                 prefac = ex_to<numeric>(ufaclst.op(2*(i-1)+1).lcoeff(x));
2059                         }
2060                         for ( size_t j=(vnlst.nops()-1)/2+1; j>0; --j ) {
2061                                 if ( abs(ftilde[j-1]) == 1 ) {
2062                                         used_flag[j-1] = true;
2063                                         continue;
2064                                 }
2065                                 numeric g = gcd(prefac, ftilde[j-1]);
2066                                 if ( g != 1 ) {
2067                                         prefac = prefac / g;
2068                                         numeric count = abs(iquo(g, ftilde[j-1]));
2069                                         used_flag[j-1] = true;
2070                                         if ( i > 0 ) {
2071                                                 if ( j == 1 ) {
2072                                                         D[i-1] = D[i-1] * pow(vnlst.op(0), count);
2073                                                 }
2074                                                 else {
2075                                                         D[i-1] = D[i-1] * pow(vnlst.op(2*(j-2)+1), count);
2076                                                 }
2077                                         }
2078                                         else {
2079                                                 ftilde[j-1] = ftilde[j-1] / prefac;
2080                                                 break;
2081                                         }
2082                                         ++j;
2083                                 }
2084                         }
2085                 }
2086
2087                 bool some_factor_unused = false;
2088                 for ( size_t i=0; i<used_flag.size(); ++i ) {
2089                         if ( !used_flag[i] ) {
2090                                 some_factor_unused = true;
2091                                 break;
2092                         }
2093                 }
2094                 if ( some_factor_unused ) {
2095                         continue;
2096                 }
2097
2098                 vector<ex> C(factor_count);
2099                 if ( delta == 1 ) {
2100                         for ( size_t i=0; i<D.size(); ++i ) {
2101                                 ex Dtilde = D[i];
2102                                 s = syms.begin();
2103                                 ++s;
2104                                 for ( size_t j=0; j<a.size(); ++j ) {
2105                                         Dtilde = Dtilde.subs(*s == a[j]);
2106                                         ++s;
2107                                 }
2108                                 C[i] = D[i] * (ufaclst.op(2*i+1).lcoeff(x) / Dtilde);
2109                         }
2110                 }
2111                 else {
2112                         for ( size_t i=0; i<D.size(); ++i ) {
2113                                 ex Dtilde = D[i];
2114                                 s = syms.begin();
2115                                 ++s;
2116                                 for ( size_t j=0; j<a.size(); ++j ) {
2117                                         Dtilde = Dtilde.subs(*s == a[j]);
2118                                         ++s;
2119                                 }
2120                                 ex ui;
2121                                 if ( i == 0 ) {
2122                                         ui = ufaclst.op(0);
2123                                 }
2124                                 else {
2125                                         ui = ufaclst.op(2*(i-1)+1);
2126                                 }
2127                                 while ( true ) {
2128                                         ex d = gcd(ui.lcoeff(x), Dtilde);
2129                                         C[i] = D[i] * ( ui.lcoeff(x) / d );
2130                                         ui = ui * ( Dtilde[i] / d );
2131                                         delta = delta / ( Dtilde[i] / d );
2132                                         if ( delta == 1 ) break;
2133                                         ui = delta * ui;
2134                                         C[i] = delta * C[i];
2135                                         pp = pp * pow(delta, D.size()-1);
2136                                 }
2137                         }
2138                 }
2139
2140                 EvalPoint ep;
2141                 vector<EvalPoint> epv;
2142                 s = syms.begin();
2143                 ++s;
2144                 for ( size_t i=0; i<a.size(); ++i ) {
2145                         ep.x = *s++;
2146                         ep.evalpoint = a[i].to_int();
2147                         epv.push_back(ep);
2148                 }
2149
2150                 // calc bound B
2151                 ex maxcoeff;
2152                 for ( int i=u.degree(x); i>=u.ldegree(x); --i ) {
2153                         maxcoeff += pow(abs(u.coeff(x, i)),2);
2154                 }
2155                 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
2156                 unsigned int maxdegree = 0;
2157                 for ( size_t i=0; i<factor_count; ++i ) {
2158                         if ( ufaclst[2*i+1].degree(x) > (int)maxdegree ) {
2159                                 maxdegree = ufaclst[2*i+1].degree(x);
2160                         }
2161                 }
2162                 cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
2163                 cl_I l = 1;
2164                 cl_I pl = prime;
2165                 while ( pl < B ) {
2166                         l = l + 1;
2167                         pl = pl * prime;
2168                 }
2169
2170                 upvec uvec;
2171                 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
2172                 for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) {
2173                         umodpoly newu;
2174                         umodpoly_from_ex(newu, ufaclst.op(i*2+1), x, R);
2175                         uvec.push_back(newu);
2176                 }
2177
2178                 ex res = hensel_multivar(ufaclst.op(0)*pp, x, epv, prime, l, uvec, C);
2179                 if ( res != lst() ) {
2180                         ex result = cont * ufaclst.op(0);
2181                         for ( size_t i=0; i<res.nops(); ++i ) {
2182                                 result *= res.op(i).content(x) * res.op(i).unit(x);
2183                                 result *= res.op(i).primpart(x);
2184                         }
2185                         return result;
2186                 }
2187         }
2188 }
2189
2190 struct find_symbols_map : public map_function {
2191         exset syms;
2192         ex operator()(const ex& e)
2193         {
2194                 if ( is_a<symbol>(e) ) {
2195                         syms.insert(e);
2196                         return e;
2197                 }
2198                 return e.map(*this);
2199         }
2200 };
2201
2202 static ex factor_sqrfree(const ex& poly)
2203 {
2204         // determine all symbols in poly
2205         find_symbols_map findsymbols;
2206         findsymbols(poly);
2207         if ( findsymbols.syms.size() == 0 ) {
2208                 return poly;
2209         }
2210
2211         if ( findsymbols.syms.size() == 1 ) {
2212                 // univariate case
2213                 const ex& x = *(findsymbols.syms.begin());
2214                 if ( poly.ldegree(x) > 0 ) {
2215                         int ld = poly.ldegree(x);
2216                         ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2217                         return res * pow(x,ld);
2218                 }
2219                 else {
2220                         ex res = factor_univariate(poly, x);
2221                         return res;
2222                 }
2223         }
2224
2225         // multivariate case
2226         ex res = factor_multivariate(poly, findsymbols.syms);
2227         return res;
2228 }
2229
2230 struct apply_factor_map : public map_function {
2231         unsigned options;
2232         apply_factor_map(unsigned options_) : options(options_) { }
2233         ex operator()(const ex& e)
2234         {
2235                 if ( e.info(info_flags::polynomial) ) {
2236                         return factor(e, options);
2237                 }
2238                 if ( is_a<add>(e) ) {
2239                         ex s1, s2;
2240                         for ( size_t i=0; i<e.nops(); ++i ) {
2241                                 if ( e.op(i).info(info_flags::polynomial) ) {
2242                                         s1 += e.op(i);
2243                                 }
2244                                 else {
2245                                         s2 += e.op(i);
2246                                 }
2247                         }
2248                         s1 = s1.eval();
2249                         s2 = s2.eval();
2250                         return factor(s1, options) + s2.map(*this);
2251                 }
2252                 return e.map(*this);
2253         }
2254 };
2255
2256 } // anonymous namespace
2257
2258 ex factor(const ex& poly, unsigned options)
2259 {
2260         // check arguments
2261         if ( !poly.info(info_flags::polynomial) ) {
2262                 if ( options & factor_options::all ) {
2263                         options &= ~factor_options::all;
2264                         apply_factor_map factor_map(options);
2265                         return factor_map(poly);
2266                 }
2267                 return poly;
2268         }
2269
2270         // determine all symbols in poly
2271         find_symbols_map findsymbols;
2272         findsymbols(poly);
2273         if ( findsymbols.syms.size() == 0 ) {
2274                 return poly;
2275         }
2276         lst syms;
2277         exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
2278         for ( ; i!=end; ++i ) {
2279                 syms.append(*i);
2280         }
2281
2282         // make poly square free
2283         ex sfpoly = sqrfree(poly, syms);
2284
2285         // factorize the square free components
2286         if ( is_a<power>(sfpoly) ) {
2287                 // case: (polynomial)^exponent
2288                 const ex& base = sfpoly.op(0);
2289                 if ( !is_a<add>(base) ) {
2290                         // simple case: (monomial)^exponent
2291                         return sfpoly;
2292                 }
2293                 ex f = factor_sqrfree(base);
2294                 return pow(f, sfpoly.op(1));
2295         }
2296         if ( is_a<mul>(sfpoly) ) {
2297                 // case: multiple factors
2298                 ex res = 1;
2299                 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
2300                         const ex& t = sfpoly.op(i);
2301                         if ( is_a<power>(t) ) {
2302                                 const ex& base = t.op(0);
2303                                 if ( !is_a<add>(base) ) {
2304                                         res *= t;
2305                                 }
2306                                 else {
2307                                         ex f = factor_sqrfree(base);
2308                                         res *= pow(f, t.op(1));
2309                                 }
2310                         }
2311                         else if ( is_a<add>(t) ) {
2312                                 ex f = factor_sqrfree(t);
2313                                 res *= f;
2314                         }
2315                         else {
2316                                 res *= t;
2317                         }
2318                 }
2319                 return res;
2320         }
2321         if ( is_a<symbol>(sfpoly) ) {
2322                 return poly;
2323         }
2324         // case: (polynomial)
2325         ex f = factor_sqrfree(sfpoly);
2326         return f;
2327 }
2328
2329 } // namespace GiNaC
2330
2331 #ifdef DEBUGFACTOR
2332 #include "test.h"
2333 #endif