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