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