typedef vector<cl_MI> mvec;
#ifdef DEBUGFACTOR
+ostream& operator<<(ostream& o, const vector<int>& v)
+{
+ vector<int>::const_iterator i = v.begin(), end = v.end();
+ while ( i != end ) {
+ o << *i++ << " ";
+ }
+ return o;
+}
+ostream& operator<<(ostream& o, const vector<cl_I>& v)
+{
+ vector<cl_I>::const_iterator i = v.begin(), end = v.end();
+ while ( i != end ) {
+ o << *i << "[" << i-v.begin() << "]" << " ";
+ ++i;
+ }
+ return o;
+}
ostream& operator<<(ostream& o, const vector<cl_MI>& v)
{
vector<cl_MI>::const_iterator i = v.begin(), end = v.end();
while ( i != end ) {
- o << *i++ << " ";
+ o << *i << "[" << i-v.begin() << "]" << " ";
+ ++i;
}
return o;
}
{
vector< vector<cl_MI> >::const_iterator i = v.begin(), end = v.end();
while ( i != end ) {
- o << *i++ << endl;
+ o << i-v.begin() << ": " << *i << endl;
+ ++i;
}
return o;
}
////////////////////////////////////////////////////////////////////////////////
// modular univariate polynomial code
-//typedef cl_UP_MI umod;
typedef std::vector<cln::cl_MI> umodpoly;
-//typedef vector<umod> umodvec;
+typedef std::vector<cln::cl_I> upoly;
typedef vector<umodpoly> upvec;
// COPY FROM UPOLY.HPP
// END COPY FROM UPOLY.HPP
-static void expt_pos(const umodpoly& a, unsigned int q, umodpoly& b)
-{
- throw runtime_error("expt_pos: not implemented!");
- // code below is not correct!
-// b.clear();
-// if ( a.empty() ) return;
-// b.resize(degree(a)*q+1, a[0].ring()->zero());
-// cl_MI norm = recip(a[0]);
-// umodpoly an = a;
-// for ( size_t i=0; i<an.size(); ++i ) {
-// an[i] = an[i] * norm;
-// }
-// b[0] = a[0].ring()->one();
-// for ( size_t i=1; i<b.size(); ++i ) {
-// for ( size_t j=1; j<i; ++j ) {
-// b[i] = b[i] + ((i-j+1)*q-i-1) * a[i-j] * b[j-1];
-// }
-// b[i] = b[i] / i;
-// }
-// cl_MI corr = expt_pos(a[0], q);
-// for ( size_t i=0; i<b.size(); ++i ) {
-// b[i] = b[i] * corr;
-// }
-}
-
-static umodpoly operator+(const umodpoly& a, const umodpoly& b)
+static void expt_pos(umodpoly& a, unsigned int q)
+{
+ if ( a.empty() ) return;
+ cl_MI zero = a[0].ring()->zero();
+ int deg = degree(a);
+ a.resize(degree(a)*q+1, zero);
+ for ( int i=deg; i>0; --i ) {
+ a[i*q] = a[i];
+ a[i] = zero;
+ }
+}
+
+template<typename T>
+static T operator+(const T& a, const T& b)
{
int sa = a.size();
int sb = b.size();
if ( sa >= sb ) {
- umodpoly r(sa);
+ T r(sa);
int i = 0;
for ( ; i<sb; ++i ) {
r[i] = a[i] + b[i];
return r;
}
else {
- umodpoly r(sb);
+ T r(sb);
int i = 0;
for ( ; i<sa; ++i ) {
r[i] = a[i] + b[i];
}
}
-static umodpoly operator-(const umodpoly& a, const umodpoly& b)
+template<typename T>
+static T operator-(const T& a, const T& b)
{
int sa = a.size();
int sb = b.size();
if ( sa >= sb ) {
- umodpoly r(sa);
+ T r(sa);
int i = 0;
for ( ; i<sb; ++i ) {
r[i] = a[i] - b[i];
return r;
}
else {
- umodpoly r(sb);
+ T r(sb);
int i = 0;
for ( ; i<sa; ++i ) {
r[i] = a[i] - b[i];
}
}
+static upoly operator*(const upoly& a, const upoly& b)
+{
+ upoly c;
+ if ( a.empty() || b.empty() ) return c;
+
+ int n = degree(a) + degree(b);
+ c.resize(n+1, 0);
+ for ( int i=0 ; i<=n; ++i ) {
+ for ( int j=0 ; j<=i; ++j ) {
+ if ( j > degree(a) || (i-j) > degree(b) ) continue;
+ c[i] = c[i] + a[j] * b[i-j];
+ }
+ }
+ canonicalize(c);
+ return c;
+}
+
static umodpoly operator*(const umodpoly& a, const umodpoly& b)
{
umodpoly c;
c.resize(n+1, a[0].ring()->zero());
for ( int i=0 ; i<=n; ++i ) {
for ( int j=0 ; j<=i; ++j ) {
- if ( j > degree(a) || (i-j) > degree(b) ) continue; // TODO optimize!
+ if ( j > degree(a) || (i-j) > degree(b) ) continue;
c[i] = c[i] + a[j] * b[i-j];
}
}
return c;
}
+static upoly operator*(const upoly& a, const cl_I& x)
+{
+ if ( zerop(x) ) {
+ upoly r;
+ return r;
+ }
+ upoly r(a.size());
+ for ( size_t i=0; i<a.size(); ++i ) {
+ r[i] = a[i] * x;
+ }
+ return r;
+}
+
+static upoly operator/(const upoly& a, const cl_I& x)
+{
+ if ( zerop(x) ) {
+ upoly r;
+ return r;
+ }
+ upoly r(a.size());
+ for ( size_t i=0; i<a.size(); ++i ) {
+ r[i] = exquo(a[i],x);
+ }
+ return r;
+}
+
static umodpoly operator*(const umodpoly& a, const cl_MI& x)
{
umodpoly r(a.size());
return r;
}
+static void upoly_from_ex(upoly& up, const ex& e, const ex& x)
+{
+ // assert: e is in Z[x]
+ int deg = e.degree(x);
+ up.resize(deg+1);
+ int ldeg = e.ldegree(x);
+ for ( ; deg>=ldeg; --deg ) {
+ up[deg] = the<cl_I>(ex_to<numeric>(e.coeff(x, deg)).to_cl_N());
+ }
+ for ( ; deg>=0; --deg ) {
+ up[deg] = 0;
+ }
+ canonicalize(up);
+}
+
+static void umodpoly_from_upoly(umodpoly& ump, const upoly& e, const cl_modint_ring& R)
+{
+ int deg = degree(e);
+ ump.resize(deg+1);
+ for ( ; deg>=0; --deg ) {
+ ump[deg] = R->canonhom(e[deg]);
+ }
+ canonicalize(ump);
+}
+
static void umodpoly_from_ex(umodpoly& ump, const ex& e, const ex& x, const cl_modint_ring& R)
{
// assert: e is in Z[x]
umodpoly_from_ex(ump, e, x, find_modint_ring(modulus));
}
-static ex umod_to_ex(const umodpoly& a, const ex& x)
+static ex upoly_to_ex(const upoly& a, const ex& x)
+{
+ if ( a.empty() ) return 0;
+ ex e;
+ for ( int i=degree(a); i>=0; --i ) {
+ e += numeric(a[i]) * pow(x, i);
+ }
+ return e;
+}
+
+static ex umodpoly_to_ex(const umodpoly& a, const ex& x)
{
if ( a.empty() ) return 0;
cl_modint_ring R = a[0].ring();
return e;
}
+static upoly umodpoly_to_upoly(const umodpoly& a)
+{
+ upoly e(a.size());
+ if ( a.empty() ) return e;
+ cl_modint_ring R = a[0].ring();
+ cl_I mod = R->modulus;
+ cl_I halfmod = (mod-1) >> 1;
+ for ( int i=degree(a); i>=0; --i ) {
+ cl_I n = R->retract(a[i]);
+ if ( n > halfmod ) {
+ e[i] = n-mod;
+ } else {
+ e[i] = n;
+ }
+ }
+ return e;
+}
+
+static umodpoly umodpoly_to_umodpoly(const umodpoly& a, const cl_modint_ring& R, unsigned int m)
+{
+ umodpoly e;
+ if ( a.empty() ) return e;
+ cl_modint_ring oldR = a[0].ring();
+ size_t sa = a.size();
+ e.resize(sa+m, R->zero());
+ for ( size_t i=0; i<sa; ++i ) {
+ e[i+m] = R->canonhom(oldR->retract(a[i]));
+ }
+ canonicalize(e);
+ return e;
+}
+
/** Divides all coefficients of the polynomial a by the integer x.
* All coefficients are supposed to be divisible by x. If they are not, the
* the<cl_I> cast will raise an exception.
umodpoly b;
deriv(a, b);
if ( b.empty() ) {
- return true;
+ return false;
}
umodpoly c;
gcd(a, b, c);
ostream& operator<<(ostream& o, const modular_matrix& m)
{
- vector<cl_MI>::const_iterator i = m.m.begin(), end = m.m.end();
- size_t wrap = 1;
- for ( ; i != end; ++i ) {
- o << *i << " ";
- if ( !(wrap++ % m.c) ) {
- o << endl;
+ cl_modint_ring R = m(0,0).ring();
+ o << "{";
+ for ( size_t i=0; i<m.rowsize(); ++i ) {
+ o << "{";
+ for ( size_t j=0; j<m.colsize()-1; ++j ) {
+ o << R->retract(m(i,j)) << ",";
+ }
+ o << R->retract(m(i,m.colsize()-1)) << "}";
+ if ( i != m.rowsize()-1 ) {
+ o << ",";
}
}
- o << endl;
+ o << "}";
return o;
}
#endif // def DEBUGFACTOR
// END modular matrix
////////////////////////////////////////////////////////////////////////////////
-static void q_matrix(const umodpoly& a, modular_matrix& Q)
+static void q_matrix(const umodpoly& a_, modular_matrix& Q)
{
+ umodpoly a = a_;
+ normalize_in_field(a);
+
int n = degree(a);
unsigned int q = cl_I_to_uint(a[0].ring()->modulus);
-// fast and buggy
-// vector<cl_MI> r(n, a.R->zero());
-// r[0] = a.R->one();
-// Q.set_row(0, r);
-// unsigned int max = (n-1) * q;
-// for ( size_t m=1; m<=max; ++m ) {
-// cl_MI rn_1 = r.back();
-// for ( size_t i=n-1; i>0; --i ) {
-// r[i] = r[i-1] - rn_1 * a[i];
-// }
-// r[0] = -rn_1 * a[0];
-// if ( (m % q) == 0 ) {
-// Q.set_row(m/q, r);
-// }
-// }
-// slow and (hopefully) correct
- cl_MI one = a[0].ring()->one();
- cl_MI zero = a[0].ring()->zero();
- for ( int i=0; i<n; ++i ) {
- umodpoly qk(i*q+1, zero);
- qk[i*q] = one;
- umodpoly r;
- rem(qk, a, r);
- mvec rvec(n, zero);
- for ( int j=0; j<=degree(r); ++j ) {
- rvec[j] = r[j];
- }
- Q.set_row(i, rvec);
+ umodpoly r(n, a[0].ring()->zero());
+ r[0] = a[0].ring()->one();
+ Q.set_row(0, r);
+ unsigned int max = (n-1) * q;
+ for ( size_t m=1; m<=max; ++m ) {
+ cl_MI rn_1 = r.back();
+ for ( size_t i=n-1; i>0; --i ) {
+ r[i] = r[i-1] - (rn_1 * a[i]);
+ }
+ r[0] = -rn_1 * a[0];
+ if ( (m % q) == 0 ) {
+ Q.set_row(m/q, r);
+ }
}
}
q_matrix(a, Q);
vector<mvec> nu;
nullspace(Q, nu);
+
const unsigned int k = nu.size();
if ( k == 1 ) {
return;
}
}
-static void rem_xq(int q, const umodpoly& b, umodpoly& c)
+static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
{
- cl_modint_ring R = b[0].ring();
-
- int n = degree(b);
- if ( n > q ) {
- c.resize(q+1, R->zero());
- c[q] = R->one();
- return;
+ size_t newdeg = degree(a)/prime;
+ ap.resize(newdeg+1);
+ ap[0] = a[0];
+ for ( size_t i=1; i<=newdeg; ++i ) {
+ ap[i] = a[i*prime];
}
+}
- c.clear();
- c.resize(n+1, R->zero());
- int k = q-n;
- c[n] = R->one();
-
- int ofs = 0;
- do {
- cl_MI qk = div(c[n-ofs], b[n]);
- if ( !zerop(qk) ) {
- for ( int i=1; i<=n; ++i ) {
- c[n-i+ofs] = c[n-i] - qk * b[n-i];
+static void modsqrfree(const umodpoly& a, upvec& factors, vector<int>& mult)
+{
+ const unsigned int prime = cl_I_to_uint(a[0].ring()->modulus);
+ int i = 1;
+ umodpoly b;
+ deriv(a, b);
+ if ( b.size() ) {
+ umodpoly c;
+ gcd(a, b, c);
+ umodpoly w;
+ div(a, c, w);
+ while ( unequal_one(w) ) {
+ umodpoly y;
+ gcd(w, c, y);
+ umodpoly z;
+ div(w, y, z);
+ factors.push_back(z);
+ mult.push_back(i);
+ ++i;
+ w = y;
+ umodpoly buf;
+ div(c, y, buf);
+ c = buf;
+ }
+ if ( unequal_one(c) ) {
+ umodpoly cp;
+ expt_1_over_p(c, prime, cp);
+ size_t previ = mult.size();
+ modsqrfree(cp, factors, mult);
+ for ( size_t i=previ; i<mult.size(); ++i ) {
+ mult[i] *= prime;
}
- ofs = ofs ? 0 : 1;
}
- } while ( k-- );
-
- if ( ofs ) {
- c.pop_back();
}
else {
- c.erase(c.begin());
+ umodpoly ap;
+ expt_1_over_p(a, prime, ap);
+ size_t previ = mult.size();
+ modsqrfree(ap, factors, mult);
+ for ( size_t i=previ; i<mult.size(); ++i ) {
+ mult[i] *= prime;
+ }
}
- canonicalize(c);
}
-static void distinct_degree_factor(const umodpoly& a_, upvec& result)
+static void distinct_degree_factor(const umodpoly& a_, vector<int>& degrees, upvec& ddfactors)
{
umodpoly a = a_;
cl_modint_ring R = a[0].ring();
int q = cl_I_to_int(R->modulus);
- int n = degree(a);
- size_t nhalf = n/2;
+ int nhalf = degree(a)/2;
- size_t i = 1;
- umodpoly w(1, R->one());
+ int i = 1;
+ umodpoly w(2);
+ w[0] = R->zero();
+ w[1] = R->one();
umodpoly x = w;
- upvec ai;
-
while ( i <= nhalf ) {
- expt_pos(w, q, w);
- rem(w, a, w);
-
+ expt_pos(w, q);
umodpoly buf;
- gcd(a, w-x, buf);
- ai.push_back(buf);
-
- if ( unequal_one(ai.back()) ) {
- div(a, ai.back(), a);
- rem(w, a, w);
+ rem(w, a, buf);
+ w = buf;
+ umodpoly wx = w - x;
+ gcd(a, wx, buf);
+ if ( unequal_one(buf) ) {
+ degrees.push_back(i);
+ ddfactors.push_back(buf);
+ }
+ if ( unequal_one(buf) ) {
+ umodpoly buf2;
+ div(a, buf, buf2);
+ a = buf2;
+ nhalf = degree(a)/2;
+ rem(w, a, buf);
+ w = buf;
}
-
++i;
}
-
- result = ai;
-}
-
-static void same_degree_factor(const umodpoly& a, upvec& result)
-{
- cl_modint_ring R = a[0].ring();
- int deg = degree(a);
-
- upvec buf;
- distinct_degree_factor(a, buf);
- int degsum = 0;
-
- for ( size_t i=0; i<buf.size(); ++i ) {
- if ( unequal_one(buf[i]) ) {
- degsum += degree(buf[i]);
- upvec upv;
- berlekamp(buf[i], upv);
- for ( size_t j=0; j<upv.size(); ++j ) {
- result.push_back(upv[j]);
- }
- }
- }
-
- if ( degsum < deg ) {
- result.push_back(a);
+ if ( unequal_one(a) ) {
+ degrees.push_back(degree(a));
+ ddfactors.push_back(a);
}
}
-static void distinct_degree_factor_BSGS(const umodpoly& a, upvec& result)
+static void same_degree_factor(const umodpoly& a, upvec& upv)
{
cl_modint_ring R = a[0].ring();
- int q = cl_I_to_int(R->modulus);
- int n = degree(a);
- cl_N pm = 0.3;
- int l = cl_I_to_int(ceiling1(the<cl_F>(expt(n, pm))));
- upvec h(l+1);
- umodpoly qk(1, R->one());
- h[0] = qk;
- for ( int i=1; i<=l; ++i ) {
- expt_pos(h[i-1], q, qk);
- rem(qk, a, h[i]);
- }
+ vector<int> degrees;
+ upvec ddfactors;
+ distinct_degree_factor(a, degrees, ddfactors);
- int m = std::ceil(((double)n)/2/l);
- upvec H(m);
- int ql = std::pow(q, l);
- H[0] = h[l];
- for ( int i=1; i<m; ++i ) {
- expt_pos(H[i-1], ql, qk);
- rem(qk, a, H[i]);
- }
-
- upvec I(m);
- umodpoly one(1, R->one());
- for ( int i=0; i<m; ++i ) {
- I[i] = one;
- for ( int j=0; j<l; ++j ) {
- I[i] = I[i] * (H[i] - h[j]);
+ for ( size_t i=0; i<degrees.size(); ++i ) {
+ if ( degrees[i] == degree(ddfactors[i]) ) {
+ upv.push_back(ddfactors[i]);
}
- rem(I[i], a, I[i]);
- }
-
- upvec F(m, one);
- umodpoly f = a;
- for ( int i=0; i<m; ++i ) {
- umodpoly g;
- gcd(f, I[i], g);
- if ( g == one ) continue;
- F[i] = g;
- div(f, g, f);
- }
-
- result.resize(n, one);
- if ( unequal_one(f) ) {
- result[n] = f;
- }
- for ( int i=0; i<m; ++i ) {
- umodpoly f = F[i];
- for ( int j=l-1; j>=0; --j ) {
- umodpoly g;
- gcd(f, H[i]-h[j], g);
- result[l*(i+1)-j-1] = g;
- div(f, g, f);
+ else {
+ berlekamp(ddfactors[i], upv);
}
}
}
-static void cantor_zassenhaus(const umodpoly& a, upvec& result)
-{
-}
+#define USE_SAME_DEGREE_FACTOR
static void factor_modular(const umodpoly& p, upvec& upv)
{
- //same_degree_factor(p, upv);
+#ifdef USE_SAME_DEGREE_FACTOR
+ same_degree_factor(p, upv);
+#else
berlekamp(p, upv);
- return;
+#endif
}
-static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& g, umodpoly& s, umodpoly& t)
+/** Calculates polynomials s and t such that a*s+b*t==1.
+ * Assertion: a and b are relatively prime and not zero.
+ *
+ * @param[in] a polynomial
+ * @param[in] b polynomial
+ * @param[out] s polynomial
+ * @param[out] t polynomial
+ */
+static void exteuclid(const umodpoly& a, const umodpoly& b, umodpoly& s, umodpoly& t)
{
if ( degree(a) < degree(b) ) {
- exteuclid(b, a, g, t, s);
+ exteuclid(b, a, t, s);
return;
}
+
umodpoly one(1, a[0].ring()->one());
umodpoly c = a; normalize_in_field(c);
umodpoly d = b; normalize_in_field(d);
- umodpoly c1 = one;
- umodpoly c2;
+ s = one;
+ t.clear();
umodpoly d1;
umodpoly d2 = one;
- while ( !d.empty() ) {
- umodpoly q;
+ umodpoly q;
+ while ( true ) {
div(c, d, q);
umodpoly r = c - q * d;
- umodpoly r1 = c1 - q * d1;
- umodpoly r2 = c2 - q * d2;
+ umodpoly r1 = s - q * d1;
+ umodpoly r2 = t - q * d2;
c = d;
- c1 = d1;
- c2 = d2;
+ s = d1;
+ t = d2;
+ if ( r.empty() ) break;
d = r;
d1 = r1;
d2 = r2;
}
- g = c; normalize_in_field(g);
- s = c1;
- for ( int i=0; i<=degree(s); ++i ) {
- s[i] = s[i] * recip(a[degree(a)] * c[degree(c)]);
+ cl_MI fac = recip(lcoeff(a) * lcoeff(c));
+ umodpoly::iterator i = s.begin(), end = s.end();
+ for ( ; i!=end; ++i ) {
+ *i = *i * fac;
}
canonicalize(s);
- s = s * g;
- t = c2;
- for ( int i=0; i<=degree(t); ++i ) {
- t[i] = t[i] * recip(b[degree(b)] * c[degree(c)]);
+ fac = recip(lcoeff(b) * lcoeff(c));
+ i = t.begin(), end = t.end();
+ for ( ; i!=end; ++i ) {
+ *i = *i * fac;
}
canonicalize(t);
- t = t * g;
}
-static ex replace_lc(const ex& poly, const ex& x, const ex& lc)
+static upoly replace_lc(const upoly& poly, const cl_I& lc)
{
- ex r = expand(poly + (lc - poly.lcoeff(x)) * pow(x, poly.degree(x)));
+ if ( poly.empty() ) return poly;
+ upoly r = poly;
+ r.back() = lc;
return r;
}
-static ex hensel_univar(const ex& a_, const ex& x, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, const ex& gamma_ = 0)
+static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
{
- ex a = a_;
+ cl_I maxcoeff = 0;
+ cl_R coeff = 0;
+ for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
+ cl_I aa = abs(the<cl_I>(ex_to<numeric>(a.coeff(x, i)).to_cl_N()));
+ if ( aa > maxcoeff ) maxcoeff = aa;
+ coeff = coeff + square(aa);
+ }
+ cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
+ cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
+ return ( B > maxcoeff ) ? B : maxcoeff;
+}
+
+static inline cl_I calc_bound(const upoly& a, int maxdeg)
+{
+ cl_I maxcoeff = 0;
+ cl_R coeff = 0;
+ for ( int i=degree(a); i>=0; --i ) {
+ cl_I aa = abs(a[i]);
+ if ( aa > maxcoeff ) maxcoeff = aa;
+ coeff = coeff + square(aa);
+ }
+ cl_I coeffnorm = ceiling1(the<cl_R>(cln::sqrt(coeff)));
+ cl_I B = coeffnorm * expt_pos(cl_I(2), cl_I(maxdeg));
+ return ( B > maxcoeff ) ? B : maxcoeff;
+}
+
+static void hensel_univar(const upoly& a_, unsigned int p, const umodpoly& u1_, const umodpoly& w1_, upoly& u, upoly& w)
+{
+ upoly a = a_;
const cl_modint_ring& R = u1_[0].ring();
// calc bound B
- ex maxcoeff;
- for ( int i=a.degree(x); i>=a.ldegree(x); --i ) {
- maxcoeff += pow(abs(a.coeff(x, i)),2);
- }
- cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
- cl_I maxdegree = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
- cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
+ int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
+ cl_I maxmodulus = 2*calc_bound(a, maxdeg);
// step 1
- ex alpha = a.lcoeff(x);
- ex gamma = gamma_;
- if ( gamma == 0 ) {
- gamma = alpha;
- }
- numeric gamma_ui = ex_to<numeric>(abs(gamma));
- a = a * gamma;
+ cl_I alpha = lcoeff(a);
+ a = a * alpha;
umodpoly nu1 = u1_;
normalize_in_field(nu1);
umodpoly nw1 = w1_;
normalize_in_field(nw1);
- ex phi;
- phi = gamma * umod_to_ex(nu1, x);
+ upoly phi;
+ phi = umodpoly_to_upoly(nu1) * alpha;
umodpoly u1;
- umodpoly_from_ex(u1, phi, x, R);
- phi = alpha * umod_to_ex(nw1, x);
+ umodpoly_from_upoly(u1, phi, R);
+ phi = umodpoly_to_upoly(nw1) * alpha;
umodpoly w1;
- umodpoly_from_ex(w1, phi, x, R);
+ umodpoly_from_upoly(w1, phi, R);
// step 2
- umodpoly g;
umodpoly s;
umodpoly t;
- exteuclid(u1, w1, g, s, t);
- if ( unequal_one(g) ) {
- throw logic_error("gcd(u1,w1) != 1");
- }
+ exteuclid(u1, w1, s, t);
// step 3
- ex u = replace_lc(umod_to_ex(u1, x), x, gamma);
- ex w = replace_lc(umod_to_ex(w1, x), x, alpha);
- ex e = expand(a - u * w);
- numeric modulus = p;
- const numeric maxmodulus = 2*numeric(B)*gamma_ui;
+ u = replace_lc(umodpoly_to_upoly(u1), alpha);
+ w = replace_lc(umodpoly_to_upoly(w1), alpha);
+ upoly e = a - u * w;
+ cl_I modulus = p;
// step 4
- while ( !e.is_zero() && modulus < maxmodulus ) {
- ex c = e / modulus;
- phi = expand(umod_to_ex(s, x) * c);
+ while ( !e.empty() && modulus < maxmodulus ) {
+ upoly c = e / modulus;
+ phi = umodpoly_to_upoly(s) * c;
umodpoly sigmatilde;
- umodpoly_from_ex(sigmatilde, phi, x, R);
- phi = expand(umod_to_ex(t, x) * c);
+ umodpoly_from_upoly(sigmatilde, phi, R);
+ phi = umodpoly_to_upoly(t) * c;
umodpoly tautilde;
- umodpoly_from_ex(tautilde, phi, x, R);
+ umodpoly_from_upoly(tautilde, phi, R);
umodpoly r, q;
remdiv(sigmatilde, w1, r, q);
umodpoly sigma = r;
- phi = expand(umod_to_ex(tautilde, x) + umod_to_ex(q, x) * umod_to_ex(u1, x));
+ phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1);
umodpoly tau;
- umodpoly_from_ex(tau, phi, x, R);
- u = expand(u + umod_to_ex(tau, x) * modulus);
- w = expand(w + umod_to_ex(sigma, x) * modulus);
- e = expand(a - u * w);
+ umodpoly_from_upoly(tau, phi, R);
+ u = u + umodpoly_to_upoly(tau) * modulus;
+ w = w + umodpoly_to_upoly(sigma) * modulus;
+ e = a - u * w;
modulus = modulus * p;
}
// step 5
- if ( e.is_zero() ) {
- ex delta = u.content(x);
- u = u / delta;
- w = w / gamma * delta;
- return lst(u, w);
+ if ( e.empty() ) {
+ cl_I g = u[0];
+ for ( size_t i=1; i<u.size(); ++i ) {
+ g = gcd(g, u[i]);
+ if ( g == 1 ) break;
+ }
+ if ( g != 1 ) {
+ u = u / g;
+ w = w * g;
+ }
+ if ( alpha != 1 ) {
+ w = w / alpha;
+ }
}
else {
- return lst();
+ u.clear();
}
}
throw logic_error("next_prime: should not reach this point!");
}
-class Partition
+class factor_partition
{
public:
- Partition(size_t n_) : n(n_)
+ factor_partition(const upvec& factors_) : factors(factors_)
{
- k.resize(n, 1);
- k[0] = 0;
- sum = n-1;
+ n = factors.size();
+ k.resize(n, 0);
+ k[0] = 1;
+ cache.resize(n-1);
+ one.resize(1, factors.front()[0].ring()->one());
+ len = 1;
+ last = 0;
+ split();
}
int operator[](size_t i) const { return k[i]; }
size_t size() const { return n; }
- size_t size_first() const { return n-sum; }
- size_t size_second() const { return sum; }
+ size_t size_left() const { return n-len; }
+ size_t size_right() const { return len; }
#ifdef DEBUGFACTOR
- void get() const
- {
- for ( size_t i=0; i<k.size(); ++i ) {
- cout << k[i] << " ";
- }
- cout << endl;
- }
+ void get() const { DCOUTVAR(k); }
#endif
bool next()
{
- for ( size_t i=n-1; i>=1; --i ) {
- if ( k[i] ) {
- --k[i];
- --sum;
- return sum > 0;
+ if ( last == n-1 ) {
+ int rem = len - 1;
+ int p = last - 1;
+ while ( rem ) {
+ if ( k[p] ) {
+ --rem;
+ --p;
+ continue;
+ }
+ last = p - 1;
+ while ( k[last] == 0 ) { --last; }
+ if ( last == 0 && n == 2*len ) return false;
+ k[last++] = 0;
+ for ( size_t i=0; i<=len-rem; ++i ) {
+ k[last] = 1;
+ ++last;
+ }
+ fill(k.begin()+last, k.end(), 0);
+ --last;
+ split();
+ return true;
}
- ++k[i];
- ++sum;
+ last = len;
+ ++len;
+ if ( len > n/2 ) return false;
+ fill(k.begin(), k.begin()+len, 1);
+ fill(k.begin()+len+1, k.end(), 0);
}
- return false;
+ else {
+ k[last++] = 0;
+ k[last] = 1;
+ }
+ split();
+ return true;
}
+ umodpoly& left() { return lr[0]; }
+ umodpoly& right() { return lr[1]; }
private:
- size_t n, sum;
- vector<int> k;
-};
-
-static void split(const upvec& factors, const Partition& part, umodpoly& a, umodpoly& b)
-{
- umodpoly one(1, factors.front()[0].ring()->one());
- a = one;
- b = one;
- for ( size_t i=0; i<part.size(); ++i ) {
- if ( part[i] ) {
- b = b * factors[i];
+ void split_cached()
+ {
+ size_t i = 0;
+ do {
+ size_t pos = i;
+ int group = k[i++];
+ size_t d = 0;
+ while ( i < n && k[i] == group ) { ++d; ++i; }
+ if ( d ) {
+ if ( cache[pos].size() >= d ) {
+ lr[group] = lr[group] * cache[pos][d-1];
+ }
+ else {
+ if ( cache[pos].size() == 0 ) {
+ cache[pos].push_back(factors[pos] * factors[pos+1]);
+ }
+ size_t j = pos + cache[pos].size() + 1;
+ d -= cache[pos].size();
+ while ( d ) {
+ umodpoly buf = cache[pos].back() * factors[j];
+ cache[pos].push_back(buf);
+ --d;
+ ++j;
+ }
+ lr[group] = lr[group] * cache[pos].back();
+ }
+ }
+ else {
+ lr[group] = lr[group] * factors[pos];
+ }
+ } while ( i < n );
+ }
+ void split()
+ {
+ lr[0] = one;
+ lr[1] = one;
+ if ( n > 6 ) {
+ split_cached();
}
else {
- a = a * factors[i];
+ for ( size_t i=0; i<n; ++i ) {
+ lr[k[i]] = lr[k[i]] * factors[i];
+ }
}
}
-}
+private:
+ umodpoly lr[2];
+ vector< vector<umodpoly> > cache;
+ upvec factors;
+ umodpoly one;
+ size_t n;
+ size_t len;
+ size_t last;
+ vector<int> k;
+};
struct ModFactors
{
- ex poly;
+ upoly poly;
upvec factors;
};
static ex factor_univariate(const ex& poly, const ex& x)
{
- ex unit, cont, prim;
- poly.unitcontprim(x, unit, cont, prim);
+ ex unit, cont, prim_ex;
+ poly.unitcontprim(x, unit, cont, prim_ex);
+ upoly prim;
+ upoly_from_ex(prim, prim_ex, x);
// determine proper prime and minimize number of modular factors
unsigned int p = 3, lastp = 3;
cl_modint_ring R;
unsigned int trials = 0;
unsigned int minfactors = 0;
- numeric lcoeff = ex_to<numeric>(prim.lcoeff(x));
+ cl_I lc = lcoeff(prim);
upvec factors;
while ( trials < 2 ) {
+ umodpoly modpoly;
while ( true ) {
p = next_prime(p);
- if ( irem(lcoeff, p) != 0 ) {
+ if ( !zerop(rem(lc, p)) ) {
R = find_modint_ring(p);
- umodpoly modpoly;
- umodpoly_from_ex(modpoly, prim, x, R);
+ umodpoly_from_upoly(modpoly, prim, R);
if ( squarefree(modpoly) ) break;
}
}
// do modular factorization
- umodpoly modpoly;
- umodpoly_from_ex(modpoly, prim, x, R);
upvec trialfactors;
factor_modular(modpoly, trialfactors);
if ( trialfactors.size() <= 1 ) {
if ( minfactors == 0 || trialfactors.size() < minfactors ) {
factors = trialfactors;
- minfactors = factors.size();
+ minfactors = trialfactors.size();
lastp = p;
trials = 1;
}
}
p = lastp;
R = find_modint_ring(p);
- cl_univpoly_modint_ring UPR = find_univpoly_ring(R);
// lift all factor combinations
stack<ModFactors> tocheck;
mf.poly = prim;
mf.factors = factors;
tocheck.push(mf);
+ upoly f1, f2;
ex result = 1;
while ( tocheck.size() ) {
const size_t n = tocheck.top().factors.size();
- Partition part(n);
+ factor_partition part(tocheck.top().factors);
while ( true ) {
- umodpoly a, b;
- split(tocheck.top().factors, part, a, b);
-
- ex answer = hensel_univar(tocheck.top().poly, x, p, a, b);
- if ( answer != lst() ) {
- if ( part.size_first() == 1 ) {
- if ( part.size_second() == 1 ) {
- result *= answer.op(0) * answer.op(1);
+ hensel_univar(tocheck.top().poly, p, part.left(), part.right(), f1, f2);
+ if ( !f1.empty() ) {
+ if ( part.size_left() == 1 ) {
+ if ( part.size_right() == 1 ) {
+ result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
tocheck.pop();
break;
}
- result *= answer.op(0);
- tocheck.top().poly = answer.op(1);
+ result *= upoly_to_ex(f1, x);
+ tocheck.top().poly = f2;
for ( size_t i=0; i<n; ++i ) {
if ( part[i] == 0 ) {
tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
}
break;
}
- else if ( part.size_second() == 1 ) {
- if ( part.size_first() == 1 ) {
- result *= answer.op(0) * answer.op(1);
+ else if ( part.size_right() == 1 ) {
+ if ( part.size_left() == 1 ) {
+ result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
tocheck.pop();
break;
}
- result *= answer.op(1);
- tocheck.top().poly = answer.op(0);
+ result *= upoly_to_ex(f2, x);
+ tocheck.top().poly = f1;
for ( size_t i=0; i<n; ++i ) {
if ( part[i] == 1 ) {
tocheck.top().factors.erase(tocheck.top().factors.begin()+i);
break;
}
else {
- upvec newfactors1(part.size_first()), newfactors2(part.size_second());
+ upvec newfactors1(part.size_left()), newfactors2(part.size_right());
upvec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
for ( size_t i=0; i<n; ++i ) {
if ( part[i] ) {
}
}
tocheck.top().factors = newfactors1;
- tocheck.top().poly = answer.op(0);
+ tocheck.top().poly = f1;
ModFactors mf;
mf.factors = newfactors2;
- mf.poly = answer.op(1);
+ mf.poly = f2;
tocheck.push(mf);
break;
}
}
else {
if ( !part.next() ) {
- result *= tocheck.top().poly;
+ result *= upoly_to_ex(tocheck.top().poly, x);
tocheck.pop();
break;
}
upvec s;
for ( size_t j=1; j<r; ++j ) {
vector<ex> mdarg(2);
- mdarg[0] = umod_to_ex(q[j-1], x);
- mdarg[1] = umod_to_ex(a[j-1], x);
+ mdarg[0] = umodpoly_to_ex(q[j-1], x);
+ mdarg[1] = umodpoly_to_ex(a[j-1], x);
vector<EvalPoint> empty;
- vector<ex> exsigma = multivar_diophant(mdarg, x, umod_to_ex(beta, x), empty, 0, p, k);
+ vector<ex> exsigma = multivar_diophant(mdarg, x, umodpoly_to_ex(beta, x), empty, 0, p, k);
umodpoly sigma1;
umodpoly_from_ex(sigma1, exsigma[0], x, R);
umodpoly sigma2;
umodpoly bmod = b;
change_modulus(R, bmod);
- umodpoly g;
umodpoly smod;
umodpoly tmod;
- exteuclid(amod, bmod, g, smod, tmod);
- if ( unequal_one(g) ) {
- throw logic_error("gcd(amod,bmod) != 1");
- }
+ exteuclid(amod, bmod, smod, tmod);
cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
umodpoly s = smod;
if ( r > 2 ) {
upvec s = multiterm_eea_lift(a, x, p, k);
for ( size_t j=0; j<r; ++j ) {
- ex phi = expand(pow(x,m) * umod_to_ex(s[j], x));
- umodpoly bmod;
- umodpoly_from_ex(bmod, phi, x, R);
+ umodpoly bmod = umodpoly_to_umodpoly(s[j], R, m);
umodpoly buf;
rem(bmod, a[j], buf);
result.push_back(buf);
}
}
else {
- umodpoly s;
- umodpoly t;
+ umodpoly s, t;
eea_lift(a[1], a[0], x, p, k, s, t);
- ex phi = expand(pow(x,m) * umod_to_ex(s, x));
- umodpoly bmod;
- umodpoly_from_ex(bmod, phi, x, R);
+ umodpoly bmod = umodpoly_to_umodpoly(s, R, m);
umodpoly buf, q;
remdiv(bmod, a[0], buf, q);
result.push_back(buf);
- phi = expand(pow(x,m) * umod_to_ex(t, x));
- umodpoly t1mod;
- umodpoly_from_ex(t1mod, phi, x, R);
- umodpoly buf2 = t1mod + q * a[1];
- result.push_back(buf2);
+ umodpoly t1mod = umodpoly_to_umodpoly(t, R, m);
+ buf = t1mod + q * a[1];
+ result.push_back(buf);
}
return result;
ex e = make_modular(buf, R);
ex monomial = 1;
- for ( size_t m=1; m<=d; ++m ) {
- while ( !e.is_zero() && e.has(xnu) ) {
- monomial *= (xnu - alphanu);
- monomial = expand(monomial);
- ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
- cm = make_modular(cm, R);
- if ( !cm.is_zero() ) {
- vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
- ex buf = e;
- for ( size_t j=0; j<delta_s.size(); ++j ) {
- delta_s[j] *= monomial;
- sigma[j] += delta_s[j];
- buf -= delta_s[j] * b[j];
- }
- e = make_modular(buf, R);
+ for ( size_t m=1; !e.is_zero() && e.has(xnu) && m<=d; ++m ) {
+ monomial *= (xnu - alphanu);
+ monomial = expand(monomial);
+ ex cm = e.diff(ex_to<symbol>(xnu), m).subs(xnu==alphanu) / factorial(m);
+ cm = make_modular(cm, R);
+ if ( !cm.is_zero() ) {
+ vector<ex> delta_s = multivar_diophant(anew, x, cm, Inew, d, p, k);
+ ex buf = e;
+ for ( size_t j=0; j<delta_s.size(); ++j ) {
+ delta_s[j] *= monomial;
+ sigma[j] += delta_s[j];
+ buf -= delta_s[j] * b[j];
}
+ e = make_modular(buf, R);
}
}
}
modcm = cl_MI(R, poscm);
for ( size_t j=0; j<delta_s.size(); ++j ) {
delta_s[j] = delta_s[j] * modcm;
- sigma[j] = sigma[j] + umod_to_ex(delta_s[j], x);
+ sigma[j] = sigma[j] + umodpoly_to_ex(delta_s[j], x);
}
if ( nterms > 1 ) {
z = c.op(i+1);
const size_t n = u.size();
vector<ex> U(n);
for ( size_t i=0; i<n; ++i ) {
- U[i] = umod_to_ex(u[i], x);
+ U[i] = umodpoly_to_ex(u[i], x);
}
for ( size_t j=2; j<=nu; ++j ) {
epv.push_back(ep);
}
- // calc bound B
- ex maxcoeff;
- for ( int i=u.degree(x); i>=u.ldegree(x); --i ) {
- maxcoeff += pow(abs(u.coeff(x, i)),2);
- }
- cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(ex_to<numeric>(maxcoeff).to_cl_N())));
- unsigned int maxdegree = 0;
+ // calc bound p^l
+ int maxdeg = 0;
for ( size_t i=0; i<factor_count; ++i ) {
- if ( ufaclst[2*i+1].degree(x) > (int)maxdegree ) {
- maxdegree = ufaclst[2*i+1].degree(x);
+ if ( ufaclst[2*i+1].degree(x) > maxdeg ) {
+ maxdeg = ufaclst[2*i+1].degree(x);
}
}
- cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
+ cl_I B = 2*calc_bound(u, x, maxdeg);
cl_I l = 1;
cl_I pl = prime;
while ( pl < B ) {
git://www.ginac.de / ginac.git / blobdiff