3 * Polynomial factorization routines.
4 * Only univariate at the moment and completely non-optimized!
8 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
10 * This program is free software; you can redistribute it and/or modify
11 * it under the terms of the GNU General Public License as published by
12 * the Free Software Foundation; either version 2 of the License, or
13 * (at your option) any later version.
15 * This program is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
18 * GNU General Public License for more details.
20 * You should have received a copy of the GNU General Public License
21 * along with this program; if not, write to the Free Software
22 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
29 #include "operators.h"
32 #include "relational.h"
50 #endif // def DEBUGFACTOR
56 typedef vector<cl_MI> Vec;
57 typedef vector<Vec> VecVec;
60 ostream& operator<<(ostream& o, const Vec& v)
62 Vec::const_iterator i = v.begin(), end = v.end();
68 #endif // def DEBUGFACTOR
71 ostream& operator<<(ostream& o, const VecVec& v)
73 VecVec::const_iterator i = v.begin(), end = v.end();
79 #endif // def DEBUGFACTOR
83 cl_MI c; // coefficient
84 unsigned int exp; // exponent >=0
88 ostream& operator<<(ostream& o, const Term& t)
91 o << "(" << t.c << ")x^" << t.exp;
94 o << "(" << t.c << ")";
98 #endif // def DEBUGFACTOR
103 list<Term> terms; // highest exponent first
105 UniPoly(const cl_modint_ring& ring) : R(ring) { }
106 UniPoly(const cl_modint_ring& ring, const ex& poly, const ex& x) : R(ring)
108 // assert: poly is in Z[x]
110 for ( int i=poly.degree(x); i>=poly.ldegree(x); --i ) {
111 int coeff = ex_to<numeric>(poly.coeff(x,i)).to_int();
113 t.c = R->canonhom(coeff);
121 UniPoly(const cl_modint_ring& ring, const Vec& v) : R(ring)
124 for ( unsigned int i=0; i<v.size(); ++i ) {
125 if ( !zerop(v[i]) ) {
132 unsigned int degree() const
134 if ( terms.size() ) {
135 return terms.front().exp;
141 bool zero() const { return (terms.size() == 0); }
142 const cl_MI operator[](unsigned int deg) const
144 list<Term>::const_iterator i = terms.begin(), end = terms.end();
145 for ( ; i != end; ++i ) {
146 if ( i->exp == deg ) {
149 if ( i->exp < deg ) {
155 void set(unsigned int deg, const cl_MI& c)
157 list<Term>::iterator i = terms.begin(), end = terms.end();
159 if ( i->exp == deg ) {
168 if ( i->exp < deg ) {
180 ex to_ex(const ex& x, bool symmetric = true) const
183 list<Term>::const_iterator i = terms.begin(), end = terms.end();
185 numeric mod(R->modulus);
186 numeric halfmod = (mod-1)/2;
187 for ( ; i != end; ++i ) {
188 numeric n(R->retract(i->c));
190 r += pow(x, i->exp) * (n-mod);
193 r += pow(x, i->exp) * n;
198 for ( ; i != end; ++i ) {
199 r += pow(x, i->exp) * numeric(R->retract(i->c));
206 if ( terms.size() ) {
207 if ( terms.front().c != R->one() ) {
208 list<Term>::iterator i = terms.begin(), end = terms.end();
211 while ( ++i != end ) {
212 i->c = div(i->c, cont);
222 return terms.front().c;
224 void divide(const cl_MI& x)
226 list<Term>::iterator i = terms.begin(), end = terms.end();
227 for ( ; i != end; ++i ) {
234 void reduce_exponents(unsigned int prime)
236 list<Term>::iterator i = terms.begin(), end = terms.end();
239 // assert: i->exp is multiple of prime
245 void deriv(UniPoly& d) const
247 list<Term>::const_iterator i = terms.begin(), end = terms.end();
250 cl_MI newc = i->c * i->exp;
251 if ( !zerop(newc) ) {
255 d.terms.push_back(t);
261 bool operator<(const UniPoly& o) const
263 if ( terms.size() != o.terms.size() ) {
264 return terms.size() < o.terms.size();
266 list<Term>::const_iterator i1 = terms.begin(), end = terms.end();
267 list<Term>::const_iterator i2 = o.terms.begin();
268 while ( i1 != end ) {
269 if ( i1->exp != i2->exp ) {
270 return i1->exp < i2->exp;
272 if ( i1->c != i2->c ) {
273 return R->retract(i1->c) < R->retract(i2->c);
279 bool operator==(const UniPoly& o) const
281 if ( terms.size() != o.terms.size() ) {
284 list<Term>::const_iterator i1 = terms.begin(), end = terms.end();
285 list<Term>::const_iterator i2 = o.terms.begin();
286 while ( i1 != end ) {
287 if ( i1->exp != i2->exp ) {
290 if ( i1->c != i2->c ) {
297 bool operator!=(const UniPoly& o) const
299 bool res = !(*this == o);
304 static UniPoly operator*(const UniPoly& a, const UniPoly& b)
306 unsigned int n = a.degree()+b.degree();
309 for ( unsigned int i=0 ; i<=n; ++i ) {
311 for ( unsigned int j=0 ; j<=i; ++j ) {
312 t.c = t.c + a[j] * b[i-j];
316 c.terms.push_front(t);
322 static UniPoly operator-(const UniPoly& a, const UniPoly& b)
324 list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
325 list<Term>::const_iterator ib = b.terms.begin(), bend = b.terms.end();
327 while ( ia != aend && ib != bend ) {
328 if ( ia->exp > ib->exp ) {
329 c.terms.push_back(*ia);
332 else if ( ia->exp < ib->exp ) {
333 c.terms.push_back(*ib);
334 c.terms.back().c = -c.terms.back().c;
342 c.terms.push_back(t);
347 while ( ia != aend ) {
348 c.terms.push_back(*ia);
351 while ( ib != bend ) {
352 c.terms.push_back(*ib);
353 c.terms.back().c = -c.terms.back().c;
359 static UniPoly operator-(const UniPoly& a)
361 list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
363 while ( ia != aend ) {
364 c.terms.push_back(*ia);
365 c.terms.back().c = -c.terms.back().c;
372 ostream& operator<<(ostream& o, const UniPoly& t)
374 list<Term>::const_iterator i = t.terms.begin(), end = t.terms.end();
379 for ( ; i != end; ) {
387 #endif // def DEBUGFACTOR
390 ostream& operator<<(ostream& o, const list<UniPoly>& t)
392 list<UniPoly>::const_iterator i = t.begin(), end = t.end();
394 for ( ; i != end; ) {
400 #endif // def DEBUGFACTOR
402 typedef vector<UniPoly> UniPolyVec;
409 UniFactor(const cl_modint_ring& ring) : p(ring) { }
410 UniFactor(const UniPoly& p_, unsigned int exp_) : p(p_), exp(exp_) { }
411 bool operator<(const UniFactor& o) const
419 vector<UniFactor> factors;
423 sort(factors.begin(), factors.end());
424 if ( factors.size() > 1 ) {
425 vector<UniFactor>::iterator i = factors.begin();
426 vector<UniFactor>::const_iterator cmp = factors.begin()+1;
427 vector<UniFactor>::iterator end = factors.end();
428 while ( cmp != end ) {
429 if ( i->p != cmp->p ) {
439 factors.erase(i+1, end);
446 ostream& operator<<(ostream& o, const UniFactorVec& ufv)
448 for ( size_t i=0; i<ufv.factors.size(); ++i ) {
449 if ( i != ufv.factors.size()-1 ) {
455 o << "[ " << ufv.factors[i].p << " ]^" << ufv.factors[i].exp << endl;
459 #endif // def DEBUGFACTOR
461 static void rem(const UniPoly& a_, const UniPoly& b, UniPoly& c)
463 if ( a_.degree() < b.degree() ) {
481 cl_MI qk = div(c[n+k], b[n]);
484 for ( unsigned int i=0; i<n; ++i ) {
486 c.set(j, c[j] - qk*b[j-k]);
492 list<Term>::iterator i = c.terms.begin(), end = c.terms.end();
494 if ( i->exp <= n-1 ) {
499 c.terms.erase(c.terms.begin(), i);
502 static void div(const UniPoly& a_, const UniPoly& b, UniPoly& q)
504 if ( a_.degree() < b.degree() ) {
517 cl_MI qk = div(c[n+k], b[n]);
522 q.terms.push_back(t);
524 for ( unsigned int i=0; i<n; ++i ) {
526 c.set(j, c[j] - qk*b[j-k]);
534 static void gcd(const UniPoly& a, const UniPoly& b, UniPoly& c)
541 if ( c.degree() < d.degree() ) {
546 while ( !d.zero() ) {
555 static bool is_one(const UniPoly& w)
557 if ( w.terms.size() == 1 && w[0] == w.R->one() ) {
563 static void sqrfree_main(const UniPoly& a, UniFactorVec& fvec)
569 UniPoly c(a.R), w(a.R);
572 while ( !is_one(w) ) {
573 UniPoly y(a.R), z(a.R);
578 fvec.factors.push_back(uf);
587 unsigned int prime = cl_I_to_uint(c.R->modulus);
588 c.reduce_exponents(prime);
589 unsigned int pos = fvec.factors.size();
590 sqrfree_main(c, fvec);
591 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
592 fvec.factors[p].exp *= prime;
598 unsigned int prime = cl_I_to_uint(a.R->modulus);
600 amod.reduce_exponents(prime);
601 unsigned int pos = fvec.factors.size();
602 sqrfree_main(amod, fvec);
603 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
604 fvec.factors[p].exp *= prime;
610 static void squarefree(const UniPoly& a, UniFactorVec& fvec)
612 sqrfree_main(a, fvec);
618 friend ostream& operator<<(ostream& o, const Matrix& m);
620 Matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
624 size_t rowsize() const { return r; }
625 size_t colsize() const { return c; }
626 cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
627 cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
628 void mul_col(size_t col, const cl_MI x)
630 Vec::iterator i = m.begin() + col;
631 for ( size_t rc=0; rc<r; ++rc ) {
636 void sub_col(size_t col1, size_t col2, const cl_MI fac)
638 Vec::iterator i1 = m.begin() + col1;
639 Vec::iterator i2 = m.begin() + col2;
640 for ( size_t rc=0; rc<r; ++rc ) {
641 *i1 = *i1 - *i2 * fac;
646 void switch_col(size_t col1, size_t col2)
649 Vec::iterator i1 = m.begin() + col1;
650 Vec::iterator i2 = m.begin() + col2;
651 for ( size_t rc=0; rc<r; ++rc ) {
652 buf = *i1; *i1 = *i2; *i2 = buf;
657 bool is_row_zero(size_t row) const
659 Vec::const_iterator i = m.begin() + row*c;
660 for ( size_t cc=0; cc<c; ++cc ) {
668 void set_row(size_t row, const vector<cl_MI>& newrow)
670 Vec::iterator i1 = m.begin() + row*c;
671 Vec::const_iterator i2 = newrow.begin(), end = newrow.end();
672 for ( ; i2 != end; ++i1, ++i2 ) {
676 Vec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
677 Vec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
684 ostream& operator<<(ostream& o, const Matrix& m)
686 vector<cl_MI>::const_iterator i = m.m.begin(), end = m.m.end();
688 for ( ; i != end; ++i ) {
690 if ( !(wrap++ % m.c) ) {
697 #endif // def DEBUGFACTOR
699 static void q_matrix(const UniPoly& a, Matrix& Q)
701 unsigned int n = a.degree();
702 unsigned int q = cl_I_to_uint(a.R->modulus);
704 // vector<cl_MI> r(n, a.R->zero());
705 // r[0] = a.R->one();
707 // unsigned int max = (n-1) * q;
708 // for ( size_t m=1; m<=max; ++m ) {
709 // cl_MI rn_1 = r.back();
710 // for ( size_t i=n-1; i>0; --i ) {
711 // r[i] = r[i-1] - rn_1 * a[i];
713 // r[0] = -rn_1 * a[0];
714 // if ( (m % q) == 0 ) {
715 // Q.set_row(m/q, r);
718 // slow and (hopefully) correct
719 for ( size_t i=0; i<n; ++i ) {
721 qk.set(i*q, a.R->one());
725 for ( size_t j=0; j<n; ++j ) {
726 rvec.push_back(r[j]);
732 static void nullspace(Matrix& M, vector<Vec>& basis)
734 const size_t n = M.rowsize();
735 const cl_MI one = M(0,0).ring()->one();
736 for ( size_t i=0; i<n; ++i ) {
737 M(i,i) = M(i,i) - one;
739 for ( size_t r=0; r<n; ++r ) {
741 for ( ; cc<n; ++cc ) {
742 if ( !zerop(M(r,cc)) ) {
744 if ( !zerop(M(cc,cc)) ) {
756 M.mul_col(r, recip(M(r,r)));
757 for ( cc=0; cc<n; ++cc ) {
759 M.sub_col(cc, r, M(r,cc));
765 for ( size_t i=0; i<n; ++i ) {
766 M(i,i) = M(i,i) - one;
768 for ( size_t i=0; i<n; ++i ) {
769 if ( !M.is_row_zero(i) ) {
770 Vec nu(M.row_begin(i), M.row_end(i));
776 static void berlekamp(const UniPoly& a, UniPolyVec& upv)
778 Matrix Q(a.degree(), a.degree(), a.R->zero());
782 const unsigned int k = nu.size();
787 list<UniPoly> factors;
788 factors.push_back(a);
789 unsigned int size = 1;
791 unsigned int q = cl_I_to_uint(a.R->modulus);
793 list<UniPoly>::iterator u = factors.begin();
796 for ( unsigned int s=0; s<q; ++s ) {
798 UniPoly nur(a.R, nu[r]);
799 nur.set(0, nur[0] - cl_MI(a.R, s));
801 if ( !is_one(g) && g != *u ) {
805 throw logic_error("berlekamp: unexpected divisor.");
810 factors.push_back(g);
812 list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
814 if ( i->degree() ) ++size;
818 list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
824 // if ( u->degree() < nur.degree() ) {
836 static void factor_modular(const UniPoly& p, UniPolyVec& upv)
842 static void exteuclid(const UniPoly& a, const UniPoly& b, UniPoly& g, UniPoly& s, UniPoly& t)
844 if ( a.degree() < b.degree() ) {
845 exteuclid(b, a, g, t, s);
848 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);
849 UniPoly c = a; c.unit_normal();
850 UniPoly d = b; d.unit_normal();
851 c1.set(0, a.R->one());
852 d2.set(0, a.R->one());
853 while ( !d.zero() ) {
866 g = c; g.unit_normal();
875 static ex replace_lc(const ex& poly, const ex& x, const ex& lc)
877 ex r = expand(poly + (lc - poly.lcoeff(x)) * pow(x, poly.degree(x)));
881 static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const UniPoly& u1_, const UniPoly& w1_, const ex& gamma_ = 0)
884 const cl_modint_ring& R = u1_.R;
888 for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
889 maxcoeff += pow(abs(a.coeff(x, i)),2);
891 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
892 unsigned int maxdegree = (u1_.degree() > w1_.degree()) ? u1_.degree() : w1_.degree();
893 unsigned int B = cl_I_to_uint(normmc * expt_pos(cl_I(2), maxdegree));
896 ex alpha = a.lcoeff(x);
901 unsigned int gamma_ui = ex_to<numeric>(abs(gamma)).to_int();
908 phi = expand(gamma * nu1.to_ex(x));
909 UniPoly u1(R, phi, x);
910 phi = expand(alpha * nw1.to_ex(x));
911 UniPoly w1(R, phi, x);
914 UniPoly s(R), t(R), g(R);
915 exteuclid(u1, w1, g, s, t);
918 ex u = replace_lc(u1.to_ex(x), x, gamma);
919 ex w = replace_lc(w1.to_ex(x), x, alpha);
920 ex e = expand(a - u * w);
921 unsigned int modulus = p;
924 while ( !e.is_zero() && modulus < 2*B*gamma_ui ) {
926 phi = expand(s.to_ex(x)*c);
927 UniPoly sigmatilde(R, phi, x);
928 phi = expand(t.to_ex(x)*c);
929 UniPoly tautilde(R, phi, x);
931 div(sigmatilde, w1, q);
932 rem(sigmatilde, w1, r);
934 phi = expand(tautilde.to_ex(x) + q.to_ex(x) * u1.to_ex(x));
935 UniPoly tau(R, phi, x);
936 u = expand(u + tau.to_ex(x) * modulus);
937 w = expand(w + sigma.to_ex(x) * modulus);
938 e = expand(a - u * w);
939 modulus = modulus * p;
944 ex delta = u.content(x);
946 w = w / gamma * delta;
954 static unsigned int next_prime(unsigned int p)
956 static vector<unsigned int> primes;
957 if ( primes.size() == 0 ) {
958 primes.push_back(3); primes.push_back(5); primes.push_back(7);
960 vector<unsigned int>::const_iterator it = primes.begin();
961 if ( p >= primes.back() ) {
962 unsigned int candidate = primes.back() + 2;
964 size_t n = primes.size()/2;
965 for ( size_t i=0; i<n; ++i ) {
966 if ( candidate % primes[i] ) continue;
970 primes.push_back(candidate);
971 if ( candidate > p ) break;
975 vector<unsigned int>::const_iterator end = primes.end();
976 for ( ; it!=end; ++it ) {
981 throw logic_error("next_prime: should not reach this point!");
987 Partition(size_t n_) : n(n_)
993 int operator[](size_t i) const { return k[i]; }
994 size_t size() const { return n; }
995 size_t size_first() const { return n-sum; }
996 size_t size_second() const { return sum; }
999 for ( size_t i=n-1; i>=1; --i ) {
1015 static void split(const UniPolyVec& factors, const Partition& part, UniPoly& a, UniPoly& b)
1017 a.set(0, a.R->one());
1018 b.set(0, a.R->one());
1019 for ( size_t i=0; i<part.size(); ++i ) {
1035 static ex factor_univariate(const ex& poly, const ex& x)
1037 ex unit, cont, prim;
1038 poly.unitcontprim(x, unit, cont, prim);
1040 // determine proper prime
1042 cl_modint_ring R = find_modint_ring(p);
1044 if ( irem(ex_to<numeric>(prim.lcoeff(x)), p) != 0 ) {
1045 UniPoly modpoly(R, prim, x);
1046 UniFactorVec sqrfree_ufv;
1047 squarefree(modpoly, sqrfree_ufv);
1048 if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break;
1051 R = find_modint_ring(p);
1054 // do modular factorization
1055 UniPoly modpoly(R, prim, x);
1057 factor_modular(modpoly, factors);
1058 if ( factors.size() <= 1 ) {
1059 // irreducible for sure
1063 // lift all factor combinations
1064 stack<ModFactors> tocheck;
1067 mf.factors = factors;
1070 while ( tocheck.size() ) {
1071 const size_t n = tocheck.top().factors.size();
1075 split(tocheck.top().factors, part, a, b);
1077 ex answer = hensel_univar(tocheck.top().poly, x, p, a, b);
1078 if ( answer != lst() ) {
1079 if ( part.size_first() == 1 ) {
1080 if ( part.size_second() == 1 ) {
1081 result *= answer.op(0) * answer.op(1);
1085 result *= answer.op(0);
1086 tocheck.top().poly = answer.op(1);
1087 for ( size_t i=0; i<n; ++i ) {
1088 if ( part[i] == 0 ) {
1089 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1095 else if ( part.size_second() == 1 ) {
1096 if ( part.size_first() == 1 ) {
1097 result *= answer.op(0) * answer.op(1);
1101 result *= answer.op(1);
1102 tocheck.top().poly = answer.op(0);
1103 for ( size_t i=0; i<n; ++i ) {
1104 if ( part[i] == 1 ) {
1105 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1112 UniPolyVec newfactors1(part.size_first(), R), newfactors2(part.size_second(), R);
1113 UniPolyVec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1114 for ( size_t i=0; i<n; ++i ) {
1116 *i2++ = tocheck.top().factors[i];
1119 *i1++ = tocheck.top().factors[i];
1122 tocheck.top().factors = newfactors1;
1123 tocheck.top().poly = answer.op(0);
1125 mf.factors = newfactors2;
1126 mf.poly = answer.op(1);
1132 if ( !part.next() ) {
1133 result *= tocheck.top().poly;
1141 return unit * cont * result;
1144 struct FindSymbolsMap : public map_function {
1146 ex operator()(const ex& e)
1148 if ( is_a<symbol>(e) ) {
1152 return e.map(*this);
1156 static ex factor_sqrfree(const ex& poly)
1158 // determine all symbols in poly
1159 FindSymbolsMap findsymbols;
1161 if ( findsymbols.syms.size() == 0 ) {
1165 if ( findsymbols.syms.size() == 1 ) {
1166 const ex& x = *(findsymbols.syms.begin());
1167 if ( poly.ldegree(x) > 0 ) {
1168 int ld = poly.ldegree(x);
1169 ex res = factor_univariate(expand(poly/pow(x, ld)), x);
1170 return res * pow(x,ld);
1173 ex res = factor_univariate(poly, x);
1178 // multivariate case not yet implemented!
1179 throw runtime_error("multivariate case not yet implemented!");
1182 } // anonymous namespace
1184 ex factor(const ex& poly)
1186 // determine all symbols in poly
1187 FindSymbolsMap findsymbols;
1189 if ( findsymbols.syms.size() == 0 ) {
1193 exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
1194 for ( ; i!=end; ++i ) {
1198 // make poly square free
1199 ex sfpoly = sqrfree(poly, syms);
1201 // factorize the square free components
1202 if ( is_a<power>(sfpoly) ) {
1203 // case: (polynomial)^exponent
1204 const ex& base = sfpoly.op(0);
1205 if ( !is_a<add>(base) ) {
1206 // simple case: (monomial)^exponent
1209 ex f = factor_sqrfree(base);
1210 return pow(f, sfpoly.op(1));
1212 if ( is_a<mul>(sfpoly) ) {
1214 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
1215 const ex& t = sfpoly.op(i);
1216 if ( is_a<power>(t) ) {
1217 const ex& base = t.op(0);
1218 if ( !is_a<add>(base) ) {
1222 ex f = factor_sqrfree(base);
1223 res *= pow(f, t.op(1));
1226 else if ( is_a<add>(t) ) {
1227 ex f = factor_sqrfree(t);
1236 if ( is_a<symbol>(sfpoly) ) {
1239 // case: (polynomial)
1240 ex f = factor_sqrfree(sfpoly);
1244 } // namespace GiNaC