>From 18e216b27de09a652626bbf579458add1e5e83bc Mon Sep 17 00:00:00 2001 From: "Vladimir V. Kisil" Date: Thu, 4 Jun 2020 15:45:04 +0100 Subject: [PATCH 3/3] Stronger normalisation method for powers. Now normalisation method try to express powers with the same base as monomials of some dummy variables. This allows to make cancellations of the sort: (x-1)/(sqrt(x)-1) -> sqrt(x)+1 Signed-off-by: Vladimir V. Kisil --- check/exam_normalization.cpp | 56 +++++++++++++++++++++++++++++++ ginac/normal.cpp | 64 ++++++++++++++++++++++++++++++++++++ 2 files changed, 120 insertions(+) diff --git a/check/exam_normalization.cpp b/check/exam_normalization.cpp index 180fe056..70921c86 100644 --- a/check/exam_normalization.cpp +++ b/check/exam_normalization.cpp @@ -255,6 +255,61 @@ static unsigned exam_exponent_law() return result; } +static unsigned exam_power_law() +{ + unsigned result = 0; + ex e, d; + + lst bases = {x, pow(x, numeric(1,3)), exp(x), sin(x)}; // We run all check for power base of different kinds + + for ( auto b : bases ) { + + // simple case + e = 4*b-9; + e /= 2*sqrt(b)-3; + d = 2*sqrt(b)+3; + result += check_normal(e, d); + + // Fractional powers + e = 4*pow(b, numeric(2,3))-9; + e /= 2*pow(b, numeric(1,3))-3; + d = 2*pow(b, numeric(1,3))+3; + result += check_normal(e, d); + + // Different powers with the same base + e = 4*b-9*sqrt(b); + e /= 2*sqrt(b)-3*pow(b, numeric(1,4)); + d = 2*sqrt(b)+3*pow(b, numeric(1,4)); + result += check_normal(e, d); + + // Non-numeric powers + e = 4*pow(b, 2*y)-9; + e /= 2*pow(b, y)-3; + d = 2*pow(b, y)+3; + result += check_normal(e, d); + + // Non-numeric fractional powers + e = 4*pow(b, y)-9; + e /= 2*pow(b, y/2)-3; + d = 2*pow(b, y/2)+3; + result += check_normal(e, d); + + // Different non-numeric powers + e = 4*pow(b, 2*y)-9*pow(b, 2*z); + e /= 2*pow(b, y)-3*pow(b, z); + d = 2*pow(b, y)+3*pow(b, z); + result += check_normal(e, d); + + // Different non-numeric fractional powers + e = 4*pow(b, y)-9*pow(b, z); + e /= 2*pow(b, y/2)-3*pow(b, z/2); + d = 2*pow(b, y/2)+3*pow(b, z/2); + result += check_normal(e, d); + } + + return result; +} + unsigned exam_normalization() { unsigned result = 0; @@ -267,6 +322,7 @@ unsigned exam_normalization() result += exam_normal4(); cout << '.' << flush; result += exam_content(); cout << '.' << flush; result += exam_exponent_law(); cout << '.' << flush; + result += exam_power_law(); cout << '.' << flush; return result; } diff --git a/ginac/normal.cpp b/ginac/normal.cpp index 7a039714..042b658e 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -2059,6 +2059,15 @@ static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, ls if (it != rev_lookup.end()) return it->second; + // The expression can be the base of substituted power, which requires a more careful search + if (! is_a(e_replaced)) + for (auto & it : repl) + if (is_a(it.second) && e_replaced.is_equal(it.second.op(0))) { + ex degree = pow(it.second.op(1), _ex_1); + if (is_a(degree) && ex_to(degree).is_integer()) + return pow(it.first, degree); + } + // We treat powers and the exponent functions differently because // they can be rationalised more efficiently if (is_a(e_replaced) && ex_to(e_replaced).get_name() == "exp") { @@ -2090,6 +2099,53 @@ static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, ls } } } + } else if (is_a(e_replaced) && !is_a(e_replaced.op(0)) // We do not replace simple monomials like x^3 or sqrt(2) + && ! (is_a(e_replaced.op(0)) + && is_a(e_replaced.op(1)) && ex_to(e_replaced.op(1)).is_integer())) { + for (auto & it : repl) { + if (e_replaced.op(0).is_equal(it.second) // The base is an allocated symbol or base of power + || (is_a(it.second) && e_replaced.op(0).is_equal(it.second.op(0)))) { + ex ratio; // We bind together two above cases + if (is_a(it.second)) + ratio = normal(e_replaced.op(1) / it.second.op(1)); + else + ratio = e_replaced.op(1); + if (is_a(ratio) && ex_to(ratio).is_rational()) { + // Different powers can be treated as powers of the same basic equation + if (ex_to(ratio).is_integer()) { + // If ratio is an integer then this is simply the power of the existing symbol. + //std::clog << e_replaced << " is a " << ratio << " power of " << it.first << std::endl; + return pow(it.first, ratio); + } else { + // otherwise we need to give the replacement pattern to change + // the previous expression... + ex es = dynallocate(); + ex Num = numer(ratio); + modifier.append(it.first == power(es, denom(ratio))); + //std::clog << e_replaced << " is power " << Num << " and " + // << it.first << " is power " << denom(ratio) << " of the common base " + // << pow(e_replaced.op(0), e_replaced.op(1)/Num) << std::endl; + // ... and modify the replacement tables + rev_lookup.erase(it.second); + rev_lookup.insert(std::make_pair(pow(e_replaced.op(0), e_replaced.op(1)/Num), es)); + repl.erase(it.first); + repl.insert(std::make_pair(es, pow(e_replaced.op(0), e_replaced.op(1)/Num))); + return power(es, Num); + } + } + } + } + // There is no an existing substitution, thus we are creating a new one. + // This needs to be done separately to treat possible occurrences of + // b = e_replaced.op(0) elsewhere in the expression as pow(b, 1). + ex degree = pow(e_replaced.op(1), _ex_1); + if (is_a(degree) && ex_to(degree).is_integer()) { + ex es = dynallocate(); + modifier.append(e_replaced.op(0) == power(es, degree)); + repl.insert(std::make_pair(es, e_replaced)); + rev_lookup.insert(std::make_pair(e_replaced, es)); + return es; + } } // Otherwise create new symbol and add to list, taking care that the @@ -2352,8 +2408,15 @@ ex mul::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const ex power::normal(exmap & repl, exmap & rev_lookup, lst & modifier) const { // Normalize basis and exponent (exponent gets reassembled) + int nmod = modifier.nops(); // To watch new modifiers to the replacement list ex n_basis = ex_to(basis).normal(repl, rev_lookup, modifier); + for (int imod = nmod; imod < modifier.nops(); ++imod) + n_basis = n_basis.subs(modifier.op(imod), subs_options::no_pattern); + + nmod = modifier.nops(); ex n_exponent = ex_to(exponent).normal(repl, rev_lookup, modifier); + for (int imod = nmod; imod < modifier.nops(); ++imod) + n_exponent = n_exponent.subs(modifier.op(imod), subs_options::no_pattern); n_exponent = n_exponent.op(0) / n_exponent.op(1); if (n_exponent.info(info_flags::integer)) { @@ -2430,6 +2493,7 @@ ex ex::normal() const ex e = bp->normal(repl, rev_lookup, modifier); GINAC_ASSERT(is_a(e)); + //if (modifier.nops() > 0) std::clog << "Numer: " << e.op(0) << "; Denom: " << e.op(1) << std::endl; // Re-insert replaced symbols if (!repl.empty()) { -- 2.28.0