3 * Polynomial factorization code (implementation).
5 * Algorithms used can be found in
6 * [W1] An Improved Multivariate Polynomial Factoring Algorithm,
7 * P.S.Wang, Mathematics of Computation, Vol. 32, No. 144 (1978) 1215--1231.
8 * [GCL] Algorithms for Computer Algebra,
9 * K.O.Geddes, S.R.Czapor, G.Labahn, Springer Verlag, 1992.
13 * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
15 * This program is free software; you can redistribute it and/or modify
16 * it under the terms of the GNU General Public License as published by
17 * the Free Software Foundation; either version 2 of the License, or
18 * (at your option) any later version.
20 * This program is distributed in the hope that it will be useful,
21 * but WITHOUT ANY WARRANTY; without even the implied warranty of
22 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 * GNU General Public License for more details.
25 * You should have received a copy of the GNU General Public License
26 * along with this program; if not, write to the Free Software
27 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
34 #include <ginac/ginac.h>
35 using namespace GiNaC;
41 #include "operators.h"
44 #include "relational.h"
67 #define DCOUT(str) cout << #str << endl
68 #define DCOUTVAR(var) cout << #var << ": " << var << endl
69 #define DCOUT2(str,var) cout << #str << ": " << var << endl
73 #define DCOUT2(str,var)
76 // forward declaration
77 ex factor(const ex& poly, unsigned options);
79 // anonymous namespace to hide all utility functions
82 typedef vector<cl_MI> mvec;
85 ostream& operator<<(ostream& o, const mvec& v)
87 mvec::const_iterator i = v.begin(), end = v.end();
93 #endif // def DEBUGFACTOR
96 ostream& operator<<(ostream& o, const vector<mvec>& v)
98 vector<mvec>::const_iterator i = v.begin(), end = v.end();
104 #endif // def DEBUGFACTOR
106 ////////////////////////////////////////////////////////////////////////////////
107 // modular univariate polynomial code
109 typedef cl_UP_MI umod;
110 typedef vector<umod> umodvec;
112 #define COPY(to,from) from.ring()->create(degree(from)); \
113 for ( int II=0; II<=degree(from); ++II ) to.set_coeff(II, coeff(from, II)); \
117 ostream& operator<<(ostream& o, const umodvec& v)
119 umodvec::const_iterator i = v.begin(), end = v.end();
121 o << *i++ << " , " << endl;
125 #endif // def DEBUGFACTOR
127 static umod umod_from_ex(const ex& e, const ex& x, const cl_univpoly_modint_ring& UPR)
129 // assert: e is in Z[x]
130 int deg = e.degree(x);
131 umod p = UPR->create(deg);
132 int ldeg = e.ldegree(x);
133 for ( ; deg>=ldeg; --deg ) {
134 cl_I coeff = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
135 p.set_coeff(deg, UPR->basering()->canonhom(coeff));
137 for ( ; deg>=0; --deg ) {
138 p.set_coeff(deg, UPR->basering()->zero());
144 static umod umod_from_ex(const ex& e, const ex& x, const cl_modint_ring& R)
146 return umod_from_ex(e, x, find_univpoly_ring(R));
149 static umod umod_from_ex(const ex& e, const ex& x, const cl_I& modulus)
151 return umod_from_ex(e, x, find_modint_ring(modulus));
154 static umod umod_from_modvec(const mvec& mv)
156 size_t n = mv.size(); // assert: n>0
157 while ( n && zerop(mv[n-1]) ) --n;
158 cl_univpoly_modint_ring UPR = find_univpoly_ring(mv.front().ring());
160 umod p = UPR->create(-1);
164 umod p = UPR->create(n-1);
165 for ( size_t i=0; i<n; ++i ) {
166 p.set_coeff(i, mv[i]);
172 static umod divide(const umod& a, const cl_I& x)
176 cl_univpoly_modint_ring UPR = a.ring();
177 cl_modint_ring R = UPR->basering();
179 umod newa = UPR->create(deg);
180 for ( int i=0; i<=deg; ++i ) {
181 cl_I c = R->retract(coeff(a, i));
182 newa.set_coeff(i, cl_MI(R, the<cl_I>(c / x)));
189 static ex umod_to_ex(const umod& a, const ex& x)
192 cl_modint_ring R = a.ring()->basering();
193 cl_I mod = R->modulus;
194 cl_I halfmod = (mod-1) >> 1;
195 for ( int i=degree(a); i>=0; --i ) {
196 cl_I n = R->retract(coeff(a, i));
198 e += numeric(n-mod) * pow(x, i);
200 e += numeric(n) * pow(x, i);
206 static void unit_normal(umod& a)
210 cl_MI lc = coeff(a, deg);
211 cl_MI one = a.ring()->basering()->one();
213 umod newa = a.ring()->create(deg);
214 newa.set_coeff(deg, one);
215 for ( --deg; deg>=0; --deg ) {
216 cl_MI nc = div(coeff(a, deg), lc);
217 newa.set_coeff(deg, nc);
225 static umod rem(const umod& a, const umod& b)
237 cl_MI qk = div(coeff(c, n+k), coeff(b, n));
240 for ( int i=0; i<n; ++i ) {
242 c.set_coeff(j, coeff(c, j) - qk * coeff(b, j-k));
247 cl_MI zero = a.ring()->basering()->zero();
248 for ( int i=degree(a); i>=n; --i ) {
249 c.set_coeff(i, zero);
256 static umod div(const umod& a, const umod& b)
262 umod q = a.ring()->create(-1);
268 umod q = a.ring()->create(k);
270 cl_MI qk = div(coeff(c, n+k), coeff(b, n));
274 for ( int i=0; i<n; ++i ) {
276 c.set_coeff(j, coeff(c, j) - qk * coeff(b, j-k));
285 static umod remdiv(const umod& a, const umod& b, umod& q)
291 q = a.ring()->create(-1);
298 q = a.ring()->create(k);
300 cl_MI qk = div(coeff(c, n+k), coeff(b, n));
304 for ( int i=0; i<n; ++i ) {
306 c.set_coeff(j, coeff(c, j) - qk * coeff(b, j-k));
311 cl_MI zero = a.ring()->basering()->zero();
312 for ( int i=degree(a); i>=n; --i ) {
313 c.set_coeff(i, zero);
321 static umod gcd(const umod& a, const umod& b)
323 if ( degree(a) < degree(b) ) return gcd(b, a);
329 while ( !zerop(d) ) {
338 static bool squarefree(const umod& a)
344 umod one = a.ring()->one();
349 // END modular univariate polynomial code
350 ////////////////////////////////////////////////////////////////////////////////
352 ////////////////////////////////////////////////////////////////////////////////
357 friend ostream& operator<<(ostream& o, const modular_matrix& m);
359 modular_matrix(size_t r_, size_t c_, const cl_MI& init) : r(r_), c(c_)
363 size_t rowsize() const { return r; }
364 size_t colsize() const { return c; }
365 cl_MI& operator()(size_t row, size_t col) { return m[row*c + col]; }
366 cl_MI operator()(size_t row, size_t col) const { return m[row*c + col]; }
367 void mul_col(size_t col, const cl_MI x)
369 mvec::iterator i = m.begin() + col;
370 for ( size_t rc=0; rc<r; ++rc ) {
375 void sub_col(size_t col1, size_t col2, const cl_MI fac)
377 mvec::iterator i1 = m.begin() + col1;
378 mvec::iterator i2 = m.begin() + col2;
379 for ( size_t rc=0; rc<r; ++rc ) {
380 *i1 = *i1 - *i2 * fac;
385 void switch_col(size_t col1, size_t col2)
388 mvec::iterator i1 = m.begin() + col1;
389 mvec::iterator i2 = m.begin() + col2;
390 for ( size_t rc=0; rc<r; ++rc ) {
391 buf = *i1; *i1 = *i2; *i2 = buf;
396 void mul_row(size_t row, const cl_MI x)
398 vector<cl_MI>::iterator i = m.begin() + row*c;
399 for ( size_t cc=0; cc<c; ++cc ) {
404 void sub_row(size_t row1, size_t row2, const cl_MI fac)
406 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
407 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
408 for ( size_t cc=0; cc<c; ++cc ) {
409 *i1 = *i1 - *i2 * fac;
414 void switch_row(size_t row1, size_t row2)
417 vector<cl_MI>::iterator i1 = m.begin() + row1*c;
418 vector<cl_MI>::iterator i2 = m.begin() + row2*c;
419 for ( size_t cc=0; cc<c; ++cc ) {
420 buf = *i1; *i1 = *i2; *i2 = buf;
425 bool is_col_zero(size_t col) const
427 mvec::const_iterator i = m.begin() + col;
428 for ( size_t rr=0; rr<r; ++rr ) {
436 bool is_row_zero(size_t row) const
438 mvec::const_iterator i = m.begin() + row*c;
439 for ( size_t cc=0; cc<c; ++cc ) {
447 void set_row(size_t row, const vector<cl_MI>& newrow)
449 mvec::iterator i1 = m.begin() + row*c;
450 mvec::const_iterator i2 = newrow.begin(), end = newrow.end();
451 for ( ; i2 != end; ++i1, ++i2 ) {
455 mvec::const_iterator row_begin(size_t row) const { return m.begin()+row*c; }
456 mvec::const_iterator row_end(size_t row) const { return m.begin()+row*c+r; }
463 modular_matrix operator*(const modular_matrix& m1, const modular_matrix& m2)
465 const unsigned int r = m1.rowsize();
466 const unsigned int c = m2.colsize();
467 modular_matrix o(r,c,m1(0,0));
469 for ( size_t i=0; i<r; ++i ) {
470 for ( size_t j=0; j<c; ++j ) {
472 buf = m1(i,0) * m2(0,j);
473 for ( size_t k=1; k<c; ++k ) {
474 buf = buf + m1(i,k)*m2(k,j);
482 ostream& operator<<(ostream& o, const modular_matrix& m)
484 vector<cl_MI>::const_iterator i = m.m.begin(), end = m.m.end();
486 for ( ; i != end; ++i ) {
488 if ( !(wrap++ % m.c) ) {
495 #endif // def DEBUGFACTOR
497 // END modular matrix
498 ////////////////////////////////////////////////////////////////////////////////
500 static void q_matrix(const umod& a, modular_matrix& Q)
503 unsigned int q = cl_I_to_uint(a.ring()->basering()->modulus);
505 // vector<cl_MI> r(n, a.R->zero());
506 // r[0] = a.R->one();
508 // unsigned int max = (n-1) * q;
509 // for ( size_t m=1; m<=max; ++m ) {
510 // cl_MI rn_1 = r.back();
511 // for ( size_t i=n-1; i>0; --i ) {
512 // r[i] = r[i-1] - rn_1 * a[i];
514 // r[0] = -rn_1 * a[0];
515 // if ( (m % q) == 0 ) {
516 // Q.set_row(m/q, r);
519 // slow and (hopefully) correct
520 cl_MI one = a.ring()->basering()->one();
521 for ( int i=0; i<n; ++i ) {
522 umod qk = a.ring()->create(i*q);
523 qk.set_coeff(i*q, one);
527 for ( int j=0; j<n; ++j ) {
528 rvec.push_back(coeff(r, j));
534 static void nullspace(modular_matrix& M, vector<mvec>& basis)
536 const size_t n = M.rowsize();
537 const cl_MI one = M(0,0).ring()->one();
538 for ( size_t i=0; i<n; ++i ) {
539 M(i,i) = M(i,i) - one;
541 for ( size_t r=0; r<n; ++r ) {
543 for ( ; cc<n; ++cc ) {
544 if ( !zerop(M(r,cc)) ) {
546 if ( !zerop(M(cc,cc)) ) {
558 M.mul_col(r, recip(M(r,r)));
559 for ( cc=0; cc<n; ++cc ) {
561 M.sub_col(cc, r, M(r,cc));
567 for ( size_t i=0; i<n; ++i ) {
568 M(i,i) = M(i,i) - one;
570 for ( size_t i=0; i<n; ++i ) {
571 if ( !M.is_row_zero(i) ) {
572 mvec nu(M.row_begin(i), M.row_end(i));
578 static void berlekamp(const umod& a, umodvec& upv)
580 cl_modint_ring R = a.ring()->basering();
581 const umod one = a.ring()->one();
583 modular_matrix Q(degree(a), degree(a), R->zero());
587 const unsigned int k = nu.size();
593 factors.push_back(a);
594 unsigned int size = 1;
596 unsigned int q = cl_I_to_uint(R->modulus);
598 list<umod>::iterator u = factors.begin();
601 for ( unsigned int s=0; s<q; ++s ) {
602 umod nur = umod_from_modvec(nu[r]);
603 cl_MI buf = coeff(nur, 0) - cl_MI(R, s);
604 nur.set_coeff(0, buf);
606 umod g = gcd(nur, *u);
607 if ( g != one && g != *u ) {
608 umod uo = div(*u, g);
610 throw logic_error("berlekamp: unexpected divisor.");
615 factors.push_back(g);
617 list<umod>::const_iterator i = factors.begin(), end = factors.end();
619 if ( degree(*i) ) ++size;
623 list<umod>::const_iterator i = factors.begin(), end = factors.end();
638 static umod rem_xq(int q, const umod& b)
640 cl_univpoly_modint_ring UPR = b.ring();
641 cl_modint_ring R = UPR->basering();
645 umod c = UPR->create(q);
646 c.set_coeff(q, R->one());
651 mvec c(n+1, R->zero());
658 cl_MI qk = div(c[n-ofs], coeff(b, n));
660 for ( int i=1; i<=n; ++i ) {
661 c[n-i+ofs] = c[n-i] - qk * coeff(b, n-i);
675 umod res = umod_from_modvec(c);
679 static void distinct_degree_factor(const umod& a_, umodvec& result)
681 umod a = COPY(a, a_);
683 DCOUT(distinct_degree_factor);
686 cl_univpoly_modint_ring UPR = a.ring();
687 cl_modint_ring R = UPR->basering();
688 int q = cl_I_to_int(R->modulus);
694 umod w = UPR->create(1);
695 w.set_coeff(1, R->one());
701 while ( i <= nhalf ) {
705 ai.push_back(gcd(a, w-x));
707 if ( ai.back() != UPR->one() ) {
708 a = div(a, ai.back());
717 DCOUT(END distinct_degree_factor);
720 static void same_degree_factor(const umod& a, umodvec& result)
722 DCOUT(same_degree_factor);
724 cl_univpoly_modint_ring UPR = a.ring();
725 cl_modint_ring R = UPR->basering();
729 distinct_degree_factor(a, buf);
732 for ( size_t i=0; i<buf.size(); ++i ) {
733 if ( buf[i] != UPR->one() ) {
734 degsum += degree(buf[i]);
736 berlekamp(buf[i], upv);
737 for ( size_t j=0; j<upv.size(); ++j ) {
738 result.push_back(upv[j]);
743 if ( degsum < deg ) {
748 DCOUT(END same_degree_factor);
751 static void distinct_degree_factor_BSGS(const umod& a, umodvec& result)
753 DCOUT(distinct_degree_factor_BSGS);
756 cl_univpoly_modint_ring UPR = a.ring();
757 cl_modint_ring R = UPR->basering();
758 int q = cl_I_to_int(R->modulus);
762 int l = cl_I_to_int(ceiling1(the<cl_F>(expt(n, pm))));
764 umodvec h(l+1, UPR->create(-1));
765 umod qk = UPR->create(1);
766 qk.set_coeff(1, R->one());
770 for ( int i=1; i<=l; ++i ) {
771 qk = expt_pos(h[i-1], q);
777 int m = std::ceil(((double)n)/2/l);
779 umodvec H(m, UPR->create(-1));
780 int ql = std::pow(q, l);
781 H[0] = COPY(H[0], h[l]);
783 for ( int i=1; i<m; ++i ) {
784 qk = expt_pos(H[i-1], ql);
790 umodvec I(m, UPR->create(-1));
791 for ( int i=0; i<m; ++i ) {
793 for ( int j=0; j<l; ++j ) {
794 I[i] = I[i] * (H[i] - h[j]);
802 umodvec F(m, UPR->one());
804 for ( int i=0; i<m; ++i ) {
806 umod g = gcd(f, I[i]);
807 if ( g == UPR->one() ) continue;
813 result.resize(n, UPR->one());
814 if ( f != UPR->one() ) {
817 for ( int i=0; i<m; ++i ) {
819 umod f = COPY(f, F[i]);
820 for ( int j=l-1; j>=0; --j ) {
821 umod g = gcd(f, H[i]-h[j]);
822 result[l*(i+1)-j-1] = g;
828 DCOUT(END distinct_degree_factor_BSGS);
831 static void cantor_zassenhaus(const umod& a, umodvec& result)
835 static void factor_modular(const umod& p, umodvec& upv)
837 //same_degree_factor(p, upv);
842 static void exteuclid(const umod& a, const umod& b, umod& g, umod& s, umod& t)
844 if ( degree(a) < degree(b) ) {
845 exteuclid(b, a, g, t, s);
848 umod c = COPY(c, a); unit_normal(c);
849 umod d = COPY(d, b); unit_normal(d);
850 umod c1 = a.ring()->one();
851 umod c2 = a.ring()->create(-1);
852 umod d1 = a.ring()->create(-1);
853 umod d2 = a.ring()->one();
854 while ( !zerop(d) ) {
857 umod r1 = c1 - q * d1;
858 umod r2 = c2 - q * d2;
866 g = COPY(g, c); unit_normal(g);
868 for ( int i=0; i<=degree(s); ++i ) {
869 s.set_coeff(i, coeff(s, i) * recip(coeff(a, degree(a)) * coeff(c, degree(c))));
873 for ( int i=0; i<=degree(t); ++i ) {
874 t.set_coeff(i, coeff(t, i) * recip(coeff(b, degree(b)) * coeff(c, degree(c))));
879 static ex replace_lc(const ex& poly, const ex& x, const ex& lc)
881 ex r = expand(poly + (lc - poly.lcoeff(x)) * pow(x, poly.degree(x)));
885 static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const umod& u1_, const umod& w1_, const ex& gamma_ = 0)
888 const cl_univpoly_modint_ring& UPR = u1_.ring();
889 const cl_modint_ring& R = UPR->basering();
893 for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
894 maxcoeff += pow(abs(a.coeff(x, i)),2);
896 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
897 cl_I maxdegree = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
898 cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
901 ex alpha = a.lcoeff(x);
906 numeric gamma_ui = ex_to<numeric>(abs(gamma));
908 umod nu1 = COPY(nu1, u1_);
910 umod nw1 = COPY(nw1, w1_);
913 phi = gamma * umod_to_ex(nu1, x);
914 umod u1 = umod_from_ex(phi, x, R);
915 phi = alpha * umod_to_ex(nw1, x);
916 umod w1 = umod_from_ex(phi, x, R);
919 umod g = UPR->create(-1);
920 umod s = UPR->create(-1);
921 umod t = UPR->create(-1);
922 exteuclid(u1, w1, g, s, t);
925 ex u = replace_lc(umod_to_ex(u1, x), x, gamma);
926 ex w = replace_lc(umod_to_ex(w1, x), x, alpha);
927 ex e = expand(a - u * w);
929 const numeric maxmodulus = 2*numeric(B)*gamma_ui;
932 while ( !e.is_zero() && modulus < maxmodulus ) {
934 phi = expand(umod_to_ex(s, x) * c);
935 umod sigmatilde = umod_from_ex(phi, x, R);
936 phi = expand(umod_to_ex(t, x) * c);
937 umod tautilde = umod_from_ex(phi, x, R);
938 umod q = UPR->create(-1);
939 umod r = remdiv(sigmatilde, w1, q);
940 umod sigma = COPY(sigma, r);
941 phi = expand(umod_to_ex(tautilde, x) + umod_to_ex(q, x) * umod_to_ex(u1, x));
942 umod tau = umod_from_ex(phi, x, R);
943 u = expand(u + umod_to_ex(tau, x) * modulus);
944 w = expand(w + umod_to_ex(sigma, x) * modulus);
945 e = expand(a - u * w);
946 modulus = modulus * p;
951 ex delta = u.content(x);
953 w = w / gamma * delta;
961 static unsigned int next_prime(unsigned int p)
963 static vector<unsigned int> primes;
964 if ( primes.size() == 0 ) {
965 primes.push_back(3); primes.push_back(5); primes.push_back(7);
967 vector<unsigned int>::const_iterator it = primes.begin();
968 if ( p >= primes.back() ) {
969 unsigned int candidate = primes.back() + 2;
971 size_t n = primes.size()/2;
972 for ( size_t i=0; i<n; ++i ) {
973 if ( candidate % primes[i] ) continue;
977 primes.push_back(candidate);
978 if ( candidate > p ) break;
982 vector<unsigned int>::const_iterator end = primes.end();
983 for ( ; it!=end; ++it ) {
988 throw logic_error("next_prime: should not reach this point!");
994 Partition(size_t n_) : n(n_)
1000 int operator[](size_t i) const { return k[i]; }
1001 size_t size() const { return n; }
1002 size_t size_first() const { return n-sum; }
1003 size_t size_second() const { return sum; }
1007 for ( size_t i=0; i<k.size(); ++i ) {
1008 cout << k[i] << " ";
1015 for ( size_t i=n-1; i>=1; --i ) {
1031 static void split(const umodvec& factors, const Partition& part, umod& a, umod& b)
1033 a = factors.front().ring()->one();
1034 b = factors.front().ring()->one();
1035 for ( size_t i=0; i<part.size(); ++i ) {
1051 static ex factor_univariate(const ex& poly, const ex& x)
1053 DCOUT(factor_univariate);
1056 ex unit, cont, prim;
1057 poly.unitcontprim(x, unit, cont, prim);
1059 // determine proper prime and minimize number of modular factors
1060 unsigned int p = 3, lastp = 3;
1062 unsigned int trials = 0;
1063 unsigned int minfactors = 0;
1064 numeric lcoeff = ex_to<numeric>(prim.lcoeff(x));
1066 while ( trials < 2 ) {
1069 if ( irem(lcoeff, p) != 0 ) {
1070 R = find_modint_ring(p);
1071 umod modpoly = umod_from_ex(prim, x, R);
1072 if ( squarefree(modpoly) ) break;
1076 // do modular factorization
1077 umod modpoly = umod_from_ex(prim, x, R);
1078 umodvec trialfactors;
1079 factor_modular(modpoly, trialfactors);
1080 if ( trialfactors.size() <= 1 ) {
1081 // irreducible for sure
1085 if ( minfactors == 0 || trialfactors.size() < minfactors ) {
1086 factors = trialfactors;
1087 minfactors = factors.size();
1096 R = find_modint_ring(p);
1097 cl_univpoly_modint_ring UPR = find_univpoly_ring(R);
1099 // lift all factor combinations
1100 stack<ModFactors> tocheck;
1103 mf.factors = factors;
1106 while ( tocheck.size() ) {
1107 const size_t n = tocheck.top().factors.size();
1110 umod a = UPR->create(-1);
1111 umod b = UPR->create(-1);
1112 split(tocheck.top().factors, part, a, b);
1114 ex answer = hensel_univar(tocheck.top().poly, x, p, a, b);
1115 if ( answer != lst() ) {
1116 if ( part.size_first() == 1 ) {
1117 if ( part.size_second() == 1 ) {
1118 result *= answer.op(0) * answer.op(1);
1122 result *= answer.op(0);
1123 tocheck.top().poly = answer.op(1);
1124 for ( size_t i=0; i<n; ++i ) {
1125 if ( part[i] == 0 ) {
1126 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1132 else if ( part.size_second() == 1 ) {
1133 if ( part.size_first() == 1 ) {
1134 result *= answer.op(0) * answer.op(1);
1138 result *= answer.op(1);
1139 tocheck.top().poly = answer.op(0);
1140 for ( size_t i=0; i<n; ++i ) {
1141 if ( part[i] == 1 ) {
1142 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1149 umodvec newfactors1(part.size_first(), UPR->create(-1)), newfactors2(part.size_second(), UPR->create(-1));
1150 umodvec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
1151 for ( size_t i=0; i<n; ++i ) {
1153 *i2++ = tocheck.top().factors[i];
1156 *i1++ = tocheck.top().factors[i];
1159 tocheck.top().factors = newfactors1;
1160 tocheck.top().poly = answer.op(0);
1162 mf.factors = newfactors2;
1163 mf.poly = answer.op(1);
1169 if ( !part.next() ) {
1170 result *= tocheck.top().poly;
1178 DCOUT(END factor_univariate);
1179 return unit * cont * result;
1190 // forward declaration
1191 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);
1193 umodvec multiterm_eea_lift(const umodvec& a, const ex& x, unsigned int p, unsigned int k)
1195 DCOUT(multiterm_eea_lift);
1200 const size_t r = a.size();
1202 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1203 cl_univpoly_modint_ring UPR = find_univpoly_ring(R);
1204 umodvec q(r-1, UPR->create(-1));
1206 for ( size_t j=r-2; j>=1; --j ) {
1207 q[j-1] = a[j] * q[j];
1210 umod beta = UPR->one();
1212 for ( size_t j=1; j<r; ++j ) {
1215 vector<ex> mdarg(2);
1216 mdarg[0] = umod_to_ex(q[j-1], x);
1217 mdarg[1] = umod_to_ex(a[j-1], x);
1218 vector<EvalPoint> empty;
1219 vector<ex> exsigma = multivar_diophant(mdarg, x, umod_to_ex(beta, x), empty, 0, p, k);
1220 umod sigma1 = umod_from_ex(exsigma[0], x, R);
1221 umod sigma2 = umod_from_ex(exsigma[1], x, R);
1222 beta = COPY(beta, sigma1);
1223 s.push_back(sigma2);
1228 DCOUT(END multiterm_eea_lift);
1232 void change_modulus(umod& out, const umod& in)
1234 // ASSERT: out and in have same degree
1235 if ( out.ring() == in.ring() ) {
1236 out = COPY(out, in);
1239 for ( int i=0; i<=degree(in); ++i ) {
1240 out.set_coeff(i, out.ring()->basering()->canonhom(in.ring()->basering()->retract(coeff(in, i))));
1246 void eea_lift(const umod& a, const umod& b, const ex& x, unsigned int p, unsigned int k, umod& s_, umod& t_)
1250 cl_modint_ring R = find_modint_ring(p);
1251 cl_univpoly_modint_ring UPR = find_univpoly_ring(R);
1252 umod amod = UPR->create(degree(a));
1253 change_modulus(amod, a);
1254 umod bmod = UPR->create(degree(b));
1255 change_modulus(bmod, b);
1257 umod g = UPR->create(-1);
1258 umod smod = UPR->create(-1);
1259 umod tmod = UPR->create(-1);
1260 exteuclid(amod, bmod, g, smod, tmod);
1262 cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
1263 cl_univpoly_modint_ring UPRpk = find_univpoly_ring(Rpk);
1264 umod s = UPRpk->create(degree(smod));
1265 change_modulus(s, smod);
1266 umod t = UPRpk->create(degree(tmod));
1267 change_modulus(t, tmod);
1270 umod one = UPRpk->one();
1271 for ( size_t j=1; j<k; ++j ) {
1272 umod e = one - a * s - b * t;
1273 e = divide(e, modulus);
1274 umod c = UPR->create(degree(e));
1275 change_modulus(c, e);
1276 umod sigmabar = smod * c;
1277 umod taubar = tmod * c;
1278 umod q = UPR->create(-1);
1279 umod sigma = remdiv(sigmabar, bmod, q);
1280 umod tau = taubar + q * amod;
1281 umod sadd = UPRpk->create(degree(sigma));
1282 change_modulus(sadd, sigma);
1283 cl_MI modmodulus(Rpk, modulus);
1284 s = s + sadd * modmodulus;
1285 umod tadd = UPRpk->create(degree(tau));
1286 change_modulus(tadd, tau);
1287 t = t + tadd * modmodulus;
1288 modulus = modulus * p;
1293 DCOUT2(check, a*s + b*t);
1294 DCOUT(END eea_lift);
1297 umodvec univar_diophant(const umodvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
1299 DCOUT(univar_diophant);
1306 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1307 cl_univpoly_modint_ring UPR = find_univpoly_ring(R);
1309 const size_t r = a.size();
1312 umodvec s = multiterm_eea_lift(a, x, p, k);
1313 for ( size_t j=0; j<r; ++j ) {
1314 ex phi = expand(pow(x,m) * umod_to_ex(s[j], x));
1315 umod bmod = umod_from_ex(phi, x, R);
1316 umod buf = rem(bmod, a[j]);
1317 result.push_back(buf);
1321 umod s = UPR->create(-1);
1322 umod t = UPR->create(-1);
1323 eea_lift(a[1], a[0], x, p, k, s, t);
1324 ex phi = expand(pow(x,m) * umod_to_ex(s, x));
1325 umod bmod = umod_from_ex(phi, x, R);
1326 umod q = UPR->create(-1);
1327 umod buf = remdiv(bmod, a[0], q);
1328 result.push_back(buf);
1329 phi = expand(pow(x,m) * umod_to_ex(t, x));
1330 umod t1mod = umod_from_ex(phi, x, R);
1331 umod buf2 = t1mod + q * a[1];
1332 result.push_back(buf2);
1336 DCOUT(END univar_diophant);
1340 struct make_modular_map : public map_function {
1342 make_modular_map(const cl_modint_ring& R_) : R(R_) { }
1343 ex operator()(const ex& e)
1345 if ( is_a<add>(e) || is_a<mul>(e) ) {
1346 return e.map(*this);
1348 else if ( is_a<numeric>(e) ) {
1349 numeric mod(R->modulus);
1350 numeric halfmod = (mod-1)/2;
1351 cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1352 numeric n(R->retract(emod));
1353 if ( n > halfmod ) {
1364 static ex make_modular(const ex& e, const cl_modint_ring& R)
1366 make_modular_map map(R);
1367 return map(e.expand());
1370 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)
1374 DCOUT(multivar_diophant);
1377 for ( size_t i=0; i<a.size(); ++i ) {
1378 cout << a[i] << " ";
1386 for ( size_t i=0; i<I.size(); ++i ) {
1387 cout << I[i].x << "=" << I[i].evalpoint << " ";
1395 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
1396 const size_t r = a.size();
1397 const size_t nu = I.size() + 1;
1403 ex xnu = I.back().x;
1404 int alphanu = I.back().evalpoint;
1407 for ( size_t i=0; i<r; ++i ) {
1411 for ( size_t i=0; i<r; ++i ) {
1412 b[i] = normal(A / a[i]);
1415 vector<ex> anew = a;
1416 for ( size_t i=0; i<r; ++i ) {
1417 anew[i] = anew[i].subs(xnu == alphanu);
1419 ex cnew = c.subs(xnu == alphanu);
1420 vector<EvalPoint> Inew = I;
1422 sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
1426 for ( size_t i=0; i<r; ++i ) {
1427 buf -= sigma[i] * b[i];
1429 ex e = make_modular(buf, R);
1434 for ( size_t m=1; m<=d; ++m ) {
1436 while ( !e.is_zero() && e.has(xnu) ) {
1437 monomial *= (xnu - alphanu);
1438 monomial = expand(monomial);
1442 ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
1443 cm = make_modular(cm, R);
1445 if ( !cm.is_zero() ) {
1446 vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
1449 for ( size_t j=0; j<delta_s.size(); ++j ) {
1450 delta_s[j] *= monomial;
1451 sigma[j] += delta_s[j];
1452 buf -= delta_s[j] * b[j];
1454 e = make_modular(buf, R);
1461 DCOUT(uniterm left);
1463 for ( size_t i=0; i<a.size(); ++i ) {
1464 umod up = umod_from_ex(a[i], x, R);
1468 sigma.insert(sigma.begin(), r, 0);
1471 if ( is_a<add>(c) ) {
1480 for ( size_t i=0; i<nterms; ++i ) {
1482 int m = z.degree(x);
1484 cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
1486 umodvec delta_s = univar_diophant(amod, x, m, p, k);
1489 while ( poscm < 0 ) {
1490 poscm = poscm + expt_pos(cl_I(p),k);
1492 modcm = cl_MI(R, poscm);
1494 for ( size_t j=0; j<delta_s.size(); ++j ) {
1495 delta_s[j] = delta_s[j] * modcm;
1496 sigma[j] = sigma[j] + umod_to_ex(delta_s[j], x);
1500 cout << "STEP " << i << " sigma ";
1501 for ( size_t p=0; p<sigma.size(); ++p ) {
1502 cout << sigma[p] << " ";
1513 for ( size_t i=0; i<sigma.size(); ++i ) {
1514 cout << sigma[i] << " ";
1519 for ( size_t i=0; i<sigma.size(); ++i ) {
1520 sigma[i] = make_modular(sigma[i], R);
1525 for ( size_t i=0; i<sigma.size(); ++i ) {
1526 cout << sigma[i] << " ";
1530 DCOUT(END multivar_diophant);
1535 ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
1537 for ( size_t i=0; i<v.size(); ++i ) {
1538 o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
1542 #endif // def DEBUGFACTOR
1545 ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I, unsigned int p, const cl_I& l, const umodvec& u, const vector<ex>& lcU)
1547 DCOUT(hensel_multivar);
1555 const size_t nu = I.size() + 1;
1556 const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
1563 for ( size_t j=nu; j>=2; --j ) {
1565 int alpha = I[j-2].evalpoint;
1566 A[j-2] = A[j-1].subs(x==alpha);
1567 A[j-2] = make_modular(A[j-2], R);
1572 for ( size_t i=0; i<A.size(); ++i) cout << A[i] << " ";
1576 int maxdeg = a.degree(I.front().x);
1577 for ( size_t i=1; i<I.size(); ++i ) {
1578 int maxdeg2 = a.degree(I[i].x);
1579 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
1583 const size_t n = u.size();
1586 for ( size_t i=0; i<n; ++i ) {
1587 U[i] = umod_to_ex(u[i], x);
1591 for ( size_t i=0; i<U.size(); ++i) cout << U[i] << " ";
1595 for ( size_t j=2; j<=nu; ++j ) {
1600 for ( size_t m=0; m<n; ++m) {
1601 if ( lcU[m] != 1 ) {
1603 for ( size_t i=j-1; i<nu-1; ++i ) {
1604 coef = coef.subs(I[i].x == I[i].evalpoint);
1606 coef = make_modular(coef, R);
1607 int deg = U[m].degree(x);
1608 U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
1613 for ( size_t i=0; i<n; ++i ) {
1616 ex e = expand(A[j-1] - Uprod);
1619 vector<EvalPoint> newI;
1620 for ( size_t i=1; i<=j-2; ++i ) {
1621 newI.push_back(I[i-1]);
1626 int alphaj = I[j-2].evalpoint;
1627 size_t deg = A[j-1].degree(xj);
1629 for ( size_t k=1; k<=deg; ++k ) {
1631 if ( !e.is_zero() ) {
1634 monomial *= (xj - alphaj);
1635 monomial = expand(monomial);
1637 ex dif = e.diff(ex_to<symbol>(xj), k);
1639 ex c = dif.subs(xj==alphaj) / factorial(k);
1641 if ( !c.is_zero() ) {
1642 vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
1643 for ( size_t i=0; i<n; ++i ) {
1645 DCOUTVAR(deltaU[i]);
1646 deltaU[i] *= monomial;
1648 U[i] = make_modular(U[i], R);
1652 for ( size_t i=0; i<n; ++i ) {
1655 DCOUTVAR(Uprod.expand());
1658 e = make_modular(e, R);
1666 for ( size_t i=0; i<U.size(); ++i ) {
1670 if ( expand(a-acand).is_zero() ) {
1672 for ( size_t i=0; i<U.size(); ++i ) {
1676 DCOUT(END hensel_multivar);
1682 DCOUT(END hensel_multivar);
1687 static ex put_factors_into_lst(const ex& e)
1689 DCOUT(put_factors_into_lst);
1694 if ( is_a<numeric>(e) ) {
1696 DCOUT(END put_factors_into_lst);
1700 if ( is_a<power>(e) ) {
1702 result.append(e.op(0));
1703 result.append(e.op(1));
1704 DCOUT(END put_factors_into_lst);
1708 if ( is_a<symbol>(e) || is_a<add>(e) ) {
1712 DCOUT(END put_factors_into_lst);
1716 if ( is_a<mul>(e) ) {
1718 for ( size_t i=0; i<e.nops(); ++i ) {
1720 if ( is_a<numeric>(op) ) {
1723 if ( is_a<power>(op) ) {
1724 result.append(op.op(0));
1725 result.append(op.op(1));
1727 if ( is_a<symbol>(op) || is_a<add>(op) ) {
1732 result.prepend(nfac);
1733 DCOUT(END put_factors_into_lst);
1737 throw runtime_error("put_factors_into_lst: bad term.");
1741 ostream& operator<<(ostream& o, const vector<numeric>& v)
1743 for ( size_t i=0; i<v.size(); ++i ) {
1748 #endif // def DEBUGFACTOR
1750 static bool checkdivisors(const lst& f, vector<numeric>& d)
1752 DCOUT(checkdivisors);
1753 const int k = f.nops()-2;
1757 d[0] = ex_to<numeric>(f.op(0) * f.op(f.nops()-1));
1758 if ( d[0] == 1 && k == 1 && abs(f.op(1)) != 1 ) {
1760 DCOUT(END checkdivisors);
1764 for ( int i=1; i<=k; ++i ) {
1766 DCOUTVAR(abs(f.op(i)));
1767 q = ex_to<numeric>(abs(f.op(i)));
1769 for ( int j=i-1; j>=0; --j ) {
1780 DCOUT(END checkdivisors);
1787 DCOUT(END checkdivisors);
1791 static bool generate_set(const ex& u, const ex& vn, const exset& syms, const ex& f, const numeric& modulus, vector<numeric>& a, vector<numeric>& d)
1793 // computation of d is actually not necessary
1794 DCOUT(generate_set);
1799 const ex& x = *syms.begin();
1805 exset::const_iterator s = syms.begin();
1807 for ( size_t i=0; i<a.size(); ++i ) {
1810 a[i] = mod(numeric(rand()), 2*modulus) - modulus;
1811 vnatry = vna.subs(*s == a[i]);
1812 } while ( vnatry == 0 );
1814 u0 = u0.subs(*s == a[i]);
1819 if ( gcd(u0,u0.diff(ex_to<symbol>(x))) != 1 ) {
1822 if ( is_a<numeric>(vn) ) {
1826 DCOUT(do substitution);
1828 lst::const_iterator i = ex_to<lst>(f).begin();
1830 bool problem = false;
1831 while ( i!=ex_to<lst>(f).end() ) {
1833 if ( !is_a<numeric>(fs) ) {
1836 for ( size_t j=0; j<a.size(); ++j ) {
1837 fs = fs.subs(*s == a[j]);
1840 if ( abs(fs) == 1 ) {
1851 ex con = u0.content(x);
1854 trying = checkdivisors(fnum, d);
1857 DCOUT(END generate_set);
1861 static ex factor_multivariate(const ex& poly, const exset& syms)
1863 DCOUT(factor_multivariate);
1866 exset::const_iterator s;
1867 const ex& x = *syms.begin();
1870 /* make polynomial primitive */
1871 ex p = poly.expand().collect(x);
1873 ex cont = p.lcoeff(x);
1874 for ( numeric i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
1875 cont = gcd(cont, p.coeff(x,ex_to<numeric>(i).to_int()));
1876 if ( cont == 1 ) break;
1879 ex pp = expand(normal(p / cont));
1881 if ( !is_a<numeric>(cont) ) {
1883 return ::factor(cont) * ::factor(pp);
1885 return factor(cont) * factor(pp);
1889 /* factor leading coefficient */
1891 ex vn = pp.lcoeff(x);
1894 if ( is_a<numeric>(vn) ) {
1899 ex vnfactors = ::factor(vn);
1901 ex vnfactors = factor(vn);
1903 vnlst = put_factors_into_lst(vnfactors);
1907 const numeric maxtrials = 3;
1908 numeric modulus = (vnlst.nops()-1 > 3) ? vnlst.nops()-1 : 3;
1910 numeric minimalr = -1;
1911 vector<numeric> a(syms.size()-1, 0);
1912 vector<numeric> d((vnlst.nops()-1)/2+1, 0);
1915 numeric trialcount = 0;
1917 unsigned int prime = 3;
1918 size_t factor_count = 0;
1921 while ( trialcount < maxtrials ) {
1922 bool problem = generate_set(pp, vn, syms, vnlst, modulus, a, d);
1933 for ( size_t i=0; i<a.size(); ++i ) {
1934 u = u.subs(*s == a[i]);
1937 delta = u.content(x);
1940 // determine proper prime
1943 cl_modint_ring R = find_modint_ring(prime);
1944 DCOUTVAR(u.lcoeff(x));
1946 if ( irem(ex_to<numeric>(u.lcoeff(x)), prime) != 0 ) {
1947 umod modpoly = umod_from_ex(u, x, R);
1948 if ( squarefree(modpoly) ) break;
1950 prime = next_prime(prime);
1952 R = find_modint_ring(prime);
1961 ufaclst = put_factors_into_lst(ufac);
1963 factor_count = (ufaclst.nops()-1)/2;
1964 DCOUTVAR(factor_count);
1966 if ( factor_count <= 1 ) {
1968 DCOUT(END factor_multivariate);
1972 if ( minimalr < 0 ) {
1973 minimalr = factor_count;
1975 else if ( minimalr == factor_count ) {
1979 else if ( minimalr > factor_count ) {
1980 minimalr = factor_count;
1983 DCOUTVAR(trialcount);
1985 if ( minimalr <= 1 ) {
1987 DCOUT(END factor_multivariate);
1992 vector<numeric> ftilde((vnlst.nops()-1)/2+1);
1993 ftilde[0] = ex_to<numeric>(vnlst.op(0));
1994 for ( size_t i=1; i<ftilde.size(); ++i ) {
1995 ex ft = vnlst.op((i-1)*2+1);
1998 for ( size_t j=0; j<a.size(); ++j ) {
1999 ft = ft.subs(*s == a[j]);
2002 ftilde[i] = ex_to<numeric>(ft);
2006 vector<bool> used_flag((vnlst.nops()-1)/2+1, false);
2007 vector<ex> D(factor_count, 1);
2008 for ( size_t i=0; i<=factor_count; ++i ) {
2012 prefac = ex_to<numeric>(ufaclst.op(0));
2013 ftilde[0] = ftilde[0] / prefac;
2014 vnlst.let_op(0) = vnlst.op(0) / prefac;
2018 prefac = ex_to<numeric>(ufaclst.op(2*(i-1)+1).lcoeff(x));
2021 for ( size_t j=(vnlst.nops()-1)/2+1; j>0; --j ) {
2024 DCOUTVAR(ftilde[j-1]);
2025 if ( abs(ftilde[j-1]) == 1 ) {
2026 used_flag[j-1] = true;
2029 numeric g = gcd(prefac, ftilde[j-1]);
2032 DCOUT(has_common_prime);
2033 prefac = prefac / g;
2034 numeric count = abs(iquo(g, ftilde[j-1]));
2036 used_flag[j-1] = true;
2039 D[i-1] = D[i-1] * pow(vnlst.op(0), count);
2042 D[i-1] = D[i-1] * pow(vnlst.op(2*(j-2)+1), count);
2046 ftilde[j-1] = ftilde[j-1] / prefac;
2048 DCOUTVAR(ftilde[j-1]);
2057 bool some_factor_unused = false;
2058 for ( size_t i=0; i<used_flag.size(); ++i ) {
2059 if ( !used_flag[i] ) {
2060 some_factor_unused = true;
2064 if ( some_factor_unused ) {
2065 DCOUT(some factor unused!);
2069 vector<ex> C(factor_count);
2073 for ( size_t i=0; i<D.size(); ++i ) {
2077 for ( size_t j=0; j<a.size(); ++j ) {
2078 Dtilde = Dtilde.subs(*s == a[j]);
2082 C[i] = D[i] * (ufaclst.op(2*i+1).lcoeff(x) / Dtilde);
2086 for ( size_t i=0; i<D.size(); ++i ) {
2090 for ( size_t j=0; j<a.size(); ++j ) {
2091 Dtilde = Dtilde.subs(*s == a[j]);
2099 ui = ufaclst.op(2*(i-1)+1);
2102 ex d = gcd(ui.lcoeff(x), Dtilde);
2103 C[i] = D[i] * ( ui.lcoeff(x) / d );
2104 ui = ui * ( Dtilde[i] / d );
2105 delta = delta / ( Dtilde[i] / d );
2106 if ( delta == 1 ) break;
2108 C[i] = delta * C[i];
2109 pp = pp * pow(delta, D.size()-1);
2116 vector<EvalPoint> epv;
2119 for ( size_t i=0; i<a.size(); ++i ) {
2121 ep.evalpoint = a[i].to_int();
2128 for ( int i=u.degree(x); i>=u.ldegree(x); --i ) {
2129 maxcoeff += pow(abs(u.coeff(x, i)),2);
2131 cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
2132 unsigned int maxdegree = 0;
2133 for ( size_t i=0; i<factor_count; ++i ) {
2134 if ( ufaclst[2*i+1].degree(x) > (int)maxdegree ) {
2135 maxdegree = ufaclst[2*i+1].degree(x);
2138 cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
2147 cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
2148 for ( size_t i=0; i<(ufaclst.nops()-1)/2; ++i ) {
2149 umod newu = umod_from_ex(ufaclst.op(i*2+1), x, R);
2150 uvec.push_back(newu);
2154 ex res = hensel_multivar(ufaclst.op(0)*pp, x, epv, prime, l, uvec, C);
2155 if ( res != lst() ) {
2156 ex result = cont * ufaclst.op(0);
2157 for ( size_t i=0; i<res.nops(); ++i ) {
2158 result *= res.op(i).content(x) * res.op(i).unit(x);
2159 result *= res.op(i).primpart(x);
2162 DCOUT(END factor_multivariate);
2168 struct find_symbols_map : public map_function {
2170 ex operator()(const ex& e)
2172 if ( is_a<symbol>(e) ) {
2176 return e.map(*this);
2180 static ex factor_sqrfree(const ex& poly)
2182 DCOUT(factor_sqrfree);
2184 // determine all symbols in poly
2185 find_symbols_map findsymbols;
2187 if ( findsymbols.syms.size() == 0 ) {
2188 DCOUT(END factor_sqrfree);
2192 if ( findsymbols.syms.size() == 1 ) {
2194 const ex& x = *(findsymbols.syms.begin());
2195 if ( poly.ldegree(x) > 0 ) {
2196 int ld = poly.ldegree(x);
2197 ex res = factor_univariate(expand(poly/pow(x, ld)), x);
2198 DCOUT(END factor_sqrfree);
2199 return res * pow(x,ld);
2202 ex res = factor_univariate(poly, x);
2203 DCOUT(END factor_sqrfree);
2208 // multivariate case
2209 ex res = factor_multivariate(poly, findsymbols.syms);
2210 DCOUT(END factor_sqrfree);
2214 struct apply_factor_map : public map_function {
2216 apply_factor_map(unsigned options_) : options(options_) { }
2217 ex operator()(const ex& e)
2219 if ( e.info(info_flags::polynomial) ) {
2221 return ::factor(e, options);
2223 return factor(e, options);
2226 if ( is_a<add>(e) ) {
2228 for ( size_t i=0; i<e.nops(); ++i ) {
2229 if ( e.op(i).info(info_flags::polynomial) ) {
2239 return ::factor(s1, options) + s2.map(*this);
2241 return factor(s1, options) + s2.map(*this);
2244 return e.map(*this);
2248 } // anonymous namespace
2251 ex factor(const ex& poly, unsigned options = 0)
2253 ex factor(const ex& poly, unsigned options)
2257 if ( !poly.info(info_flags::polynomial) ) {
2258 if ( options & factor_options::all ) {
2259 options &= ~factor_options::all;
2260 apply_factor_map factor_map(options);
2261 return factor_map(poly);
2266 // determine all symbols in poly
2267 find_symbols_map findsymbols;
2269 if ( findsymbols.syms.size() == 0 ) {
2273 exset::const_iterator i=findsymbols.syms.begin(), end=findsymbols.syms.end();
2274 for ( ; i!=end; ++i ) {
2278 // make poly square free
2279 ex sfpoly = sqrfree(poly, syms);
2281 // factorize the square free components
2282 if ( is_a<power>(sfpoly) ) {
2283 // case: (polynomial)^exponent
2284 const ex& base = sfpoly.op(0);
2285 if ( !is_a<add>(base) ) {
2286 // simple case: (monomial)^exponent
2289 ex f = factor_sqrfree(base);
2290 return pow(f, sfpoly.op(1));
2292 if ( is_a<mul>(sfpoly) ) {
2293 // case: multiple factors
2295 for ( size_t i=0; i<sfpoly.nops(); ++i ) {
2296 const ex& t = sfpoly.op(i);
2297 if ( is_a<power>(t) ) {
2298 const ex& base = t.op(0);
2299 if ( !is_a<add>(base) ) {
2303 ex f = factor_sqrfree(base);
2304 res *= pow(f, t.op(1));
2307 else if ( is_a<add>(t) ) {
2308 ex f = factor_sqrfree(t);
2317 if ( is_a<symbol>(sfpoly) ) {
2320 // case: (polynomial)
2321 ex f = factor_sqrfree(sfpoly);
2325 } // namespace GiNaC