From 8ba901b532b844d49ffbd11f906f9124f357db1f Mon Sep 17 00:00:00 2001 From: Jens Vollinga Date: Wed, 11 Jul 2007 20:59:19 +0000 Subject: [PATCH] Additional transformations for mul and power [Sheplyakov]. --- check/exam_pseries.cpp | 3 ++- ginac/mul.cpp | 57 ++++++++++++++++++++++++++++++++++++++++-- ginac/power.cpp | 29 +++++++++++++++++++-- 3 files changed, 84 insertions(+), 5 deletions(-) diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp index 1c8acd93..2e5d2d6e 100644 --- a/check/exam_pseries.cpp +++ b/check/exam_pseries.cpp @@ -354,7 +354,8 @@ static unsigned exam_series13() { unsigned result = 0; - ex e = pow(2,x)*( 1/x*(-(1+x)/(1-x)) + (1+x)/x/(1-x)); + ex e = (new mul(pow(2,x), (1/x*(-(1+x)/(1-x)) + (1+x)/x/(1-x))) + )->setflag(status_flags::evaluated); ex d = Order(x); result += check_series(e,0,d,1); diff --git a/ginac/mul.cpp b/ginac/mul.cpp index 9bbcd611..b4ef8972 100644 --- a/ginac/mul.cpp +++ b/ginac/mul.cpp @@ -446,9 +446,62 @@ ex mul::eval(int level) const } return (new add(distrseq, ex_to(addref.overall_coeff). - mul_dyn(ex_to(overall_coeff)))) - ->setflag(status_flags::dynallocated | status_flags::evaluated); + mul_dyn(ex_to(overall_coeff))) + )->setflag(status_flags::dynallocated | status_flags::evaluated); + } else if (seq_size >= 2) { + // Strip the content and the unit part from each term. Thus + // things like (-x+a)*(3*x-3*a) automagically turn into - 3*(x-a)2 + + epvector::const_iterator last = seq.end(); + epvector::const_iterator i = seq.begin(); + while (i!=last) { + if (! (is_a(i->rest) && i->coeff.is_equal(_ex1))) { + // power::eval has such a rule, no need to handle powers here + ++i; + continue; + } + + // XXX: What is the best way to check if the polynomial is a primitive? + numeric c = i->rest.integer_content(); + const numeric& lead_coeff = + ex_to(ex_to(i->rest).seq.begin()->coeff).div_dyn(c); + const bool canonicalizable = lead_coeff.is_integer(); + + // XXX: The main variable is chosen in a random way, so this code + // does NOT transform the term into the canonical form (thus, in some + // very unlucky event it can even loop forever). Hopefully the main + // variable will be the same for all terms in *this + const bool unit_normal = lead_coeff.is_pos_integer(); + if ((c == *_num1_p) && ((! canonicalizable) || unit_normal)) { + ++i; + continue; + } + + std::auto_ptr s(new epvector); + s->reserve(seq.size()); + + epvector::const_iterator j=seq.begin(); + while (j!=i) { + s->push_back(*j); + ++j; + } + + if (! unit_normal) { + c = c.mul(*_num_1_p); + } + const ex primitive = (i->rest)/c; + s->push_back(expair(primitive, _ex1)); + ++j; + + while (j!=last) { + s->push_back(*j); + ++j; + } + return (new mul(s, ex_to(overall_coeff).mul_dyn(c)) + )->setflag(status_flags::dynallocated); + } } + return this->hold(); } diff --git a/ginac/power.cpp b/ginac/power.cpp index fb5e351e..8829bbae 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -467,8 +467,9 @@ ex power::eval(int level) const if (is_exactly_a(sub_exponent)) { const numeric & num_sub_exponent = ex_to(sub_exponent); GINAC_ASSERT(num_sub_exponent!=numeric(1)); - if (num_exponent->is_integer() || (abs(num_sub_exponent) - (*_num1_p)).is_negative()) + if (num_exponent->is_integer() || (abs(num_sub_exponent) - (*_num1_p)).is_negative()) { return power(sub_basis,num_sub_exponent.mul(*num_exponent)); + } } } @@ -476,7 +477,31 @@ ex power::eval(int level) const if (num_exponent->is_integer() && is_exactly_a(ebasis)) { return expand_mul(ex_to(ebasis), *num_exponent, 0); } - + + // (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4) + if (num_exponent->is_integer() && is_exactly_a(ebasis)) { + const numeric icont = ebasis.integer_content(); + const numeric& lead_coeff = + ex_to(ex_to(ebasis).seq.begin()->coeff).div_dyn(icont); + + const bool canonicalizable = lead_coeff.is_integer(); + const bool unit_normal = lead_coeff.is_pos_integer(); + + if (icont != *_num1_p) { + return (new mul(power(ebasis/icont, *num_exponent), power(icont, *num_exponent)) + )->setflag(status_flags::dynallocated); + } + + if (canonicalizable && (! unit_normal)) { + if (num_exponent->is_even()) { + return power(-ebasis, *num_exponent); + } else { + return (new mul(power(-ebasis, *num_exponent), *_num_1_p) + )->setflag(status_flags::dynallocated); + } + } + } + // ^(*(...,x;c1),c2) -> *(^(*(...,x;1),c2),c1^c2) (c1, c2 numeric(), c1>0) // ^(*(...,x;c1),c2) -> *(^(*(...,x;-1),c2),(-c1)^c2) (c1, c2 numeric(), c1<0) if (is_exactly_a(ebasis)) { -- 2.44.0