From 1b0be38f3eb98912a80f83d73cbec551d5b2013d Mon Sep 17 00:00:00 2001 From: "Vladimir V. Kisil" Date: Sun, 28 Jul 2013 14:48:56 +0100 Subject: [PATCH] Add an expansion law in power::expand() (x*p)^c -> x^c*p^c 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 | 24 ++++++++++++++++++++++++ ginac/power.cpp | 35 +++++++++++++++++++++++++++++++---- 2 files changed, 55 insertions(+), 4 deletions(-) diff --git a/check/exam_powerlaws.cpp b/check/exam_powerlaws.cpp index cab16df..f0e8b67 100644 --- a/check/exam_powerlaws.cpp +++ b/check/exam_powerlaws.cpp @@ -288,6 +288,29 @@ 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; + } + + return 0; +} + unsigned exam_powerlaws() { unsigned result = 0; @@ -299,6 +322,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/power.cpp b/ginac/power.cpp index 9a8f7d5..d01a88f 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -824,11 +824,38 @@ ex power::expand(unsigned options) const if (!is_exactly_a(expanded_exponent) || !ex_to(expanded_exponent).is_integer()) { - if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent)) { - return this->hold(); - } else { + if (is_exactly_a(expanded_basis) && expanded_exponent.info(info_flags::real)) { + const mul &m = ex_to(expanded_basis); + exvector prodseq, 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; + + while (cit!=last) { + ex e=m.recombine_pair_to_ex(*cit); + if (e.info(info_flags::positive)) + prodseq.push_back(power(e, expanded_exponent)); + else if (e.info(info_flags::negative)) { + prodseq.push_back(power(-e, expanded_exponent)); + possign = !possign; + } else + powseq.push_back(m.recombine_pair_to_ex(*cit)); + ++cit; + } + + if (m.overall_coeff.info(info_flags::positive)) + prodseq.push_back(power(m.overall_coeff, expanded_exponent)); + else if (m.overall_coeff.info(info_flags::negative)) { + prodseq.push_back(power(-m.overall_coeff, expanded_exponent)); + possign = !possign; + } else + powseq.push_back(m.overall_coeff); + + return (new mul(prodseq))->setflag(status_flags::dynallocated)*(new power((possign? _ex1 : _ex_1)*mul(powseq), expanded_exponent))->setflag(status_flags::dynallocated); + } else return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)); - } } // integer numeric exponent -- 1.7.10.4