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