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 <ginac/ginac.h>
30 using namespace GiNaC;
36 #include "operators.h"
39 #include "relational.h"
61 #define DCOUT(str) cout << #str << endl
62 #define DCOUTVAR(var) cout << #var << ": " << var << endl
63 #define DCOUT2(str,var) cout << #str << ": " << var << endl
67 #define DCOUT2(str,var)
72 typedef vector<cl_MI> Vec;
73 typedef vector<Vec> VecVec;
76 ostream& operator<<(ostream& o, const Vec& v)
78 Vec::const_iterator i = v.begin(), end = v.end();
84 #endif // def DEBUGFACTOR
87 ostream& operator<<(ostream& o, const VecVec& v)
89 VecVec::const_iterator i = v.begin(), end = v.end();
95 #endif // def DEBUGFACTOR
99 cl_MI c; // coefficient
100 unsigned int exp; // exponent >=0
104 ostream& operator<<(ostream& o, const Term& t)
107 o << "(" << t.c << ")x^" << t.exp;
110 o << "(" << t.c << ")";
114 #endif // def DEBUGFACTOR
119 list<Term> terms; // highest exponent first
121 UniPoly(const cl_modint_ring& ring) : R(ring) { }
122 UniPoly(const cl_modint_ring& ring, const ex& poly, const ex& x) : R(ring)
124 // assert: poly is in Z[x]
126 for ( int i=poly.degree(x); i>=poly.ldegree(x); --i ) {
127 int coeff = ex_to<numeric>(poly.coeff(x,i)).to_int();
129 t.c = R->canonhom(coeff);
137 UniPoly(const cl_modint_ring& ring, const UniPoly& poly) : R(ring)
139 if ( R->modulus == poly.R->modulus ) {
143 list<Term>::const_iterator i=poly.terms.begin(), end=poly.terms.end();
144 for ( ; i!=end; ++i ) {
146 terms.back().c = R->canonhom(poly.R->retract(i->c));
147 if ( zerop(terms.back().c) ) {
153 UniPoly(const cl_modint_ring& ring, const Vec& v) : R(ring)
156 for ( unsigned int i=0; i<v.size(); ++i ) {
157 if ( !zerop(v[i]) ) {
164 unsigned int degree() const
166 if ( terms.size() ) {
167 return terms.front().exp;
173 bool zero() const { return (terms.size() == 0); }
174 const cl_MI operator[](unsigned int deg) const
176 list<Term>::const_iterator i = terms.begin(), end = terms.end();
177 for ( ; i != end; ++i ) {
178 if ( i->exp == deg ) {
181 if ( i->exp < deg ) {
187 void set(unsigned int deg, const cl_MI& c)
189 list<Term>::iterator i = terms.begin(), end = terms.end();
191 if ( i->exp == deg ) {
200 if ( i->exp < deg ) {
212 ex to_ex(const ex& x, bool symmetric = true) const
215 list<Term>::const_iterator i = terms.begin(), end = terms.end();
217 numeric mod(R->modulus);
218 numeric halfmod = (mod-1)/2;
219 for ( ; i != end; ++i ) {
220 numeric n(R->retract(i->c));
222 r += pow(x, i->exp) * (n-mod);
225 r += pow(x, i->exp) * n;
230 for ( ; i != end; ++i ) {
231 r += pow(x, i->exp) * numeric(R->retract(i->c));
238 if ( terms.size() ) {
239 if ( terms.front().c != R->one() ) {
240 list<Term>::iterator i = terms.begin(), end = terms.end();
243 while ( ++i != end ) {
244 i->c = div(i->c, cont);
254 return terms.front().c;
256 void divide(const cl_MI& x)
258 list<Term>::iterator i = terms.begin(), end = terms.end();
259 for ( ; i != end; ++i ) {
266 void divide(const cl_I& x)
268 list<Term>::iterator i = terms.begin(), end = terms.end();
269 for ( ; i != end; ++i ) {
270 i->c = cl_MI(R, the<cl_I>(R->retract(i->c) / x));
273 void reduce_exponents(unsigned int prime)
275 list<Term>::iterator i = terms.begin(), end = terms.end();
278 // assert: i->exp is multiple of prime
284 void deriv(UniPoly& d) const
286 list<Term>::const_iterator i = terms.begin(), end = terms.end();
289 cl_MI newc = i->c * i->exp;
290 if ( !zerop(newc) ) {
294 d.terms.push_back(t);
300 bool operator<(const UniPoly& o) const
302 if ( terms.size() != o.terms.size() ) {
303 return terms.size() < o.terms.size();
305 list<Term>::const_iterator i1 = terms.begin(), end = terms.end();
306 list<Term>::const_iterator i2 = o.terms.begin();
307 while ( i1 != end ) {
308 if ( i1->exp != i2->exp ) {
309 return i1->exp < i2->exp;
311 if ( i1->c != i2->c ) {
312 return R->retract(i1->c) < R->retract(i2->c);
318 bool operator==(const UniPoly& o) const
320 if ( terms.size() != o.terms.size() ) {
323 list<Term>::const_iterator i1 = terms.begin(), end = terms.end();
324 list<Term>::const_iterator i2 = o.terms.begin();
325 while ( i1 != end ) {
326 if ( i1->exp != i2->exp ) {
329 if ( i1->c != i2->c ) {
336 bool operator!=(const UniPoly& o) const
338 bool res = !(*this == o);
343 static UniPoly operator*(const UniPoly& a, const UniPoly& b)
345 unsigned int n = a.degree()+b.degree();
348 for ( unsigned int i=0 ; i<=n; ++i ) {
350 for ( unsigned int j=0 ; j<=i; ++j ) {
351 t.c = t.c + a[j] * b[i-j];
355 c.terms.push_front(t);
361 static UniPoly operator-(const UniPoly& a, const UniPoly& b)
363 list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
364 list<Term>::const_iterator ib = b.terms.begin(), bend = b.terms.end();
366 while ( ia != aend && ib != bend ) {
367 if ( ia->exp > ib->exp ) {
368 c.terms.push_back(*ia);
371 else if ( ia->exp < ib->exp ) {
372 c.terms.push_back(*ib);
373 c.terms.back().c = -c.terms.back().c;
381 c.terms.push_back(t);
386 while ( ia != aend ) {
387 c.terms.push_back(*ia);
390 while ( ib != bend ) {
391 c.terms.push_back(*ib);
392 c.terms.back().c = -c.terms.back().c;
398 static UniPoly operator*(const UniPoly& a, const cl_MI& fac)
400 unsigned int n = a.degree();
403 for ( unsigned int i=0 ; i<=n; ++i ) {
407 c.terms.push_front(t);
413 static UniPoly operator+(const UniPoly& a, const UniPoly& b)
415 list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
416 list<Term>::const_iterator ib = b.terms.begin(), bend = b.terms.end();
418 while ( ia != aend && ib != bend ) {
419 if ( ia->exp > ib->exp ) {
420 c.terms.push_back(*ia);
423 else if ( ia->exp < ib->exp ) {
424 c.terms.push_back(*ib);
432 c.terms.push_back(t);
437 while ( ia != aend ) {
438 c.terms.push_back(*ia);
441 while ( ib != bend ) {
442 c.terms.push_back(*ib);
448 // static UniPoly operator-(const UniPoly& a)
450 // list<Term>::const_iterator ia = a.terms.begin(), aend = a.terms.end();
452 // while ( ia != aend ) {
453 // c.terms.push_back(*ia);
454 // c.terms.back().c = -c.terms.back().c;
461 ostream& operator<<(ostream& o, const UniPoly& t)
463 list<Term>::const_iterator i = t.terms.begin(), end = t.terms.end();
468 for ( ; i != end; ) {
476 #endif // def DEBUGFACTOR
479 ostream& operator<<(ostream& o, const list<UniPoly>& t)
481 list<UniPoly>::const_iterator i = t.begin(), end = t.end();
483 for ( ; i != end; ) {
489 #endif // def DEBUGFACTOR
491 typedef vector<UniPoly> UniPolyVec;
494 ostream& operator<<(ostream& o, const UniPolyVec& v)
496 UniPolyVec::const_iterator i = v.begin(), end = v.end();
498 o << *i++ << " , " << endl;
502 #endif // def DEBUGFACTOR
509 UniFactor(const cl_modint_ring& ring) : p(ring) { }
510 UniFactor(const UniPoly& p_, unsigned int exp_) : p(p_), exp(exp_) { }
511 bool operator<(const UniFactor& o) const
519 vector<UniFactor> factors;
523 sort(factors.begin(), factors.end());
524 if ( factors.size() > 1 ) {
525 vector<UniFactor>::iterator i = factors.begin();
526 vector<UniFactor>::const_iterator cmp = factors.begin()+1;
527 vector<UniFactor>::iterator end = factors.end();
528 while ( cmp != end ) {
529 if ( i->p != cmp->p ) {
539 factors.erase(i+1, end);
546 ostream& operator<<(ostream& o, const UniFactorVec& ufv)
548 for ( size_t i=0; i<ufv.factors.size(); ++i ) {
549 if ( i != ufv.factors.size()-1 ) {
555 o << "[ " << ufv.factors[i].p << " ]^" << ufv.factors[i].exp << endl;
559 #endif // def DEBUGFACTOR
561 static void rem(const UniPoly& a_, const UniPoly& b, UniPoly& c)
563 if ( a_.degree() < b.degree() ) {
581 cl_MI qk = div(c[n+k], b[n]);
584 for ( unsigned int i=0; i<n; ++i ) {
586 c.set(j, c[j] - qk*b[j-k]);
592 list<Term>::iterator i = c.terms.begin(), end = c.terms.end();
594 if ( i->exp <= n-1 ) {
599 c.terms.erase(c.terms.begin(), i);
602 static void div(const UniPoly& a_, const UniPoly& b, UniPoly& q)
604 if ( a_.degree() < b.degree() ) {
617 cl_MI qk = div(c[n+k], b[n]);
622 q.terms.push_back(t);
624 for ( unsigned int i=0; i<n; ++i ) {
626 c.set(j, c[j] - qk*b[j-k]);
634 static void gcd(const UniPoly& a, const UniPoly& b, UniPoly& c)
641 if ( c.degree() < d.degree() ) {
646 while ( !d.zero() ) {
655 static bool is_one(const UniPoly& w)
657 if ( w.terms.size() == 1 && w[0] == w.R->one() ) {
663 static void sqrfree_main(const UniPoly& a, UniFactorVec& fvec)
669 UniPoly c(a.R), w(a.R);
672 while ( !is_one(w) ) {
673 UniPoly y(a.R), z(a.R);
678 fvec.factors.push_back(uf);
687 unsigned int prime = cl_I_to_uint(c.R->modulus);
688 c.reduce_exponents(prime);
689 unsigned int pos = fvec.factors.size();
690 sqrfree_main(c, fvec);
691 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
692 fvec.factors[p].exp *= prime;
698 unsigned int prime = cl_I_to_uint(a.R->modulus);
700 amod.reduce_exponents(prime);
701 unsigned int pos = fvec.factors.size();
702 sqrfree_main(amod, fvec);
703 for ( unsigned int p=pos; p<fvec.factors.size(); ++p ) {
704 fvec.factors[p].exp *= prime;
710 static void squarefree(const UniPoly& a, UniFactorVec& fvec)
712 sqrfree_main(a, fvec);
718 friend ostream& operator<<(ostream& o, const Matrix& m);
720 Matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
724 size_t rowsize() const { return r; }
725 size_t colsize() const { return c; }
726 cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
727 cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
728 void mul_col(size_t col, const cl_MI x)
730 Vec::iterator i = m.begin() + col;
731 for ( size_t rc=0; rc<r; ++rc ) {
736 void sub_col(size_t col1, size_t col2, const cl_MI fac)
738 Vec::iterator i1 = m.begin() + col1;
739 Vec::iterator i2 = m.begin() + col2;
740 for ( size_t rc=0; rc<r; ++rc ) {
741 *i1 = *i1 - *i2 * fac;
746 void switch_col(size_t col1, size_t col2)
749 Vec::iterator i1 = m.begin() + col1;
750 Vec::iterator i2 = m.begin() + col2;
751 for ( size_t rc=0; rc<r; ++rc ) {
752 buf = *i1; *i1 = *i2; *i2 = buf;
757 void mul_row(size_t row, const cl_MI x)
759 vector<cl_MI>::iterator i = m.begin() + row*c;
760 for ( size_t cc=0; cc<c; ++cc ) {
765 void sub_row(size_t row1, size_t row2, const cl_MI fac)
767 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
768 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
769 for ( size_t cc=0; cc<c; ++cc ) {
770 *i1 = *i1 - *i2 * fac;
775 void switch_row(size_t row1, size_t row2)
778 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
779 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
780 for ( size_t cc=0; cc<c; ++cc ) {
781 buf = *i1; *i1 = *i2; *i2 = buf;
786 bool is_col_zero(size_t col) const
788 Vec::const_iterator i = m.begin() + col;
789 for ( size_t rr=0; rr<r; ++rr ) {
797 bool is_row_zero(size_t row) const
799 Vec::const_iterator i = m.begin() + row*c;
800 for ( size_t cc=0; cc<c; ++cc ) {
808 void set_row(size_t row, const vector<cl_MI>& newrow)
810 Vec::iterator i1 = m.begin() + row*c;
811 Vec::const_iterator i2 = newrow.begin(), end = newrow.end();
812 for ( ; i2 != end; ++i1, ++i2 ) {
816 Vec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
817 Vec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
824 Matrix operator*(const Matrix& m1, const Matrix& m2)
826 const unsigned int r = m1.rowsize();
827 const unsigned int c = m2.colsize();
828 Matrix o(r,c,m1(0,0));
830 for ( size_t i=0; i<r; ++i ) {
831 for ( size_t j=0; j<c; ++j ) {
833 buf = m1(i,0) * m2(0,j);
834 for ( size_t k=1; k<c; ++k ) {
835 buf = buf + m1(i,k)*m2(k,j);
843 ostream& operator<<(ostream& o, const Matrix& m)
845 vector<cl_MI>::const_iterator i = m.m.begin(), end = m.m.end();
847 for ( ; i != end; ++i ) {
849 if ( !(wrap++ % m.c) ) {
856 #endif // def DEBUGFACTOR
858 static void q_matrix(const UniPoly& a, Matrix& Q)
860 unsigned int n = a.degree();
861 unsigned int q = cl_I_to_uint(a.R->modulus);
863 // vector<cl_MI> r(n, a.R->zero());
864 // r[0] = a.R->one();
866 // unsigned int max = (n-1) * q;
867 // for ( size_t m=1; m<=max; ++m ) {
868 // cl_MI rn_1 = r.back();
869 // for ( size_t i=n-1; i>0; --i ) {
870 // r[i] = r[i-1] - rn_1 * a[i];
872 // r[0] = -rn_1 * a[0];
873 // if ( (m % q) == 0 ) {
874 // Q.set_row(m/q, r);
877 // slow and (hopefully) correct
878 for ( size_t i=0; i<n; ++i ) {
880 qk.set(i*q, a.R->one());
884 for ( size_t j=0; j<n; ++j ) {
885 rvec.push_back(r[j]);
891 static void nullspace(Matrix& M, vector<Vec>& basis)
893 const size_t n = M.rowsize();
894 const cl_MI one = M(0,0).ring()->one();
895 for ( size_t i=0; i<n; ++i ) {
896 M(i,i) = M(i,i) - one;
898 for ( size_t r=0; r<n; ++r ) {
900 for ( ; cc<n; ++cc ) {
901 if ( !zerop(M(r,cc)) ) {
903 if ( !zerop(M(cc,cc)) ) {
915 M.mul_col(r, recip(M(r,r)));
916 for ( cc=0; cc<n; ++cc ) {
918 M.sub_col(cc, r, M(r,cc));
924 for ( size_t i=0; i<n; ++i ) {
925 M(i,i) = M(i,i) - one;
927 for ( size_t i=0; i<n; ++i ) {
928 if ( !M.is_row_zero(i) ) {
929 Vec nu(M.row_begin(i), M.row_end(i));
935 static void berlekamp(const UniPoly& a, UniPolyVec& upv)
937 Matrix Q(a.degree(), a.degree(), a.R->zero());
941 const unsigned int k = nu.size();
946 list<UniPoly> factors;
947 factors.push_back(a);
948 unsigned int size = 1;
950 unsigned int q = cl_I_to_uint(a.R->modulus);
952 list<UniPoly>::iterator u = factors.begin();
955 for ( unsigned int s=0; s<q; ++s ) {
957 UniPoly nur(a.R, nu[r]);
958 nur.set(0, nur[0] - cl_MI(a.R, s));
960 if ( !is_one(g) && g != *u ) {
964 throw logic_error("berlekamp: unexpected divisor.");
969 factors.push_back(g);
971 list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
973 if ( i->degree() ) ++size;
977 list<UniPoly>::const_iterator i = factors.begin(), end = factors.end();
983 // if ( u->degree() < nur.degree() ) {
995 static void factor_modular(const UniPoly& p, UniPolyVec& upv)
1001 static void exteuclid(const UniPoly& a, const UniPoly& b, UniPoly& g, UniPoly& s, UniPoly& t)
1003 if ( a.degree() < b.degree() ) {
1004 exteuclid(b, a, g, t, s);
1007 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);
1008 UniPoly c = a; c.unit_normal();
1009 UniPoly d = b; d.unit_normal();
1010 c1.set(0, a.R->one());
1011 d2.set(0, a.R->one());
1012 while ( !d.zero() ) {
1025 g = c; g.unit_normal();
1034 static ex replace_lc(const ex& poly, const ex& x, const ex& lc)
1036 ex r = expand(poly + (lc - poly.lcoeff(x)) * pow(x, poly.degree(x)));
1040 static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const UniPoly& u1_, const UniPoly& w1_, const ex& gamma_ = 0)
1043 const cl_modint_ring& R = u1_.R;
1047 for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
1048 maxcoeff += pow(abs(a.coeff(x, i)),2);
1050 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
1051 cl_I maxdegree = (u1_.degree() > w1_.degree()) ? u1_.degree() : w1_.degree();
1052 cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
1055 ex alpha = a.lcoeff(x);
1060 unsigned int gamma_ui = ex_to<numeric>(abs(gamma)).to_int();
1067 phi = expand(gamma * nu1.to_ex(x));
1068 UniPoly u1(R, phi, x);
1069 phi = expand(alpha * nw1.to_ex(x));
1070 UniPoly w1(R, phi, x);
1073 UniPoly s(R), t(R), g(R);
1074 exteuclid(u1, w1, g, s, t);
1077 ex u = replace_lc(u1.to_ex(x), x, gamma);
1078 ex w = replace_lc(w1.to_ex(x), x, alpha);
1079 ex e = expand(a - u * w);
1080 numeric modulus = p;
1081 const numeric maxmodulus(2*B*gamma_ui);
1084 while ( !e.is_zero() && modulus < maxmodulus ) {
1086 phi = expand(s.to_ex(x)*c);
1087 UniPoly sigmatilde(R, phi, x);
1088 phi = expand(t.to_ex(x)*c);
1089 UniPoly tautilde(R, phi, x);
1091 div(sigmatilde, w1, q);
1092 rem(sigmatilde, w1, r);
1094 phi = expand(tautilde.to_ex(x) + q.to_ex(x) * u1.to_ex(x));
1095 UniPoly tau(R, phi, x);
1096 u = expand(u + tau.to_ex(x) * modulus);
1097 w = expand(w + sigma.to_ex(x) * modulus);
1098 e = expand(a - u * w);
1099 modulus = modulus * p;
1103 if ( e.is_zero() ) {
1104 ex delta = u.content(x);
1106 w = w / gamma * delta;
1114 static unsigned int next_prime(unsigned int p)
1116 static vector<unsigned int> primes;
1117 if ( primes.size() == 0 ) {
1118 primes.push_back(3); primes.push_back(5); primes.push_back(7);
1120 vector<unsigned int>::const_iterator it = primes.begin();
1121 if ( p >= primes.back() ) {
1122 unsigned int candidate = primes.back() + 2;
1124 size_t n = primes.size()/2;
1125 for ( size_t i=0; i<n; ++i ) {
1126 if ( candidate % primes[i] ) continue;
1130 primes.push_back(candidate);
1131 if ( candidate > p ) break;
1135 vector<unsigned int>::const_iterator end = primes.end();
1136 for ( ; it!=end; ++it ) {
1141 throw logic_error("next_prime: should not reach this point!");
1147 Partition(size_t n_) : n(n_)
1153 int operator[](size_t i) const { return k[i]; }
1154 size_t size() const { return n; }
1155 size_t size_first() const { return n-sum; }
1156 size_t size_second() const { return sum; }
1160 for ( size_t i=0; i<k.size(); ++i ) {
1161 cout << k[i] << " ";
1168 for ( size_t i=n-1; i>=1; --i ) {
1184 static void split(const UniPolyVec& factors, const Partition& part, UniPoly& a, UniPoly& b)
1186 a.set(0, a.R->one());
1187 b.set(0, a.R->one());
1188 for ( size_t i=0; i<part.size(); ++i ) {
1204 static ex factor_univariate(const ex& poly, const ex& x)
1206 ex unit, cont, prim;
1207 poly.unitcontprim(x, unit, cont, prim);
1209 // determine proper prime
1211 cl_modint_ring R = find_modint_ring(p);
1213 if ( irem(ex_to<numeric>(prim.lcoeff(x)), p) != 0 ) {
1214 UniPoly modpoly(R, prim, x);
1215 UniFactorVec sqrfree_ufv;
1216 squarefree(modpoly, sqrfree_ufv);
1217 if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break;
1220 R = find_modint_ring(p);
1223 // do modular factorization
1224 UniPoly modpoly(R, prim, x);
1226 factor_modular(modpoly, factors);
1227 if ( factors.size() <= 1 ) {
1228 // irreducible for sure
1232 // lift all factor combinations
1233 stack<ModFactors> tocheck;
1236 mf.factors = factors;
1239 while ( tocheck.size() ) {
1240 const size_t n = tocheck.top().factors.size();
1244 split(tocheck.top().factors, part, a, b);
1246 ex answer = hensel_univar(tocheck.top().poly, x, p, a, b);
1247 if ( answer != lst() ) {
1248 if ( part.size_first() == 1 ) {
1249 if ( part.size_second() == 1 ) {
1250 result *= answer.op(0) * answer.op(1);
1254 result *= answer.op(0);
1255 tocheck.top().poly = answer.op(1);
1256 for ( size_t i=0; i<n; ++i ) {
1257 if ( part[i] == 0 ) {
1258 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1264 else if ( part.size_second() == 1 ) {
1265 if ( part.size_first() == 1 ) {
1266 result *= answer.op(0) * answer.op(1);
1270 result *= answer.op(1);
1271 tocheck.top().poly = answer.op(0);
1272 for ( size_t i=0; i<n; ++i ) {
1273 if ( part[i] == 1 ) {
1274 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1281 UniPolyVec newfactors1(part.size_first(), R), newfactors2(part.size_second(), R);
1282 UniPolyVec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1283 for ( size_t i=0; i<n; ++i ) {
1285 *i2++ = tocheck.top().factors[i];
1288 *i1++ = tocheck.top().factors[i];
1291 tocheck.top().factors = newfactors1;
1292 tocheck.top().poly = answer.op(0);
1294 mf.factors = newfactors2;
1295 mf.poly = answer.op(1);
1301 if ( !part.next() ) {
1302 result *= tocheck.top().poly;
1310 return unit * cont * result;
1313 struct FindSymbolsMap : public map_function {
1315 ex operator()(const ex& e)
1317 if ( is_a<symbol>(e) ) {
1321 return e.map(*this);
1331 // forward declaration
1332 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);
1334 UniPolyVec multiterm_eea_lift(const UniPolyVec& a, const ex& x, unsigned int p, unsigned int k)
1336 DCOUT(multiterm_eea_lift);
1341 const size_t r = a.size();
1343 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1345 UniPolyVec q(r-1, fill);
1347 for ( size_t j=r-2; j>=1; --j ) {
1348 q[j-1] = a[j] * q[j];
1352 beta.set(0, R->one());
1354 for ( size_t j=1; j<r; ++j ) {
1357 vector<ex> mdarg(2);
1358 mdarg[0] = q[j-1].to_ex(x);
1359 mdarg[1] = a[j-1].to_ex(x);
1360 vector<EvalPoint> empty;
1361 vector<ex> exsigma = multivar_diophant(mdarg, x, beta.to_ex(x), empty, 0, p, k);
1362 UniPoly sigma1(R, exsigma[0], x);
1363 UniPoly sigma2(R, exsigma[1], x);
1365 s.push_back(sigma2);
1370 DCOUT(END multiterm_eea_lift);
1374 void eea_lift(const UniPoly& a, const UniPoly& b, const ex& x, unsigned int p, unsigned int k, UniPoly& s_, UniPoly& t_)
1383 cl_modint_ring R = find_modint_ring(p);
1389 UniPoly smod(R), tmod(R), g(R);
1390 exteuclid(amod, bmod, g, smod, tmod);
1396 cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1397 UniPoly s(Rpk, smod);
1398 UniPoly t(Rpk, tmod);
1405 one.set(0, Rpk->one());
1406 for ( size_t j=1; j<k; ++j ) {
1407 UniPoly e = one - a * s - b * t;
1410 UniPoly sigmabar(R);
1411 sigmabar = smod * c;
1415 div(sigmabar, bmod, q);
1417 rem(sigmabar, bmod, sigma);
1419 tau = taubar + q * amod;
1420 UniPoly sadd(Rpk, sigma);
1421 cl_MI modmodulus(Rpk, modulus);
1422 s = s + sadd * modmodulus;
1423 UniPoly tadd(Rpk, tau);
1424 t = t + tadd * modmodulus;
1425 modulus = modulus * p;
1432 DCOUT2(check, a*s + b*t);
1433 DCOUT(END eea_lift);
1436 UniPolyVec univar_diophant(const UniPolyVec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1438 DCOUT(univar_diophant);
1445 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1447 const size_t r = a.size();
1450 UniPolyVec s = multiterm_eea_lift(a, x, p, k);
1451 for ( size_t j=0; j<r; ++j ) {
1452 ex phi = expand(pow(x,m)*s[j].to_ex(x));
1453 UniPoly bmod(R, phi, x);
1455 rem(bmod, a[j], buf);
1456 result.push_back(buf);
1461 eea_lift(a[1], a[0], x, p, k, s, t);
1462 ex phi = expand(pow(x,m)*s.to_ex(x));
1463 UniPoly bmod(R, phi, x);
1465 rem(bmod, a[0], buf);
1466 result.push_back(buf);
1469 phi = expand(pow(x,m)*t.to_ex(x));
1470 UniPoly t1mod(R, phi, x);
1471 buf = t1mod + q * a[1];
1472 result.push_back(buf);
1476 DCOUT(END univar_diophant);
1480 struct make_modular_map : public map_function {
1482 make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1483 ex operator()(const ex& e)
1485 if ( is_a<add>(e) || is_a<mul>(e) ) {
1486 return e.map(*this);
1488 else if ( is_a<numeric>(e) ) {
1489 numeric mod(R->modulus);
1490 numeric halfmod = (mod-1)/2;
1491 cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1492 numeric n(R->retract(emod));
1493 if ( n > halfmod ) {
1504 static ex make_modular(const ex& e, const cl_modint_ring& R)
1506 make_modular_map map(R);
1510 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)
1514 DCOUT(multivar_diophant);
1517 for ( size_t i=0; i<a.size(); ++i ) {
1518 cout << a[i] << " ";
1526 for ( size_t i=0; i<I.size(); ++i ) {
1527 cout << I[i].x << "=" << I[i].evalpoint << " ";
1535 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1536 const size_t r = a.size();
1537 const size_t nu = I.size() + 1;
1543 ex xnu = I.back().x;
1544 int alphanu = I.back().evalpoint;
1547 for ( size_t i=0; i<r; ++i ) {
1551 for ( size_t i=0; i<r; ++i ) {
1552 b[i] = normal(A / a[i]);
1555 vector<ex> anew = a;
1556 for ( size_t i=0; i<r; ++i ) {
1557 a[i] = a[i].subs(xnu == alphanu);
1559 ex cnew = c.subs(xnu == alphanu);
1560 vector<EvalPoint> Inew = I;
1562 vector<ex> sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1565 for ( size_t i=0; i<r; ++i ) {
1566 buf -= sigma[i] * b[i];
1569 e = make_modular(e, R);
1572 for ( size_t m=1; m<=d; ++m ) {
1573 while ( !e.is_zero() ) {
1574 monomial *= (xnu - alphanu);
1575 monomial = expand(monomial);
1576 ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1577 if ( !cm.is_zero() ) {
1578 vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1580 for ( size_t j=0; j<delta_s.size(); ++j ) {
1581 delta_s[j] *= monomial;
1582 sigma[j] += delta_s[j];
1583 buf -= delta_s[j] * b[j];
1586 e = make_modular(e, R);
1593 for ( size_t i=0; i<a.size(); ++i ) {
1594 UniPoly up(R, a[i], x);
1598 sigma.insert(sigma.begin(), r, 0);
1601 if ( is_a<add>(c) ) {
1610 for ( size_t i=0; i<nterms; ++i ) {
1612 int m = z.degree(x);
1614 cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1616 UniPolyVec delta_s = univar_diophant(amod, x, m, p, k);
1619 while ( poscm < 0 ) {
1620 poscm = poscm + expt_pos(cl_I(p),k);
1622 modcm = cl_MI(R, poscm);
1624 for ( size_t j=0; j<delta_s.size(); ++j ) {
1625 delta_s[j] = delta_s[j] * modcm;
1626 sigma[j] = sigma[j] + delta_s[j].to_ex(x);
1630 cout << "STEP " << i << " sigma ";
1631 for ( size_t p=0; p<sigma.size(); ++p ) {
1632 cout << sigma[p] << " ";
1643 for ( size_t i=0; i<sigma.size(); ++i ) {
1644 cout << sigma[i] << " ";
1649 for ( size_t i=0; i<sigma.size(); ++i ) {
1650 sigma[i] = make_modular(sigma[i], R);
1655 for ( size_t i=0; i<sigma.size(); ++i ) {
1656 cout << sigma[i] << " ";
1660 DCOUT(END multivar_diophant);
1665 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1667 for ( size_t i=0; i<v.size(); ++i ) {
1668 o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1672 #endif // def DEBUGFACTOR
1675 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)
1677 DCOUT(hensel_multivar);
1685 const size_t nu = I.size() + 1;
1686 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1693 for ( size_t j=nu; j>=2; --j ) {
1695 int alpha = I[j-2].evalpoint;
1696 A[j-2] = A[j-1].subs(x==alpha);
1697 A[j-2] = make_modular(A[j-2], R);
1702 for ( size_t i=0; i<A.size(); ++i) cout << A[i] << " ";
1706 int maxdeg = a.degree(I.front().x);
1707 for ( size_t i=1; i<I.size(); ++i ) {
1708 int maxdeg2 = a.degree(I[i].x);
1709 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
1713 const size_t n = u.size();
1716 for ( size_t i=0; i<n; ++i ) {
1717 U[i] = u[i].to_ex(x);
1721 for ( size_t i=0; i<U.size(); ++i) cout << U[i] << " ";
1725 for ( size_t j=2; j<=nu; ++j ) {
1729 for ( size_t m=0; m<n; ++m) {
1730 if ( lcU[m] != 1 ) {
1732 for ( size_t i=j-1; i<nu-1; ++i ) {
1733 coef = coef.subs(I[i].x == I[i].evalpoint);
1735 coef = expand(coef);
1736 coef = make_modular(coef, R);
1737 int deg = U[m].degree(x);
1738 U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
1742 for ( size_t i=0; i<n; ++i ) {
1745 ex e = expand(A[j-1] - Uprod);
1749 int alphaj = I[j-2].evalpoint;
1750 size_t deg = A[j-1].degree(xj);
1752 for ( size_t k=1; k<=deg; ++k ) {
1754 if ( !e.is_zero() ) {
1757 monomial *= (xj - alphaj);
1758 monomial = expand(monomial);
1760 ex dif = e.diff(ex_to<symbol>(xj), k);
1762 ex c = dif.subs(xj==alphaj) / factorial(k);
1764 if ( !c.is_zero() ) {
1765 vector<EvalPoint> newI = I;
1767 vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
1768 for ( size_t i=0; i<n; ++i ) {
1770 DCOUTVAR(deltaU[i]);
1771 deltaU[i] *= monomial;
1773 U[i] = make_modular(U[i], R);
1776 for ( size_t i=0; i<n; ++i ) {
1779 e = expand(A[j-1] - Uprod);
1780 e = make_modular(e, R);
1791 for ( size_t i=0; i<U.size(); ++i ) {
1795 if ( expand(a-acand).is_zero() ) {
1797 for ( size_t i=0; i<U.size(); ++i ) {
1801 DCOUT(END hensel_multivar);
1807 DCOUT(END hensel_multivar);
1812 static ex put_factors_into_lst(const ex& e)
1814 DCOUT(put_factors_into_lst);
1819 if ( is_a<numeric>(e) ) {
1821 DCOUT(END put_factors_into_lst);
1825 if ( is_a<power>(e) ) {
1827 result.append(e.op(0));
1828 result.append(e.op(1));
1829 DCOUT(END put_factors_into_lst);
1833 if ( is_a<symbol>(e) ) {
1837 DCOUT(END put_factors_into_lst);
1841 if ( is_a<mul>(e) ) {
1843 for ( size_t i=0; i<e.nops(); ++i ) {
1845 if ( is_a<numeric>(op) ) {
1848 if ( is_a<power>(op) ) {
1849 result.append(op.op(0));
1850 result.append(op.op(1));
1852 if ( is_a<symbol>(op) || is_a<add>(op) ) {
1857 result.prepend(nfac);
1858 DCOUT(END put_factors_into_lst);
1862 throw runtime_error("put_factors_into_lst: bad term.");
1865 static bool checkdivisors(const lst& f, vector<numeric>& d)
1867 const int k = f.nops()-2;
1869 d[0] = ex_to<numeric>(f.op(0) * f.op(f.nops()-1));
1870 for ( int i=1; i<=k; ++i ) {
1871 q = ex_to<numeric>(abs(f.op(i-1)));
1872 for ( int j=i-1; j>=0; --j ) {
1887 static void generate_set(const ex& u, const ex& vn, const exset& syms, const ex& f, const numeric& modulus, vector<numeric>& a, vector<numeric>& d)
1889 const ex& x = *syms.begin();
1895 exset::const_iterator s = syms.begin();
1897 for ( size_t i=0; i<a.size(); ++i ) {
1899 a[i] = mod(numeric(rand()), 2*modulus) - modulus;
1900 vnatry = vna.subs(*s == a[i]);
1901 } while ( vnatry == 0 );
1903 u0 = u0.subs(*s == a[i]);
1905 if ( gcd(u0,u0.diff(ex_to<symbol>(x))) != 1 ) {
1908 if ( is_a<numeric>(vn) ) {
1914 lst::const_iterator i = ex_to<lst>(f).begin();
1916 while ( i!=ex_to<lst>(f).end() ) {
1920 for ( size_t j=0; j<a.size(); ++j ) {
1921 fs = fs.subs(*s == a[j]);
1926 ex con = u0.content(x);
1928 trying = checkdivisors(fnum, d);
1934 ostream& operator<<(ostream& o, const vector<numeric>& v)
1936 for ( size_t i=0; i<v.size(); ++i ) {
1941 #endif // def DEBUGFACTOR
1943 static ex factor_multivariate(const ex& poly, const exset& syms)
1945 DCOUT(factor_multivariate);
1948 exset::const_iterator s;
1949 const ex& x = *syms.begin();
1952 /* make polynomial primitive */
1953 ex p = poly.expand().collect(x);
1955 ex cont = p.lcoeff(x);
1956 for ( numeric i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
1957 cont = gcd(cont, p.coeff(x,ex_to<numeric>(i).to_int()));
1958 if ( cont == 1 ) break;
1961 ex pp = expand(normal(p / cont));
1963 if ( !is_a<numeric>(cont) ) {
1964 return factor(cont) * factor(pp);
1967 /* factor leading coefficient */
1969 ex vn = p.lcoeff(x);
1971 if ( is_a<numeric>(vn) ) {
1975 ex vnfactors = factor(vn);
1976 vnlst = put_factors_into_lst(vnfactors);
1980 const numeric maxtrials = 3;
1981 numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3;
1982 numeric minimalr = -1;
1983 vector<numeric> a(syms.size()-1);
1984 vector<numeric> d(syms.size()-1);
1987 numeric trialcount = 0;
1991 while ( trialcount < maxtrials ) {
1993 generate_set(pp, vn, syms, vnlst, modulus, a, d);
1999 for ( size_t i=0; i<a.size(); ++i ) {
2000 u = u.subs(*s == a[i]);
2003 delta = u.content(x);
2005 // determine proper prime
2007 cl_modint_ring R = find_modint_ring(prime);
2009 if ( irem(ex_to<numeric>(u.lcoeff(x)), prime) != 0 ) {
2010 UniPoly modpoly(R, u, x);
2011 UniFactorVec sqrfree_ufv;
2012 squarefree(modpoly, sqrfree_ufv);
2013 if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break;
2015 prime = next_prime(prime);
2016 R = find_modint_ring(prime);
2019 UniPoly umod(R, u, x);
2021 factor_modular(umod, uvec);
2024 if ( uvec.size() == 1 ) {
2026 DCOUT(END factor_multivariate);
2030 if ( minimalr < 0 ) {
2031 minimalr = uvec.size();
2033 else if ( minimalr == uvec.size() ) {
2037 else if ( minimalr > uvec.size() ) {
2038 minimalr = uvec.size();
2041 DCOUTVAR(trialcount);
2043 if ( minimalr == 0 ) {
2045 DCOUT(END factor_multivariate);
2051 if ( vnlst.nops() == 1 ) {
2052 C.resize(uvec.size(), 1);
2056 vector<numeric> ftilde((vnlst.nops()-1)/2);
2057 for ( size_t i=0; i<ftilde.size(); ++i ) {
2058 ex ft = vnlst.op(i*2+1);
2061 for ( size_t j=0; j<a.size(); ++j ) {
2062 ft = ft.subs(*s == a[j]);
2065 ftilde[i] = ex_to<numeric>(ft);
2070 vector<bool> fflag(ftilde.size(), false);
2071 for ( size_t i=0; i<uvec.size(); ++i ) {
2072 ex ui = uvec[i].to_ex(x);
2074 numeric coeff = ex_to<numeric>(ui.lcoeff(x));
2075 for ( size_t j=0; j<ftilde.size(); ++j ) {
2076 if ( numeric(coeff / ftilde[j]).is_integer() ) {
2077 coeff = coeff / ftilde[j];
2083 D.push_back(Di.expand());
2085 for ( size_t i=0; i<fflag.size(); ++i ) {
2095 for ( size_t i=0; i<D.size(); ++i ) {
2099 for ( size_t j=0; j<a.size(); ++j ) {
2100 Dtilde = Dtilde.subs(*s == a[j]);
2103 ex Ci = D[i] * (uvec[i].to_ex(x).lcoeff(x) / Dtilde);
2108 for ( size_t i=0; i<D.size(); ++i ) {
2112 for ( size_t j=0; j<a.size(); ++j ) {
2113 Dtilde = Dtilde.subs(*s == a[j]);
2116 ex ui = uvec[i].to_ex(x);
2119 ex d = gcd(ui.lcoeff(x), Dtilde);
2120 Ci = D[i] * ( ui.lcoeff(x) / d );
2121 ui = ui * ( Dtilde[i] / d );
2122 delta = delta / ( Dtilde[i] / d );
2123 if ( delta == 1 ) break;
2126 pp = pp * pow(delta, D.size()-1);
2134 vector<EvalPoint> epv;
2137 for ( size_t i=0; i<a.size(); ++i ) {
2139 ep.evalpoint = a[i].to_int();
2145 for ( int i=u.degree(x); i>=u.ldegree(x); --i ) {
2146 maxcoeff += pow(abs(u.coeff(x, i)),2);
2148 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
2149 unsigned int maxdegree = 0;
2150 for ( size_t i=0; i<uvec.size(); ++i ) {
2151 if ( uvec[i].degree() > maxdegree ) {
2152 maxdegree = uvec[i].degree();
2155 unsigned int B = cl_I_to_uint(normmc * expt_pos(cl_I(2), maxdegree));
2157 ex res = hensel_multivar(poly, x, epv, prime, B, uvec, C);
2158 if ( res != lst() ) {
2160 for ( size_t i=0; i<res.nops(); ++i ) {
2161 result *= res.op(i).content(x) * res.op(i).unit(x);
2162 result *= res.op(i).primpart(x);
2165 DCOUT(END factor_multivariate);
2171 static ex factor_sqrfree(const ex& poly)
2173 // determine all symbols in poly
2174 FindSymbolsMap findsymbols;
2176 if ( findsymbols.syms.size() == 0 ) {
2180 if ( findsymbols.syms.size() == 1 ) {
2182 const ex& x = *(findsymbols.syms.begin());
2183 if ( poly.ldegree(x) > 0 ) {
2184 int ld = poly.ldegree(x);
2185 ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2186 return res * pow(x,ld);
2189 ex res = factor_univariate(poly, x);
2194 // multivariate case
2195 ex res = factor_multivariate(poly, findsymbols.syms);
2199 } // anonymous namespace
2201 ex factor(const ex& poly)
2203 // determine all symbols in poly
2204 FindSymbolsMap findsymbols;
2206 if ( findsymbols.syms.size() == 0 ) {
2210 exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
2211 for ( ; i!=end; ++i ) {
2215 // make poly square free
2216 ex sfpoly = sqrfree(poly, syms);
2218 // factorize the square free components
2219 if ( is_a<power>(sfpoly) ) {
2220 // case: (polynomial)^exponent
2221 const ex& base = sfpoly.op(0);
2222 if ( !is_a<add>(base) ) {
2223 // simple case: (monomial)^exponent
2226 ex f = factor_sqrfree(base);
2227 return pow(f, sfpoly.op(1));
2229 if ( is_a<mul>(sfpoly) ) {
2231 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
2232 const ex& t = sfpoly.op(i);
2233 if ( is_a<power>(t) ) {
2234 const ex& base = t.op(0);
2235 if ( !is_a<add>(base) ) {
2239 ex f = factor_sqrfree(base);
2240 res *= pow(f, t.op(1));
2243 else if ( is_a<add>(t) ) {
2244 ex f = factor_sqrfree(t);
2253 if ( is_a<symbol>(sfpoly) ) {
2256 // case: (polynomial)
2257 ex f = factor_sqrfree(sfpoly);
2261 } // namespace GiNaC