3 * Polynomial factorization code (Implementation).
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.
13 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
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.
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.
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
34 #include <ginac/ginac.h>
35 using namespace GiNaC;
41 #include "operators.h"
44 #include "relational.h"
66 #define DCOUT(str) cout << #str << endl
67 #define DCOUTVAR(var) cout << #var << ": " << var << endl
68 #define DCOUT2(str,var) cout << #str << ": " << var << endl
72 #define DCOUT2(str,var)
75 // forward declaration
76 ex factor(const ex& poly, unsigned options);
78 // anonymous namespace to hide all utility functions
81 typedef vector<cl_MI> Vec;
82 typedef vector<Vec> VecVec;
85 ostream& operator<<(ostream& o, const Vec& v)
87 Vec::const_iterator i = v.begin(), end = v.end();
93 #endif // def DEBUGFACTOR
96 ostream& operator<<(ostream& o, const VecVec& v)
98 VecVec::const_iterator i = v.begin(), end = v.end();
104 #endif // def DEBUGFACTOR
108 cl_MI c; // coefficient
109 unsigned int exp; // exponent >=0
113 ostream& operator<<(ostream& o, const Term& t)
116 o << "(" << t.c << ")x^" << t.exp;
119 o << "(" << t.c << ")";
123 #endif // def DEBUGFACTOR
128 list<Term> terms; // highest exponent first
130 UniPoly(const cl_modint_ring& ring) : R(ring) { }
131 UniPoly(const cl_modint_ring& ring, const ex& poly, const ex& x) : R(ring)
133 // assert: poly is in Z[x]
135 for ( int i=poly.degree(x); i>=poly.ldegree(x); --i ) {
136 cl_I coeff = the<cl_I>(ex_to<numeric>(poly.coeff(x,i)).to_cl_N());
137 if ( !zerop(coeff) ) {
138 t.c = R->canonhom(coeff);
146 UniPoly(const cl_modint_ring& ring, const UniPoly& poly) : R(ring)
148 if ( R->modulus == poly.R->modulus ) {
152 list<Term>::const_iterator i=poly.terms.begin(), end=poly.terms.end();
153 for ( ; i!=end; ++i ) {
155 terms.back().c = R->canonhom(poly.R->retract(i->c));
156 if ( zerop(terms.back().c) ) {
162 UniPoly(const cl_modint_ring& ring, const Vec& v) : R(ring)
165 for ( unsigned int i=0; i<v.size(); ++i ) {
166 if ( !zerop(v[i]) ) {
173 unsigned int degree() const
175 if ( terms.size() ) {
176 return terms.front().exp;
182 bool zero() const { return (terms.size() == 0); }
183 const cl_MI operator[](unsigned int deg) const
185 list<Term>::const_iterator i = terms.begin(), end = terms.end();
186 for ( ; i != end; ++i ) {
187 if ( i->exp == deg ) {
190 if ( i->exp < deg ) {
196 void set(unsigned int deg, const cl_MI& c)
198 list<Term>::iterator i = terms.begin(), end = terms.end();
200 if ( i->exp == deg ) {
209 if ( i->exp < deg ) {
221 ex to_ex(const ex& x, bool symmetric = true) const
224 list<Term>::const_iterator i = terms.begin(), end = terms.end();
226 numeric mod(R->modulus);
227 numeric halfmod = (mod-1)/2;
228 for ( ; i != end; ++i ) {
229 numeric n(R->retract(i->c));
231 r += pow(x, i->exp) * (n-mod);
234 r += pow(x, i->exp) * n;
239 for ( ; i != end; ++i ) {
240 r += pow(x, i->exp) * numeric(R->retract(i->c));
247 if ( terms.size() ) {
248 if ( terms.front().c != R->one() ) {
249 list<Term>::iterator i = terms.begin(), end = terms.end();
252 while ( ++i != end ) {
253 i->c = div(i->c, cont);
263 return terms.front().c;
265 void divide(const cl_MI& x)
267 list<Term>::iterator i = terms.begin(), end = terms.end();
268 for ( ; i != end; ++i ) {
275 void divide(const cl_I& x)
277 list<Term>::iterator i = terms.begin(), end = terms.end();
278 for ( ; i != end; ++i ) {
279 i->c = cl_MI(R, the<cl_I>(R->retract(i->c) / x));
282 void reduce_exponents(unsigned int prime)
284 list<Term>::iterator i = terms.begin(), end = terms.end();
287 // assert: i->exp is multiple of prime
293 void deriv(UniPoly& d) const
295 list<Term>::const_iterator i = terms.begin(), end = terms.end();
298 cl_MI newc = i->c * i->exp;
299 if ( !zerop(newc) ) {
303 d.terms.push_back(t);
309 bool operator<(const UniPoly& o) const
311 if ( terms.size() != o.terms.size() ) {
312 return terms.size() < o.terms.size();
314 list<Term>::const_iterator i1 = terms.begin(), end = terms.end();
315 list<Term>::const_iterator i2 = o.terms.begin();
316 while ( i1 != end ) {
317 if ( i1->exp != i2->exp ) {
318 return i1->exp < i2->exp;
320 if ( i1->c != i2->c ) {
321 return R->retract(i1->c) < R->retract(i2->c);
327 bool operator==(const UniPoly& o) const
329 if ( terms.size() != o.terms.size() ) {
332 list<Term>::const_iterator i1 = terms.begin(), end = terms.end();
333 list<Term>::const_iterator i2 = o.terms.begin();
334 while ( i1 != end ) {
335 if ( i1->exp != i2->exp ) {
338 if ( i1->c != i2->c ) {
345 bool operator!=(const UniPoly& o) const
347 bool res = !(*this == o);
352 static UniPoly operator*(const UniPoly& a, const UniPoly& b)
354 unsigned int n = a.degree()+b.degree();
357 for ( unsigned int i=0 ; i<=n; ++i ) {
359 for ( unsigned int j=0 ; j<=i; ++j ) {
360 t.c = t.c + a[j] * b[i-j];
364 c.terms.push_front(t);
370 static UniPoly operator-(const UniPoly& a, const UniPoly& b)
372 list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
373 list<Term>::const_iterator ib = b.terms.begin(), bend = b.terms.end();
375 while ( ia != aend && ib != bend ) {
376 if ( ia->exp > ib->exp ) {
377 c.terms.push_back(*ia);
380 else if ( ia->exp < ib->exp ) {
381 c.terms.push_back(*ib);
382 c.terms.back().c = -c.terms.back().c;
390 c.terms.push_back(t);
395 while ( ia != aend ) {
396 c.terms.push_back(*ia);
399 while ( ib != bend ) {
400 c.terms.push_back(*ib);
401 c.terms.back().c = -c.terms.back().c;
407 static UniPoly operator*(const UniPoly& a, const cl_MI& fac)
409 unsigned int n = a.degree();
412 for ( unsigned int i=0 ; i<=n; ++i ) {
416 c.terms.push_front(t);
422 static UniPoly operator+(const UniPoly& a, const UniPoly& b)
424 list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
425 list<Term>::const_iterator ib = b.terms.begin(), bend = b.terms.end();
427 while ( ia != aend && ib != bend ) {
428 if ( ia->exp > ib->exp ) {
429 c.terms.push_back(*ia);
432 else if ( ia->exp < ib->exp ) {
433 c.terms.push_back(*ib);
441 c.terms.push_back(t);
446 while ( ia != aend ) {
447 c.terms.push_back(*ia);
450 while ( ib != bend ) {
451 c.terms.push_back(*ib);
457 // static UniPoly operator-(const UniPoly& a)
459 // list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
461 // while ( ia != aend ) {
462 // c.terms.push_back(*ia);
463 // c.terms.back().c = -c.terms.back().c;
470 ostream& operator<<(ostream& o, const UniPoly& t)
472 list<Term>::const_iterator i = t.terms.begin(), end = t.terms.end();
477 for ( ; i != end; ) {
485 #endif // def DEBUGFACTOR
488 ostream& operator<<(ostream& o, const list<UniPoly>& t)
490 list<UniPoly>::const_iterator i = t.begin(), end = t.end();
492 for ( ; i != end; ) {
498 #endif // def DEBUGFACTOR
500 typedef vector<UniPoly> UniPolyVec;
503 ostream& operator<<(ostream& o, const UniPolyVec& v)
505 UniPolyVec::const_iterator i = v.begin(), end = v.end();
507 o << *i++ << " , " << endl;
511 #endif // def DEBUGFACTOR
518 UniFactor(const cl_modint_ring& ring) : p(ring) { }
519 UniFactor(const UniPoly& p_, unsigned int exp_) : p(p_), exp(exp_) { }
520 bool operator<(const UniFactor& o) const
528 vector<UniFactor> factors;
532 sort(factors.begin(), factors.end());
533 if ( factors.size() > 1 ) {
534 vector<UniFactor>::iterator i = factors.begin();
535 vector<UniFactor>::const_iterator cmp = factors.begin()+1;
536 vector<UniFactor>::iterator end = factors.end();
537 while ( cmp != end ) {
538 if ( i->p != cmp->p ) {
548 factors.erase(i+1, end);
555 ostream& operator<<(ostream& o, const UniFactorVec& ufv)
557 for ( size_t i=0; i<ufv.factors.size(); ++i ) {
558 if ( i != ufv.factors.size()-1 ) {
564 o << "[ " << ufv.factors[i].p << " ]^" << ufv.factors[i].exp << endl;
568 #endif // def DEBUGFACTOR
570 static void rem(const UniPoly& a_, const UniPoly& b, UniPoly& c)
572 if ( a_.degree() < b.degree() ) {
590 cl_MI qk = div(c[n+k], b[n]);
593 for ( unsigned int i=0; i<n; ++i ) {
595 c.set(j, c[j] - qk*b[j-k]);
601 list<Term>::iterator i = c.terms.begin(), end = c.terms.end();
603 if ( i->exp <= n-1 ) {
608 c.terms.erase(c.terms.begin(), i);
611 static void div(const UniPoly& a_, const UniPoly& b, UniPoly& q)
613 if ( a_.degree() < b.degree() ) {
626 cl_MI qk = div(c[n+k], b[n]);
631 q.terms.push_back(t);
633 for ( unsigned int i=0; i<n; ++i ) {
635 c.set(j, c[j] - qk*b[j-k]);
643 static void gcd(const UniPoly& a, const UniPoly& b, UniPoly& c)
650 if ( c.degree() < d.degree() ) {
655 while ( !d.zero() ) {
664 static bool is_one(const UniPoly& w)
666 if ( w.terms.size() == 1 && w[0] == w.R->one() ) {
672 static void sqrfree_main(const UniPoly& a, UniFactorVec& fvec)
678 UniPoly c(a.R), w(a.R);
681 while ( !is_one(w) ) {
682 UniPoly y(a.R), z(a.R);
687 fvec.factors.push_back(uf);
696 unsigned int prime = cl_I_to_uint(c.R->modulus);
697 c.reduce_exponents(prime);
698 unsigned int pos = fvec.factors.size();
699 sqrfree_main(c, fvec);
700 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
701 fvec.factors[p].exp *= prime;
707 unsigned int prime = cl_I_to_uint(a.R->modulus);
709 amod.reduce_exponents(prime);
710 unsigned int pos = fvec.factors.size();
711 sqrfree_main(amod, fvec);
712 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
713 fvec.factors[p].exp *= prime;
719 static void squarefree(const UniPoly& a, UniFactorVec& fvec)
721 sqrfree_main(a, fvec);
727 friend ostream& operator<<(ostream& o, const Matrix& m);
729 Matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
733 size_t rowsize() const { return r; }
734 size_t colsize() const { return c; }
735 cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
736 cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
737 void mul_col(size_t col, const cl_MI x)
739 Vec::iterator i = m.begin() + col;
740 for ( size_t rc=0; rc<r; ++rc ) {
745 void sub_col(size_t col1, size_t col2, const cl_MI fac)
747 Vec::iterator i1 = m.begin() + col1;
748 Vec::iterator i2 = m.begin() + col2;
749 for ( size_t rc=0; rc<r; ++rc ) {
750 *i1 = *i1 - *i2 * fac;
755 void switch_col(size_t col1, size_t col2)
758 Vec::iterator i1 = m.begin() + col1;
759 Vec::iterator i2 = m.begin() + col2;
760 for ( size_t rc=0; rc<r; ++rc ) {
761 buf = *i1; *i1 = *i2; *i2 = buf;
766 void mul_row(size_t row, const cl_MI x)
768 vector<cl_MI>::iterator i = m.begin() + row*c;
769 for ( size_t cc=0; cc<c; ++cc ) {
774 void sub_row(size_t row1, size_t row2, const cl_MI fac)
776 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
777 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
778 for ( size_t cc=0; cc<c; ++cc ) {
779 *i1 = *i1 - *i2 * fac;
784 void switch_row(size_t row1, size_t row2)
787 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
788 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
789 for ( size_t cc=0; cc<c; ++cc ) {
790 buf = *i1; *i1 = *i2; *i2 = buf;
795 bool is_col_zero(size_t col) const
797 Vec::const_iterator i = m.begin() + col;
798 for ( size_t rr=0; rr<r; ++rr ) {
806 bool is_row_zero(size_t row) const
808 Vec::const_iterator i = m.begin() + row*c;
809 for ( size_t cc=0; cc<c; ++cc ) {
817 void set_row(size_t row, const vector<cl_MI>& newrow)
819 Vec::iterator i1 = m.begin() + row*c;
820 Vec::const_iterator i2 = newrow.begin(), end = newrow.end();
821 for ( ; i2 != end; ++i1, ++i2 ) {
825 Vec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
826 Vec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
833 Matrix operator*(const Matrix& m1, const Matrix& m2)
835 const unsigned int r = m1.rowsize();
836 const unsigned int c = m2.colsize();
837 Matrix o(r,c,m1(0,0));
839 for ( size_t i=0; i<r; ++i ) {
840 for ( size_t j=0; j<c; ++j ) {
842 buf = m1(i,0) * m2(0,j);
843 for ( size_t k=1; k<c; ++k ) {
844 buf = buf + m1(i,k)*m2(k,j);
852 ostream& operator<<(ostream& o, const Matrix& m)
854 vector<cl_MI>::const_iterator i = m.m.begin(), end = m.m.end();
856 for ( ; i != end; ++i ) {
858 if ( !(wrap++ % m.c) ) {
865 #endif // def DEBUGFACTOR
867 static void q_matrix(const UniPoly& a, Matrix& Q)
869 unsigned int n = a.degree();
870 unsigned int q = cl_I_to_uint(a.R->modulus);
872 // vector<cl_MI> r(n, a.R->zero());
873 // r[0] = a.R->one();
875 // unsigned int max = (n-1) * q;
876 // for ( size_t m=1; m<=max; ++m ) {
877 // cl_MI rn_1 = r.back();
878 // for ( size_t i=n-1; i>0; --i ) {
879 // r[i] = r[i-1] - rn_1 * a[i];
881 // r[0] = -rn_1 * a[0];
882 // if ( (m % q) == 0 ) {
883 // Q.set_row(m/q, r);
886 // slow and (hopefully) correct
887 for ( size_t i=0; i<n; ++i ) {
889 qk.set(i*q, a.R->one());
893 for ( size_t j=0; j<n; ++j ) {
894 rvec.push_back(r[j]);
900 static void nullspace(Matrix& M, vector<Vec>& basis)
902 const size_t n = M.rowsize();
903 const cl_MI one = M(0,0).ring()->one();
904 for ( size_t i=0; i<n; ++i ) {
905 M(i,i) = M(i,i) - one;
907 for ( size_t r=0; r<n; ++r ) {
909 for ( ; cc<n; ++cc ) {
910 if ( !zerop(M(r,cc)) ) {
912 if ( !zerop(M(cc,cc)) ) {
924 M.mul_col(r, recip(M(r,r)));
925 for ( cc=0; cc<n; ++cc ) {
927 M.sub_col(cc, r, M(r,cc));
933 for ( size_t i=0; i<n; ++i ) {
934 M(i,i) = M(i,i) - one;
936 for ( size_t i=0; i<n; ++i ) {
937 if ( !M.is_row_zero(i) ) {
938 Vec nu(M.row_begin(i), M.row_end(i));
944 static void berlekamp(const UniPoly& a, UniPolyVec& upv)
946 Matrix Q(a.degree(), a.degree(), a.R->zero());
950 const unsigned int k = nu.size();
955 list<UniPoly> factors;
956 factors.push_back(a);
957 unsigned int size = 1;
959 unsigned int q = cl_I_to_uint(a.R->modulus);
961 list<UniPoly>::iterator u = factors.begin();
964 for ( unsigned int s=0; s<q; ++s ) {
966 UniPoly nur(a.R, nu[r]);
967 nur.set(0, nur[0] - cl_MI(a.R, s));
969 if ( !is_one(g) && g != *u ) {
973 throw logic_error("berlekamp: unexpected divisor.");
978 factors.push_back(g);
980 list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
982 if ( i->degree() ) ++size;
986 list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
992 // if ( u->degree() < nur.degree() ) {
1004 static void factor_modular(const UniPoly& p, UniPolyVec& upv)
1010 static void exteuclid(const UniPoly& a, const UniPoly& b, UniPoly& g, UniPoly& s, UniPoly& t)
1012 if ( a.degree() < b.degree() ) {
1013 exteuclid(b, a, g, t, s);
1016 UniPoly c1(a.R), c2(a.R), d1(a.R), d2(a.R), q(a.R), r(a.R), r1(a.R), r2(a.R);
1017 UniPoly c = a; c.unit_normal();
1018 UniPoly d = b; d.unit_normal();
1019 c1.set(0, a.R->one());
1020 d2.set(0, a.R->one());
1021 while ( !d.zero() ) {
1034 g = c; g.unit_normal();
1043 static ex replace_lc(const ex& poly, const ex& x, const ex& lc)
1045 ex r = expand(poly + (lc - poly.lcoeff(x)) * pow(x, poly.degree(x)));
1049 static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const UniPoly& u1_, const UniPoly& w1_, const ex& gamma_ = 0)
1052 const cl_modint_ring& R = u1_.R;
1056 for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
1057 maxcoeff += pow(abs(a.coeff(x, i)),2);
1059 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
1060 cl_I maxdegree = (u1_.degree() > w1_.degree()) ? u1_.degree() : w1_.degree();
1061 cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
1064 ex alpha = a.lcoeff(x);
1069 numeric gamma_ui = ex_to<numeric>(abs(gamma));
1076 phi = expand(gamma * nu1.to_ex(x));
1077 UniPoly u1(R, phi, x);
1078 phi = expand(alpha * nw1.to_ex(x));
1079 UniPoly w1(R, phi, x);
1082 UniPoly s(R), t(R), g(R);
1083 exteuclid(u1, w1, g, s, t);
1086 ex u = replace_lc(u1.to_ex(x), x, gamma);
1087 ex w = replace_lc(w1.to_ex(x), x, alpha);
1088 ex e = expand(a - u * w);
1089 numeric modulus = p;
1090 const numeric maxmodulus = 2*numeric(B)*gamma_ui;
1093 while ( !e.is_zero() && modulus < maxmodulus ) {
1095 phi = expand(s.to_ex(x)*c);
1096 UniPoly sigmatilde(R, phi, x);
1097 phi = expand(t.to_ex(x)*c);
1098 UniPoly tautilde(R, phi, x);
1100 div(sigmatilde, w1, q);
1101 rem(sigmatilde, w1, r);
1103 phi = expand(tautilde.to_ex(x) + q.to_ex(x) * u1.to_ex(x));
1104 UniPoly tau(R, phi, x);
1105 u = expand(u + tau.to_ex(x) * modulus);
1106 w = expand(w + sigma.to_ex(x) * modulus);
1107 e = expand(a - u * w);
1108 modulus = modulus * p;
1112 if ( e.is_zero() ) {
1113 ex delta = u.content(x);
1115 w = w / gamma * delta;
1123 static unsigned int next_prime(unsigned int p)
1125 static vector<unsigned int> primes;
1126 if ( primes.size() == 0 ) {
1127 primes.push_back(3); primes.push_back(5); primes.push_back(7);
1129 vector<unsigned int>::const_iterator it = primes.begin();
1130 if ( p >= primes.back() ) {
1131 unsigned int candidate = primes.back() + 2;
1133 size_t n = primes.size()/2;
1134 for ( size_t i=0; i<n; ++i ) {
1135 if ( candidate % primes[i] ) continue;
1139 primes.push_back(candidate);
1140 if ( candidate > p ) break;
1144 vector<unsigned int>::const_iterator end = primes.end();
1145 for ( ; it!=end; ++it ) {
1150 throw logic_error("next_prime: should not reach this point!");
1156 Partition(size_t n_) : n(n_)
1162 int operator[](size_t i) const { return k[i]; }
1163 size_t size() const { return n; }
1164 size_t size_first() const { return n-sum; }
1165 size_t size_second() const { return sum; }
1169 for ( size_t i=0; i<k.size(); ++i ) {
1170 cout << k[i] << " ";
1177 for ( size_t i=n-1; i>=1; --i ) {
1193 static void split(const UniPolyVec& factors, const Partition& part, UniPoly& a, UniPoly& b)
1195 a.set(0, a.R->one());
1196 b.set(0, a.R->one());
1197 for ( size_t i=0; i<part.size(); ++i ) {
1213 static ex factor_univariate(const ex& poly, const ex& x)
1215 ex unit, cont, prim;
1216 poly.unitcontprim(x, unit, cont, prim);
1218 // determine proper prime
1220 cl_modint_ring R = find_modint_ring(p);
1222 if ( irem(ex_to<numeric>(prim.lcoeff(x)), p) != 0 ) {
1223 UniPoly modpoly(R, prim, x);
1224 UniFactorVec sqrfree_ufv;
1225 squarefree(modpoly, sqrfree_ufv);
1226 if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break;
1229 R = find_modint_ring(p);
1232 // do modular factorization
1233 UniPoly modpoly(R, prim, x);
1235 factor_modular(modpoly, factors);
1236 if ( factors.size() <= 1 ) {
1237 // irreducible for sure
1241 // lift all factor combinations
1242 stack<ModFactors> tocheck;
1245 mf.factors = factors;
1248 while ( tocheck.size() ) {
1249 const size_t n = tocheck.top().factors.size();
1253 split(tocheck.top().factors, part, a, b);
1255 ex answer = hensel_univar(tocheck.top().poly, x, p, a, b);
1256 if ( answer != lst() ) {
1257 if ( part.size_first() == 1 ) {
1258 if ( part.size_second() == 1 ) {
1259 result *= answer.op(0) * answer.op(1);
1263 result *= answer.op(0);
1264 tocheck.top().poly = answer.op(1);
1265 for ( size_t i=0; i<n; ++i ) {
1266 if ( part[i] == 0 ) {
1267 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1273 else if ( part.size_second() == 1 ) {
1274 if ( part.size_first() == 1 ) {
1275 result *= answer.op(0) * answer.op(1);
1279 result *= answer.op(1);
1280 tocheck.top().poly = answer.op(0);
1281 for ( size_t i=0; i<n; ++i ) {
1282 if ( part[i] == 1 ) {
1283 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1290 UniPolyVec newfactors1(part.size_first(), R), newfactors2(part.size_second(), R);
1291 UniPolyVec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1292 for ( size_t i=0; i<n; ++i ) {
1294 *i2++ = tocheck.top().factors[i];
1297 *i1++ = tocheck.top().factors[i];
1300 tocheck.top().factors = newfactors1;
1301 tocheck.top().poly = answer.op(0);
1303 mf.factors = newfactors2;
1304 mf.poly = answer.op(1);
1310 if ( !part.next() ) {
1311 result *= tocheck.top().poly;
1319 return unit * cont * result;
1328 // forward declaration
1329 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);
1331 UniPolyVec multiterm_eea_lift(const UniPolyVec& a, const ex& x, unsigned int p, unsigned int k)
1333 DCOUT(multiterm_eea_lift);
1338 const size_t r = a.size();
1340 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1342 UniPolyVec q(r-1, fill);
1344 for ( size_t j=r-2; j>=1; --j ) {
1345 q[j-1] = a[j] * q[j];
1349 beta.set(0, R->one());
1351 for ( size_t j=1; j<r; ++j ) {
1354 vector<ex> mdarg(2);
1355 mdarg[0] = q[j-1].to_ex(x);
1356 mdarg[1] = a[j-1].to_ex(x);
1357 vector<EvalPoint> empty;
1358 vector<ex> exsigma = multivar_diophant(mdarg, x, beta.to_ex(x), empty, 0, p, k);
1359 UniPoly sigma1(R, exsigma[0], x);
1360 UniPoly sigma2(R, exsigma[1], x);
1362 s.push_back(sigma2);
1367 DCOUT(END multiterm_eea_lift);
1371 void eea_lift(const UniPoly& a, const UniPoly& b, const ex& x, unsigned int p, unsigned int k, UniPoly& s_, UniPoly& t_)
1380 cl_modint_ring R = find_modint_ring(p);
1386 UniPoly smod(R), tmod(R), g(R);
1387 exteuclid(amod, bmod, g, smod, tmod);
1393 cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1394 UniPoly s(Rpk, smod);
1395 UniPoly t(Rpk, tmod);
1402 one.set(0, Rpk->one());
1403 for ( size_t j=1; j<k; ++j ) {
1404 UniPoly e = one - a * s - b * t;
1407 UniPoly sigmabar(R);
1408 sigmabar = smod * c;
1412 div(sigmabar, bmod, q);
1414 rem(sigmabar, bmod, sigma);
1416 tau = taubar + q * amod;
1417 UniPoly sadd(Rpk, sigma);
1418 cl_MI modmodulus(Rpk, modulus);
1419 s = s + sadd * modmodulus;
1420 UniPoly tadd(Rpk, tau);
1421 t = t + tadd * modmodulus;
1422 modulus = modulus * p;
1429 DCOUT2(check, a*s + b*t);
1430 DCOUT(END eea_lift);
1433 UniPolyVec univar_diophant(const UniPolyVec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1435 DCOUT(univar_diophant);
1442 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1444 const size_t r = a.size();
1447 UniPolyVec s = multiterm_eea_lift(a, x, p, k);
1448 for ( size_t j=0; j<r; ++j ) {
1449 ex phi = expand(pow(x,m)*s[j].to_ex(x));
1450 UniPoly bmod(R, phi, x);
1452 rem(bmod, a[j], buf);
1453 result.push_back(buf);
1458 eea_lift(a[1], a[0], x, p, k, s, t);
1459 ex phi = expand(pow(x,m)*s.to_ex(x));
1460 UniPoly bmod(R, phi, x);
1462 rem(bmod, a[0], buf);
1463 result.push_back(buf);
1466 phi = expand(pow(x,m)*t.to_ex(x));
1467 UniPoly t1mod(R, phi, x);
1468 buf = t1mod + q * a[1];
1469 result.push_back(buf);
1473 DCOUT(END univar_diophant);
1477 struct make_modular_map : public map_function {
1479 make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1480 ex operator()(const ex& e)
1482 if ( is_a<add>(e) || is_a<mul>(e) ) {
1483 return e.map(*this);
1485 else if ( is_a<numeric>(e) ) {
1486 numeric mod(R->modulus);
1487 numeric halfmod = (mod-1)/2;
1488 cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1489 numeric n(R->retract(emod));
1490 if ( n > halfmod ) {
1501 static ex make_modular(const ex& e, const cl_modint_ring& R)
1503 make_modular_map map(R);
1507 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)
1511 DCOUT(multivar_diophant);
1514 for ( size_t i=0; i<a.size(); ++i ) {
1515 cout << a[i] << " ";
1523 for ( size_t i=0; i<I.size(); ++i ) {
1524 cout << I[i].x << "=" << I[i].evalpoint << " ";
1532 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1533 const size_t r = a.size();
1534 const size_t nu = I.size() + 1;
1540 ex xnu = I.back().x;
1541 int alphanu = I.back().evalpoint;
1544 for ( size_t i=0; i<r; ++i ) {
1548 for ( size_t i=0; i<r; ++i ) {
1549 b[i] = normal(A / a[i]);
1552 vector<ex> anew = a;
1553 for ( size_t i=0; i<r; ++i ) {
1554 anew[i] = anew[i].subs(xnu == alphanu);
1556 ex cnew = c.subs(xnu == alphanu);
1557 vector<EvalPoint> Inew = I;
1559 sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1563 for ( size_t i=0; i<r; ++i ) {
1564 buf -= sigma[i] * b[i];
1567 e = make_modular(e, R);
1573 for ( size_t m=1; m<=d; ++m ) {
1575 while ( !e.is_zero() && e.has(xnu) ) {
1576 monomial *= (xnu - alphanu);
1577 monomial = expand(monomial);
1580 ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1582 if ( !cm.is_zero() ) {
1583 vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1586 for ( size_t j=0; j<delta_s.size(); ++j ) {
1587 delta_s[j] *= monomial;
1588 sigma[j] += delta_s[j];
1589 buf -= delta_s[j] * b[j];
1592 e = make_modular(e, R);
1598 DCOUT(uniterm left);
1600 for ( size_t i=0; i<a.size(); ++i ) {
1601 UniPoly up(R, a[i], x);
1605 sigma.insert(sigma.begin(), r, 0);
1608 if ( is_a<add>(c) ) {
1617 for ( size_t i=0; i<nterms; ++i ) {
1619 int m = z.degree(x);
1621 cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1623 UniPolyVec delta_s = univar_diophant(amod, x, m, p, k);
1626 while ( poscm < 0 ) {
1627 poscm = poscm + expt_pos(cl_I(p),k);
1629 modcm = cl_MI(R, poscm);
1631 for ( size_t j=0; j<delta_s.size(); ++j ) {
1632 delta_s[j] = delta_s[j] * modcm;
1633 sigma[j] = sigma[j] + delta_s[j].to_ex(x);
1637 cout << "STEP " << i << " sigma ";
1638 for ( size_t p=0; p<sigma.size(); ++p ) {
1639 cout << sigma[p] << " ";
1650 for ( size_t i=0; i<sigma.size(); ++i ) {
1651 cout << sigma[i] << " ";
1656 for ( size_t i=0; i<sigma.size(); ++i ) {
1657 sigma[i] = make_modular(sigma[i], R);
1662 for ( size_t i=0; i<sigma.size(); ++i ) {
1663 cout << sigma[i] << " ";
1667 DCOUT(END multivar_diophant);
1672 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1674 for ( size_t i=0; i<v.size(); ++i ) {
1675 o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1679 #endif // def DEBUGFACTOR
1682 ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I, unsigned int p, const cl_I& l, const UniPolyVec& u, const vector<ex>& lcU)
1684 DCOUT(hensel_multivar);
1692 const size_t nu = I.size() + 1;
1693 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1700 for ( size_t j=nu; j>=2; --j ) {
1702 int alpha = I[j-2].evalpoint;
1703 A[j-2] = A[j-1].subs(x==alpha);
1704 A[j-2] = make_modular(A[j-2], R);
1709 for ( size_t i=0; i<A.size(); ++i) cout << A[i] << " ";
1713 int maxdeg = a.degree(I.front().x);
1714 for ( size_t i=1; i<I.size(); ++i ) {
1715 int maxdeg2 = a.degree(I[i].x);
1716 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
1720 const size_t n = u.size();
1723 for ( size_t i=0; i<n; ++i ) {
1724 U[i] = u[i].to_ex(x);
1728 for ( size_t i=0; i<U.size(); ++i) cout << U[i] << " ";
1732 for ( size_t j=2; j<=nu; ++j ) {
1737 for ( size_t m=0; m<n; ++m) {
1738 if ( lcU[m] != 1 ) {
1740 for ( size_t i=j-1; i<nu-1; ++i ) {
1741 coef = coef.subs(I[i].x == I[i].evalpoint);
1743 coef = expand(coef);
1744 coef = make_modular(coef, R);
1745 int deg = U[m].degree(x);
1746 U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
1751 for ( size_t i=0; i<n; ++i ) {
1754 ex e = expand(A[j-1] - Uprod);
1757 vector<EvalPoint> newI;
1758 for ( size_t i=1; i<=j-2; ++i ) {
1759 newI.push_back(I[i-1]);
1764 int alphaj = I[j-2].evalpoint;
1765 size_t deg = A[j-1].degree(xj);
1767 for ( size_t k=1; k<=deg; ++k ) {
1769 if ( !e.is_zero() ) {
1772 monomial *= (xj - alphaj);
1773 monomial = expand(monomial);
1775 ex dif = e.diff(ex_to<symbol>(xj), k);
1777 ex c = dif.subs(xj==alphaj) / factorial(k);
1779 if ( !c.is_zero() ) {
1780 vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
1781 for ( size_t i=0; i<n; ++i ) {
1783 DCOUTVAR(deltaU[i]);
1784 deltaU[i] *= monomial;
1786 U[i] = make_modular(U[i], R);
1787 U[i] = U[i].expand();
1791 for ( size_t i=0; i<n; ++i ) {
1794 DCOUTVAR(Uprod.expand());
1796 e = expand(A[j-1] - Uprod);
1797 e = make_modular(e, R);
1808 for ( size_t i=0; i<U.size(); ++i ) {
1812 if ( expand(a-acand).is_zero() ) {
1814 for ( size_t i=0; i<U.size(); ++i ) {
1818 DCOUT(END hensel_multivar);
1824 DCOUT(END hensel_multivar);
1829 static ex put_factors_into_lst(const ex& e)
1831 DCOUT(put_factors_into_lst);
1836 if ( is_a<numeric>(e) ) {
1838 DCOUT(END put_factors_into_lst);
1842 if ( is_a<power>(e) ) {
1844 result.append(e.op(0));
1845 result.append(e.op(1));
1846 DCOUT(END put_factors_into_lst);
1850 if ( is_a<symbol>(e) || is_a<add>(e) ) {
1854 DCOUT(END put_factors_into_lst);
1858 if ( is_a<mul>(e) ) {
1860 for ( size_t i=0; i<e.nops(); ++i ) {
1862 if ( is_a<numeric>(op) ) {
1865 if ( is_a<power>(op) ) {
1866 result.append(op.op(0));
1867 result.append(op.op(1));
1869 if ( is_a<symbol>(op) || is_a<add>(op) ) {
1874 result.prepend(nfac);
1875 DCOUT(END put_factors_into_lst);
1879 throw runtime_error("put_factors_into_lst: bad term.");
1883 ostream& operator<<(ostream& o, const vector<numeric>& v)
1885 for ( size_t i=0; i<v.size(); ++i ) {
1890 #endif // def DEBUGFACTOR
1892 static bool checkdivisors(const lst& f, vector<numeric>& d)
1894 DCOUT(checkdivisors);
1895 const int k = f.nops()-2;
1899 d[0] = ex_to<numeric>(f.op(0) * f.op(f.nops()-1));
1900 if ( d[0] == 1 && k == 1 && abs(f.op(1)) != 1 ) {
1902 DCOUT(END checkdivisors);
1906 for ( int i=1; i<=k; ++i ) {
1908 DCOUTVAR(abs(f.op(i)));
1909 q = ex_to<numeric>(abs(f.op(i)));
1911 for ( int j=i-1; j>=0; --j ) {
1922 DCOUT(END checkdivisors);
1929 DCOUT(END checkdivisors);
1933 static bool generate_set(const ex& u, const ex& vn, const exset& syms, const ex& f, const numeric& modulus, vector<numeric>& a, vector<numeric>& d)
1935 // computation of d is actually not necessary
1936 DCOUT(generate_set);
1941 const ex& x = *syms.begin();
1947 exset::const_iterator s = syms.begin();
1949 for ( size_t i=0; i<a.size(); ++i ) {
1952 a[i] = mod(numeric(rand()), 2*modulus) - modulus;
1953 vnatry = vna.subs(*s == a[i]);
1954 } while ( vnatry == 0 );
1956 u0 = u0.subs(*s == a[i]);
1961 if ( gcd(u0,u0.diff(ex_to<symbol>(x))) != 1 ) {
1964 if ( is_a<numeric>(vn) ) {
1968 DCOUT(do substitution);
1970 lst::const_iterator i = ex_to<lst>(f).begin();
1972 bool problem = false;
1973 while ( i!=ex_to<lst>(f).end() ) {
1975 if ( !is_a<numeric>(fs) ) {
1978 for ( size_t j=0; j<a.size(); ++j ) {
1979 fs = fs.subs(*s == a[j]);
1982 if ( abs(fs) == 1 ) {
1993 ex con = u0.content(x);
1996 trying = checkdivisors(fnum, d);
1999 DCOUT(END generate_set);
2003 static ex factor_multivariate(const ex& poly, const exset& syms)
2005 DCOUT(factor_multivariate);
2008 exset::const_iterator s;
2009 const ex& x = *syms.begin();
2012 /* make polynomial primitive */
2013 ex p = poly.expand().collect(x);
2015 ex cont = p.lcoeff(x);
2016 for ( numeric i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
2017 cont = gcd(cont, p.coeff(x,ex_to<numeric>(i).to_int()));
2018 if ( cont == 1 ) break;
2021 ex pp = expand(normal(p / cont));
2023 if ( !is_a<numeric>(cont) ) {
2025 return ::factor(cont) * ::factor(pp);
2027 return factor(cont) * factor(pp);
2031 /* factor leading coefficient */
2033 ex vn = pp.lcoeff(x);
2036 if ( is_a<numeric>(vn) ) {
2041 ex vnfactors = ::factor(vn);
2043 ex vnfactors = factor(vn);
2045 vnlst = put_factors_into_lst(vnfactors);
2049 const numeric maxtrials = 3;
2050 numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3;
2052 numeric minimalr = -1;
2053 vector<numeric> a(syms.size()-1, 0);
2054 vector<numeric> d((vnlst.nops()-1)/2+1, 0);
2057 numeric trialcount = 0;
2060 size_t factor_count;
2063 while ( trialcount < maxtrials ) {
2064 bool problem = generate_set(pp, vn, syms, vnlst, modulus, a, d);
2075 for ( size_t i=0; i<a.size(); ++i ) {
2076 u = u.subs(*s == a[i]);
2079 delta = u.content(x);
2082 // determine proper prime
2085 cl_modint_ring R = find_modint_ring(prime);
2086 DCOUTVAR(u.lcoeff(x));
2088 if ( irem(ex_to<numeric>(u.lcoeff(x)), prime) != 0 ) {
2089 UniPoly modpoly(R, u, x);
2090 UniFactorVec sqrfree_ufv;
2091 squarefree(modpoly, sqrfree_ufv);
2092 DCOUTVAR(sqrfree_ufv);
2093 if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break;
2095 prime = next_prime(prime);
2097 R = find_modint_ring(prime);
2106 ufaclst = put_factors_into_lst(ufac);
2108 factor_count = (ufaclst.nops()-1)/2;
2109 DCOUTVAR(factor_count);
2111 if ( factor_count <= 1 ) {
2113 DCOUT(END factor_multivariate);
2117 if ( minimalr < 0 ) {
2118 minimalr = factor_count;
2120 else if ( minimalr == factor_count ) {
2124 else if ( minimalr > factor_count ) {
2125 minimalr = factor_count;
2128 DCOUTVAR(trialcount);
2130 if ( minimalr <= 1 ) {
2132 DCOUT(END factor_multivariate);
2137 vector<numeric> ftilde((vnlst.nops()-1)/2+1);
2138 ftilde[0] = ex_to<numeric>(vnlst.op(0));
2139 for ( size_t i=1; i<ftilde.size(); ++i ) {
2140 ex ft = vnlst.op((i-1)*2+1);
2143 for ( size_t j=0; j<a.size(); ++j ) {
2144 ft = ft.subs(*s == a[j]);
2147 ftilde[i] = ex_to<numeric>(ft);
2151 vector<bool> used_flag((vnlst.nops()-1)/2+1, false);
2152 vector<ex> D(factor_count, 1);
2153 for ( size_t i=0; i<=factor_count; ++i ) {
2157 prefac = ex_to<numeric>(ufaclst.op(0));
2158 ftilde[0] = ftilde[0] / prefac;
2159 vnlst.let_op(0) = vnlst.op(0) / prefac;
2163 prefac = ex_to<numeric>(ufaclst.op(2*(i-1)+1).lcoeff(x));
2166 for ( size_t j=(vnlst.nops()-1)/2+1; j>0; --j ) {
2169 DCOUTVAR(ftilde[j-1]);
2170 if ( abs(ftilde[j-1]) == 1 ) {
2171 used_flag[j-1] = true;
2174 numeric g = gcd(prefac, ftilde[j-1]);
2177 DCOUT(has_common_prime);
2178 prefac = prefac / g;
2179 numeric count = abs(iquo(g, ftilde[j-1]));
2181 used_flag[j-1] = true;
2184 D[i-1] = D[i-1] * pow(vnlst.op(0), count);
2187 D[i-1] = D[i-1] * pow(vnlst.op(2*(j-2)+1), count);
2191 ftilde[j-1] = ftilde[j-1] / prefac;
2193 DCOUTVAR(ftilde[j-1]);
2202 bool some_factor_unused = false;
2203 for ( size_t i=0; i<used_flag.size(); ++i ) {
2204 if ( !used_flag[i] ) {
2205 some_factor_unused = true;
2209 if ( some_factor_unused ) {
2210 DCOUT(some factor unused!);
2214 vector<ex> C(factor_count);
2218 for ( size_t i=0; i<D.size(); ++i ) {
2222 for ( size_t j=0; j<a.size(); ++j ) {
2223 Dtilde = Dtilde.subs(*s == a[j]);
2227 C[i] = D[i] * (ufaclst.op(2*i+1).lcoeff(x) / Dtilde);
2231 for ( size_t i=0; i<D.size(); ++i ) {
2235 for ( size_t j=0; j<a.size(); ++j ) {
2236 Dtilde = Dtilde.subs(*s == a[j]);
2244 ui = ufaclst.op(2*(i-1)+1);
2247 ex d = gcd(ui.lcoeff(x), Dtilde);
2248 C[i] = D[i] * ( ui.lcoeff(x) / d );
2249 ui = ui * ( Dtilde[i] / d );
2250 delta = delta / ( Dtilde[i] / d );
2251 if ( delta == 1 ) break;
2253 C[i] = delta * C[i];
2254 pp = pp * pow(delta, D.size()-1);
2261 vector<EvalPoint> epv;
2264 for ( size_t i=0; i<a.size(); ++i ) {
2266 ep.evalpoint = a[i].to_int();
2273 for ( int i=u.degree(x); i>=u.ldegree(x); --i ) {
2274 maxcoeff += pow(abs(u.coeff(x, i)),2);
2276 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
2277 unsigned int maxdegree = 0;
2278 for ( size_t i=0; i<factor_count; ++i ) {
2279 if ( ufaclst[2*i+1].degree(x) > (int)maxdegree ) {
2280 maxdegree = ufaclst[2*i+1].degree(x);
2283 cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
2292 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
2293 for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) {
2294 UniPoly newu(R, ufaclst.op(i*2+1), x);
2295 uvec.push_back(newu);
2299 ex res = hensel_multivar(ufaclst.op(0)*pp, x, epv, prime, l, uvec, C);
2300 if ( res != lst() ) {
2301 ex result = cont * ufaclst.op(0);
2302 for ( size_t i=0; i<res.nops(); ++i ) {
2303 result *= res.op(i).content(x) * res.op(i).unit(x);
2304 result *= res.op(i).primpart(x);
2307 DCOUT(END factor_multivariate);
2313 struct find_symbols_map : public map_function {
2315 ex operator()(const ex& e)
2317 if ( is_a<symbol>(e) ) {
2321 return e.map(*this);
2325 static ex factor_sqrfree(const ex& poly)
2327 // determine all symbols in poly
2328 find_symbols_map findsymbols;
2330 if ( findsymbols.syms.size() == 0 ) {
2334 if ( findsymbols.syms.size() == 1 ) {
2336 const ex& x = *(findsymbols.syms.begin());
2337 if ( poly.ldegree(x) > 0 ) {
2338 int ld = poly.ldegree(x);
2339 ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2340 return res * pow(x,ld);
2343 ex res = factor_univariate(poly, x);
2348 // multivariate case
2349 ex res = factor_multivariate(poly, findsymbols.syms);
2353 struct apply_factor_map : public map_function {
2355 apply_factor_map(unsigned options_) : options(options_) { }
2356 ex operator()(const ex& e)
2358 if ( e.info(info_flags::polynomial) ) {
2360 return ::factor(e, options);
2362 return factor(e, options);
2365 if ( is_a<add>(e) ) {
2367 for ( size_t i=0; i<e.nops(); ++i ) {
2368 if ( e.op(i).info(info_flags::polynomial) ) {
2378 return ::factor(s1, options) + s2.map(*this);
2380 return factor(s1, options) + s2.map(*this);
2383 return e.map(*this);
2387 } // anonymous namespace
2390 ex factor(const ex& poly, unsigned options = 0)
2392 ex factor(const ex& poly, unsigned options)
2396 if ( !poly.info(info_flags::polynomial) ) {
2397 if ( options & factor_options::all ) {
2398 options &= ~factor_options::all;
2399 apply_factor_map factor_map(options);
2400 return factor_map(poly);
2405 // determine all symbols in poly
2406 find_symbols_map findsymbols;
2408 if ( findsymbols.syms.size() == 0 ) {
2412 exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
2413 for ( ; i!=end; ++i ) {
2417 // make poly square free
2418 ex sfpoly = sqrfree(poly, syms);
2420 // factorize the square free components
2421 if ( is_a<power>(sfpoly) ) {
2422 // case: (polynomial)^exponent
2423 const ex& base = sfpoly.op(0);
2424 if ( !is_a<add>(base) ) {
2425 // simple case: (monomial)^exponent
2428 ex f = factor_sqrfree(base);
2429 return pow(f, sfpoly.op(1));
2431 if ( is_a<mul>(sfpoly) ) {
2432 // case: multiple factors
2434 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
2435 const ex& t = sfpoly.op(i);
2436 if ( is_a<power>(t) ) {
2437 const ex& base = t.op(0);
2438 if ( !is_a<add>(base) ) {
2442 ex f = factor_sqrfree(base);
2443 res *= pow(f, t.op(1));
2446 else if ( is_a<add>(t) ) {
2447 ex f = factor_sqrfree(t);
2456 if ( is_a<symbol>(sfpoly) ) {
2459 // case: (polynomial)
2460 ex f = factor_sqrfree(sfpoly);
2464 } // namespace GiNaC