84 #define DCOUT(str) cout << #str << endl
85 #define DCOUTVAR(var) cout << #var << ": " << var << endl
86 #define DCOUT2(str,var) cout << #str << ": " << var << endl
87 ostream&
operator<<(ostream& o,
const vector<int>& v)
89 auto i = v.begin(), end = v.end();
96 static ostream&
operator<<(ostream& o,
const vector<cl_I>& v)
98 auto i = v.begin(), end = v.end();
100 o << *i <<
"[" << i-v.begin() <<
"]" <<
" ";
105 static ostream&
operator<<(ostream& o,
const vector<cl_MI>& v)
107 auto i = v.begin(), end = v.end();
109 o << *i <<
"[" << i-v.begin() <<
"]" <<
" ";
114 ostream&
operator<<(ostream& o,
const vector<numeric>& v)
116 for (
size_t i=0; i<v.size(); ++i ) {
121 ostream&
operator<<(ostream& o,
const vector<vector<cl_MI>>& v)
123 auto i = v.begin(), end = v.end();
125 o << i-v.begin() <<
": " << *i << endl;
132 #define DCOUTVAR(var)
133 #define DCOUT2(str,var)
142 typedef std::vector<cln::cl_MI> umodpoly;
143 typedef std::vector<cln::cl_I> upoly;
144 typedef vector<umodpoly> upvec;
149 template<
typename T>
static int degree(
const T& p)
154 template<
typename T>
static typename T::value_type lcoeff(
const T& p)
156 return p[p.size() - 1];
159 static bool normalize_in_field(umodpoly& a)
163 if ( lcoeff(a) == a[0].ring()->
one() ) {
167 const cln::cl_MI lc_1 = recip(lcoeff(a));
168 for (std::size_t
k = a.size();
k-- != 0; )
173 template<
typename T>
static void
174 canonicalize(T& p,
const typename T::size_type hint = std::numeric_limits<typename T::size_type>::max())
179 std::size_t i = p.size() - 1;
205 p.erase(p.begin() + i, p.end());
210 static void expt_pos(umodpoly& a,
unsigned int q)
212 if ( a.empty() )
return;
213 cl_MI zero = a[0].ring()->zero();
215 a.resize(
degree(a)*q+1, zero);
216 for (
int i=deg; i>0; --i ) {
222 template<
bool COND,
typename T =
void>
struct enable_if
227 template<
typename T>
struct enable_if<false, T> { };
229 template<
typename T>
struct uvar_poly_p
234 template<>
struct uvar_poly_p<upoly>
236 static const bool value =
true;
239 template<>
struct uvar_poly_p<umodpoly>
241 static const bool value =
true;
254 for ( ; i<sb; ++i ) {
257 for ( ; i<sa; ++i ) {
266 for ( ; i<sa; ++i ) {
269 for ( ; i<sb; ++i ) {
289 for ( ; i<sb; ++i ) {
292 for ( ; i<sa; ++i ) {
301 for ( ; i<sa; ++i ) {
304 for ( ; i<sb; ++i ) {
312 static upoly
operator*(
const upoly& a,
const upoly& b)
315 if ( a.empty() || b.empty() )
return c;
319 for (
int i=0 ; i<=
n; ++i ) {
320 for (
int j=0 ; j<=i; ++j ) {
322 c[i] =
c[i] + a[j] * b[i-j];
329 static umodpoly
operator*(
const umodpoly& a,
const umodpoly& b)
332 if ( a.empty() || b.empty() )
return c;
335 c.resize(
n+1, a[0].ring()->zero());
336 for (
int i=0 ; i<=
n; ++i ) {
337 for (
int j=0 ; j<=i; ++j ) {
339 c[i] =
c[i] + a[j] * b[i-j];
346 static upoly
operator*(
const upoly& a,
const cl_I&
x)
353 for (
size_t i=0; i<a.size(); ++i ) {
359 static upoly
operator/(
const upoly& a,
const cl_I&
x)
366 for (
size_t i=0; i<a.size(); ++i ) {
367 r[i] = exquo(a[i],
x);
372 static umodpoly
operator*(
const umodpoly& a,
const cl_MI&
x)
374 umodpoly
r(a.size());
375 for (
size_t i=0; i<a.size(); ++i ) {
382 static void upoly_from_ex(upoly& up,
const ex& e,
const ex&
x)
385 int deg = e.degree(
x);
387 int ldeg = e.ldegree(
x);
388 for ( ; deg>=ldeg; --deg ) {
389 up[deg] = the<cl_I>(ex_to<numeric>(e.coeff(
x, deg)).to_cl_N());
391 for ( ; deg>=0; --deg ) {
397 static void umodpoly_from_upoly(umodpoly& ump,
const upoly& e,
const cl_modint_ring&
R)
401 for ( ; deg>=0; --deg ) {
402 ump[deg] =
R->canonhom(e[deg]);
407 static void umodpoly_from_ex(umodpoly& ump,
const ex& e,
const ex&
x,
const cl_modint_ring&
R)
410 int deg = e.degree(
x);
412 int ldeg = e.ldegree(
x);
413 for ( ; deg>=ldeg; --deg ) {
414 cl_I
coeff = the<cl_I>(ex_to<numeric>(e.coeff(
x, deg)).to_cl_N());
415 ump[deg] =
R->canonhom(
coeff);
417 for ( ; deg>=0; --deg ) {
418 ump[deg] =
R->zero();
424 static void umodpoly_from_ex(umodpoly& ump,
const ex& e,
const ex&
x,
const cl_I& modulus)
426 umodpoly_from_ex(ump, e,
x, find_modint_ring(modulus));
430 static ex upoly_to_ex(
const upoly& a,
const ex&
x)
432 if ( a.empty() )
return 0;
434 for (
int i=
degree(a); i>=0; --i ) {
435 e += numeric(a[i]) *
pow(
x, i);
440 static ex umodpoly_to_ex(
const umodpoly& a,
const ex&
x)
442 if ( a.empty() )
return 0;
443 cl_modint_ring
R = a[0].ring();
444 cl_I
mod =
R->modulus;
445 cl_I halfmod = (
mod-1) >> 1;
447 for (
int i=
degree(a); i>=0; --i ) {
448 cl_I
n =
R->retract(a[i]);
452 e += numeric(
n) *
pow(
x, i);
458 static upoly umodpoly_to_upoly(
const umodpoly& a)
461 if ( a.empty() )
return e;
462 cl_modint_ring
R = a[0].ring();
463 cl_I
mod =
R->modulus;
464 cl_I halfmod = (
mod-1) >> 1;
465 for (
int i=
degree(a); i>=0; --i ) {
466 cl_I
n =
R->retract(a[i]);
476 static umodpoly umodpoly_to_umodpoly(
const umodpoly& a,
const cl_modint_ring&
R,
unsigned int m)
479 if ( a.empty() )
return e;
480 cl_modint_ring oldR = a[0].ring();
481 size_t sa = a.size();
482 e.resize(sa+
m,
R->zero());
483 for (
size_t i=0; i<sa; ++i ) {
484 e[i+
m] =
R->canonhom(oldR->retract(a[i]));
497 static void reduce_coeff(umodpoly& a,
const cl_I&
x)
499 if ( a.empty() )
return;
501 cl_modint_ring
R = a[0].ring();
504 cl_I
c =
R->retract(i);
505 i = cl_MI(
R, the<cl_I>(
c /
x));
516 static void rem(
const umodpoly& a,
const umodpoly& b, umodpoly&
r)
525 cl_MI qk = div(
r[
n+
k], b[
n]);
527 for (
int i=0; i<
n; ++i ) {
528 unsigned int j =
n +
k - 1 - i;
529 r[j] =
r[j] - qk * b[j-
k];
534 fill(
r.begin()+
n,
r.end(), a[0].ring()->zero());
545 static void div(
const umodpoly& a,
const umodpoly& b, umodpoly& q)
554 q.resize(
k+1, a[0].ring()->zero());
556 cl_MI qk = div(
r[
n+
k], b[
n]);
559 for (
int i=0; i<
n; ++i ) {
560 unsigned int j =
n +
k - 1 - i;
561 r[j] =
r[j] - qk * b[j-
k];
577 static void remdiv(
const umodpoly& a,
const umodpoly& b, umodpoly&
r, umodpoly& q)
586 q.resize(
k+1, a[0].ring()->zero());
588 cl_MI qk = div(
r[
n+
k], b[
n]);
591 for (
int i=0; i<
n; ++i ) {
592 unsigned int j =
n +
k - 1 - i;
593 r[j] =
r[j] - qk * b[j-
k];
598 fill(
r.begin()+
n,
r.end(), a[0].ring()->zero());
609 static void gcd(
const umodpoly& a,
const umodpoly& b, umodpoly&
c)
614 normalize_in_field(
c);
616 normalize_in_field(d);
618 while ( !d.empty() ) {
623 normalize_in_field(
c);
631 static void deriv(
const umodpoly& a, umodpoly& d)
634 if ( a.size() <= 1 )
return;
636 d.insert(d.begin(), a.begin()+1, a.end());
638 for (
int i=1; i<max; ++i ) {
644 static bool unequal_one(
const umodpoly& a)
646 if ( a.empty() )
return true;
647 return ( a.size() != 1 || a[0] != a[0].ring()->one() );
650 static bool equal_one(
const umodpoly& a)
652 return ( a.size() == 1 && a[0] == a[0].ring()->one() );
660 static bool squarefree(
const umodpoly& a)
678 typedef vector<cl_MI> mvec;
683 friend ostream&
operator<<(ostream& o,
const modular_matrix&
m);
686 modular_matrix(
size_t r_,
size_t c_,
const cl_MI& init) :
r(r_),
c(c_)
690 size_t rowsize()
const {
return r; }
691 size_t colsize()
const {
return c; }
692 cl_MI& operator()(
size_t row,
size_t col) {
return m[row*
c + col]; }
693 cl_MI operator()(
size_t row,
size_t col)
const {
return m[row*
c + col]; }
694 void mul_col(
size_t col,
const cl_MI
x)
696 for (
size_t rc=0; rc<
r; ++rc ) {
697 std::size_t i =
c*rc + col;
701 void sub_col(
size_t col1,
size_t col2,
const cl_MI fac)
703 for (
size_t rc=0; rc<
r; ++rc ) {
704 std::size_t i1 = col1 +
c*rc;
705 std::size_t i2 = col2 +
c*rc;
706 m[i1] =
m[i1] -
m[i2]*fac;
709 void switch_col(
size_t col1,
size_t col2)
711 for (
size_t rc=0; rc<
r; ++rc ) {
712 std::size_t i1 = col1 + rc*
c;
713 std::size_t i2 = col2 + rc*
c;
717 void mul_row(
size_t row,
const cl_MI
x)
719 for (
size_t cc=0; cc<
c; ++cc ) {
720 std::size_t i = row*
c + cc;
724 void sub_row(
size_t row1,
size_t row2,
const cl_MI fac)
726 for (
size_t cc=0; cc<
c; ++cc ) {
727 std::size_t i1 = row1*
c + cc;
728 std::size_t i2 = row2*
c + cc;
729 m[i1] =
m[i1] -
m[i2]*fac;
732 void switch_row(
size_t row1,
size_t row2)
734 for (
size_t cc=0; cc<
c; ++cc ) {
735 std::size_t i1 = row1*
c + cc;
736 std::size_t i2 = row2*
c + cc;
740 bool is_col_zero(
size_t col)
const
742 for (
size_t rr=0; rr<
r; ++rr ) {
743 std::size_t i = col + rr*
c;
744 if ( !zerop(
m[i]) ) {
750 bool is_row_zero(
size_t row)
const
752 for (
size_t cc=0; cc<
c; ++cc ) {
753 std::size_t i = row*
c + cc;
754 if ( !zerop(
m[i]) ) {
760 void set_row(
size_t row,
const vector<cl_MI>& newrow)
762 for (std::size_t i2 = 0; i2 < newrow.size(); ++i2) {
763 std::size_t i1 = row*
c + i2;
767 mvec::const_iterator row_begin(
size_t row)
const {
return m.begin()+row*
c; }
768 mvec::const_iterator row_end(
size_t row)
const {
return m.begin()+row*
c+
r; }
775 modular_matrix
operator*(
const modular_matrix& m1,
const modular_matrix& m2)
777 const unsigned int r = m1.rowsize();
778 const unsigned int c = m2.colsize();
779 modular_matrix o(
r,
c,m1(0,0));
781 for (
size_t i=0; i<
r; ++i ) {
782 for (
size_t j=0; j<
c; ++j ) {
784 buf = m1(i,0) * m2(0,j);
785 for (
size_t k=1;
k<
c; ++
k ) {
786 buf = buf + m1(i,
k)*m2(
k,j);
794 ostream&
operator<<(ostream& o,
const modular_matrix&
m)
796 cl_modint_ring
R =
m(0,0).ring();
798 for (
size_t i=0; i<
m.rowsize(); ++i ) {
800 for (
size_t j=0; j<
m.colsize()-1; ++j ) {
801 o <<
R->retract(
m(i,j)) <<
",";
803 o <<
R->retract(
m(i,
m.colsize()-1)) <<
"}";
804 if ( i !=
m.rowsize()-1 ) {
821 static void q_matrix(
const umodpoly& a_, modular_matrix& Q)
824 normalize_in_field(a);
827 unsigned int q = cl_I_to_uint(a[0].ring()->modulus);
828 umodpoly
r(
n, a[0].ring()->zero());
829 r[0] = a[0].ring()->one();
831 unsigned int max = (
n-1) * q;
832 for (
size_t m=1;
m<=max; ++
m ) {
833 cl_MI rn_1 =
r.back();
834 for (
size_t i=
n-1; i>0; --i ) {
835 r[i] =
r[i-1] - (rn_1 * a[i]);
838 if ( (
m % q) == 0 ) {
849 static void nullspace(modular_matrix& M, vector<mvec>& basis)
851 const size_t n = M.rowsize();
852 const cl_MI
one = M(0,0).ring()->one();
853 for (
size_t i=0; i<
n; ++i ) {
854 M(i,i) = M(i,i) -
one;
856 for (
size_t r=0;
r<
n; ++
r ) {
858 for ( ; cc<
n; ++cc ) {
859 if ( !zerop(M(
r,cc)) ) {
861 if ( !zerop(M(cc,cc)) ) {
873 M.mul_col(
r, recip(M(
r,
r)));
874 for ( cc=0; cc<
n; ++cc ) {
876 M.sub_col(cc,
r, M(
r,cc));
882 for (
size_t i=0; i<
n; ++i ) {
883 M(i,i) = M(i,i) -
one;
885 for (
size_t i=0; i<
n; ++i ) {
886 if ( !M.is_row_zero(i) ) {
887 mvec nu(M.row_begin(i), M.row_end(i));
901 static void berlekamp(
const umodpoly& a, upvec& upv)
903 cl_modint_ring
R = a[0].ring();
904 umodpoly
one(1,
R->one());
912 const unsigned int k = nu.size();
919 unsigned int size = 1;
921 unsigned int q = cl_I_to_uint(
R->modulus);
923 list<umodpoly>::iterator u =
factors.begin();
927 for (
unsigned int s=0; s<q; ++s ) {
928 umodpoly nur = nu[
r];
929 nur[0] = nur[0] - cl_MI(
R, s);
933 if ( unequal_one(g) && g != *u ) {
936 if ( equal_one(uo) ) {
937 throw logic_error(
"berlekamp: unexpected divisor.");
972 static void expt_1_over_p(
const umodpoly& a,
unsigned int prime, umodpoly& ap)
974 size_t newdeg =
degree(a)/prime;
977 for (
size_t i=1; i<=newdeg; ++i ) {
988 static void modsqrfree(
const umodpoly& a, upvec&
factors, vector<int>& mult)
990 const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus);
999 while ( unequal_one(w) ) {
1012 if ( unequal_one(
c) ) {
1014 expt_1_over_p(
c, prime, cp);
1015 size_t previ = mult.size();
1016 modsqrfree(cp,
factors, mult);
1017 for (
size_t i=previ; i<mult.size(); ++i ) {
1023 expt_1_over_p(a, prime, ap);
1024 size_t previ = mult.size();
1025 modsqrfree(ap,
factors, mult);
1026 for (
size_t i=previ; i<mult.size(); ++i ) {
1044 static void distinct_degree_factor(
const umodpoly& a_, vector<int>& degrees, upvec& ddfactors)
1048 cl_modint_ring
R = a[0].ring();
1049 int q = cl_I_to_int(
R->modulus);
1058 while ( i <= nhalf ) {
1063 umodpoly wx = w -
x;
1065 if ( unequal_one(buf) ) {
1066 degrees.push_back(i);
1067 ddfactors.push_back(buf);
1069 if ( unequal_one(buf) ) {
1079 if ( unequal_one(a) ) {
1080 degrees.push_back(
degree(a));
1081 ddfactors.push_back(a);
1095 static void same_degree_factor(
const umodpoly& a, upvec& upv)
1097 cl_modint_ring
R = a[0].ring();
1099 vector<int> degrees;
1101 distinct_degree_factor(a, degrees, ddfactors);
1103 for (
size_t i=0; i<degrees.size(); ++i ) {
1104 if ( degrees[i] ==
degree(ddfactors[i]) ) {
1105 upv.push_back(ddfactors[i]);
1107 berlekamp(ddfactors[i], upv);
1113 #define USE_SAME_DEGREE_FACTOR
1125 static void factor_modular(
const umodpoly& p, upvec& upv)
1127 #ifdef USE_SAME_DEGREE_FACTOR
1128 same_degree_factor(p, upv);
1142 static void exteuclid(
const umodpoly& a,
const umodpoly& b, umodpoly& s, umodpoly& t)
1145 exteuclid(b, a, t, s);
1149 umodpoly
one(1, a[0].ring()->
one());
1150 umodpoly
c = a; normalize_in_field(
c);
1151 umodpoly d = b; normalize_in_field(d);
1159 umodpoly
r =
c - q * d;
1160 umodpoly r1 = s - q * d1;
1161 umodpoly r2 = t - q * d2;
1165 if (
r.empty() )
break;
1170 cl_MI fac = recip(lcoeff(a) * lcoeff(
c));
1171 for (
auto & i : s) {
1175 fac = recip(lcoeff(b) * lcoeff(
c));
1176 for (
auto & i : t) {
1188 static upoly replace_lc(
const upoly&
poly,
const cl_I& lc)
1199 static inline cl_I calc_bound(
const ex& a,
const ex&
x,
int maxdeg)
1203 for (
int i=a.degree(
x); i>=a.ldegree(
x); --i ) {
1204 cl_I aa =
abs(the<cl_I>(ex_to<numeric>(a.coeff(
x, i)).to_cl_N()));
1205 if ( aa > maxcoeff ) maxcoeff = aa;
1209 cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1210 return ( B > maxcoeff ) ? B : maxcoeff;
1216 static inline cl_I calc_bound(
const upoly& a,
int maxdeg)
1220 for (
int i=
degree(a); i>=0; --i ) {
1221 cl_I aa =
abs(a[i]);
1222 if ( aa > maxcoeff ) maxcoeff = aa;
1226 cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
1227 return ( B > maxcoeff ) ? B : maxcoeff;
1242 static void hensel_univar(
const upoly& a_,
unsigned int p,
const umodpoly& u1_,
const umodpoly& w1_, upoly& u, upoly& w)
1245 const cl_modint_ring&
R = u1_[0].ring();
1249 cl_I maxmodulus = 2*calc_bound(a, maxdeg);
1252 cl_I alpha = lcoeff(a);
1255 normalize_in_field(nu1);
1257 normalize_in_field(nw1);
1259 phi = umodpoly_to_upoly(nu1) * alpha;
1261 umodpoly_from_upoly(u1, phi,
R);
1262 phi = umodpoly_to_upoly(nw1) * alpha;
1264 umodpoly_from_upoly(w1, phi,
R);
1269 exteuclid(u1, w1, s, t);
1272 u = replace_lc(umodpoly_to_upoly(u1), alpha);
1273 w = replace_lc(umodpoly_to_upoly(w1), alpha);
1274 upoly e = a - u * w;
1278 while ( !e.empty() && modulus < maxmodulus ) {
1279 upoly
c = e / modulus;
1280 phi = umodpoly_to_upoly(s) *
c;
1281 umodpoly sigmatilde;
1282 umodpoly_from_upoly(sigmatilde, phi,
R);
1283 phi = umodpoly_to_upoly(t) *
c;
1285 umodpoly_from_upoly(tautilde, phi,
R);
1287 remdiv(sigmatilde, w1,
r, q);
1289 phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1);
1291 umodpoly_from_upoly(tau, phi,
R);
1292 u = u + umodpoly_to_upoly(tau) * modulus;
1293 w = w + umodpoly_to_upoly(sigma) * modulus;
1295 modulus = modulus * p;
1301 for (
size_t i=1; i<u.size(); ++i ) {
1303 if ( g == 1 )
break;
1322 static unsigned int next_prime(
unsigned int p)
1324 static vector<unsigned int> primes;
1325 if (primes.empty()) {
1328 if ( p >= primes.back() ) {
1329 unsigned int candidate = primes.back() + 2;
1331 size_t n = primes.size()/2;
1332 for (
size_t i=0; i<
n; ++i ) {
1333 if (candidate % primes[i])
1338 primes.push_back(candidate);
1344 for (
auto & it : primes) {
1349 throw logic_error(
"next_prime: should not reach this point!");
1354 class factor_partition
1358 factor_partition(
const upvec& factors_) :
factors(factors_)
1364 one.resize(1,
factors.front()[0].ring()->one());
1369 int operator[](
size_t i)
const {
return k[i]; }
1370 size_t size()
const {
return n; }
1371 size_t size_left()
const {
return n-
len; }
1372 size_t size_right()
const {
return len; }
1377 if (
last ==
n-1 ) {
1388 if (
last == 0 &&
n == 2*
len )
return false;
1390 for (
size_t i=0; i<=
len-
rem; ++i ) {
1394 fill(
k.begin()+
last,
k.end(), 0);
1401 if (
len >
n/2 )
return false;
1402 fill(
k.begin(),
k.begin()+
len, 1);
1403 fill(
k.begin()+
len+1,
k.end(), 0);
1412 umodpoly& left() {
return lr[0]; }
1414 umodpoly& right() {
return lr[1]; }
1423 while ( i <
n &&
k[i] == group ) { ++d; ++i; }
1425 if (
cache[pos].size() >= d ) {
1426 lr[group] =
lr[group] *
cache[pos][d-1];
1428 if (
cache[pos].size() == 0 ) {
1431 size_t j = pos +
cache[pos].size() + 1;
1432 d -=
cache[pos].size();
1435 cache[pos].push_back(buf);
1439 lr[group] =
lr[group] *
cache[pos].back();
1453 for (
size_t i=0; i<
n; ++i ) {
1488 static ex factor_univariate(
const ex&
poly,
const ex&
x,
unsigned int& prime)
1490 ex unit, cont, prim_ex;
1491 poly.unitcontprim(
x, unit, cont, prim_ex);
1493 upoly_from_ex(prim, prim_ex,
x);
1494 if (prim_ex.is_equal(1)) {
1500 unsigned int lastp = prime;
1502 unsigned int trials = 0;
1503 unsigned int minfactors = 0;
1505 const numeric& cont_n = ex_to<numeric>(cont);
1507 if (cont_n.is_integer()) {
1508 i_cont = the<cl_I>(cont_n.to_cl_N());
1514 cl_I lc = lcoeff(prim)*i_cont;
1516 while ( trials < 2 ) {
1519 prime = next_prime(prime);
1520 if ( !zerop(
rem(lc, prime)) ) {
1521 R = find_modint_ring(prime);
1522 umodpoly_from_upoly(modpoly, prim,
R);
1523 if ( squarefree(modpoly) )
break;
1529 factor_modular(modpoly, trialfactors);
1530 if ( trialfactors.size() <= 1 ) {
1535 if ( minfactors == 0 || trialfactors.size() < minfactors ) {
1537 minfactors = trialfactors.size();
1545 R = find_modint_ring(prime);
1548 stack<ModFactors> tocheck;
1555 while ( tocheck.size() ) {
1556 const size_t n = tocheck.top().factors.size();
1557 factor_partition part(tocheck.top().factors);
1560 hensel_univar(tocheck.top().poly, prime, part.left(), part.right(), f1, f2);
1561 if ( !f1.empty() ) {
1563 if ( part.size_left() == 1 ) {
1564 if ( part.size_right() == 1 ) {
1565 result *= upoly_to_ex(f1,
x) * upoly_to_ex(f2,
x);
1569 result *= upoly_to_ex(f1,
x);
1570 tocheck.top().poly = f2;
1571 for (
size_t i=0; i<
n; ++i ) {
1572 if ( part[i] == 0 ) {
1573 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1579 else if ( part.size_right() == 1 ) {
1580 if ( part.size_left() == 1 ) {
1581 result *= upoly_to_ex(f1,
x) * upoly_to_ex(f2,
x);
1585 result *= upoly_to_ex(f2,
x);
1586 tocheck.top().poly = f1;
1587 for (
size_t i=0; i<
n; ++i ) {
1588 if ( part[i] == 1 ) {
1589 tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
1595 upvec newfactors1(part.size_left()), newfactors2(part.size_right());
1596 auto i1 = newfactors1.begin(), i2 = newfactors2.begin();
1597 for (
size_t i=0; i<
n; ++i ) {
1599 *i2++ = tocheck.top().factors[i];
1601 *i1++ = tocheck.top().factors[i];
1604 tocheck.top().factors = newfactors1;
1605 tocheck.top().poly = f1;
1607 mf.factors = newfactors2;
1614 if ( !part.next() ) {
1617 result *= upoly_to_ex(tocheck.top().poly,
x);
1625 return unit * cont * result;
1631 static inline ex factor_univariate(
const ex&
poly,
const ex&
x)
1634 return factor_univariate(
poly,
x, prime);
1646 ostream&
operator<<(ostream& o,
const vector<EvalPoint>& v)
1648 for (
size_t i=0; i<v.size(); ++i ) {
1649 o <<
"(" << v[i].x <<
"==" << v[i].evalpoint <<
") ";
1656 static 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);
1673 static upvec multiterm_eea_lift(
const upvec& a,
const ex&
x,
unsigned int p,
unsigned int k)
1675 const size_t r = a.size();
1676 cl_modint_ring
R = find_modint_ring(expt_pos(cl_I(p),
k));
1679 for (
size_t j=
r-2; j>=1; --j ) {
1680 q[j-1] = a[j] * q[j];
1682 umodpoly beta(1,
R->one());
1684 for (
size_t j=1; j<
r; ++j ) {
1685 vector<ex> mdarg(2);
1686 mdarg[0] = umodpoly_to_ex(q[j-1],
x);
1687 mdarg[1] = umodpoly_to_ex(a[j-1],
x);
1688 vector<EvalPoint> empty;
1689 vector<ex> exsigma = multivar_diophant(mdarg,
x, umodpoly_to_ex(beta,
x), empty, 0, p,
k);
1691 umodpoly_from_ex(sigma1, exsigma[0],
x,
R);
1693 umodpoly_from_ex(sigma2, exsigma[1],
x,
R);
1695 s.push_back(sigma2);
1706 static void change_modulus(
const cl_modint_ring&
R, umodpoly& a)
1708 if ( a.empty() )
return;
1709 cl_modint_ring oldR = a[0].ring();
1710 for (
auto & i : a) {
1711 i =
R->canonhom(oldR->retract(i));
1730 static void eea_lift(
const umodpoly& a,
const umodpoly& b,
const ex&
x,
unsigned int p,
unsigned int k, umodpoly& s_, umodpoly& t_)
1732 cl_modint_ring
R = find_modint_ring(p);
1734 change_modulus(
R, amod);
1736 change_modulus(
R, bmod);
1740 exteuclid(amod, bmod,
smod, tmod);
1742 cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),
k));
1744 change_modulus(Rpk, s);
1746 change_modulus(Rpk, t);
1749 umodpoly
one(1, Rpk->one());
1750 for (
size_t j=1; j<
k; ++j ) {
1751 umodpoly e =
one - a * s - b * t;
1752 reduce_coeff(e, modulus);
1754 change_modulus(
R,
c);
1755 umodpoly sigmabar =
smod *
c;
1756 umodpoly taubar = tmod *
c;
1758 remdiv(sigmabar, bmod, sigma, q);
1759 umodpoly tau = taubar + q * amod;
1760 umodpoly sadd = sigma;
1761 change_modulus(Rpk, sadd);
1762 cl_MI modmodulus(Rpk, modulus);
1763 s = s + sadd * modmodulus;
1764 umodpoly tadd = tau;
1765 change_modulus(Rpk, tadd);
1766 t = t + tadd * modmodulus;
1767 modulus = modulus * p;
1788 static upvec univar_diophant(
const upvec& a,
const ex&
x,
unsigned int m,
unsigned int p,
unsigned int k)
1790 cl_modint_ring
R = find_modint_ring(expt_pos(cl_I(p),
k));
1792 const size_t r = a.size();
1795 upvec s = multiterm_eea_lift(a,
x, p,
k);
1796 for (
size_t j=0; j<
r; ++j ) {
1797 umodpoly bmod = umodpoly_to_umodpoly(s[j],
R,
m);
1799 rem(bmod, a[j], buf);
1800 result.push_back(buf);
1804 eea_lift(a[1], a[0],
x, p,
k, s, t);
1805 umodpoly bmod = umodpoly_to_umodpoly(s,
R,
m);
1807 remdiv(bmod, a[0], buf, q);
1808 result.push_back(buf);
1809 umodpoly t1mod = umodpoly_to_umodpoly(t,
R,
m);
1810 buf = t1mod + q * a[1];
1811 result.push_back(buf);
1821 struct make_modular_map :
public map_function {
1823 make_modular_map(
const cl_modint_ring& R_) :
R(R_) { }
1824 ex operator()(
const ex& e)
override
1826 if ( is_a<add>(e) || is_a<mul>(e) ) {
1827 return e.map(*
this);
1829 else if ( is_a<numeric>(e) ) {
1830 numeric
mod(
R->modulus);
1831 numeric halfmod = (
mod-1)/2;
1832 cl_MI emod =
R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
1833 numeric
n(
R->retract(emod));
1834 if (
n > halfmod ) {
1851 static ex make_modular(
const ex& e,
const cl_modint_ring&
R)
1853 make_modular_map map(
R);
1854 return map(e.expand());
1874 static vector<ex> multivar_diophant(
const vector<ex>& a_,
const ex&
x,
const ex&
c,
const vector<EvalPoint>&
I,
1875 unsigned int d,
unsigned int p,
unsigned int k)
1879 const cl_I modulus = expt_pos(cl_I(p),
k);
1880 const cl_modint_ring
R = find_modint_ring(modulus);
1881 const size_t r = a.size();
1882 const size_t nu =
I.size() + 1;
1886 ex xnu =
I.back().x;
1887 int alphanu =
I.back().evalpoint;
1890 for (
size_t i=0; i<
r; ++i ) {
1894 for (
size_t i=0; i<
r; ++i ) {
1898 vector<ex> anew = a;
1899 for (
size_t i=0; i<
r; ++i ) {
1900 anew[i] = anew[i].subs(xnu == alphanu);
1902 ex cnew =
c.subs(xnu == alphanu);
1903 vector<EvalPoint> Inew =
I;
1905 sigma = multivar_diophant(anew,
x, cnew, Inew, d, p,
k);
1908 for (
size_t i=0; i<
r; ++i ) {
1909 buf -= sigma[i] * b[i];
1911 ex e = make_modular(buf,
R);
1914 for (
size_t m=1; !e.is_zero() && e.has(xnu) &&
m<=d; ++
m ) {
1915 monomial *= (xnu - alphanu);
1916 monomial =
expand(monomial);
1917 ex cm = e.diff(ex_to<symbol>(xnu),
m).subs(xnu==alphanu) /
factorial(
m);
1918 cm = make_modular(cm,
R);
1919 if ( !cm.is_zero() ) {
1920 vector<ex> delta_s = multivar_diophant(anew,
x, cm, Inew, d, p,
k);
1922 for (
size_t j=0; j<delta_s.size(); ++j ) {
1923 delta_s[j] *= monomial;
1924 sigma[j] += delta_s[j];
1925 buf -= delta_s[j] * b[j];
1927 e = make_modular(buf,
R);
1932 for (
size_t i=0; i<a.size(); ++i ) {
1934 umodpoly_from_ex(up, a[i],
x,
R);
1938 sigma.insert(sigma.begin(),
r, 0);
1941 if ( is_a<add>(
c) ) {
1948 for (
size_t i=0; i<nterms; ++i ) {
1949 int m = z.degree(
x);
1950 cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(
x)).to_cl_N());
1951 upvec delta_s = univar_diophant(amod,
x,
m, p,
k);
1953 cl_I poscm = plusp(cm) ? cm :
mod(cm, modulus);
1954 modcm = cl_MI(
R, poscm);
1955 for (
size_t j=0; j<delta_s.size(); ++j ) {
1956 delta_s[j] = delta_s[j] * modcm;
1957 sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j],
x);
1959 if ( nterms > 1 && i+1 != nterms ) {
1965 for (
size_t i=0; i<sigma.size(); ++i ) {
1966 sigma[i] = make_modular(sigma[i],
R);
1988 static ex hensel_multivar(
const ex& a,
const ex&
x,
const vector<EvalPoint>&
I,
1989 unsigned int p,
const cl_I& l,
const upvec& u,
const vector<ex>& lcU)
1991 const size_t nu =
I.size() + 1;
1992 const cl_modint_ring
R = find_modint_ring(expt_pos(cl_I(p),l));
1997 for (
size_t j=nu; j>=2; --j ) {
1999 int alpha =
I[j-2].evalpoint;
2000 A[j-2] = A[j-1].
subs(
x==alpha);
2001 A[j-2] = make_modular(A[j-2],
R);
2004 int maxdeg = a.degree(
I.front().x);
2005 for (
size_t i=1; i<
I.size(); ++i ) {
2006 int maxdeg2 = a.degree(
I[i].
x);
2007 if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
2010 const size_t n = u.size();
2012 for (
size_t i=0; i<
n; ++i ) {
2013 U[i] = umodpoly_to_ex(u[i],
x);
2016 for (
size_t j=2; j<=nu; ++j ) {
2019 for (
size_t m=0;
m<
n; ++
m) {
2020 if ( lcU[
m] != 1 ) {
2022 for (
size_t i=j-1; i<nu-1; ++i ) {
2025 coef = make_modular(coef,
R);
2026 int deg = U[
m].degree(
x);
2027 U[
m] = U[
m] - U[
m].lcoeff(
x) *
pow(
x,deg) + coef *
pow(
x,deg);
2031 for (
size_t i=0; i<
n; ++i ) {
2034 ex e =
expand(A[j-1] - Uprod);
2036 vector<EvalPoint> newI;
2037 for (
size_t i=1; i<=j-2; ++i ) {
2038 newI.push_back(
I[i-1]);
2042 int alphaj =
I[j-2].evalpoint;
2043 size_t deg = A[j-1].
degree(xj);
2044 for (
size_t k=1;
k<=deg; ++
k ) {
2045 if ( !e.is_zero() ) {
2046 monomial *= (xj - alphaj);
2047 monomial =
expand(monomial);
2048 ex dif = e.diff(ex_to<symbol>(xj),
k);
2050 if ( !
c.is_zero() ) {
2051 vector<ex> deltaU = multivar_diophant(U1,
x,
c, newI, maxdeg, p, cl_I_to_uint(l));
2052 for (
size_t i=0; i<
n; ++i ) {
2053 deltaU[i] *= monomial;
2055 U[i] = make_modular(U[i],
R);
2058 for (
size_t i=0; i<
n; ++i ) {
2062 e = make_modular(e,
R);
2069 for (
size_t i=0; i<U.size(); ++i ) {
2073 return lst(U.begin(), U.end());
2083 static ex put_factors_into_lst(
const ex& e)
2086 if ( is_a<numeric>(e) ) {
2090 if ( is_a<power>(e) ) {
2092 result.append(e.op(0));
2095 if ( is_a<symbol>(e) || is_a<add>(e) ) {
2096 ex icont(e.integer_content());
2097 result.append(icont);
2098 result.append(e/icont);
2101 if ( is_a<mul>(e) ) {
2103 for (
size_t i=0; i<e.nops(); ++i ) {
2105 if ( is_a<numeric>(
op) ) {
2108 if ( is_a<power>(
op) ) {
2109 result.append(
op.
op(0));
2111 if ( is_a<symbol>(
op) || is_a<add>(
op) ) {
2115 result.prepend(nfac);
2118 throw runtime_error(
"put_factors_into_lst: bad term.");
2127 static bool checkdivisors(
const lst& f)
2129 const int k = f.nops();
2131 vector<numeric> d(
k);
2132 d[0] = ex_to<numeric>(
abs(f.op(0)));
2133 for (
int i=1; i<
k; ++i ) {
2134 q = ex_to<numeric>(
abs(f.op(i)));
2135 for (
int j=i-1; j>=0; --j ) {
2166 static void generate_set(
const ex& u,
const ex& vn,
const exset&
syms,
const lst& f,
2167 numeric& modulus, ex& u0, vector<numeric>& a)
2169 const ex&
x = *
syms.begin();
2176 exset::const_iterator s =
syms.begin();
2178 for (
size_t i=0; i<a.size(); ++i ) {
2180 a[i] =
mod(numeric(rand()), 2*modulus) - modulus;
2181 vnatry = vna.
subs(*s == a[i]);
2183 }
while ( vnatry == 0 );
2185 u0 = u0.
subs(*s == a[i]);
2189 ex g =
gcd(u0, u0.diff(ex_to<symbol>(
x)));
2190 if ( !is_a<numeric>(g) ) {
2193 if ( !is_a<numeric>(vn) ) {
2196 fnum.let_op(0) = fnum.op(0) * u0.content(
x);
2197 for (
size_t i=1; i<fnum.nops(); ++i ) {
2198 if ( !is_a<numeric>(fnum.op(i)) ) {
2201 for (
size_t j=0; j<a.size(); ++j, ++s ) {
2202 fnum.let_op(i) = fnum.op(i).subs(*s == a[j]);
2206 if ( checkdivisors(fnum) ) {
2216 static ex factor_sqrfree(
const ex&
poly);
2233 static ex factor_multivariate(
const ex&
poly,
const exset&
syms)
2235 const ex&
x = *
syms.begin();
2239 poly.unitcontprim(
x, unit, cont, pp);
2240 if ( !is_a<numeric>(cont) ) {
2241 return unit * factor_sqrfree(cont) * factor_sqrfree(pp);
2245 ex vn = pp.collect(
x).lcoeff(
x);
2247 if ( is_a<numeric>(vn) ) {
2251 ex vnfactors =
factor(vn);
2252 vnlst = put_factors_into_lst(vnfactors);
2255 const unsigned int maxtrials = 3;
2256 numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3;
2257 vector<numeric> a(
syms.size()-1, 0);
2262 unsigned int trialcount = 0;
2264 int factor_count = 0;
2265 int min_factor_count = -1;
2270 while ( trialcount < maxtrials ) {
2273 generate_set(pp, vn,
syms, ex_to<lst>(vnlst), modulus, u, a);
2275 ufac = factor_univariate(u,
x, prime);
2276 ufaclst = put_factors_into_lst(ufac);
2277 factor_count = ufaclst.nops()-1;
2278 delta = ufaclst.op(0);
2280 if ( factor_count <= 1 ) {
2284 if ( min_factor_count < 0 ) {
2286 min_factor_count = factor_count;
2288 else if ( min_factor_count == factor_count ) {
2292 else if ( min_factor_count > factor_count ) {
2294 min_factor_count = factor_count;
2300 vector<ex> C(factor_count);
2301 if ( is_a<numeric>(vn) ) {
2303 for (
size_t i=1; i<ufaclst.nops(); ++i ) {
2304 C[i-1] = ufaclst.op(i).lcoeff(
x);
2311 vector<numeric> ftilde(vnlst.nops()-1);
2312 for (
size_t i=0; i<ftilde.size(); ++i ) {
2313 ex ft = vnlst.op(i+1);
2314 auto s =
syms.begin();
2316 for (
size_t j=0; j<a.size(); ++j ) {
2317 ft = ft.subs(*s == a[j]);
2320 ftilde[i] = ex_to<numeric>(ft);
2323 vector<bool> used_flag(ftilde.size(),
false);
2324 vector<ex> D(factor_count, 1);
2326 for (
int i=0; i<factor_count; ++i ) {
2327 numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(
x));
2328 for (
int j=ftilde.size()-1; j>=0; --j ) {
2330 while (
irem(prefac, ftilde[j]) == 0 ) {
2331 prefac =
iquo(prefac, ftilde[j]);
2335 used_flag[j] =
true;
2336 D[i] = D[i] *
pow(vnlst.op(j+1), count);
2339 C[i] = D[i] * prefac;
2342 for (
int i=0; i<factor_count; ++i ) {
2343 numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(
x));
2344 for (
int j=ftilde.size()-1; j>=0; --j ) {
2346 while (
irem(prefac, ftilde[j]) == 0 ) {
2347 prefac =
iquo(prefac, ftilde[j]);
2350 while (
irem(ex_to<numeric>(delta)*prefac, ftilde[j]) == 0 ) {
2351 numeric g =
gcd(prefac, ex_to<numeric>(ftilde[j]));
2352 prefac =
iquo(prefac, g);
2353 delta = delta / (ftilde[j]/g);
2354 ufaclst.let_op(i+1) = ufaclst.op(i+1) * (ftilde[j]/g);
2358 used_flag[j] =
true;
2359 D[i] = D[i] *
pow(vnlst.op(j+1), count);
2362 C[i] = D[i] * prefac;
2366 bool some_factor_unused =
false;
2367 for (
size_t i=0; i<used_flag.size(); ++i ) {
2368 if ( !used_flag[i] ) {
2369 some_factor_unused =
true;
2373 if ( some_factor_unused ) {
2381 C[0] = C[0] * delta;
2382 ufaclst.let_op(1) = ufaclst.op(1) * delta;
2387 vector<EvalPoint> epv;
2388 auto s =
syms.begin();
2390 for (
size_t i=0; i<a.size(); ++i ) {
2392 ep.evalpoint = a[i].to_int();
2398 for (
int i=1; i<=factor_count; ++i ) {
2399 if ( ufaclst.op(i).degree(
x) > maxdeg ) {
2400 maxdeg = ufaclst[i].degree(
x);
2403 cl_I B = 2*calc_bound(u,
x, maxdeg);
2412 cl_modint_ring
R = find_modint_ring(expt_pos(cl_I(prime),l));
2413 upvec modfactors(ufaclst.nops()-1);
2414 for (
size_t i=1; i<ufaclst.nops(); ++i ) {
2415 umodpoly_from_ex(modfactors[i-1], ufaclst.op(i),
x,
R);
2419 ex res = hensel_multivar(pp,
x, epv, prime, l, modfactors, C);
2420 if ( res !=
lst{} ) {
2421 ex result = cont * unit;
2422 for (
size_t i=0; i<res.nops(); ++i ) {
2423 result *= res.op(i).content(
x) * res.op(i).unit(
x);
2424 result *= res.op(i).primpart(
x);
2433 struct find_symbols_map :
public map_function {
2435 ex operator()(
const ex& e)
override
2437 if ( is_a<symbol>(e) ) {
2441 return e.map(*
this);
2448 static ex factor_sqrfree(
const ex&
poly)
2451 find_symbols_map findsymbols;
2453 if ( findsymbols.syms.size() == 0 ) {
2457 if ( findsymbols.syms.size() == 1 ) {
2459 const ex&
x = *(findsymbols.syms.begin());
2460 int ld =
poly.ldegree(
x);
2464 return res *
pow(
x,ld);
2466 ex res = factor_univariate(
poly,
x);
2472 ex res = factor_multivariate(
poly, findsymbols.syms);
2479 struct apply_factor_map :
public map_function {
2481 apply_factor_map(
unsigned options_) :
options(options_) { }
2482 ex operator()(
const ex& e)
override
2484 if ( e.info(info_flags::polynomial) ) {
2487 if ( is_a<add>(e) ) {
2489 for (
size_t i=0; i<e.nops(); ++i ) {
2490 if ( e.op(i).info(info_flags::polynomial) ) {
2498 return e.
map(*
this);
2508 template <
typename F>
void
2509 factor_iter(
const ex &e, F yield)
2512 for (
const auto &f : e) {
2513 if (is_a<power>(f)) {
2514 yield(f.op(0), f.op(1));
2520 if (is_a<power>(e)) {
2521 yield(e.op(0), e.op(1));
2536 static ex factor1(
const ex&
poly,
unsigned options)
2539 if ( !
poly.info(info_flags::polynomial) ) {
2540 if (
options & factor_options::all ) {
2541 options &= ~factor_options::all;
2542 apply_factor_map factor_map(
options);
2543 return factor_map(
poly);
2549 find_symbols_map findsymbols;
2551 if ( findsymbols.syms.size() == 0 ) {
2555 for (
auto & i : findsymbols.syms ) {
2565 [&](
const ex &f,
const ex &
k) {
2566 if ( is_a<add>(f) ) {
2567 res *=
pow(factor_sqrfree(f),
k);
2585 [&](
const ex &f1,
const ex &k1) {
2586 factor_iter(factor1(f1,
options),
2587 [&](
const ex &f2,
const ex &k2) {
2588 result *=
pow(f2, k1*k2);
Interface to GiNaC's sums of expressions.
Lightweight wrapper for GiNaC's symbolic objects.
ex map(map_function &f) const
ex subs(const exmap &m, unsigned options=0) const
ex subs(const exmap &m, unsigned options=0) const override
Substitute a set of objects by arbitrary expressions.
int degree(const ex &s) const override
Return degree of highest power in object s.
Interface to GiNaC's light-weight expression handles.
vector< vector< umodpoly > > cache
Polynomial factorization.
Interface to GiNaC's initially known functions.
Interface to GiNaC's products of expressions.
bool is_zero(const ex &thisex)
const numeric I
Imaginary unit.
const numeric pow(const numeric &x, const numeric &y)
container< std::list > lst
ex sqrfree(const ex &a, const lst &l)
Compute a square-free factorization of a multivariate polynomial in Q[X].
std::set< ex, ex_is_less > exset
const numeric mod(const numeric &a, const numeric &b)
Modulus (in positive representation).
const ex operator/(const ex &lh, const ex &rh)
const numeric abs(const numeric &x)
Absolute value.
ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned options)
Compute GCD (Greatest Common Divisor) of multivariate polynomials a(X) and b(X) in Z[X].
const numeric sqrt(const numeric &x)
Numeric square root.
const numeric irem(const numeric &a, const numeric &b)
Numeric integer remainder.
const ex operator+(const ex &lh, const ex &rh)
const ex operator*(const ex &lh, const ex &rh)
const numeric factorial(const numeric &n)
Factorial combinatorial function.
const numeric smod(const numeric &a_, const numeric &b_)
Modulus (in symmetric representation).
int degree(const ex &thisex, const ex &s)
std::ostream & operator<<(std::ostream &os, const archive_node &n)
Write archive_node to binary data stream.
const numeric iquo(const numeric &a, const numeric &b)
Numeric integer quotient.
ex normal(const ex &thisex)
ex op(const ex &thisex, size_t i)
ex coeff(const ex &thisex, const ex &s, int n=1)
ex factor(const ex &poly, unsigned options)
Interface function to the outside world.
int canonicalize(exvector::iterator v, const symmetry &symm)
Canonicalize the order of elements of an expression vector, according to the symmetry properties defi...
ex rem(const ex &a, const ex &b, const ex &x, bool check_args)
Remainder r(x) of polynomials a(x) and b(x) in Q[x].
const ex operator-(const ex &lh, const ex &rh)
ex expand(const ex &thisex, unsigned options=0)
void swap(GiNaC::ex &a, GiNaC::ex &b)
Specialization of std::swap() for ex objects.
This file defines several functions that work on univariate and multivariate polynomials and rational...
Makes the interface to the underlying bignum package available.
Interface to GiNaC's overloaded operators.
Interface to GiNaC's symbolic exponentiation (basis^exponent).
Interface to relations between expressions.
Interface to GiNaC's symbolic objects.