+/** Berlekamp's modular factorization.
+ *
+ * The implementation follows the algorithm in chapter 8 of [GCL].
+ *
+ * @param[in] a modular polynomial
+ * @param[out] upv vector containing modular factors. if upv was not empty the
+ * new elements are added at the end
+ */
+static void berlekamp(const umodpoly& a, upvec& upv)
+{
+ cl_modint_ring R = a[0].ring();
+ umodpoly one(1, R->one());
+
+ // find nullspace of Q matrix
+ modular_matrix Q(degree(a), degree(a), R->zero());
+ q_matrix(a, Q);
+ vector<mvec> nu;
+ nullspace(Q, nu);
+
+ const unsigned int k = nu.size();
+ if ( k == 1 ) {
+ // irreducible
+ return;
+ }
+
+ list<umodpoly> factors;
+ factors.push_back(a);
+ unsigned int size = 1;
+ unsigned int r = 1;
+ unsigned int q = cl_I_to_uint(R->modulus);
+
+ list<umodpoly>::iterator u = factors.begin();
+
+ // calculate all gcd's
+ while ( true ) {
+ for ( unsigned int s=0; s<q; ++s ) {
+ umodpoly nur = nu[r];
+ nur[0] = nur[0] - cl_MI(R, s);
+ canonicalize(nur);
+ umodpoly g;
+ gcd(nur, *u, g);
+ if ( unequal_one(g) && g != *u ) {
+ umodpoly uo;
+ div(*u, g, uo);
+ if ( equal_one(uo) ) {
+ throw logic_error("berlekamp: unexpected divisor.");
+ }
+ else {
+ *u = uo;
+ }
+ factors.push_back(g);
+ size = 0;
+ list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
+ while ( i != end ) {
+ if ( degree(*i) ) ++size;
+ ++i;
+ }
+ if ( size == k ) {
+ list<umodpoly>::const_iterator i = factors.begin(), end = factors.end();
+ while ( i != end ) {
+ upv.push_back(*i++);
+ }
+ return;
+ }
+ }
+ }
+ if ( ++r == k ) {
+ r = 1;
+ ++u;
+ }
+ }
+}
+
+// modular square free factorization is not used at the moment so we deactivate
+// the code
+#if 0
+
+/** Calculates a^(1/prime).
+ *
+ * @param[in] a polynomial
+ * @param[in] prime prime number -> exponent 1/prime
+ * @param[in] ap resulting polynomial
+ */
+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];
+ }
+}
+
+/** Modular square free factorization.
+ *
+ * @param[in] a polynomial
+ * @param[out] factors modular factors
+ * @param[out] mult corresponding multiplicities (exponents)
+ */
+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;
+ }
+ }
+}
+
+#endif // deactivation of square free factorization
+
+/** Distinct degree factorization (DDF).
+ *
+ * The implementation follows the algorithm in chapter 8 of [GCL].
+ *
+ * @param[in] a_ modular polynomial
+ * @param[out] degrees vector containing the degrees of the factors of the
+ * corresponding polynomials in ddfactors.
+ * @param[out] ddfactors vector containing polynomials which factors have the
+ * degree given in degrees.
+ */
+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;
+
+ 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);
+ }
+}
+
+/** Modular same degree factorization.
+ * Same degree factorization is a kind of misnomer. It performs distinct degree
+ * factorization, but instead of using the Cantor-Zassenhaus algorithm it
+ * (sub-optimally) uses Berlekamp's algorithm for the factors of the same
+ * degree.
+ *
+ * @param[in] a modular polynomial
+ * @param[out] upv vector containing modular factors. if upv was not empty the
+ * new elements are added at the end
+ */
+static void same_degree_factor(const umodpoly& a, upvec& upv)
+{
+ cl_modint_ring R = a[0].ring();
+
+ 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);
+ }
+ }
+}
+
+// Yes, we can (choose).
+#define USE_SAME_DEGREE_FACTOR
+
+/** Modular univariate factorization.
+ *
+ * In principle, we have two algorithms at our disposal: Berlekamp's algorithm
+ * and same degree factorization (SDF). SDF seems to be slightly faster in
+ * almost all cases so it is activated as default.
+ *
+ * @param[in] p modular polynomial
+ * @param[out] upv vector containing modular factors. if upv was not empty the
+ * new elements are added at the end
+ */
+static void factor_modular(const umodpoly& p, upvec& upv)
+{
+#ifdef USE_SAME_DEGREE_FACTOR
+ same_degree_factor(p, upv);
+#else
+ berlekamp(p, upv);
+#endif
+}
+
+/** Calculates modular 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);
+}
+
+/** Replaces the leading coefficient in a polynomial by a given number.
+ *
+ * @param[in] poly polynomial to change
+ * @param[in] lc new leading coefficient
+ * @return changed polynomial
+ */
+static upoly replace_lc(const upoly& poly, const cl_I& lc)
+{
+ if ( poly.empty() ) return poly;
+ upoly r = poly;
+ r.back() = lc;
+ return r;
+}
+
+/** Calculates the bound for the modulus.
+ * See [Mig].
+ */
+static inline cl_I calc_bound(const ex& a, const ex& x, int maxdeg)
+{
+ 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;
+}
+
+/** Calculates the bound for the modulus.
+ * See [Mig].
+ */
+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;
+}
+
+/** Hensel lifting as used by factor_univariate().
+ *
+ * The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ * @param[in] a_ primitive univariate polynomials
+ * @param[in] p prime number that does not divide lcoeff(a)
+ * @param[in] u1_ modular factor of a (mod p)
+ * @param[in] w1_ modular factor of a (mod p), relatively prime to u1_,
+ * fulfilling u1_*w1_ == a mod p
+ * @param[out] u lifted factor
+ * @param[out] w lifted factor, u*w = a
+ */
+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
+ int maxdeg = (degree(u1_) > degree(w1_)) ? degree(u1_) : degree(w1_);
+ cl_I maxmodulus = 2*calc_bound(a, maxdeg);
+
+ // 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;
+
+ // step 4
+ while ( !e.empty() && modulus < maxmodulus ) {
+ 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;
+ }
+
+ // 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();
+ }
+}
+
+/** Returns a new prime number.
+ *
+ * @param[in] p prime number
+ * @return next prime number after p
+ */
+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!");
+}
+
+/** Manages the splitting a vector of of modular factors into two partitions.
+ */
+class factor_partition
+{
+public:
+ /** Takes the vector of modular factors and initializes the first partition */
+ factor_partition(const upvec& factors_) : factors(factors_)
+ {
+ 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_left() const { return n-len; }
+ size_t size_right() const { return len; }
+ /** Initializes the next partition.
+ Returns true, if there is one, false otherwise. */
+ bool next()
+ {
+ 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;
+ }
+ last = len;
+ ++len;
+ if ( len > n/2 ) return false;
+ fill(k.begin(), k.begin()+len, 1);
+ fill(k.begin()+len+1, k.end(), 0);
+ }
+ else {
+ k[last++] = 0;
+ k[last] = 1;
+ }
+ split();
+ return true;
+ }
+ /** Get first partition */
+ umodpoly& left() { return lr[0]; }
+ /** Get second partition */
+ umodpoly& right() { return lr[1]; }
+private:
+ 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 {
+ 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;
+};
+
+/** Contains a pair of univariate polynomial and its modular factors.
+ * Used by factor_univariate().
+ */
+struct ModFactors
+{
+ upoly poly;
+ upvec factors;
+};
+
+/** Univariate polynomial factorization.
+ *
+ * Modular factorization is tried for several primes to minimize the number of
+ * modular factors. Then, Hensel lifting is performed.
+ *
+ * @param[in] poly expanded square free univariate polynomial
+ * @param[in] x symbol
+ * @param[in,out] prime prime number to start trying modular factorization with,
+ * output value is the prime number actually used
+ */
+static ex factor_univariate(const ex& poly, const ex& x, unsigned int& prime)
+{
+ 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
+ prime = 3;
+ unsigned int lastp = prime;
+ cl_modint_ring R;
+ unsigned int trials = 0;
+ unsigned int minfactors = 0;
+
+ const numeric& cont_n = ex_to<numeric>(cont);
+ cl_I i_cont;
+ if (cont_n.is_integer()) {
+ i_cont = the<cl_I>(cont_n.to_cl_N());
+ } else {
+ // poly \in Q[x] => poly = q ipoly, ipoly \in Z[x], q \in Q
+ // factor(poly) \equiv q factor(ipoly)
+ i_cont = cl_I(1);
+ }
+ cl_I lc = lcoeff(prim)*i_cont;
+ upvec factors;
+ while ( trials < 2 ) {
+ umodpoly modpoly;
+ while ( true ) {
+ prime = next_prime(prime);
+ if ( !zerop(rem(lc, prime)) ) {
+ R = find_modint_ring(prime);
+ 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 = trialfactors.size();
+ lastp = prime;
+ trials = 1;
+ }
+ else {
+ ++trials;
+ }
+ }
+ prime = lastp;
+ R = find_modint_ring(prime);
+
+ // 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 ) {
+ // call Hensel lifting
+ hensel_univar(tocheck.top().poly, prime, part.left(), part.right(), f1, f2);
+ if ( !f1.empty() ) {
+ // successful, update the stack and the result
+ 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 *= 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_right() == 1 ) {
+ if ( part.size_left() == 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_left()), newfactors2(part.size_right());
+ 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 {
+ // not successful
+ if ( !part.next() ) {
+ // if no more combinations left, return polynomial as
+ // irreducible
+ result *= upoly_to_ex(tocheck.top().poly, x);
+ tocheck.pop();
+ break;
+ }
+ }
+ }
+ }
+
+ return unit * cont * result;
+}
+
+/** Second interface to factor_univariate() to be used if the information about
+ * the prime is not needed.
+ */
+static inline ex factor_univariate(const ex& poly, const ex& x)
+{
+ unsigned int prime;
+ return factor_univariate(poly, x, prime);
+}
+
+/** Represents an evaluation point (<symbol>==<integer>).
+ */
+struct EvalPoint
+{
+ ex x;
+ int evalpoint;
+};
+
+#ifdef DEBUGFACTOR
+ostream& operator<<(ostream& o, const vector<EvalPoint>& v)
+{
+ for ( size_t i=0; i<v.size(); ++i ) {
+ o << "(" << v[i].x << "==" << v[i].evalpoint << ") ";
+ }
+ return o;
+}
+#endif // def DEBUGFACTOR
+
+// forward declaration
+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);
+
+/** Utility function for multivariate Hensel lifting.
+ *
+ * Solves the equation
+ * s_1*b_1 + ... + s_r*b_r == 1 mod p^k
+ * with deg(s_i) < deg(a_i)
+ * and with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
+ *
+ * The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ * @param[in] a vector of modular univariate polynomials
+ * @param[in] x symbol
+ * @param[in] p prime number
+ * @param[in] k p^k is modulus
+ * @return vector of polynomials (s_i)
+ */
+static 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;
+}
+
+/** Changes the modulus of a modular polynomial. Used by eea_lift().
+ *
+ * @param[in] R new modular ring
+ * @param[in,out] a polynomial to change (in situ)
+ */
+static 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);
+}
+
+/** Utility function for multivariate Hensel lifting.
+ *
+ * Solves s*a + t*b == 1 mod p^k given a,b.
+ *
+ * The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ * @param[in] a polynomial
+ * @param[in] b polynomial
+ * @param[in] x symbol
+ * @param[in] p prime number
+ * @param[in] k p^k is modulus
+ * @param[out] s_ output polynomial
+ * @param[out] t_ output polynomial
+ */
+static 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;
+}
+
+/** Utility function for multivariate Hensel lifting.
+ *
+ * Solves the equation
+ * s_1*b_1 + ... + s_r*b_r == x^m mod p^k
+ * with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
+ *
+ * The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ * @param a vector with univariate polynomials mod p^k
+ * @param x symbol
+ * @param m exponent of x^m in the equation to solve
+ * @param p prime number
+ * @param k p^k is modulus
+ * @return vector of polynomials (s_i)
+ */
+static 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;
+}
+
+/** Map used by function make_modular().
+ * Finds every coefficient in a polynomial and replaces it by is value in the
+ * given modular ring R (symmetric representation).
+ */
+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;
+ }
+};
+
+/** Helps mimicking modular multivariate polynomial arithmetic.
+ *
+ * @param e expression of which to make the coefficients equal to their value
+ * in the modular ring R (symmetric representation)
+ * @param R modular ring
+ * @return resulting expression
+ */
+static ex make_modular(const ex& e, const cl_modint_ring& R)
+{
+ make_modular_map map(R);
+ return map(e.expand());
+}
+
+/** Utility function for multivariate Hensel lifting.
+ *
+ * Returns the polynomials s_i that fulfill
+ * s_1*b_1 + ... + s_r*b_r == c mod <I^(d+1),p^k>
+ * with given b_1 = a_1 * ... * a_{i-1} * a_{i+1} * ... * a_r
+ *
+ * The implementation follows the algorithm in chapter 6 of [GCL].
+ *
+ * @param a_ vector of multivariate factors mod p^k
+ * @param x symbol (equiv. x_1 in [GCL])
+ * @param c polynomial mod p^k
+ * @param I vector of evaluation points
+ * @param d maximum total degree of result
+ * @param p prime number
+ * @param k p^k is modulus
+ * @return vector of polynomials (s_i)
+ */
+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)