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