X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fnormal.cpp;h=9f8b7b4a406fdb5fe5ee39454e86eb6340f0d86a;hp=c6dc39e6b0ef1d4108a5dcfeb4628054f3b5a9f4;hb=4ffb3cbb3ab5f642e461bcbf8fb29743752c5d58;hpb=faf8c206bdcd01eee758103d56b83a634f787e7a diff --git a/ginac/normal.cpp b/ginac/normal.cpp index c6dc39e6..9f8b7b4a 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -6,7 +6,7 @@ * computation, square-free factorization and rational function normalization. */ /* - * GiNaC Copyright (C) 1999-2015 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2018 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -146,7 +146,7 @@ struct sym_desc { /** Maximum number of terms of leading coefficient of symbol in both polynomials */ size_t max_lcnops; - /** Commparison operator for sorting */ + /** Comparison operator for sorting */ bool operator<(const sym_desc &x) const { if (max_deg == x.max_deg) @@ -196,8 +196,8 @@ static void collect_symbols(const ex &e, sym_desc_vec &v) * @param v vector of sym_desc structs (filled in) */ static void get_symbol_stats(const ex &a, const ex &b, sym_desc_vec &v) { - collect_symbols(a.eval(), v); // eval() to expand assigned symbols - collect_symbols(b.eval(), v); + collect_symbols(a, v); + collect_symbols(b, v); for (auto & it : v) { int deg_a = a.degree(it.sym); int deg_b = b.degree(it.sym); @@ -270,9 +270,15 @@ static numeric lcm_of_coefficients_denominators(const ex &e) * @param lcm LCM to multiply in */ static ex multiply_lcm(const ex &e, const numeric &lcm) { + if (lcm.is_equal(*_num1_p)) + // e * 1 -> e; + return e; + if (is_exactly_a(e)) { + // (a*b*...)*lcm -> (a*lcma)*(b*lcmb)*...*(lcm/(lcma*lcmb*...)) size_t num = e.nops(); - exvector v; v.reserve(num + 1); + exvector v; + v.reserve(num + 1); numeric lcm_accum = *_num1_p; for (size_t i=0; isetflag(status_flags::dynallocated); + return dynallocate(v); } else if (is_exactly_a(e)) { + // (a+b+...)*lcm -> a*lcm+b*lcm+... size_t num = e.nops(); - exvector v; v.reserve(num); + exvector v; + v.reserve(num); for (size_t i=0; isetflag(status_flags::dynallocated); + return dynallocate(v); } else if (is_exactly_a(e)) { - if (is_a(e.op(0))) - return e * lcm; - else - return pow(multiply_lcm(e.op(0), lcm.power(ex_to(e.op(1)).inverse())), e.op(1)); - } else - return e * lcm; + if (!is_a(e.op(0))) { + // (b^e)*lcm -> (b*lcm^(1/e))^e if lcm^(1/e) ∈ ℚ (i.e. not a float) + // but not for symbolic b, as evaluation would undo this again + numeric root_of_lcm = lcm.power(ex_to(e.op(1)).inverse()); + if (root_of_lcm.is_rational()) + return pow(multiply_lcm(e.op(0), root_of_lcm), e.op(1)); + } + } + // can't recurse down into e + return dynallocate(e, lcm); } @@ -386,16 +398,16 @@ ex quo(const ex &a, const ex &b, const ex &x, bool check_args) term = rcoeff / blcoeff; else { if (!divide(rcoeff, blcoeff, term, false)) - return (new fail())->setflag(status_flags::dynallocated); + return dynallocate(); } - term *= power(x, rdeg - bdeg); + term *= pow(x, rdeg - bdeg); v.push_back(term); r -= (term * b).expand(); if (r.is_zero()) break; rdeg = r.degree(x); } - return (new add(v))->setflag(status_flags::dynallocated); + return dynallocate(v); } @@ -439,9 +451,9 @@ ex rem(const ex &a, const ex &b, const ex &x, bool check_args) term = rcoeff / blcoeff; else { if (!divide(rcoeff, blcoeff, term, false)) - return (new fail())->setflag(status_flags::dynallocated); + return dynallocate(); } - term *= power(x, rdeg - bdeg); + term *= pow(x, rdeg - bdeg); r -= (term * b).expand(); if (r.is_zero()) break; @@ -501,23 +513,23 @@ ex prem(const ex &a, const ex &b, const ex &x, bool check_args) if (bdeg == 0) eb = _ex0; else - eb -= blcoeff * power(x, bdeg); + eb -= blcoeff * pow(x, bdeg); } else blcoeff = _ex1; int delta = rdeg - bdeg + 1, i = 0; while (rdeg >= bdeg && !r.is_zero()) { ex rlcoeff = r.coeff(x, rdeg); - ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand(); + ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand(); if (rdeg == 0) r = _ex0; else - r -= rlcoeff * power(x, rdeg); + r -= rlcoeff * pow(x, rdeg); r = (blcoeff * r).expand() - term; rdeg = r.degree(x); i++; } - return power(blcoeff, delta - i) * r; + return pow(blcoeff, delta - i) * r; } @@ -553,17 +565,17 @@ ex sprem(const ex &a, const ex &b, const ex &x, bool check_args) if (bdeg == 0) eb = _ex0; else - eb -= blcoeff * power(x, bdeg); + eb -= blcoeff * pow(x, bdeg); } else blcoeff = _ex1; while (rdeg >= bdeg && !r.is_zero()) { ex rlcoeff = r.coeff(x, rdeg); - ex term = (power(x, rdeg - bdeg) * eb * rlcoeff).expand(); + ex term = (pow(x, rdeg - bdeg) * eb * rlcoeff).expand(); if (rdeg == 0) r = _ex0; else - r -= rlcoeff * power(x, rdeg); + r -= rlcoeff * pow(x, rdeg); r = (blcoeff * r).expand() - term; rdeg = r.degree(x); } @@ -653,7 +665,7 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) else resv.push_back(a.op(j)); } - q = (new mul(resv))->setflag(status_flags::dynallocated); + q = dynallocate(resv); return true; } } else if (is_exactly_a(a)) { @@ -663,7 +675,7 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) int a_exp = ex_to(a.op(1)).to_int(); ex rem_i; if (divide(ab, b, rem_i, false)) { - q = rem_i*power(ab, a_exp - 1); + q = rem_i * pow(ab, a_exp - 1); return true; } // code below is commented-out because it leads to a significant slowdown @@ -693,11 +705,11 @@ bool divide(const ex &a, const ex &b, ex &q, bool check_args) else if (!divide(rcoeff, blcoeff, term, false)) return false; - term *= power(x, rdeg - bdeg); + term *= pow(x, rdeg - bdeg); v.push_back(term); r -= (term * b).expand(); if (r.is_zero()) { - q = (new add(v))->setflag(status_flags::dynallocated); + q = dynallocate(v); return true; } rdeg = r.degree(x); @@ -876,11 +888,11 @@ static bool divide_in_z(const ex &a, const ex &b, ex &q, sym_desc_vec::const_ite ex term, rcoeff = r.coeff(x, rdeg); if (!divide_in_z(rcoeff, blcoeff, term, var+1)) break; - term = (term * power(x, rdeg - bdeg)).expand(); + term = (term * pow(x, rdeg - bdeg)).expand(); v.push_back(term); r -= (term * eb).expand(); if (r.is_zero()) { - q = (new add(v))->setflag(status_flags::dynallocated); + q = dynallocate(v); #if USE_REMEMBER dr_remember[ex2(a, b)] = exbool(q, true); #endif @@ -1210,7 +1222,7 @@ ex add::smod(const numeric &xi) const } GINAC_ASSERT(is_exactly_a(overall_coeff)); numeric coeff = GiNaC::smod(ex_to(overall_coeff), xi); - return (new add(newseq,coeff))->setflag(status_flags::dynallocated); + return dynallocate(std::move(newseq), coeff); } ex mul::smod(const numeric &xi) const @@ -1220,12 +1232,12 @@ ex mul::smod(const numeric &xi) const GINAC_ASSERT(!is_exactly_a(recombine_pair_to_ex(it))); } #endif // def DO_GINAC_ASSERT - mul * mulcopyp = new mul(*this); + mul & mulcopy = dynallocate(*this); GINAC_ASSERT(is_exactly_a(overall_coeff)); - mulcopyp->overall_coeff = GiNaC::smod(ex_to(overall_coeff),xi); - mulcopyp->clearflag(status_flags::evaluated); - mulcopyp->clearflag(status_flags::hash_calculated); - return mulcopyp->setflag(status_flags::dynallocated); + mulcopy.overall_coeff = GiNaC::smod(ex_to(overall_coeff),xi); + mulcopy.clearflag(status_flags::evaluated); + mulcopy.clearflag(status_flags::hash_calculated); + return mulcopy; } @@ -1237,10 +1249,10 @@ static ex interpolate(const ex &gamma, const numeric &xi, const ex &x, int degre numeric rxi = xi.inverse(); for (int i=0; !e.is_zero(); i++) { ex gi = e.smod(xi); - g.push_back(gi * power(x, i)); + g.push_back(gi * pow(x, i)); e = (e - gi) * rxi; } - return (new add(g))->setflag(status_flags::dynallocated); + return dynallocate(g); } /** Exception thrown by heur_gcd() to signal failure. */ @@ -1536,7 +1548,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio // The symbol with least degree which is contained in both polynomials // is our main variable - sym_desc_vec::iterator vari = sym_stats.begin(); + auto vari = sym_stats.begin(); while ((vari != sym_stats.end()) && (((vari->ldeg_b == 0) && (vari->deg_b == 0)) || ((vari->ldeg_a == 0) && (vari->deg_a == 0)))) @@ -1563,7 +1575,7 @@ ex gcd(const ex &a, const ex &b, ex *ca, ex *cb, bool check_args, unsigned optio int ldeg_b = var->ldeg_b; int min_ldeg = std::min(ldeg_a,ldeg_b); if (min_ldeg > 0) { - ex common = power(x, min_ldeg); + ex common = pow(x, min_ldeg); return gcd((aex / common).expand(), (bex / common).expand(), ca, cb, false) * common; } @@ -1644,14 +1656,14 @@ static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb) if (ca) *ca = _ex1; if (cb) - *cb = power(p, exp_b - exp_a); - return power(p, exp_a); + *cb = pow(p, exp_b - exp_a); + return pow(p, exp_a); } else { if (ca) - *ca = power(p, exp_a - exp_b); + *ca = pow(p, exp_a - exp_b); if (cb) *cb = _ex1; - return power(p, exp_b); + return pow(p, exp_b); } } @@ -1671,11 +1683,11 @@ static ex gcd_pf_pow_pow(const ex& a, const ex& b, ex* ca, ex* cb) // a(x) = g(x)^n A(x)^n, b(x) = g(x)^m B(x)^m ==> // gcd(a, b) = g(x)^n gcd(A(x)^n, g(x)^(n-m) B(x)^m if (exp_a < exp_b) { - ex pg = gcd(power(p_co, exp_a), power(p_gcd, exp_b-exp_a)*power(pb_co, exp_b), ca, cb, false); - return power(p_gcd, exp_a)*pg; + ex pg = gcd(pow(p_co, exp_a), pow(p_gcd, exp_b-exp_a)*pow(pb_co, exp_b), ca, cb, false); + return pow(p_gcd, exp_a)*pg; } else { - ex pg = gcd(power(p_gcd, exp_a - exp_b)*power(p_co, exp_a), power(pb_co, exp_b), ca, cb, false); - return power(p_gcd, exp_b)*pg; + ex pg = gcd(pow(p_gcd, exp_a - exp_b)*pow(p_co, exp_a), pow(pb_co, exp_b), ca, cb, false); + return pow(p_gcd, exp_b)*pg; } } @@ -1694,7 +1706,7 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb) if (p.is_equal(b)) { // a = p^n, b = p, gcd = p if (ca) - *ca = power(p, a.op(1) - 1); + *ca = pow(p, a.op(1) - 1); if (cb) *cb = _ex1; return p; @@ -1712,7 +1724,7 @@ static ex gcd_pf_pow(const ex& a, const ex& b, ex* ca, ex* cb) return _ex1; } // a(x) = g(x)^n A(x)^n, b(x) = g(x) B(x) ==> gcd(a, b) = g(x) gcd(g(x)^(n-1) A(x)^n, B(x)) - ex rg = gcd(power(p_gcd, exp_a-1)*power(p_co, exp_a), bpart_co, ca, cb, false); + ex rg = gcd(pow(p_gcd, exp_a-1)*pow(p_co, exp_a), bpart_co, ca, cb, false); return p_gcd*rg; } @@ -1737,10 +1749,10 @@ static ex gcd_pf_mul(const ex& a, const ex& b, ex* ca, ex* cb) part_b = part_cb; } if (ca) - *ca = (new mul(acc_ca))->setflag(status_flags::dynallocated); + *ca = dynallocate(acc_ca); if (cb) *cb = part_b; - return (new mul(g))->setflag(status_flags::dynallocated); + return dynallocate(g); } /** Compute LCM (Least Common Multiple) of multivariate polynomials in Z[X]. @@ -1771,34 +1783,47 @@ ex lcm(const ex &a, const ex &b, bool check_args) * Yun's algorithm. Used internally by sqrfree(). * * @param a multivariate polynomial over Z[X], treated here as univariate - * polynomial in x. + * polynomial in x (needs not be expanded). * @param x variable to factor in - * @return vector of factors sorted in ascending degree */ -static exvector sqrfree_yun(const ex &a, const symbol &x) + * @return vector of expairs (factor, exponent), sorted by exponents */ +static epvector sqrfree_yun(const ex &a, const symbol &x) { - exvector res; ex w = a; ex z = w.diff(x); ex g = gcd(w, z); + if (g.is_zero()) { + return epvector{}; + } if (g.is_equal(_ex1)) { - res.push_back(a); - return res; + return epvector{expair(a, _ex1)}; } - ex y; + epvector results; + ex exponent = _ex0; do { w = quo(w, g, x); - y = quo(z, g, x); - z = y - w.diff(x); + if (w.is_zero()) { + return res; + } + z = quo(z, g, x) - w.diff(x); + exponent = exponent + 1; + if (w.is_equal(x)) { + // shortcut for x^n with n ∈ ℕ + exponent += quo(z, w.diff(x), x); + results.push_back(expair(w, exponent)); + break; + } g = gcd(w, z); - res.push_back(g); + if (!g.is_equal(_ex1)) { + results.push_back(expair(g, exponent)); + } } while (!z.is_zero()); - return res; + return results; } /** Compute a square-free factorization of a multivariate polynomial in Q[X]. * - * @param a multivariate polynomial over Q[X] + * @param a multivariate polynomial over Q[X] (needs not be expanded) * @param l lst of variables to factor in, may be left empty for autodetection * @return a square-free factorization of \p a. * @@ -1833,8 +1858,8 @@ static exvector sqrfree_yun(const ex &a, const symbol &x) */ ex sqrfree(const ex &a, const lst &l) { - if (is_exactly_a(a) || // algorithm does not trap a==0 - is_a(a)) // shortcut + if (is_exactly_a(a) || + is_a(a)) // shortcuts return a; // If no lst of variables to factorize in was specified we have to @@ -1860,30 +1885,28 @@ ex sqrfree(const ex &a, const lst &l) const ex tmp = multiply_lcm(a,lcm); // find the factors - exvector factors = sqrfree_yun(tmp, x); + epvector factors = sqrfree_yun(tmp, x); - // construct the next list of symbols with the first element popped - lst newargs = args; - newargs.remove_first(); + // remove symbol x and proceed recursively with the remaining symbols + args.remove_first(); // recurse down the factors in remaining variables - if (newargs.nops()>0) { + if (args.nops()>0) { for (auto & it : factors) - it = sqrfree(it, newargs); + it.rest = sqrfree(it.rest, args); } // Done with recursion, now construct the final result ex result = _ex1; - int p = 1; for (auto & it : factors) - result *= power(it, p++); + result *= pow(it.rest, it.coeff); // Yun's algorithm does not account for constant factors. (For univariate // polynomials it works only in the monic case.) We can correct this by // inserting what has been lost back into the result. For completeness // we'll also have to recurse down that factor in the remaining variables. - if (newargs.nops()>0) - result *= sqrfree(quo(tmp, result, x), newargs); + if (args.nops()>0) + result *= sqrfree(quo(tmp, result, x), args); else result *= quo(tmp, result, x); @@ -1911,24 +1934,21 @@ ex sqrfree_parfrac(const ex & a, const symbol & x) //clog << "red_poly = " << red_poly << ", red_numer = " << red_numer << endl; // Factorize denominator and compute cofactors - exvector yun = sqrfree_yun(denom, x); -//clog << "yun factors: " << exprseq(yun) << endl; - size_t num_yun = yun.size(); - exvector factor; factor.reserve(num_yun); - exvector cofac; cofac.reserve(num_yun); - for (size_t i=0; i(yun.back().coeff).to_int(); + exvector factor, cofac; + for (size_t i=0; i(yun[i].coeff); + for (size_t j=0; jsetflag(status_flags::dynallocated); + ex es = dynallocate(); repl.insert(std::make_pair(es, e_replaced)); rev_lookup.insert(std::make_pair(e_replaced, es)); return es; @@ -2016,7 +2036,7 @@ static ex replace_with_symbol(const ex & e, exmap & repl) // Otherwise create new symbol and add to list, taking care that the // replacement expression doesn't itself contain symbols from repl, // because subs() is not recursive - ex es = (new symbol)->setflag(status_flags::dynallocated); + ex es = dynallocate(); repl.insert(std::make_pair(es, e_replaced)); return es; } @@ -2024,36 +2044,27 @@ static ex replace_with_symbol(const ex & e, exmap & repl) /** Function object to be applied by basic::normal(). */ struct normal_map_function : public map_function { - int level; - normal_map_function(int l) : level(l) {} - ex operator()(const ex & e) { return normal(e, level); } + ex operator()(const ex & e) override { return normal(e); } }; /** Default implementation of ex::normal(). It normalizes the children and * replaces the object with a temporary symbol. * @see ex::normal */ -ex basic::normal(exmap & repl, exmap & rev_lookup, int level) const +ex basic::normal(exmap & repl, exmap & rev_lookup) const { if (nops() == 0) - return (new lst(replace_with_symbol(*this, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); - else { - if (level == 1) - return (new lst(replace_with_symbol(*this, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); - else if (level == -max_recursion_level) - throw(std::runtime_error("max recursion level reached")); - else { - normal_map_function map_normal(level - 1); - return (new lst(replace_with_symbol(map(map_normal), repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); - } - } + return dynallocate({replace_with_symbol(*this, repl, rev_lookup), _ex1}); + + normal_map_function map_normal; + return dynallocate({replace_with_symbol(map(map_normal), repl, rev_lookup), _ex1}); } /** Implementation of ex::normal() for symbols. This returns the unmodified symbol. * @see ex::normal */ -ex symbol::normal(exmap & repl, exmap & rev_lookup, int level) const +ex symbol::normal(exmap & repl, exmap & rev_lookup) const { - return (new lst(*this, _ex1))->setflag(status_flags::dynallocated); + return dynallocate({*this, _ex1}); } @@ -2061,7 +2072,7 @@ ex symbol::normal(exmap & repl, exmap & rev_lookup, int level) const * into re+I*im and replaces I and non-rational real numbers with a temporary * symbol. * @see ex::normal */ -ex numeric::normal(exmap & repl, exmap & rev_lookup, int level) const +ex numeric::normal(exmap & repl, exmap & rev_lookup) const { numeric num = numer(); ex numex = num; @@ -2077,7 +2088,7 @@ ex numeric::normal(exmap & repl, exmap & rev_lookup, int level) const } // Denominator is always a real integer (see numeric::denom()) - return (new lst(numex, denom()))->setflag(status_flags::dynallocated); + return dynallocate({numex, denom()}); } @@ -2095,11 +2106,11 @@ static ex frac_cancel(const ex &n, const ex &d) // Handle trivial case where denominator is 1 if (den.is_equal(_ex1)) - return (new lst(num, den))->setflag(status_flags::dynallocated); + return dynallocate({num, den}); // Handle special cases where numerator or denominator is 0 if (num.is_zero()) - return (new lst(num, _ex1))->setflag(status_flags::dynallocated); + return dynallocate({num, _ex1}); if (den.expand().is_zero()) throw(std::overflow_error("frac_cancel: division by zero in frac_cancel")); @@ -2138,30 +2149,25 @@ static ex frac_cancel(const ex &n, const ex &d) // Return result as list //std::clog << " returns num = " << num << ", den = " << den << ", pre_factor = " << pre_factor << std::endl; - return (new lst(num * pre_factor.numer(), den * pre_factor.denom()))->setflag(status_flags::dynallocated); + return dynallocate({num * pre_factor.numer(), den * pre_factor.denom()}); } /** Implementation of ex::normal() for a sum. It expands terms and performs * fractional addition. * @see ex::normal */ -ex add::normal(exmap & repl, exmap & rev_lookup, int level) const +ex add::normal(exmap & repl, exmap & rev_lookup) const { - if (level == 1) - return (new lst(replace_with_symbol(*this, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); - else if (level == -max_recursion_level) - throw(std::runtime_error("max recursion level reached")); - // Normalize children and split each one into numerator and denominator exvector nums, dens; nums.reserve(seq.size()+1); dens.reserve(seq.size()+1); for (auto & it : seq) { - ex n = ex_to(recombine_pair_to_ex(it)).normal(repl, rev_lookup, level-1); + ex n = ex_to(recombine_pair_to_ex(it)).normal(repl, rev_lookup); nums.push_back(n.op(0)); dens.push_back(n.op(1)); } - ex n = ex_to(overall_coeff).normal(repl, rev_lookup, level-1); + ex n = ex_to(overall_coeff).normal(repl, rev_lookup); nums.push_back(n.op(0)); dens.push_back(n.op(1)); GINAC_ASSERT(nums.size() == dens.size()); @@ -2202,29 +2208,23 @@ ex add::normal(exmap & repl, exmap & rev_lookup, int level) const /** Implementation of ex::normal() for a product. It cancels common factors * from fractions. * @see ex::normal() */ -ex mul::normal(exmap & repl, exmap & rev_lookup, int level) const +ex mul::normal(exmap & repl, exmap & rev_lookup) const { - if (level == 1) - return (new lst(replace_with_symbol(*this, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); - else if (level == -max_recursion_level) - throw(std::runtime_error("max recursion level reached")); - // Normalize children, separate into numerator and denominator exvector num; num.reserve(seq.size()); exvector den; den.reserve(seq.size()); ex n; for (auto & it : seq) { - n = ex_to(recombine_pair_to_ex(it)).normal(repl, rev_lookup, level-1); + n = ex_to(recombine_pair_to_ex(it)).normal(repl, rev_lookup); num.push_back(n.op(0)); den.push_back(n.op(1)); } - n = ex_to(overall_coeff).normal(repl, rev_lookup, level-1); + n = ex_to(overall_coeff).normal(repl, rev_lookup); num.push_back(n.op(0)); den.push_back(n.op(1)); // Perform fraction cancellation - return frac_cancel((new mul(num))->setflag(status_flags::dynallocated), - (new mul(den))->setflag(status_flags::dynallocated)); + return frac_cancel(dynallocate(num), dynallocate(den)); } @@ -2232,16 +2232,11 @@ ex mul::normal(exmap & repl, exmap & rev_lookup, int level) const * distributes integer exponents to numerator and denominator, and replaces * non-integer powers by temporary symbols. * @see ex::normal */ -ex power::normal(exmap & repl, exmap & rev_lookup, int level) const +ex power::normal(exmap & repl, exmap & rev_lookup) const { - if (level == 1) - return (new lst(replace_with_symbol(*this, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); - else if (level == -max_recursion_level) - throw(std::runtime_error("max recursion level reached")); - // Normalize basis and exponent (exponent gets reassembled) - ex n_basis = ex_to(basis).normal(repl, rev_lookup, level-1); - ex n_exponent = ex_to(exponent).normal(repl, rev_lookup, level-1); + ex n_basis = ex_to(basis).normal(repl, rev_lookup); + ex n_exponent = ex_to(exponent).normal(repl, rev_lookup); n_exponent = n_exponent.op(0) / n_exponent.op(1); if (n_exponent.info(info_flags::integer)) { @@ -2249,12 +2244,12 @@ ex power::normal(exmap & repl, exmap & rev_lookup, int level) const if (n_exponent.info(info_flags::positive)) { // (a/b)^n -> {a^n, b^n} - return (new lst(power(n_basis.op(0), n_exponent), power(n_basis.op(1), n_exponent)))->setflag(status_flags::dynallocated); + return dynallocate({pow(n_basis.op(0), n_exponent), pow(n_basis.op(1), n_exponent)}); } else if (n_exponent.info(info_flags::negative)) { // (a/b)^-n -> {b^n, a^n} - return (new lst(power(n_basis.op(1), -n_exponent), power(n_basis.op(0), -n_exponent)))->setflag(status_flags::dynallocated); + return dynallocate({pow(n_basis.op(1), -n_exponent), pow(n_basis.op(0), -n_exponent)}); } } else { @@ -2262,32 +2257,32 @@ ex power::normal(exmap & repl, exmap & rev_lookup, int level) const if (n_exponent.info(info_flags::positive)) { // (a/b)^x -> {sym((a/b)^x), 1} - return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); + return dynallocate({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1}); } else if (n_exponent.info(info_flags::negative)) { if (n_basis.op(1).is_equal(_ex1)) { // a^-x -> {1, sym(a^x)} - return (new lst(_ex1, replace_with_symbol(power(n_basis.op(0), -n_exponent), repl, rev_lookup)))->setflag(status_flags::dynallocated); + return dynallocate({_ex1, replace_with_symbol(pow(n_basis.op(0), -n_exponent), repl, rev_lookup)}); } else { // (a/b)^-x -> {sym((b/a)^x), 1} - return (new lst(replace_with_symbol(power(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); + return dynallocate({replace_with_symbol(pow(n_basis.op(1) / n_basis.op(0), -n_exponent), repl, rev_lookup), _ex1}); } } } // (a/b)^x -> {sym((a/b)^x, 1} - return (new lst(replace_with_symbol(power(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); + return dynallocate({replace_with_symbol(pow(n_basis.op(0) / n_basis.op(1), n_exponent), repl, rev_lookup), _ex1}); } /** Implementation of ex::normal() for pseries. It normalizes each coefficient * and replaces the series by a temporary symbol. * @see ex::normal */ -ex pseries::normal(exmap & repl, exmap & rev_lookup, int level) const +ex pseries::normal(exmap & repl, exmap & rev_lookup) const { epvector newseq; for (auto & it : seq) { @@ -2296,7 +2291,7 @@ ex pseries::normal(exmap & repl, exmap & rev_lookup, int level) const newseq.push_back(expair(restexp, it.coeff)); } ex n = pseries(relational(var,point), std::move(newseq)); - return (new lst(replace_with_symbol(n, repl, rev_lookup), _ex1))->setflag(status_flags::dynallocated); + return dynallocate({replace_with_symbol(n, repl, rev_lookup), _ex1}); } @@ -2310,13 +2305,12 @@ ex pseries::normal(exmap & repl, exmap & rev_lookup, int level) const * expression can be treated as a rational function). normal() is applied * recursively to arguments of functions etc. * - * @param level maximum depth of recursion * @return normalized expression */ -ex ex::normal(int level) const +ex ex::normal() const { exmap repl, rev_lookup; - ex e = bp->normal(repl, rev_lookup, level); + ex e = bp->normal(repl, rev_lookup); GINAC_ASSERT(is_a(e)); // Re-insert replaced symbols @@ -2337,7 +2331,7 @@ ex ex::numer() const { exmap repl, rev_lookup; - ex e = bp->normal(repl, rev_lookup, 0); + ex e = bp->normal(repl, rev_lookup); GINAC_ASSERT(is_a(e)); // Re-insert replaced symbols @@ -2357,7 +2351,7 @@ ex ex::denom() const { exmap repl, rev_lookup; - ex e = bp->normal(repl, rev_lookup, 0); + ex e = bp->normal(repl, rev_lookup); GINAC_ASSERT(is_a(e)); // Re-insert replaced symbols @@ -2377,7 +2371,7 @@ ex ex::numer_denom() const { exmap repl, rev_lookup; - ex e = bp->normal(repl, rev_lookup, 0); + ex e = bp->normal(repl, rev_lookup); GINAC_ASSERT(is_a(e)); // Re-insert replaced symbols @@ -2406,47 +2400,11 @@ ex ex::to_rational(exmap & repl) const return bp->to_rational(repl); } -// GiNaC 1.1 compatibility function -ex ex::to_rational(lst & repl_lst) const -{ - // Convert lst to exmap - exmap m; - for (auto & it : repl_lst) - m.insert(std::make_pair(it.op(0), it.op(1))); - - ex ret = bp->to_rational(m); - - // Convert exmap back to lst - repl_lst.remove_all(); - for (auto & it : m) - repl_lst.append(it.first == it.second); - - return ret; -} - ex ex::to_polynomial(exmap & repl) const { return bp->to_polynomial(repl); } -// GiNaC 1.1 compatibility function -ex ex::to_polynomial(lst & repl_lst) const -{ - // Convert lst to exmap - exmap m; - for (auto & it : repl_lst) - m.insert(std::make_pair(it.op(0), it.op(1))); - - ex ret = bp->to_polynomial(m); - - // Convert exmap back to lst - repl_lst.remove_all(); - for (auto & it : m) - repl_lst.append(it.first == it.second); - - return ret; -} - /** Default implementation of ex::to_rational(). This replaces the object with * a temporary symbol. */ ex basic::to_rational(exmap & repl) const @@ -2517,7 +2475,7 @@ ex numeric::to_polynomial(exmap & repl) const ex power::to_rational(exmap & repl) const { if (exponent.info(info_flags::integer)) - return power(basis.to_rational(repl), exponent); + return pow(basis.to_rational(repl), exponent); else return replace_with_symbol(*this, repl); } @@ -2527,17 +2485,17 @@ ex power::to_rational(exmap & repl) const ex power::to_polynomial(exmap & repl) const { if (exponent.info(info_flags::posint)) - return power(basis.to_rational(repl), exponent); + return pow(basis.to_rational(repl), exponent); else if (exponent.info(info_flags::negint)) { ex basis_pref = collect_common_factors(basis); if (is_exactly_a(basis_pref) || is_exactly_a(basis_pref)) { // (A*B)^n will be automagically transformed to A^n*B^n - ex t = power(basis_pref, exponent); + ex t = pow(basis_pref, exponent); return t.to_polynomial(repl); } else - return power(replace_with_symbol(power(basis, _ex_1), repl), -exponent); + return pow(replace_with_symbol(pow(basis, _ex_1), repl), -exponent); } else return replace_with_symbol(*this, repl); @@ -2556,7 +2514,7 @@ ex expairseq::to_rational(exmap & repl) const if (oc.info(info_flags::numeric)) return thisexpairseq(std::move(s), overall_coeff); else - s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1)); + s.push_back(expair(oc, _ex1)); return thisexpairseq(std::move(s), default_overall_coeff()); } @@ -2572,7 +2530,7 @@ ex expairseq::to_polynomial(exmap & repl) const if (oc.info(info_flags::numeric)) return thisexpairseq(std::move(s), overall_coeff); else - s.push_back(combine_ex_with_coeff_to_pair(oc, _ex1)); + s.push_back(expair(oc, _ex1)); return thisexpairseq(std::move(s), default_overall_coeff()); } @@ -2628,7 +2586,7 @@ static ex find_common_factor(const ex & e, ex & factor, exmap & repl) else v.push_back(t.op(k)); } - t = (new mul(v))->setflag(status_flags::dynallocated); + t = dynallocate(v); goto term_done; } } @@ -2638,7 +2596,7 @@ static ex find_common_factor(const ex & e, ex & factor, exmap & repl) t = x; term_done: ; } - return (new add(terms))->setflag(status_flags::dynallocated); + return dynallocate(terms); } else if (is_exactly_a(e)) { @@ -2648,7 +2606,7 @@ term_done: ; for (size_t i=0; isetflag(status_flags::dynallocated); + return dynallocate(v); } else if (is_exactly_a(e)) { const ex e_exp(e.op(1)); @@ -2656,8 +2614,8 @@ term_done: ; ex eb = e.op(0).to_polynomial(repl); ex factor_local(_ex1); ex pre_res = find_common_factor(eb, factor_local, repl); - factor *= power(factor_local, e_exp); - return power(pre_res, e_exp); + factor *= pow(factor_local, e_exp); + return pow(pre_res, e_exp); } else return e.to_polynomial(repl);