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 unsigned int maxdegree = (u1_.degree() > w1_.degree()) ? u1_.degree() : w1_.degree();
1052 unsigned int B = cl_I_to_uint(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 unsigned int modulus = p;
1083 while ( !e.is_zero() && modulus < 2*B*gamma_ui ) {
1085 phi = expand(s.to_ex(x)*c);
1086 UniPoly sigmatilde(R, phi, x);
1087 phi = expand(t.to_ex(x)*c);
1088 UniPoly tautilde(R, phi, x);
1090 div(sigmatilde, w1, q);
1091 rem(sigmatilde, w1, r);
1093 phi = expand(tautilde.to_ex(x) + q.to_ex(x) * u1.to_ex(x));
1094 UniPoly tau(R, phi, x);
1095 u = expand(u + tau.to_ex(x) * modulus);
1096 w = expand(w + sigma.to_ex(x) * modulus);
1097 e = expand(a - u * w);
1098 modulus = modulus * p;
1102 if ( e.is_zero() ) {
1103 ex delta = u.content(x);
1105 w = w / gamma * delta;
1113 static unsigned int next_prime(unsigned int p)
1115 static vector<unsigned int> primes;
1116 if ( primes.size() == 0 ) {
1117 primes.push_back(3); primes.push_back(5); primes.push_back(7);
1119 vector<unsigned int>::const_iterator it = primes.begin();
1120 if ( p >= primes.back() ) {
1121 unsigned int candidate = primes.back() + 2;
1123 size_t n = primes.size()/2;
1124 for ( size_t i=0; i<n; ++i ) {
1125 if ( candidate % primes[i] ) continue;
1129 primes.push_back(candidate);
1130 if ( candidate > p ) break;
1134 vector<unsigned int>::const_iterator end = primes.end();
1135 for ( ; it!=end; ++it ) {
1140 throw logic_error("next_prime: should not reach this point!");
1146 Partition(size_t n_) : n(n_)
1152 int operator[](size_t i) const { return k[i]; }
1153 size_t size() const { return n; }
1154 size_t size_first() const { return n-sum; }
1155 size_t size_second() const { return sum; }
1159 for ( size_t i=0; i<k.size(); ++i ) {
1160 cout << k[i] << " ";
1167 for ( size_t i=n-1; i>=1; --i ) {
1183 static void split(const UniPolyVec& factors, const Partition& part, UniPoly& a, UniPoly& b)
1185 a.set(0, a.R->one());
1186 b.set(0, a.R->one());
1187 for ( size_t i=0; i<part.size(); ++i ) {
1203 static ex factor_univariate(const ex& poly, const ex& x)
1205 ex unit, cont, prim;
1206 poly.unitcontprim(x, unit, cont, prim);
1208 // determine proper prime
1210 cl_modint_ring R = find_modint_ring(p);
1212 if ( irem(ex_to<numeric>(prim.lcoeff(x)), p) != 0 ) {
1213 UniPoly modpoly(R, prim, x);
1214 UniFactorVec sqrfree_ufv;
1215 squarefree(modpoly, sqrfree_ufv);
1216 if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break;
1219 R = find_modint_ring(p);
1222 // do modular factorization
1223 UniPoly modpoly(R, prim, x);
1225 factor_modular(modpoly, factors);
1226 if ( factors.size() <= 1 ) {
1227 // irreducible for sure
1231 // lift all factor combinations
1232 stack<ModFactors> tocheck;
1235 mf.factors = factors;
1238 while ( tocheck.size() ) {
1239 const size_t n = tocheck.top().factors.size();
1243 split(tocheck.top().factors, part, a, b);
1245 ex answer = hensel_univar(tocheck.top().poly, x, p, a, b);
1246 if ( answer != lst() ) {
1247 if ( part.size_first() == 1 ) {
1248 if ( part.size_second() == 1 ) {
1249 result *= answer.op(0) * answer.op(1);
1253 result *= answer.op(0);
1254 tocheck.top().poly = answer.op(1);
1255 for ( size_t i=0; i<n; ++i ) {
1256 if ( part[i] == 0 ) {
1257 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1263 else if ( part.size_second() == 1 ) {
1264 if ( part.size_first() == 1 ) {
1265 result *= answer.op(0) * answer.op(1);
1269 result *= answer.op(1);
1270 tocheck.top().poly = answer.op(0);
1271 for ( size_t i=0; i<n; ++i ) {
1272 if ( part[i] == 1 ) {
1273 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1280 UniPolyVec newfactors1(part.size_first(), R), newfactors2(part.size_second(), R);
1281 UniPolyVec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1282 for ( size_t i=0; i<n; ++i ) {
1284 *i2++ = tocheck.top().factors[i];
1287 *i1++ = tocheck.top().factors[i];
1290 tocheck.top().factors = newfactors1;
1291 tocheck.top().poly = answer.op(0);
1293 mf.factors = newfactors2;
1294 mf.poly = answer.op(1);
1300 if ( !part.next() ) {
1301 result *= tocheck.top().poly;
1309 return unit * cont * result;
1312 struct FindSymbolsMap : public map_function {
1314 ex operator()(const ex& e)
1316 if ( is_a<symbol>(e) ) {
1320 return e.map(*this);
1330 // forward declaration
1331 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);
1333 UniPolyVec multiterm_eea_lift(const UniPolyVec& a, const ex& x, unsigned int p, unsigned int k)
1335 DCOUT(multiterm_eea_lift);
1340 const size_t r = a.size();
1342 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1344 UniPolyVec q(r-1, fill);
1346 for ( size_t j=r-2; j>=1; --j ) {
1347 q[j-1] = a[j] * q[j];
1351 beta.set(0, R->one());
1353 for ( size_t j=1; j<r; ++j ) {
1356 vector<ex> mdarg(2);
1357 mdarg[0] = q[j-1].to_ex(x);
1358 mdarg[1] = a[j-1].to_ex(x);
1359 vector<EvalPoint> empty;
1360 vector<ex> exsigma = multivar_diophant(mdarg, x, beta.to_ex(x), empty, 0, p, k);
1361 UniPoly sigma1(R, exsigma[0], x);
1362 UniPoly sigma2(R, exsigma[1], x);
1364 s.push_back(sigma2);
1369 DCOUT(END multiterm_eea_lift);
1373 void eea_lift(const UniPoly& a, const UniPoly& b, const ex& x, unsigned int p, unsigned int k, UniPoly& s_, UniPoly& t_)
1382 cl_modint_ring R = find_modint_ring(p);
1388 UniPoly smod(R), tmod(R), g(R);
1389 exteuclid(amod, bmod, g, smod, tmod);
1395 cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1396 UniPoly s(Rpk, smod);
1397 UniPoly t(Rpk, tmod);
1404 one.set(0, Rpk->one());
1405 for ( size_t j=1; j<k; ++j ) {
1406 UniPoly e = one - a * s - b * t;
1409 UniPoly sigmabar(R);
1410 sigmabar = smod * c;
1414 div(sigmabar, bmod, q);
1416 rem(sigmabar, bmod, sigma);
1418 tau = taubar + q * amod;
1419 UniPoly sadd(Rpk, sigma);
1420 cl_MI modmodulus(Rpk, modulus);
1421 s = s + sadd * modmodulus;
1422 UniPoly tadd(Rpk, tau);
1423 t = t + tadd * modmodulus;
1424 modulus = modulus * p;
1431 DCOUT2(check, a*s + b*t);
1432 DCOUT(END eea_lift);
1435 UniPolyVec univar_diophant(const UniPolyVec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1437 DCOUT(univar_diophant);
1444 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1446 const size_t r = a.size();
1449 UniPolyVec s = multiterm_eea_lift(a, x, p, k);
1450 for ( size_t j=0; j<r; ++j ) {
1451 ex phi = expand(pow(x,m)*s[j].to_ex(x));
1452 UniPoly bmod(R, phi, x);
1454 rem(bmod, a[j], buf);
1455 result.push_back(buf);
1460 eea_lift(a[1], a[0], x, p, k, s, t);
1461 ex phi = expand(pow(x,m)*s.to_ex(x));
1462 UniPoly bmod(R, phi, x);
1464 rem(bmod, a[0], buf);
1465 result.push_back(buf);
1468 phi = expand(pow(x,m)*t.to_ex(x));
1469 UniPoly t1mod(R, phi, x);
1470 buf = t1mod + q * a[1];
1471 result.push_back(buf);
1475 DCOUT(END univar_diophant);
1479 struct make_modular_map : public map_function {
1481 make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1482 ex operator()(const ex& e)
1484 if ( is_a<add>(e) || is_a<mul>(e) ) {
1485 return e.map(*this);
1487 else if ( is_a<numeric>(e) ) {
1488 numeric mod(R->modulus);
1489 numeric halfmod = (mod-1)/2;
1490 cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1491 numeric n(R->retract(emod));
1492 if ( n > halfmod ) {
1503 static ex make_modular(const ex& e, const cl_modint_ring& R)
1505 make_modular_map map(R);
1509 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)
1513 DCOUT(multivar_diophant);
1516 for ( size_t i=0; i<a.size(); ++i ) {
1517 cout << a[i] << " ";
1525 for ( size_t i=0; i<I.size(); ++i ) {
1526 cout << I[i].x << "=" << I[i].evalpoint << " ";
1534 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1535 const size_t r = a.size();
1536 const size_t nu = I.size() + 1;
1542 ex xnu = I.back().x;
1543 int alphanu = I.back().evalpoint;
1546 for ( size_t i=0; i<r; ++i ) {
1550 for ( size_t i=0; i<r; ++i ) {
1551 b[i] = normal(A / a[i]);
1554 vector<ex> anew = a;
1555 for ( size_t i=0; i<r; ++i ) {
1556 a[i] = a[i].subs(xnu == alphanu);
1558 ex cnew = c.subs(xnu == alphanu);
1559 vector<EvalPoint> Inew = I;
1561 vector<ex> sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1564 for ( size_t i=0; i<r; ++i ) {
1565 buf -= sigma[i] * b[i];
1568 e = make_modular(e, R);
1571 for ( size_t m=1; m<=d; ++m ) {
1572 while ( !e.is_zero() ) {
1573 monomial *= (xnu - alphanu);
1574 monomial = expand(monomial);
1575 ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1576 if ( !cm.is_zero() ) {
1577 vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1579 for ( size_t j=0; j<delta_s.size(); ++j ) {
1580 delta_s[j] *= monomial;
1581 sigma[j] += delta_s[j];
1582 buf -= delta_s[j] * b[j];
1585 e = make_modular(e, R);
1592 for ( size_t i=0; i<a.size(); ++i ) {
1593 UniPoly up(R, a[i], x);
1597 sigma.insert(sigma.begin(), r, 0);
1600 if ( is_a<add>(c) ) {
1609 for ( size_t i=0; i<nterms; ++i ) {
1611 int m = z.degree(x);
1613 cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1615 UniPolyVec delta_s = univar_diophant(amod, x, m, p, k);
1618 while ( poscm < 0 ) {
1619 poscm = poscm + expt_pos(cl_I(p),k);
1621 modcm = cl_MI(R, poscm);
1623 for ( size_t j=0; j<delta_s.size(); ++j ) {
1624 delta_s[j] = delta_s[j] * modcm;
1625 sigma[j] = sigma[j] + delta_s[j].to_ex(x);
1629 cout << "STEP " << i << " sigma ";
1630 for ( size_t p=0; p<sigma.size(); ++p ) {
1631 cout << sigma[p] << " ";
1642 for ( size_t i=0; i<sigma.size(); ++i ) {
1643 cout << sigma[i] << " ";
1648 for ( size_t i=0; i<sigma.size(); ++i ) {
1649 sigma[i] = make_modular(sigma[i], R);
1654 for ( size_t i=0; i<sigma.size(); ++i ) {
1655 cout << sigma[i] << " ";
1659 DCOUT(END multivar_diophant);
1664 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1666 for ( size_t i=0; i<v.size(); ++i ) {
1667 o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1671 #endif // def DEBUGFACTOR
1674 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)
1676 DCOUT(hensel_multivar);
1684 const size_t nu = I.size() + 1;
1685 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1692 for ( size_t j=nu; j>=2; --j ) {
1694 int alpha = I[j-2].evalpoint;
1695 A[j-2] = A[j-1].subs(x==alpha);
1696 A[j-2] = make_modular(A[j-2], R);
1701 for ( size_t i=0; i<A.size(); ++i) cout << A[i] << " ";
1705 int maxdeg = a.degree(I.front().x);
1706 for ( size_t i=1; i<I.size(); ++i ) {
1707 int maxdeg2 = a.degree(I[i].x);
1708 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
1712 const size_t n = u.size();
1715 for ( size_t i=0; i<n; ++i ) {
1716 U[i] = u[i].to_ex(x);
1720 for ( size_t i=0; i<U.size(); ++i) cout << U[i] << " ";
1724 for ( size_t j=2; j<=nu; ++j ) {
1728 for ( size_t m=0; m<n; ++m) {
1729 if ( lcU[m] != 1 ) {
1731 for ( size_t i=j-1; i<nu-1; ++i ) {
1732 coef = coef.subs(I[i].x == I[i].evalpoint);
1734 coef = expand(coef);
1735 coef = make_modular(coef, R);
1736 int deg = U[m].degree(x);
1737 U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
1741 for ( size_t i=0; i<n; ++i ) {
1744 ex e = expand(A[j-1] - Uprod);
1748 int alphaj = I[j-2].evalpoint;
1749 size_t deg = A[j-1].degree(xj);
1751 for ( size_t k=1; k<=deg; ++k ) {
1753 if ( !e.is_zero() ) {
1756 monomial *= (xj - alphaj);
1757 monomial = expand(monomial);
1759 ex dif = e.diff(ex_to<symbol>(xj), k);
1761 ex c = dif.subs(xj==alphaj) / factorial(k);
1763 if ( !c.is_zero() ) {
1764 vector<EvalPoint> newI = I;
1766 vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
1767 for ( size_t i=0; i<n; ++i ) {
1769 DCOUTVAR(deltaU[i]);
1770 deltaU[i] *= monomial;
1772 U[i] = make_modular(U[i], R);
1775 for ( size_t i=0; i<n; ++i ) {
1778 e = expand(A[j-1] - Uprod);
1779 e = make_modular(e, R);
1790 for ( size_t i=0; i<U.size(); ++i ) {
1794 if ( expand(a-acand).is_zero() ) {
1796 for ( size_t i=0; i<U.size(); ++i ) {
1800 DCOUT(END hensel_multivar);
1806 DCOUT(END hensel_multivar);
1811 static ex put_factors_into_lst(const ex& e)
1813 DCOUT(put_factors_into_lst);
1818 if ( is_a<numeric>(e) ) {
1820 DCOUT(END put_factors_into_lst);
1824 if ( is_a<power>(e) ) {
1826 result.append(e.op(0));
1827 result.append(e.op(1));
1828 DCOUT(END put_factors_into_lst);
1832 if ( is_a<symbol>(e) ) {
1836 DCOUT(END put_factors_into_lst);
1840 if ( is_a<mul>(e) ) {
1842 for ( size_t i=0; i<e.nops(); ++i ) {
1844 if ( is_a<numeric>(op) ) {
1847 if ( is_a<power>(op) ) {
1848 result.append(op.op(0));
1849 result.append(op.op(1));
1851 if ( is_a<symbol>(op) || is_a<add>(op) ) {
1856 result.prepend(nfac);
1857 DCOUT(END put_factors_into_lst);
1861 throw runtime_error("put_factors_into_lst: bad term.");
1864 static bool checkdivisors(const lst& f, vector<numeric>& d)
1866 const int k = f.nops()-2;
1868 d[0] = ex_to<numeric>(f.op(0) * f.op(f.nops()-1));
1869 for ( int i=1; i<=k; ++i ) {
1870 q = ex_to<numeric>(abs(f.op(i-1)));
1871 for ( int j=i-1; j>=0; --j ) {
1886 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)
1888 const ex& x = *syms.begin();
1894 exset::const_iterator s = syms.begin();
1896 for ( size_t i=0; i<a.size(); ++i ) {
1898 a[i] = mod(numeric(rand()), 2*modulus) - modulus;
1899 vnatry = vna.subs(*s == a[i]);
1900 } while ( vnatry == 0 );
1902 u0 = u0.subs(*s == a[i]);
1904 if ( gcd(u0,u0.diff(ex_to<symbol>(x))) != 1 ) {
1907 if ( is_a<numeric>(vn) ) {
1913 lst::const_iterator i = ex_to<lst>(f).begin();
1915 while ( i!=ex_to<lst>(f).end() ) {
1919 for ( size_t j=0; j<a.size(); ++j ) {
1920 fs = fs.subs(*s == a[j]);
1925 ex con = u0.content(x);
1927 trying = checkdivisors(fnum, d);
1933 ostream& operator<<(ostream& o, const vector<numeric>& v)
1935 for ( size_t i=0; i<v.size(); ++i ) {
1940 #endif // def DEBUGFACTOR
1942 static ex factor_multivariate(const ex& poly, const exset& syms)
1944 DCOUT(factor_multivariate);
1947 exset::const_iterator s;
1948 const ex& x = *syms.begin();
1951 /* make polynomial primitive */
1952 ex p = poly.expand().collect(x);
1954 ex cont = p.lcoeff(x);
1955 for ( numeric i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
1956 cont = gcd(cont, p.coeff(x,ex_to<numeric>(i).to_int()));
1957 if ( cont == 1 ) break;
1960 ex pp = expand(normal(p / cont));
1962 if ( !is_a<numeric>(cont) ) {
1963 return factor(cont) * factor(pp);
1966 /* factor leading coefficient */
1968 ex vn = p.lcoeff(x);
1970 if ( is_a<numeric>(vn) ) {
1974 ex vnfactors = factor(vn);
1975 vnlst = put_factors_into_lst(vnfactors);
1979 const numeric maxtrials = 3;
1980 numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3;
1981 numeric minimalr = -1;
1982 vector<numeric> a(syms.size()-1);
1983 vector<numeric> d(syms.size()-1);
1986 numeric trialcount = 0;
1990 while ( trialcount < maxtrials ) {
1992 generate_set(pp, vn, syms, vnlst, modulus, a, d);
1998 for ( size_t i=0; i<a.size(); ++i ) {
1999 u = u.subs(*s == a[i]);
2002 delta = u.content(x);
2004 // determine proper prime
2006 cl_modint_ring R = find_modint_ring(prime);
2008 if ( irem(ex_to<numeric>(u.lcoeff(x)), prime) != 0 ) {
2009 UniPoly modpoly(R, u, x);
2010 UniFactorVec sqrfree_ufv;
2011 squarefree(modpoly, sqrfree_ufv);
2012 if ( sqrfree_ufv.factors.size() == 1 && sqrfree_ufv.factors.front().exp == 1 ) break;
2014 prime = next_prime(prime);
2015 R = find_modint_ring(prime);
2018 UniPoly umod(R, u, x);
2020 factor_modular(umod, uvec);
2023 if ( uvec.size() == 1 ) {
2025 DCOUT(END factor_multivariate);
2029 if ( minimalr < 0 ) {
2030 minimalr = uvec.size();
2032 else if ( minimalr == uvec.size() ) {
2036 else if ( minimalr > uvec.size() ) {
2037 minimalr = uvec.size();
2040 DCOUTVAR(trialcount);
2042 if ( minimalr == 0 ) {
2044 DCOUT(END factor_multivariate);
2050 if ( vnlst.nops() == 1 ) {
2051 C.resize(uvec.size(), 1);
2055 vector<numeric> ftilde((vnlst.nops()-1)/2);
2056 for ( size_t i=0; i<ftilde.size(); ++i ) {
2057 ex ft = vnlst.op(i*2+1);
2060 for ( size_t j=0; j<a.size(); ++j ) {
2061 ft = ft.subs(*s == a[j]);
2064 ftilde[i] = ex_to<numeric>(ft);
2069 vector<bool> fflag(ftilde.size(), false);
2070 for ( size_t i=0; i<uvec.size(); ++i ) {
2071 ex ui = uvec[i].to_ex(x);
2073 numeric coeff = ex_to<numeric>(ui.lcoeff(x));
2074 for ( size_t j=0; j<ftilde.size(); ++j ) {
2075 if ( numeric(coeff / ftilde[j]).is_integer() ) {
2076 coeff = coeff / ftilde[j];
2082 D.push_back(Di.expand());
2084 for ( size_t i=0; i<fflag.size(); ++i ) {
2094 for ( size_t i=0; i<D.size(); ++i ) {
2098 for ( size_t j=0; j<a.size(); ++j ) {
2099 Dtilde = Dtilde.subs(*s == a[j]);
2102 ex Ci = D[i] * (uvec[i].to_ex(x).lcoeff(x) / Dtilde);
2107 for ( size_t i=0; i<D.size(); ++i ) {
2111 for ( size_t j=0; j<a.size(); ++j ) {
2112 Dtilde = Dtilde.subs(*s == a[j]);
2115 ex ui = uvec[i].to_ex(x);
2118 ex d = gcd(ui.lcoeff(x), Dtilde);
2119 Ci = D[i] * ( ui.lcoeff(x) / d );
2120 ui = ui * ( Dtilde[i] / d );
2121 delta = delta / ( Dtilde[i] / d );
2122 if ( delta == 1 ) break;
2125 pp = pp * pow(delta, D.size()-1);
2133 vector<EvalPoint> epv;
2136 for ( size_t i=0; i<a.size(); ++i ) {
2138 ep.evalpoint = a[i].to_int();
2144 for ( int i=u.degree(x); i>=u.ldegree(x); --i ) {
2145 maxcoeff += pow(abs(u.coeff(x, i)),2);
2147 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
2148 unsigned int maxdegree = 0;
2149 for ( size_t i=0; i<uvec.size(); ++i ) {
2150 if ( uvec[i].degree() > maxdegree ) {
2151 maxdegree = uvec[i].degree();
2154 unsigned int B = cl_I_to_uint(normmc * expt_pos(cl_I(2), maxdegree));
2156 ex res = hensel_multivar(poly, x, epv, prime, B, uvec, C);
2157 if ( res != lst() ) {
2159 for ( size_t i=0; i<res.nops(); ++i ) {
2160 result *= res.op(i).content(x) * res.op(i).unit(x);
2161 result *= res.op(i).primpart(x);
2164 DCOUT(END factor_multivariate);
2170 static ex factor_sqrfree(const ex& poly)
2172 // determine all symbols in poly
2173 FindSymbolsMap findsymbols;
2175 if ( findsymbols.syms.size() == 0 ) {
2179 if ( findsymbols.syms.size() == 1 ) {
2181 const ex& x = *(findsymbols.syms.begin());
2182 if ( poly.ldegree(x) > 0 ) {
2183 int ld = poly.ldegree(x);
2184 ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2185 return res * pow(x,ld);
2188 ex res = factor_univariate(poly, x);
2193 // multivariate case
2194 ex res = factor_multivariate(poly, findsymbols.syms);
2198 } // anonymous namespace
2200 ex factor(const ex& poly)
2202 // determine all symbols in poly
2203 FindSymbolsMap findsymbols;
2205 if ( findsymbols.syms.size() == 0 ) {
2209 exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
2210 for ( ; i!=end; ++i ) {
2214 // make poly square free
2215 ex sfpoly = sqrfree(poly, syms);
2217 // factorize the square free components
2218 if ( is_a<power>(sfpoly) ) {
2219 // case: (polynomial)^exponent
2220 const ex& base = sfpoly.op(0);
2221 if ( !is_a<add>(base) ) {
2222 // simple case: (monomial)^exponent
2225 ex f = factor_sqrfree(base);
2226 return pow(f, sfpoly.op(1));
2228 if ( is_a<mul>(sfpoly) ) {
2230 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
2231 const ex& t = sfpoly.op(i);
2232 if ( is_a<power>(t) ) {
2233 const ex& base = t.op(0);
2234 if ( !is_a<add>(base) ) {
2238 ex f = factor_sqrfree(base);
2239 res *= pow(f, t.op(1));
2242 else if ( is_a<add>(t) ) {
2243 ex f = factor_sqrfree(t);
2252 if ( is_a<symbol>(sfpoly) ) {
2255 // case: (polynomial)
2256 ex f = factor_sqrfree(sfpoly);
2260 } // namespace GiNaC