From aec32ae2e31d7fda07af3912c44a9a0ce7dd7fe4 Mon Sep 17 00:00:00 2001 From: "Vladimir V. Kisil" Date: Mon, 29 Jul 2013 22:25:58 +0100 Subject: [PATCH 2/9] Add an expansion law in power::expand(). The expansion (x*p)^c -> x^c*p^c is done for a positive expression p and a real exponent c. The corresponding check is added to exam_powerlaws. Signed-off-by: Vladimir V. Kisil --- check/exam_powerlaws.cpp | 32 ++++++++++++++++++++++++++++++++ ginac/flags.h | 6 ++++++ ginac/mul.cpp | 11 +++++++++++ ginac/power.cpp | 47 +++++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 96 insertions(+) diff --git a/check/exam_powerlaws.cpp b/check/exam_powerlaws.cpp index beb4dfb..aa5e108 100644 --- a/check/exam_powerlaws.cpp +++ b/check/exam_powerlaws.cpp @@ -288,6 +288,37 @@ static unsigned exam_powerlaws5() return 0; } +static unsigned exam_powerlaws6() +{ + // check expansion rules for positive symbols + + symbol a("a"); + symbol b("b"); + realsymbol c("c"); + realsymbol x("x"); + realsymbol y("y"); + possymbol p("p"); + possymbol q("q"); + numeric half=numeric(1,2); + + ex e1 = pow(5*pow(3*a*b*x*y*p*q,2),7*half*c).expand(); + ex e2 = pow(p,7*c)*pow(q,7*c)*pow(pow(a*b*x*y,2),numeric(7,2)*c)*pow(45,numeric(7,2)*c); + if (!e1.is_equal(e2)) { + clog << "Could not expand real exponents with positive bases in " << e1 << endl; + return 1; + } + + ex e3 = pow(-pow(-a*x*p,3)*pow(b*y*p,3),half*c).expand().normal(); + ex e4 = pow(p,3*c)*pow(pow(a*b*x*y,3),half*c); + + if (!e3.is_equal(e4)) { + clog << "Could not expand real exponents with positive bases in " << e3 << endl; + return 1; + } + + return 0; +} + unsigned exam_powerlaws() { unsigned result = 0; @@ -299,6 +330,7 @@ unsigned exam_powerlaws() result += exam_powerlaws3(); cout << '.' << flush; result += exam_powerlaws4(); cout << '.' << flush; result += exam_powerlaws5(); cout << '.' << flush; + result += exam_powerlaws6(); cout << '.' << flush; return result; } diff --git a/ginac/flags.h b/ginac/flags.h index 8f8f25d..f4f0f6b 100644 --- a/ginac/flags.h +++ b/ginac/flags.h @@ -197,6 +197,9 @@ public: not_shareable = 0x0010, ///< don't share instances of this object between different expressions unless explicitly asked to (used by ex::compare()) has_indices = 0x0020, has_no_indices = 0x0040 // ! (has_indices || has_no_indices) means "don't know" + is_positive = 0x0080, + is_negative = 0x0100, + purely_indefinite = 0x0200, // If set in a mul, then it does not contains any terms with determined signs, used in power::expand() }; }; @@ -221,6 +224,9 @@ public: odd, prime, + // is meaningful for mul only + indefinite, + // answered by class relation relation, relation_equal, diff --git a/ginac/mul.cpp b/ginac/mul.cpp index cee5cd7..3ec279d 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -307,6 +307,11 @@ bool mul::info(unsigned inf) const } case info_flags::positive: case info_flags::negative: { + if ((inf==info_flags::positive) && (flags & status_flags::is_positive)) + return true; + else if ((inf==info_flags::negative) && (flags & status_flags::is_negative)) + return true; + bool pos = true; epvector::const_iterator i = seq.begin(), end = seq.end(); while (i != end) { @@ -320,9 +325,12 @@ bool mul::info(unsigned inf) const } if (overall_coeff.info(info_flags::negative)) pos = !pos; + flags |= (pos? status_flags::is_positive : status_flags::is_negative); return (inf ==info_flags::positive? pos : !pos); } case info_flags::nonnegative: { + if (flags & status_flags::is_positive) + return true; bool pos = true; epvector::const_iterator i = seq.begin(), end = seq.end(); while (i != end) { @@ -373,6 +381,9 @@ bool mul::info(unsigned inf) const return false; return pos; } + case info_flags::indefinite: { + return (flags & status_flags::purely_indefinite); + } } return inherited::info(inf); } diff --git a/ginac/power.cpp b/ginac/power.cpp index b815d17..ec3717e 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -792,6 +792,53 @@ ex power::expand(unsigned options) const return *this; } + // (x*p)^c -> x^c * p^c, if p>0 and c is real + // make sense before the basis was expanded + if (is_exactly_a(basis) && !basis.info(info_flags::indefinite) && exponent.info(info_flags::real)) { + const mul &m = ex_to(basis); + exvector prodseq; + epvector powseq; + prodseq.reserve(m.seq.size() + 1); + powseq.reserve(m.seq.size() + 1); + epvector::const_iterator last = m.seq.end(); + epvector::const_iterator cit = m.seq.begin(); + bool possign=true; + + // searching for positive/negative factors + while (cit!=last) { + ex e=m.recombine_pair_to_ex(*cit); + if (e.info(info_flags::positive)) + prodseq.push_back(pow(e, exponent).expand(options)); + else if (e.info(info_flags::negative)) { + prodseq.push_back(pow(-e, exponent).expand(options)); + possign = !possign; + } else + powseq.push_back(*cit); + ++cit; + } + + // take care on the numeric coefficient + ex coeff=(possign? _ex1 : _ex_1); + if (m.overall_coeff.info(info_flags::positive) && m.overall_coeff != _ex1) + prodseq.push_back(power(m.overall_coeff, exponent)); + else if (m.overall_coeff.info(info_flags::negative) && m.overall_coeff != _ex_1) { + prodseq.push_back(power(-m.overall_coeff, exponent)); + possign = !possign; + } else { + coeff *=m.overall_coeff; + } + + // If positive/negative factors are found then extract them. + // In either case we set a flag to avoid the second run on a part, + // which does not have positive/negative terms + if (prodseq.size() > 0) { + ex newbasis=coeff*mul(powseq); + ex_to(newbasis).setflag(status_flags::purely_indefinite); + return ((new mul(prodseq))->setflag(status_flags::dynallocated)*(new power(newbasis, exponent))->setflag(status_flags::dynallocated).expand(options)).expand(options); + } else + ex_to(basis).setflag(status_flags::purely_indefinite); + } + const ex expanded_basis = basis.expand(options); const ex expanded_exponent = exponent.expand(options); -- 1.8.4.rc3