+ return;
+ }
+ }
+ }
+ if ( ++r == k ) {
+ r = 1;
+ ++u;
+ }
+ }
+}
+
+static void expt_1_over_p(const umodpoly& a, unsigned int prime, umodpoly& ap)
+{
+ 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];
+ }
+}
+
+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;
+ }
+ }
+ }
+ else {
+ 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;
+ }
+ }
+}
+
+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 nhalf = degree(a)/2;
+
+ int i = 1;
+ umodpoly w(2);
+ w[0] = R->zero();
+ w[1] = R->one();
+ umodpoly x = w;
+
+ bool nontrivial = false;
+ while ( i <= nhalf ) {
+ expt_pos(w, q);
+ umodpoly buf;
+ 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;
+ }
+ if ( unequal_one(a) ) {
+ degrees.push_back(degree(a));
+ ddfactors.push_back(a);
+ }
+}
+
+static void same_degree_factor(const umodpoly& a, upvec& upv)
+{
+ cl_modint_ring R = a[0].ring();
+ int deg = degree(a);
+
+ vector<int> degrees;
+ upvec ddfactors;
+ distinct_degree_factor(a, degrees, ddfactors);
+
+ for ( size_t i=0; i<degrees.size(); ++i ) {
+ if ( degrees[i] == degree(ddfactors[i]) ) {
+ upv.push_back(ddfactors[i]);
+ }
+ else {
+ berlekamp(ddfactors[i], upv);
+ }
+ }
+}
+
+static void factor_modular(const umodpoly& p, upvec& upv)
+{
+ upvec factors;
+ vector<int> mult;
+ modsqrfree(p, factors, mult);
+
+#define USE_SAME_DEGREE_FACTOR
+#ifdef USE_SAME_DEGREE_FACTOR
+ for ( size_t i=0; i<factors.size(); ++i ) {
+ upvec upvbuf;
+ same_degree_factor(factors[i], upvbuf);
+ for ( int j=mult[i]; j>0; --j ) {
+ upv.insert(upv.end(), upvbuf.begin(), upvbuf.end());
+ }
+ }
+#else
+ for ( size_t i=0; i<factors.size(); ++i ) {
+ upvec upvbuf;
+ berlekamp(factors[i], upvbuf);
+ if ( upvbuf.size() ) {
+ for ( size_t j=0; j<upvbuf.size(); ++j ) {
+ upv.insert(upv.end(), mult[i], upvbuf[j]);
+ }
+ }
+ else {
+ for ( int j=mult[i]; j>0; --j ) {
+ upv.push_back(factors[i]);
+ }
+ }
+ }
+#endif
+}
+
+/** 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, t, s);
+ return;
+ }
+
+ umodpoly one(1, a[0].ring()->one());
+ umodpoly c = a; normalize_in_field(c);
+ umodpoly d = b; normalize_in_field(d);
+ s = one;
+ t.clear();
+ umodpoly d1;
+ umodpoly d2 = one;
+ umodpoly q;
+ while ( true ) {
+ div(c, d, q);
+ umodpoly r = c - q * d;
+ umodpoly r1 = s - q * d1;
+ umodpoly r2 = t - q * d2;
+ c = d;
+ s = d1;
+ t = d2;
+ if ( r.empty() ) break;
+ d = r;
+ d1 = r1;
+ d2 = r2;
+ }
+ 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);
+ fac = recip(lcoeff(b) * lcoeff(c));
+ i = t.begin(), end = t.end();
+ for ( ; i!=end; ++i ) {
+ *i = *i * fac;
+ }
+ canonicalize(t);
+}
+
+static upoly replace_lc(const upoly& poly, const cl_I& lc)
+{
+ if ( poly.empty() ) return poly;
+ upoly r = poly;
+ r.back() = lc;
+ return r;
+}
+
+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
+ cl_R maxcoeff = 0;
+ for ( int i=degree(a); i>=0; --i ) {
+ maxcoeff = maxcoeff + square(abs(a[i]));
+ }
+ cl_I normmc = ceiling1(the<cl_R>(cln::sqrt(maxcoeff)));
+ cl_I maxdegree = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
+ cl_I B = normmc * expt_pos(cl_I(2), maxdegree);
+
+ // step 1
+ cl_I alpha = lcoeff(a);
+ a = a * alpha;
+ umodpoly nu1 = u1_;
+ normalize_in_field(nu1);
+ umodpoly nw1 = w1_;
+ normalize_in_field(nw1);
+ upoly phi;
+ phi = umodpoly_to_upoly(nu1) * alpha;
+ umodpoly u1;
+ umodpoly_from_upoly(u1, phi, R);
+ phi = umodpoly_to_upoly(nw1) * alpha;
+ umodpoly w1;
+ umodpoly_from_upoly(w1, phi, R);
+
+ // step 2
+ umodpoly s;
+ umodpoly t;
+ exteuclid(u1, w1, s, t);
+
+ // step 3
+ 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;
+ const cl_I maxmodulus = 2*B*abs(alpha);
+
+ // step 4
+ while ( !e.empty() && modulus < maxmodulus ) {
+ // ad-hoc divisablity check
+ for ( size_t k=0; k<e.size(); ++k ) {
+ if ( !zerop(mod(e[k], modulus)) ) {
+ goto quickexit;
+ }
+ }
+ upoly c = e / modulus;
+ phi = umodpoly_to_upoly(s) * c;
+ umodpoly sigmatilde;
+ umodpoly_from_upoly(sigmatilde, phi, R);
+ phi = umodpoly_to_upoly(t) * c;
+ umodpoly tautilde;
+ umodpoly_from_upoly(tautilde, phi, R);
+ umodpoly r, q;
+ remdiv(sigmatilde, w1, r, q);
+ umodpoly sigma = r;
+ phi = umodpoly_to_upoly(tautilde) + umodpoly_to_upoly(q) * umodpoly_to_upoly(u1);
+ umodpoly tau;
+ 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;
+ }
+quickexit: ;
+
+ // step 5
+ 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 {
+ u.clear();
+ }
+}
+
+static unsigned int next_prime(unsigned int p)
+{
+ static vector<unsigned int> primes;
+ if ( primes.size() == 0 ) {
+ primes.push_back(3); primes.push_back(5); primes.push_back(7);
+ }
+ vector<unsigned int>::const_iterator it = primes.begin();
+ if ( p >= primes.back() ) {
+ unsigned int candidate = primes.back() + 2;
+ while ( true ) {
+ size_t n = primes.size()/2;
+ for ( size_t i=0; i<n; ++i ) {
+ if ( candidate % primes[i] ) continue;
+ candidate += 2;
+ i=-1;
+ }
+ primes.push_back(candidate);
+ if ( candidate > p ) break;
+ }
+ return candidate;
+ }
+ vector<unsigned int>::const_iterator end = primes.end();
+ for ( ; it!=end; ++it ) {
+ if ( *it > p ) {
+ return *it;
+ }
+ }
+ throw logic_error("next_prime: should not reach this point!");
+}
+
+class factor_partition
+{
+public:
+ factor_partition(const upvec& factors_) : factors(factors_)
+ {
+ n = factors.size();
+ k.resize(n, 1);
+ k[0] = 0;
+ sum = n-1;
+ one.resize(1, factors.front()[0].ring()->one());
+ 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; }
+#ifdef DEBUGFACTOR
+ void get() const { DCOUTVAR(k); }
+#endif
+ bool next()
+ {
+ for ( size_t i=n-1; i>=1; --i ) {
+ if ( k[i] ) {
+ --k[i];
+ --sum;
+ if ( sum > 0 ) {
+ split();
+ return true;
+ }
+ else {
+ return false;
+ }
+ }
+ ++k[i];
+ ++sum;
+ }
+ return false;
+ }
+ void split()
+ {
+ left = one;
+ right = one;
+ for ( size_t i=0; i<n; ++i ) {
+ if ( k[i] ) {
+ right = right * factors[i];
+ }
+ else {
+ left = left * factors[i];
+ }
+ }
+ }
+public:
+ umodpoly left, right;
+private:
+ upvec factors;
+ umodpoly one;
+ size_t n, sum;
+ vector<int> k;
+};
+
+struct ModFactors
+{
+ upoly poly;
+ upvec factors;
+};
+
+static ex factor_univariate(const ex& poly, const ex& x)
+{
+ 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;
+ cl_I lc = lcoeff(prim);
+ upvec factors;
+ while ( trials < 2 ) {
+ umodpoly modpoly;
+ while ( true ) {
+ p = next_prime(p);
+ if ( !zerop(rem(lc, p)) ) {
+ R = find_modint_ring(p);
+ umodpoly_from_upoly(modpoly, prim, R);
+ if ( squarefree(modpoly) ) break;
+ }
+ }
+
+ // do modular factorization
+ upvec trialfactors;
+ factor_modular(modpoly, trialfactors);
+ if ( trialfactors.size() <= 1 ) {
+ // irreducible for sure
+ return poly;
+ }
+
+ if ( minfactors == 0 || trialfactors.size() < minfactors ) {
+ factors = trialfactors;
+ minfactors = factors.size();
+ lastp = p;
+ trials = 1;
+ }
+ else {
+ ++trials;
+ }
+ }
+ p = lastp;
+ R = find_modint_ring(p);
+
+ // lift all factor combinations
+ stack<ModFactors> tocheck;
+ ModFactors mf;
+ 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();
+ factor_partition part(tocheck.top().factors);
+ while ( true ) {
+ hensel_univar(tocheck.top().poly, p, part.left, part.right, f1, f2);
+ if ( !f1.empty() ) {
+ if ( part.size_first() == 1 ) {
+ if ( part.size_second() == 1 ) {
+ result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
+ tocheck.pop();
+ break;
+ }
+ 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;
+ }
+ }
+ break;
+ }
+ else if ( part.size_second() == 1 ) {
+ if ( part.size_first() == 1 ) {
+ result *= upoly_to_ex(f1, x) * upoly_to_ex(f2, x);
+ tocheck.pop();
+ break;
+ }
+ 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;
+ }
+ }
+ break;
+ }
+ else {
+ upvec newfactors1(part.size_first()), newfactors2(part.size_second());
+ upvec::iterator i1 = newfactors1.begin(), i2 = newfactors2.begin();
+ for ( size_t i=0; i<n; ++i ) {
+ if ( part[i] ) {
+ *i2++ = tocheck.top().factors[i];
+ }
+ else {
+ *i1++ = tocheck.top().factors[i];
+ }
+ }
+ tocheck.top().factors = newfactors1;
+ tocheck.top().poly = f1;
+ ModFactors mf;
+ mf.factors = newfactors2;
+ mf.poly = f2;
+ tocheck.push(mf);
+ break;
+ }
+ }
+ else {
+ if ( !part.next() ) {
+ result *= upoly_to_ex(tocheck.top().poly, x);
+ tocheck.pop();
+ break;
+ }
+ }
+ }
+ }
+
+ return unit * cont * result;
+}
+
+struct EvalPoint
+{
+ ex x;
+ int evalpoint;
+};
+
+// forward declaration
+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);
+
+upvec multiterm_eea_lift(const upvec& a, const ex& x, unsigned int p, unsigned int k)
+{
+ const size_t r = a.size();
+ cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
+ upvec q(r-1);
+ q[r-2] = a[r-1];
+ for ( size_t j=r-2; j>=1; --j ) {
+ q[j-1] = a[j] * q[j];
+ }
+ umodpoly beta(1, R->one());
+ upvec s;
+ for ( size_t j=1; j<r; ++j ) {
+ vector<ex> mdarg(2);
+ 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, umodpoly_to_ex(beta, x), empty, 0, p, k);
+ umodpoly sigma1;
+ umodpoly_from_ex(sigma1, exsigma[0], x, R);
+ umodpoly sigma2;
+ umodpoly_from_ex(sigma2, exsigma[1], x, R);
+ beta = sigma1;
+ s.push_back(sigma2);
+ }
+ s.push_back(beta);
+ return s;
+}
+
+/**
+ * Assert: a not empty.
+ */
+void change_modulus(const cl_modint_ring& R, umodpoly& a)
+{
+ if ( a.empty() ) return;
+ cl_modint_ring oldR = a[0].ring();
+ umodpoly::iterator i = a.begin(), end = a.end();
+ for ( ; i!=end; ++i ) {
+ *i = R->canonhom(oldR->retract(*i));
+ }
+ canonicalize(a);
+}
+
+void eea_lift(const umodpoly& a, const umodpoly& b, const ex& x, unsigned int p, unsigned int k, umodpoly& s_, umodpoly& t_)
+{
+ cl_modint_ring R = find_modint_ring(p);
+ umodpoly amod = a;
+ change_modulus(R, amod);
+ umodpoly bmod = b;
+ change_modulus(R, bmod);
+
+ umodpoly smod;
+ umodpoly tmod;
+ exteuclid(amod, bmod, smod, tmod);
+
+ cl_modint_ring Rpk = find_modint_ring(expt_pos(cl_I(p),k));
+ umodpoly s = smod;
+ change_modulus(Rpk, s);
+ umodpoly t = tmod;
+ change_modulus(Rpk, t);
+
+ cl_I modulus(p);
+ umodpoly one(1, Rpk->one());
+ for ( size_t j=1; j<k; ++j ) {
+ umodpoly e = one - a * s - b * t;
+ reduce_coeff(e, modulus);
+ umodpoly c = e;
+ change_modulus(R, c);
+ umodpoly sigmabar = smod * c;
+ umodpoly taubar = tmod * c;
+ umodpoly sigma, q;
+ remdiv(sigmabar, bmod, sigma, q);
+ umodpoly tau = taubar + q * amod;
+ umodpoly sadd = sigma;
+ change_modulus(Rpk, sadd);
+ cl_MI modmodulus(Rpk, modulus);
+ s = s + sadd * modmodulus;
+ umodpoly tadd = tau;
+ change_modulus(Rpk, tadd);
+ t = t + tadd * modmodulus;
+ modulus = modulus * p;
+ }
+
+ s_ = s; t_ = t;
+}
+
+upvec univar_diophant(const upvec& a, const ex& x, unsigned int m, unsigned int p, unsigned int k)
+{
+ cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
+
+ const size_t r = a.size();
+ upvec result;
+ if ( r > 2 ) {
+ upvec s = multiterm_eea_lift(a, x, p, k);
+ for ( size_t j=0; j<r; ++j ) {
+ umodpoly bmod = umodpoly_to_umodpoly(s[j], R, m);
+ umodpoly buf;
+ rem(bmod, a[j], buf);
+ result.push_back(buf);
+ }
+ }
+ else {
+ umodpoly s, t;
+ eea_lift(a[1], a[0], x, p, k, s, t);
+ umodpoly bmod = umodpoly_to_umodpoly(s, R, m);
+ umodpoly buf, q;
+ remdiv(bmod, a[0], buf, q);
+ result.push_back(buf);
+ umodpoly t1mod = umodpoly_to_umodpoly(t, R, m);
+ buf = t1mod + q * a[1];
+ result.push_back(buf);
+ }
+
+ return result;
+}
+
+struct make_modular_map : public map_function {
+ cl_modint_ring R;
+ make_modular_map(const cl_modint_ring& R_) : R(R_) { }
+ ex operator()(const ex& e)
+ {
+ if ( is_a<add>(e) || is_a<mul>(e) ) {
+ return e.map(*this);
+ }
+ else if ( is_a<numeric>(e) ) {
+ numeric mod(R->modulus);
+ numeric halfmod = (mod-1)/2;
+ cl_MI emod = R->canonhom(the<cl_I>(ex_to<numeric>(e).to_cl_N()));
+ numeric n(R->retract(emod));
+ if ( n > halfmod ) {
+ return n-mod;
+ }
+ else {
+ return n;
+ }
+ }
+ return e;
+ }
+};
+
+static ex make_modular(const ex& e, const cl_modint_ring& R)
+{
+ make_modular_map map(R);
+ return map(e.expand());
+}
+
+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)
+{
+ vector<ex> a = a_;
+
+ const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),k));
+ const size_t r = a.size();
+ const size_t nu = I.size() + 1;
+
+ vector<ex> sigma;
+ if ( nu > 1 ) {
+ ex xnu = I.back().x;
+ int alphanu = I.back().evalpoint;
+
+ ex A = 1;
+ for ( size_t i=0; i<r; ++i ) {
+ A *= a[i];
+ }
+ vector<ex> b(r);
+ for ( size_t i=0; i<r; ++i ) {
+ b[i] = normal(A / a[i]);
+ }
+
+ vector<ex> anew = a;
+ for ( size_t i=0; i<r; ++i ) {
+ anew[i] = anew[i].subs(xnu == alphanu);
+ }
+ ex cnew = c.subs(xnu == alphanu);
+ vector<EvalPoint> Inew = I;
+ Inew.pop_back();
+ sigma = multivar_diophant(anew, x, cnew, Inew, d, p, k);
+
+ ex buf = c;
+ for ( size_t i=0; i<r; ++i ) {
+ buf -= sigma[i] * b[i];
+ }
+ 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);
git://www.ginac.de / ginac.git / blobdiff