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);
577 UniFactor uf(z, i++);
578 fvec.factors.push_back(uf);
586 unsigned int prime = cl_I_to_uint(c.R->modulus);
587 c.reduce_exponents(prime);
588 unsigned int pos = fvec.factors.size();
589 sqrfree_main(c, fvec);
590 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
591 fvec.factors[p].exp *= prime;
597 unsigned int prime = cl_I_to_uint(a.R->modulus);
599 amod.reduce_exponents(prime);
600 unsigned int pos = fvec.factors.size();
601 sqrfree_main(amod, fvec);
602 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
603 fvec.factors[p].exp *= prime;
609 static void squarefree(const UniPoly& a, UniFactorVec& fvec)
611 sqrfree_main(a, fvec);
617 friend ostream& operator<<(ostream& o, const Matrix& m);
619 Matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
623 size_t rowsize() const { return r; }
624 size_t colsize() const { return c; }
625 cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
626 cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
627 void mul_col(size_t col, const cl_MI x)
629 Vec::iterator i = m.begin() + col;
630 for ( size_t rc=0; rc<r; ++rc ) {
635 void sub_col(size_t col1, size_t col2, const cl_MI fac)
637 Vec::iterator i1 = m.begin() + col1;
638 Vec::iterator i2 = m.begin() + col2;
639 for ( size_t rc=0; rc<r; ++rc ) {
640 *i1 = *i1 - *i2 * fac;
645 void switch_col(size_t col1, size_t col2)
648 Vec::iterator i1 = m.begin() + col1;
649 Vec::iterator i2 = m.begin() + col2;
650 for ( size_t rc=0; rc<r; ++rc ) {
651 buf = *i1; *i1 = *i2; *i2 = buf;
656 bool is_row_zero(size_t row) const
658 Vec::const_iterator i = m.begin() + row*c;
659 for ( size_t cc=0; cc<c; ++cc ) {
667 void set_row(size_t row, const vector<cl_MI>& newrow)
669 Vec::iterator i1 = m.begin() + row*c;
670 Vec::const_iterator i2 = newrow.begin(), end = newrow.end();
671 for ( ; i2 != end; ++i1, ++i2 ) {
675 Vec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
676 Vec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
683 ostream& operator<<(ostream& o, const Matrix& m)
685 vector<cl_MI>::const_iterator i = m.m.begin(), end = m.m.end();
687 for ( ; i != end; ++i ) {
689 if ( !(wrap++ % m.c) ) {
696 #endif // def DEBUGFACTOR
698 static void q_matrix(const UniPoly& a, Matrix& Q)
700 unsigned int n = a.degree();
701 unsigned int q = cl_I_to_uint(a.R->modulus);
702 vector<cl_MI> r(n, a.R->zero());
705 unsigned int max = (n-1) * q;
706 for ( size_t m=1; m<=max; ++m ) {
707 cl_MI rn_1 = r.back();
708 for ( size_t i=n-1; i>0; --i ) {
709 r[i] = r[i-1] - rn_1 * a[i];
712 if ( (m % q) == 0 ) {
718 static void nullspace(Matrix& M, vector<Vec>& basis)
720 const size_t n = M.rowsize();
721 const cl_MI one = M(0,0).ring()->one();
722 for ( size_t i=0; i<n; ++i ) {
723 M(i,i) = M(i,i) - one;
725 for ( size_t r=0; r<n; ++r ) {
727 for ( ; cc<n; ++cc ) {
728 if ( !zerop(M(r,cc)) ) {
730 if ( !zerop(M(cc,cc)) ) {
742 M.mul_col(r, recip(M(r,r)));
743 for ( cc=0; cc<n; ++cc ) {
745 M.sub_col(cc, r, M(r,cc));
751 for ( size_t i=0; i<n; ++i ) {
752 M(i,i) = M(i,i) - one;
754 for ( size_t i=0; i<n; ++i ) {
755 if ( !M.is_row_zero(i) ) {
756 Vec nu(M.row_begin(i), M.row_end(i));
762 static void berlekamp(const UniPoly& a, UniPolyVec& upv)
764 Matrix Q(a.degree(), a.degree(), a.R->zero());
768 const unsigned int k = nu.size();
773 list<UniPoly> factors;
774 factors.push_back(a);
775 unsigned int size = 1;
777 unsigned int q = cl_I_to_uint(a.R->modulus);
779 list<UniPoly>::iterator u = factors.begin();
782 for ( unsigned int s=0; s<q; ++s ) {
784 UniPoly nur(a.R, nu[r]);
785 nur.set(0, nur[0] - cl_MI(a.R, s));
787 if ( !is_one(g) && g != *u ) {
791 throw logic_error("berlekamp: unexpected divisor.");
796 factors.push_back(g);
799 list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
805 if ( u->degree() < nur.degree() ) {
817 static void factor_modular(const UniPoly& p, UniPolyVec& upv)
823 static void exteuclid(const UniPoly& a, const UniPoly& b, UniPoly& g, UniPoly& s, UniPoly& t)
825 if ( a.degree() < b.degree() ) {
826 exteuclid(b, a, g, t, s);
829 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);
830 UniPoly c = a; c.unit_normal();
831 UniPoly d = b; d.unit_normal();
832 c1.set(0, a.R->one());
833 d2.set(0, a.R->one());
834 while ( !d.zero() ) {
847 g = c; g.unit_normal();
856 static ex replace_lc(const ex& poly, const ex& x, const ex& lc)
858 ex r = expand(poly + (lc - poly.lcoeff(x)) * pow(x, poly.degree(x)));
862 static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const UniPoly& u1_, const UniPoly& w1_, const ex& gamma_ = 0)
865 const cl_modint_ring& R = u1_.R;
869 for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
870 maxcoeff += pow(abs(a.coeff(x, i)),2);
872 cl_I normmc = ceiling1(the<cl_F>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
873 unsigned int maxdegree = (u1_.degree() > w1_.degree()) ? u1_.degree() : w1_.degree();
874 unsigned int B = cl_I_to_uint(normmc * expt_pos(cl_I(2), maxdegree));
877 ex alpha = a.lcoeff(x);
882 unsigned int gamma_ui = ex_to<numeric>(abs(gamma)).to_int();
889 phi = expand(gamma * nu1.to_ex(x));
890 UniPoly u1(R, phi, x);
891 phi = expand(alpha * nw1.to_ex(x));
892 UniPoly w1(R, phi, x);
895 UniPoly s(R), t(R), g(R);
896 exteuclid(u1, w1, g, s, t);
899 ex u = replace_lc(u1.to_ex(x), x, gamma);
900 ex w = replace_lc(w1.to_ex(x), x, alpha);
901 ex e = expand(a - u * w);
902 unsigned int modulus = p;
905 while ( !e.is_zero() && modulus < 2*B*gamma_ui ) {
907 phi = expand(s.to_ex(x)*c);
908 UniPoly sigmatilde(R, phi, x);
909 phi = expand(t.to_ex(x)*c);
910 UniPoly tautilde(R, phi, x);
912 div(sigmatilde, w1, q);
913 rem(sigmatilde, w1, r);
915 phi = expand(tautilde.to_ex(x) + q.to_ex(x) * u1.to_ex(x));
916 UniPoly tau(R, phi, x);
917 u = expand(u + tau.to_ex(x) * modulus);
918 w = expand(w + sigma.to_ex(x) * modulus);
919 e = expand(a - u * w);
920 modulus = modulus * p;
925 ex delta = u.content(x);
927 w = w / gamma * delta;
935 static unsigned int next_prime(unsigned int p)
937 static vector<unsigned int> primes;
938 if ( primes.size() == 0 ) {
939 primes.push_back(3); primes.push_back(5); primes.push_back(7);
941 vector<unsigned int>::const_iterator it = primes.begin();
942 if ( p >= primes.back() ) {
943 unsigned int candidate = primes.back() + 2;
945 size_t n = primes.size()/2;
946 for ( size_t i=0; i<n; ++i ) {
947 if ( candidate % primes[i] ) continue;
951 primes.push_back(candidate);
952 if ( candidate > p ) break;
956 vector<unsigned int>::const_iterator end = primes.end();
957 for ( ; it!=end; ++it ) {
962 throw logic_error("next_prime: should not reach this point!");
968 Partition(size_t n_) : n(n_)
974 int operator[](size_t i) const { return k[i]; }
975 size_t size() const { return n; }
976 size_t size_first() const { return n-sum; }
977 size_t size_second() const { return sum; }
980 for ( size_t i=n-1; i>=1; --i ) {
996 static void split(const UniPolyVec& factors, const Partition& part, UniPoly& a, UniPoly& b)
998 a.set(0, a.R->one());
999 b.set(0, a.R->one());
1000 for ( size_t i=0; i<part.size(); ++i ) {
1016 static ex factor_univariate(const ex& poly, const ex& x)
1018 ex unit, cont, prim;
1019 poly.unitcontprim(x, unit, cont, prim);
1021 // determine proper prime
1023 cl_modint_ring R = find_modint_ring(p);
1025 if ( irem(ex_to<numeric>(prim.lcoeff(x)), p) != 0 ) {
1026 UniPoly modpoly(R, prim, x);
1027 UniFactorVec sqrfree_ufv;
1028 squarefree(modpoly, sqrfree_ufv);
1029 if ( sqrfree_ufv.factors.size() == 1 ) break;
1032 R = find_modint_ring(p);
1035 // do modular factorization
1036 UniPoly modpoly(R, prim, x);
1038 factor_modular(modpoly, factors);
1039 if ( factors.size() <= 1 ) {
1040 // irreducible for sure
1044 // lift all factor combinations
1045 stack<ModFactors> tocheck;
1048 mf.factors = factors;
1051 while ( tocheck.size() ) {
1052 const size_t n = tocheck.top().factors.size();
1056 split(tocheck.top().factors, part, a, b);
1058 ex answer = hensel_univar(tocheck.top().poly, x, p, a, b);
1059 if ( answer != lst() ) {
1060 if ( part.size_first() == 1 ) {
1061 if ( part.size_second() == 1 ) {
1062 result *= answer.op(0) * answer.op(1);
1066 result *= answer.op(0);
1067 tocheck.top().poly = answer.op(1);
1068 for ( size_t i=0; i<n; ++i ) {
1069 if ( part[i] == 0 ) {
1070 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1076 else if ( part.size_second() == 1 ) {
1077 if ( part.size_first() == 1 ) {
1078 result *= answer.op(0) * answer.op(1);
1082 result *= answer.op(1);
1083 tocheck.top().poly = answer.op(0);
1084 for ( size_t i=0; i<n; ++i ) {
1085 if ( part[i] == 1 ) {
1086 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1093 UniPolyVec newfactors1(part.size_first(), R), newfactors2(part.size_second(), R);
1094 UniPolyVec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1095 for ( size_t i=0; i<n; ++i ) {
1097 *i2++ = tocheck.top().factors[i];
1100 *i1++ = tocheck.top().factors[i];
1103 tocheck.top().factors = newfactors1;
1104 tocheck.top().poly = answer.op(0);
1106 mf.factors = newfactors2;
1107 mf.poly = answer.op(1);
1112 if ( !part.next() ) {
1113 result *= tocheck.top().poly;
1121 return unit * cont * result;
1124 struct FindSymbolsMap : public map_function {
1126 ex operator()(const ex& e)
1128 if ( is_a<symbol>(e) ) {
1132 return e.map(*this);
1136 static ex factor_sqrfree(const ex& poly)
1138 // determine all symbols in poly
1139 FindSymbolsMap findsymbols;
1141 if ( findsymbols.syms.size() == 0 ) {
1145 if ( findsymbols.syms.size() == 1 ) {
1146 const ex& x = *(findsymbols.syms.begin());
1147 if ( poly.ldegree(x) > 0 ) {
1148 int ld = poly.ldegree(x);
1149 ex res = factor_univariate(expand(poly/pow(x, ld)), x);
1150 return res * pow(x,ld);
1153 ex res = factor_univariate(poly, x);
1158 // multivariate case not yet implemented!
1159 throw runtime_error("multivariate case not yet implemented!");
1162 } // anonymous namespace
1164 ex factor(const ex& poly)
1166 // determine all symbols in poly
1167 FindSymbolsMap findsymbols;
1169 if ( findsymbols.syms.size() == 0 ) {
1173 exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
1174 for ( ; i!=end; ++i ) {
1178 // make poly square free
1179 ex sfpoly = sqrfree(poly, syms);
1181 // factorize the square free components
1182 if ( is_a<power>(sfpoly) ) {
1183 // case: (polynomial)^exponent
1184 const ex& base = sfpoly.op(0);
1185 if ( !is_a<add>(base) ) {
1186 // simple case: (monomial)^exponent
1189 ex f = factor_sqrfree(base);
1190 return pow(f, sfpoly.op(1));
1192 if ( is_a<mul>(sfpoly) ) {
1194 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
1195 const ex& t = sfpoly.op(i);
1196 if ( is_a<power>(t) ) {
1197 const ex& base = t.op(0);
1198 if ( !is_a<add>(base) ) {
1202 ex f = factor_sqrfree(base);
1203 res *= pow(f, t.op(1));
1206 else if ( is_a<add>(t) ) {
1207 ex f = factor_sqrfree(t);
1216 // case: (polynomial)
1217 ex f = factor_sqrfree(sfpoly);
1221 } // namespace GiNaC