+static inline ex factor_univariate(const ex& poly, const ex& x)
+{
+ unsigned int prime;
+ return factor_univariate(poly, x, prime);
+}
+
+struct EvalPoint
+{
+ ex x;
+ int evalpoint;
+};
+
+// 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);
+
+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;
+}
+
+/**
+ * Assert: a not empty.
+ */
+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);
+}
+
+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;
+}
+
+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;
+}
+
+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());
+}
+
+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)
+{
+ 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; !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);
+ }
+ }
+ }
+ else {
+ upvec amod;
+ for ( size_t i=0; i<a.size(); ++i ) {
+ umodpoly up;
+ umodpoly_from_ex(up, a[i], x, R);
+ amod.push_back(up);
+ }
+
+ sigma.insert(sigma.begin(), r, 0);
+ size_t nterms;
+ ex z;
+ if ( is_a<add>(c) ) {
+ nterms = c.nops();
+ z = c.op(0);
+ }
+ else {
+ nterms = 1;
+ z = c;
+ }
+ for ( size_t i=0; i<nterms; ++i ) {
+ int m = z.degree(x);
+ cl_I cm = the<cl_I>(ex_to<numeric>(z.lcoeff(x)).to_cl_N());
+ upvec delta_s = univar_diophant(amod, x, m, p, k);
+ cl_MI modcm;
+ cl_I poscm = cm;
+ while ( poscm < 0 ) {
+ poscm = poscm + expt_pos(cl_I(p),k);
+ }
+ 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] + umodpoly_to_ex(delta_s[j], x);
+ }
+ if ( nterms > 1 ) {
+ z = c.op(i+1);
+ }
+ }
+ }
+
+ for ( size_t i=0; i<sigma.size(); ++i ) {
+ sigma[i] = make_modular(sigma[i], R);
+ }
+
+ return sigma;
+}
+
+#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
+
+static ex hensel_multivar(const ex& a, const ex& x, const vector<EvalPoint>& I, unsigned int p, const cl_I& l, const upvec& u, const vector<ex>& lcU)
+{
+ const size_t nu = I.size() + 1;
+ const cl_modint_ring R = find_modint_ring(expt_pos(cl_I(p),l));
+
+ vector<ex> A(nu);
+ A[nu-1] = a;
+
+ for ( size_t j=nu; j>=2; --j ) {
+ ex x = I[j-2].x;
+ int alpha = I[j-2].evalpoint;
+ A[j-2] = A[j-1].subs(x==alpha);
+ A[j-2] = make_modular(A[j-2], R);
+ }
+
+ int maxdeg = a.degree(I.front().x);
+ for ( size_t i=1; i<I.size(); ++i ) {
+ int maxdeg2 = a.degree(I[i].x);
+ if ( maxdeg2 > maxdeg ) maxdeg = maxdeg2;
+ }
+
+ const size_t n = u.size();
+ vector<ex> U(n);
+ for ( size_t i=0; i<n; ++i ) {
+ U[i] = umodpoly_to_ex(u[i], x);
+ }
+
+ for ( size_t j=2; j<=nu; ++j ) {
+ vector<ex> U1 = U;
+ ex monomial = 1;
+ for ( size_t m=0; m<n; ++m) {
+ if ( lcU[m] != 1 ) {
+ ex coef = lcU[m];
+ for ( size_t i=j-1; i<nu-1; ++i ) {
+ coef = coef.subs(I[i].x == I[i].evalpoint);
+ }
+ coef = make_modular(coef, R);
+ int deg = U[m].degree(x);
+ U[m] = U[m] - U[m].lcoeff(x) * pow(x,deg) + coef * pow(x,deg);
+ }
+ }
+ ex Uprod = 1;
+ for ( size_t i=0; i<n; ++i ) {
+ Uprod *= U[i];
+ }
+ ex e = expand(A[j-1] - Uprod);
+
+ vector<EvalPoint> newI;
+ for ( size_t i=1; i<=j-2; ++i ) {
+ newI.push_back(I[i-1]);
+ }
+
+ ex xj = I[j-2].x;
+ int alphaj = I[j-2].evalpoint;
+ size_t deg = A[j-1].degree(xj);
+ for ( size_t k=1; k<=deg; ++k ) {
+ if ( !e.is_zero() ) {
+ monomial *= (xj - alphaj);
+ monomial = expand(monomial);
+ ex dif = e.diff(ex_to<symbol>(xj), k);
+ ex c = dif.subs(xj==alphaj) / factorial(k);
+ if ( !c.is_zero() ) {
+ vector<ex> deltaU = multivar_diophant(U1, x, c, newI, maxdeg, p, cl_I_to_uint(l));
+ for ( size_t i=0; i<n; ++i ) {
+ deltaU[i] *= monomial;
+ U[i] += deltaU[i];
+ U[i] = make_modular(U[i], R);
+ }
+ ex Uprod = 1;
+ for ( size_t i=0; i<n; ++i ) {
+ Uprod *= U[i];
+ }
+ e = A[j-1] - Uprod;
+ e = make_modular(e, R);
+ }
+ }
+ }
+ }
+
+ ex acand = 1;
+ for ( size_t i=0; i<U.size(); ++i ) {
+ acand *= U[i];
+ }
+ if ( expand(a-acand).is_zero() ) {
+ lst res;
+ for ( size_t i=0; i<U.size(); ++i ) {
+ res.append(U[i]);
+ }
+ return res;
+ }
+ else {
+ lst res;
+ return lst();
+ }
+}
+
+static ex put_factors_into_lst(const ex& e)
+{
+ lst result;
+
+ if ( is_a<numeric>(e) ) {
+ result.append(e);
+ return result;
+ }
+ if ( is_a<power>(e) ) {
+ result.append(1);
+ result.append(e.op(0));
+ return result;
+ }
+ if ( is_a<symbol>(e) || is_a<add>(e) ) {
+ result.append(1);
+ result.append(e);
+ return result;
+ }
+ if ( is_a<mul>(e) ) {
+ ex nfac = 1;
+ for ( size_t i=0; i<e.nops(); ++i ) {
+ ex op = e.op(i);
+ if ( is_a<numeric>(op) ) {
+ nfac = op;
+ }
+ if ( is_a<power>(op) ) {
+ result.append(op.op(0));
+ }
+ if ( is_a<symbol>(op) || is_a<add>(op) ) {
+ result.append(op);
+ }
+ }
+ result.prepend(nfac);
+ return result;
+ }
+ throw runtime_error("put_factors_into_lst: bad term.");
+}
+
+#ifdef DEBUGFACTOR
+ostream& operator<<(ostream& o, const vector<numeric>& v)
+{
+ for ( size_t i=0; i<v.size(); ++i ) {
+ o << v[i] << " ";
+ }
+ return o;
+}
+#endif // def DEBUGFACTOR
+
+/** Checks whether in a set of numbers each has a unique prime factor.
+ *
+ * @param[in] f list of numbers to check
+ * @return true: if number set is bad, false: otherwise
+ */
+static bool checkdivisors(const lst& f)
+{
+ const int k = f.nops();
+ numeric q, r;
+ vector<numeric> d(k);
+ d[0] = ex_to<numeric>(abs(f.op(0)));
+ for ( int i=1; i<k; ++i ) {
+ q = ex_to<numeric>(abs(f.op(i)));
+ for ( int j=i-1; j>=0; --j ) {
+ r = d[j];
+ do {
+ r = gcd(r, q);
+ q = q/r;
+ } while ( r != 1 );
+ if ( q == 1 ) {
+ return true;
+ }
+ }
+ d[i] = q;
+ }
+ return false;
+}
+
+/** Generates a set of evaluation points for a multivariate polynomial.
+ * The set fulfills the following conditions:
+ * 1. lcoeff(evaluated_polynomial) does not vanish
+ * 2. factors of lcoeff(evaluated_polynomial) have each a unique prime factor
+ * 3. evaluated_polynomial is square free
+ * See [W1] for more details.
+ *
+ * @param[in] u multivariate polynomial to be factored
+ * @param[in] vn leading coefficient of u in x (x==first symbol in syms)
+ * @param[in] syms set of symbols that appear in u
+ * @param[in] f lst containing the factors of the leading coefficient vn
+ * @param[in,out] modulus integer modulus for random number generation (i.e. |a_i| < modulus)
+ * @param[out] u0 returns the evaluated (univariate) polynomial
+ * @param[out] a returns the valid evaluation points. must have initial size equal
+ * number of symbols-1 before calling generate_set
+ */
+static void generate_set(const ex& u, const ex& vn, const exset& syms, const lst& f,
+ numeric& modulus, ex& u0, vector<numeric>& a)
+{
+ const ex& x = *syms.begin();
+ while ( true ) {
+ ++modulus;
+ /* generate a set of integers ... */
+ u0 = u;
+ ex vna = vn;
+ ex vnatry;
+ exset::const_iterator s = syms.begin();
+ ++s;
+ for ( size_t i=0; i<a.size(); ++i ) {
+ do {
+ a[i] = mod(numeric(rand()), 2*modulus) - modulus;
+ vnatry = vna.subs(*s == a[i]);
+ /* ... for which the leading coefficient doesn't vanish ... */
+ } while ( vnatry == 0 );
+ vna = vnatry;
+ u0 = u0.subs(*s == a[i]);
+ ++s;
+ }
+ /* ... for which u0 is square free ... */
+ ex g = gcd(u0, u0.diff(ex_to<symbol>(x)));
+ if ( !is_a<numeric>(g) ) {
+ continue;
+ }
+ if ( !is_a<numeric>(vn) ) {
+ /* ... and for which the evaluated factors have each an unique prime factor */
+ lst fnum = f;
+ fnum.let_op(0) = fnum.op(0) * u0.content(x);
+ for ( size_t i=1; i<fnum.nops(); ++i ) {
+ if ( !is_a<numeric>(fnum.op(i)) ) {
+ s = syms.begin();
+ ++s;
+ for ( size_t j=0; j<a.size(); ++j, ++s ) {
+ fnum.let_op(i) = fnum.op(i).subs(*s == a[j]);
+ }
+ }
+ }
+ if ( checkdivisors(fnum) ) {
+ continue;
+ }
+ }
+ /* ok, we have a valid set now */
+ return;
+ }
+}
+
+// forward declaration
+static ex factor_sqrfree(const ex& poly);
+
+/**
+ * ASSERT: poly is expanded
+ */
+static ex factor_multivariate(const ex& poly, const exset& syms)
+{
+ exset::const_iterator s;
+ const ex& x = *syms.begin();
+
+ /* make polynomial primitive */
+ ex p = poly.collect(x);
+ ex cont = p.lcoeff(x);
+ for ( int i=p.degree(x)-1; i>=p.ldegree(x); --i ) {
+ cont = gcd(cont, p.coeff(x,i));
+ if ( cont == 1 ) break;
+ }
+ ex pp = expand(normal(p / cont));
+ if ( !is_a<numeric>(cont) ) {
+ return factor_sqrfree(cont) * factor_sqrfree(pp);
+ }
+
+ /* factor leading coefficient */
+ ex vn = pp.collect(x).lcoeff(x);
+ ex vnlst;
+ if ( is_a<numeric>(vn) ) {
+ vnlst = lst(vn);
+ }
+ else {
+ ex vnfactors = factor(vn);
+ vnlst = put_factors_into_lst(vnfactors);
+ }
+
+ const unsigned int maxtrials = 3;
+ numeric modulus = (vnlst.nops() > 3) ? vnlst.nops() : 3;
+ vector<numeric> a(syms.size()-1, 0);
+
+ /* try now to factorize until we are successful */
+ while ( true ) {
+
+ unsigned int trialcount = 0;
+ unsigned int prime;
+ int factor_count = 0;
+ int min_factor_count = -1;
+ ex u, delta;
+ ex ufac, ufaclst;
+
+ /* try several evaluation points to reduce the number of modular factors */
+ while ( trialcount < maxtrials ) {
+
+ /* generate a set of valid evaluation points */
+ generate_set(pp, vn, syms, ex_to<lst>(vnlst), modulus, u, a);
+
+ ufac = factor_univariate(u, x, prime);
+ ufaclst = put_factors_into_lst(ufac);
+ factor_count = ufaclst.nops()-1;
+ delta = ufaclst.op(0);
+
+ if ( factor_count <= 1 ) {
+ /* irreducible */
+ return poly;
+ }
+ if ( min_factor_count < 0 ) {
+ /* first time here */
+ min_factor_count = factor_count;
+ }
+ else if ( min_factor_count == factor_count ) {
+ /* one less to try */
+ ++trialcount;
+ }
+ else if ( min_factor_count > factor_count ) {
+ /* new minimum, reset trial counter */
+ min_factor_count = factor_count;
+ trialcount = 0;
+ }
+ }
+
+ // determine true leading coefficients for the Hensel lifting
+ vector<ex> C(factor_count);
+ if ( is_a<numeric>(vn) ) {
+ for ( size_t i=1; i<ufaclst.nops(); ++i ) {
+ C[i-1] = ufaclst.op(i).lcoeff(x);
+ }
+ }
+ else {
+ vector<numeric> ftilde(vnlst.nops()-1);
+ for ( size_t i=0; i<ftilde.size(); ++i ) {
+ ex ft = vnlst.op(i+1);
+ s = syms.begin();
+ ++s;
+ for ( size_t j=0; j<a.size(); ++j ) {
+ ft = ft.subs(*s == a[j]);
+ ++s;
+ }
+ ftilde[i] = ex_to<numeric>(ft);
+ }
+
+ vector<bool> used_flag(ftilde.size(), false);
+ vector<ex> D(factor_count, 1);
+ if ( delta == 1 ) {
+ for ( int i=0; i<factor_count; ++i ) {
+ numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
+ for ( int j=ftilde.size()-1; j>=0; --j ) {
+ int count = 0;
+ while ( irem(prefac, ftilde[j]) == 0 ) {
+ prefac = iquo(prefac, ftilde[j]);
+ ++count;
+ }
+ if ( count ) {
+ used_flag[j] = true;
+ D[i] = D[i] * pow(vnlst.op(j+1), count);
+ }
+ }
+ C[i] = D[i] * prefac;
+ }
+ }
+ else {
+ for ( int i=0; i<factor_count; ++i ) {
+ numeric prefac = ex_to<numeric>(ufaclst.op(i+1).lcoeff(x));
+ for ( int j=ftilde.size()-1; j>=0; --j ) {
+ int count = 0;
+ while ( irem(prefac, ftilde[j]) == 0 ) {
+ prefac = iquo(prefac, ftilde[j]);
+ ++count;
+ }
+ while ( irem(ex_to<numeric>(delta)*prefac, ftilde[j]) == 0 ) {
+ numeric g = gcd(prefac, ex_to<numeric>(ftilde[j]));
+ prefac = iquo(prefac, g);
+ delta = delta / (ftilde[j]/g);
+ ufaclst.let_op(i+1) = ufaclst.op(i+1) * (ftilde[j]/g);
+ ++count;
+ }
+ if ( count ) {
+ used_flag[j] = true;
+ D[i] = D[i] * pow(vnlst.op(j+1), count);
+ }
+ }
+ C[i] = D[i] * prefac;
+ }
+ }
+
+ bool some_factor_unused = false;
+ for ( size_t i=0; i<used_flag.size(); ++i ) {
+ if ( !used_flag[i] ) {
+ some_factor_unused = true;
+ break;
+ }
+ }
+ if ( some_factor_unused ) {
+ continue;
+ }
+ }
+
+ if ( delta != 1 ) {
+ C[0] = C[0] * delta;
+ ufaclst.let_op(1) = ufaclst.op(1) * delta;
+ }
+
+ EvalPoint ep;
+ vector<EvalPoint> epv;
+ s = syms.begin();
+ ++s;
+ for ( size_t i=0; i<a.size(); ++i ) {
+ ep.x = *s++;
+ ep.evalpoint = a[i].to_int();
+ epv.push_back(ep);
+ }
+
+ // calc bound p^l
+ int maxdeg = 0;
+ for ( int i=1; i<=factor_count; ++i ) {
+ if ( ufaclst.op(i).degree(x) > maxdeg ) {
+ maxdeg = ufaclst[i].degree(x);
+ }
+ }
+ cl_I B = 2*calc_bound(u, x, maxdeg);
+ cl_I l = 1;
+ cl_I pl = prime;
+ while ( pl < B ) {
+ l = l + 1;
+ pl = pl * prime;
+ }
+ cl_modint_ring R = find_modint_ring(expt_pos(cl_I(prime),l));
+ upvec modfactors(ufaclst.nops()-1);
+ for ( size_t i=1; i<ufaclst.nops(); ++i ) {
+ umodpoly_from_ex(modfactors[i-1], ufaclst.op(i), x, R);
+ }
+
+ ex res = hensel_multivar(pp, x, epv, prime, l, modfactors, C);
+ if ( res != lst() ) {
+ ex result = cont;
+ for ( size_t i=0; i<res.nops(); ++i ) {
+ result *= res.op(i).content(x) * res.op(i).unit(x);
+ result *= res.op(i).primpart(x);
+ }
+ return result;
+ }
+ }
+}
+
+struct find_symbols_map : public map_function {