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