[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 <V.Kisilv@leeds.ac.uk>

diff --git a/check/exam_normalization.cpp b/check/exam_normalization.cpp
index c432102193ec2ef5534bc8923b42ebe7f0cb040e..93161cb31e04f22d124d21a862d99be95f67f6cc 100644 (file)
--- 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 342ef06b411480e4efc42f269d280065d675922e..828866b57a36fc8646943629a17aa5d6c21ccdd2 100644 (file)
--- 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<numeric>(e_replaced))
+               for (auto & it : repl)
+                       if (is_a<power>(it.second) && e_replaced.is_equal(it.second.op(0))) {
+                               ex degree = pow(it.second.op(1), _ex_1);
+                               if (is_a<numeric>(degree) && ex_to<numeric>(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<function>(e_replaced) && is_ex_the_function(e_replaced, exp)) {
@@ -2090,6 +2099,53 @@ static ex replace_with_symbol(const ex & e, exmap & repl, exmap & rev_lookup, ls
                                }
                        }
                }
+       } else if (is_a<power>(e_replaced) && !is_a<numeric>(e_replaced.op(0)) // We do not replace simple monomials like x^3 or sqrt(2)
+                  && ! (is_a<symbol>(e_replaced.op(0))
+                      && is_a<numeric>(e_replaced.op(1)) && ex_to<numeric>(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<power>(it.second) && e_replaced.op(0).is_equal(it.second.op(0)))) {
+                               ex ratio; // We bind together two above cases
+                               if (is_a<power>(it.second))
+                                       ratio = normal(e_replaced.op(1) / it.second.op(1));
+                               else
+                                       ratio = e_replaced.op(1);
+                               if (is_a<numeric>(ratio) && ex_to<numeric>(ratio).is_rational())  {
+                                       // Different powers can be treated as powers of the same basic equation
+                                       if (ex_to<numeric>(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 dynallocate<power>(it.first, ratio);
+                                       } else {
+                                               // otherwise we need to give the replacement pattern to change
+                                               // the previous expression...
+                                               ex es = dynallocate<symbol>();
+                                               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({pow(e_replaced.op(0), e_replaced.op(1)/Num), es});
+                                               repl.erase(it.first);
+                                               repl.insert({es, pow(e_replaced.op(0), e_replaced.op(1)/Num)});
+                                               return dynallocate<power>(es, Num);
+                                       }
+                               }
+                       }
+               }
+               // There is no 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<numeric>(degree) && ex_to<numeric>(degree).is_integer()) {
+                       ex es = dynallocate<symbol>();
+                       modifier.append(e_replaced.op(0) == power(es, degree));
+                       repl.insert({es, e_replaced});
+                       rev_lookup.insert({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<basic>(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<basic>(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)) {