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