+ return inherited::eval_ncmul(v);
+}
+
+bool tryfactsubs(const ex & origfactor, const ex & patternfactor, int & nummatches, lst & repls)
+{
+ ex origbase;
+ int origexponent;
+ int origexpsign;
+
+ if (is_exactly_a<power>(origfactor) && origfactor.op(1).info(info_flags::integer)) {
+ origbase = origfactor.op(0);
+ int expon = ex_to<numeric>(origfactor.op(1)).to_int();
+ origexponent = expon > 0 ? expon : -expon;
+ origexpsign = expon > 0 ? 1 : -1;
+ } else {
+ origbase = origfactor;
+ origexponent = 1;
+ origexpsign = 1;
+ }
+
+ ex patternbase;
+ int patternexponent;
+ int patternexpsign;
+
+ if (is_exactly_a<power>(patternfactor) && patternfactor.op(1).info(info_flags::integer)) {
+ patternbase = patternfactor.op(0);
+ int expon = ex_to<numeric>(patternfactor.op(1)).to_int();
+ patternexponent = expon > 0 ? expon : -expon;
+ patternexpsign = expon > 0 ? 1 : -1;
+ } else {
+ patternbase = patternfactor;
+ patternexponent = 1;
+ patternexpsign = 1;
+ }
+
+ lst saverepls = repls;
+ if (origexponent < patternexponent || origexpsign != patternexpsign || !origbase.match(patternbase,saverepls))
+ return false;
+ repls = saverepls;
+
+ int newnummatches = origexponent / patternexponent;
+ if (newnummatches < nummatches)
+ nummatches = newnummatches;
+ return true;
+}
+
+/** Checks wheter e matches to the pattern pat and the (possibly to be updated
+ * list of replacements repls. This matching is in the sense of algebraic
+ * substitutions. Matching starts with pat.op(factor) of the pattern because
+ * the factors before this one have already been matched. The (possibly
+ * updated) number of matches is in nummatches. subsed[i] is true for factors
+ * that already have been replaced by previous substitutions and matched[i]
+ * is true for factors that have been matched by the current match.
+ */
+bool algebraic_match_mul_with_mul(const mul &e, const ex &pat, lst &repls,
+ int factor, int &nummatches, const std::vector<bool> &subsed,
+ std::vector<bool> &matched)
+{
+ if (factor == pat.nops())
+ return true;
+
+ for (size_t i=0; i<e.nops(); ++i) {
+ if(subsed[i] || matched[i])
+ continue;
+ lst newrepls = repls;
+ int newnummatches = nummatches;
+ if (tryfactsubs(e.op(i), pat.op(factor), newnummatches, newrepls)) {
+ matched[i] = true;
+ if (algebraic_match_mul_with_mul(e, pat, newrepls, factor+1,
+ newnummatches, subsed, matched)) {
+ repls = newrepls;
+ nummatches = newnummatches;
+ return true;
+ }
+ else
+ matched[i] = false;
+ }
+ }
+
+ return false;
+}
+
+bool mul::has(const ex & pattern, unsigned options) const
+{
+ if(!(options&has_options::algebraic))
+ return basic::has(pattern,options);
+ if(is_a<mul>(pattern)) {
+ lst repls;
+ int nummatches = std::numeric_limits<int>::max();
+ std::vector<bool> subsed(seq.size(), false);
+ std::vector<bool> matched(seq.size(), false);
+ if(algebraic_match_mul_with_mul(*this, pattern, repls, 0, nummatches,
+ subsed, matched))
+ return true;
+ }
+ return basic::has(pattern, options);
+}
+
+ex mul::algebraic_subs_mul(const exmap & m, unsigned options) const
+{
+ std::vector<bool> subsed(seq.size(), false);
+ exvector subsresult(seq.size());
+
+ for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
+
+ if (is_exactly_a<mul>(it->first)) {
+retry1:
+ int nummatches = std::numeric_limits<int>::max();
+ std::vector<bool> currsubsed(seq.size(), false);
+ bool succeed = true;
+ lst repls;
+
+ if(!algebraic_match_mul_with_mul(*this, it->first, repls, 0, nummatches, subsed, currsubsed))
+ continue;
+
+ bool foundfirstsubsedfactor = false;
+ for (size_t j=0; j<subsed.size(); j++) {
+ if (currsubsed[j]) {
+ if (foundfirstsubsedfactor)
+ subsresult[j] = op(j);
+ else {
+ foundfirstsubsedfactor = true;
+ subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches);
+ }
+ subsed[j] = true;
+ }
+ }
+ goto retry1;
+
+ } else {
+retry2:
+ int nummatches = std::numeric_limits<int>::max();
+ lst repls;
+
+ for (size_t j=0; j<this->nops(); j++) {
+ if (!subsed[j] && tryfactsubs(op(j), it->first, nummatches, repls)) {
+ subsed[j] = true;
+ subsresult[j] = op(j) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches);
+ goto retry2;
+ }
+ }
+ }
+ }
+
+ bool subsfound = false;
+ for (size_t i=0; i<subsed.size(); i++) {
+ if (subsed[i]) {
+ subsfound = true;
+ break;
+ }
+ }
+ if (!subsfound)
+ return subs_one_level(m, options | subs_options::algebraic);
+
+ exvector ev; ev.reserve(nops());
+ for (size_t i=0; i<nops(); i++) {
+ if (subsed[i])
+ ev.push_back(subsresult[i]);
+ else
+ ev.push_back(op(i));
+ }
+
+ return (new mul(ev))->setflag(status_flags::dynallocated);