+ 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;
+ }
+ }