X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpower.cpp;h=f5ab330d450f0c76ac2207ef0a047b7f1b6a39bb;hp=ec339b9a660afc6aac21fa9307b6bec5e26144ce;hb=2c7d26f82ab2938ccc42372459a6640137a910a1;hpb=67467d256b44f5e08498ca81c946d9ffaa25d1e2 diff --git a/ginac/power.cpp b/ginac/power.cpp index ec339b9a..f5ab330d 100644 --- a/ginac/power.cpp +++ b/ginac/power.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's symbolic exponentiation (basis^exponent). */ /* - * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2014 Johannes Gutenberg University Mainz, Germany * * This program is free software; you can redistribute it and/or modify * it under the terms of the GNU General Public License as published by @@ -20,11 +20,6 @@ * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA */ -#include -#include -#include -#include - #include "power.h" #include "expairseq.h" #include "add.h" @@ -43,6 +38,11 @@ #include "relational.h" #include "compiler.h" +#include +#include +#include +#include + namespace GiNaC { GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic, @@ -59,7 +59,7 @@ typedef std::vector intvector; // default constructor ////////// -power::power() : inherited(&power::tinfo_static) { } +power::power() { } ////////// // other constructors @@ -71,8 +71,9 @@ power::power() : inherited(&power::tinfo_static) { } // archiving ////////// -power::power(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst) +void power::read_archive(const archive_node &n, lst &sym_lst) { + inherited::read_archive(n, sym_lst); n.find_ex("basis", basis, sym_lst); n.find_ex("exponent", exponent, sym_lst); } @@ -84,8 +85,6 @@ void power::archive(archive_node &n) const n.add_ex("exponent", exponent); } -DEFAULT_UNARCHIVE(power) - ////////// // functions overriding virtual functions from base classes ////////// @@ -241,6 +240,10 @@ bool power::info(unsigned inf) const basis.info(inf); case info_flags::expanded: return (flags & status_flags::expanded); + case info_flags::positive: + return basis.info(info_flags::positive) && exponent.info(info_flags::real); + case info_flags::nonnegative: + return basis.info(info_flags::real) && exponent.info(info_flags::integer) && exponent.info(info_flags::even); case info_flags::has_indices: { if (flags & status_flags::has_indices) return true; @@ -286,11 +289,16 @@ ex power::map(map_function & f) const bool power::is_polynomial(const ex & var) const { - if (exponent.has(var)) - return false; - if (!exponent.info(info_flags::nonnegint)) - return false; - return basis.is_polynomial(var); + if (basis.is_polynomial(var)) { + if (basis.has(var)) + // basis is non-constant polynomial in var + return exponent.info(info_flags::nonnegint); + else + // basis is constant in var + return !exponent.has(var); + } + // basis is a non-polynomial function of var + return false; } int power::degree(const ex & s) const @@ -361,7 +369,7 @@ ex power::coeff(const ex & s, int n) const * - ^(1,x) -> 1 * - ^(c1,c2) -> *(c1^n,c1^(c2-n)) (so that 0<(c2-n)<1, try to evaluate roots, possibly in numerator and denominator of c1) * - ^(^(x,c1),c2) -> ^(x,c1*c2) if x is positive and c1 is real. - * - ^(^(x,c1),c2) -> ^(x,c1*c2) (c2 integer or -1 < c1 <= 1, case c1=1 should not happen, see below!) + * - ^(^(x,c1),c2) -> ^(x,c1*c2) (c2 integer or -1 < c1 <= 1 or (c1=-1 and c2>0), case c1=1 should not happen, see below!) * - ^(*(x,y,z),c) -> *(x^c,y^c,z^c) (if c integer) * - ^(*(x,c1),c2) -> ^(x,c2)*c1^c2 (c1>0) * - ^(*(x,c1),c2) -> ^(-x,c2)*c1^c2 (c1<0) @@ -377,17 +385,13 @@ ex power::eval(int level) const const ex & ebasis = level==1 ? basis : basis.eval(level-1); const ex & eexponent = level==1 ? exponent : exponent.eval(level-1); - bool basis_is_numerical = false; - bool exponent_is_numerical = false; - const numeric *num_basis; - const numeric *num_exponent; + const numeric *num_basis = NULL; + const numeric *num_exponent = NULL; if (is_exactly_a(ebasis)) { - basis_is_numerical = true; num_basis = &ex_to(ebasis); } if (is_exactly_a(eexponent)) { - exponent_is_numerical = true; num_exponent = &ex_to(eexponent); } @@ -404,7 +408,7 @@ ex power::eval(int level) const return ebasis; // ^(0,c1) -> 0 or exception (depending on real value of c1) - if (ebasis.is_zero() && exponent_is_numerical) { + if ( ebasis.is_zero() && num_exponent ) { if ((num_exponent->real()).is_zero()) throw (std::domain_error("power::eval(): pow(0,I) is undefined")); else if ((num_exponent->real()).is_negative()) @@ -425,11 +429,11 @@ ex power::eval(int level) const if (is_exactly_a(ebasis) && ebasis.op(0).info(info_flags::positive) && ebasis.op(1).info(info_flags::real)) return power(ebasis.op(0), ebasis.op(1) * eexponent); - if (exponent_is_numerical) { + if ( num_exponent ) { // ^(c1,c2) -> c1^c2 (c1, c2 numeric(), // except if c1,c2 are rational, but c1^c2 is not) - if (basis_is_numerical) { + if ( num_basis ) { const bool basis_is_crational = num_basis->is_crational(); const bool exponent_is_crational = num_exponent->is_crational(); if (!basis_is_crational || !exponent_is_crational) { @@ -483,7 +487,7 @@ ex power::eval(int level) const } // ^(^(x,c1),c2) -> ^(x,c1*c2) - // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1, + // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1 or (c1=-1 and c2>0), // case c1==1 should not happen, see below!) if (is_exactly_a(ebasis)) { const power & sub_power = ex_to(ebasis); @@ -492,7 +496,8 @@ 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() + || (num_sub_exponent == *_num_1_p && num_exponent->is_positive())) { return power(sub_basis,num_sub_exponent.mul(*num_exponent)); } } @@ -542,6 +547,7 @@ ex power::eval(int level) const if (num_coeff.is_positive()) { mul *mulp = new mul(mulref); mulp->overall_coeff = _ex1; + mulp->setflag(status_flags::dynallocated); mulp->clearflag(status_flags::evaluated); mulp->clearflag(status_flags::hash_calculated); return (new mul(power(*mulp,exponent), @@ -551,6 +557,7 @@ ex power::eval(int level) const if (!num_coeff.is_equal(*_num_1_p)) { mul *mulp = new mul(mulref); mulp->overall_coeff = _ex_1; + mulp->setflag(status_flags::dynallocated); mulp->clearflag(status_flags::evaluated); mulp->clearflag(status_flags::hash_calculated); return (new mul(power(*mulp,exponent), @@ -635,7 +642,7 @@ bool power::has(const ex & other, unsigned options) const } // from mul.cpp -extern bool tryfactsubs(const ex &, const ex &, int &, lst &); +extern bool tryfactsubs(const ex &, const ex &, int &, exmap&); ex power::subs(const exmap & m, unsigned options) const { @@ -651,9 +658,13 @@ ex power::subs(const exmap & m, unsigned options) const for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) { int nummatches = std::numeric_limits::max(); - lst repls; - if (tryfactsubs(*this, it->first, nummatches, repls)) - return (ex_to((*this) * power(it->second.subs(ex(repls), subs_options::no_pattern) / it->first.subs(ex(repls), subs_options::no_pattern), nummatches))).subs_one_level(m, options); + exmap repls; + if (tryfactsubs(*this, it->first, nummatches, repls)) { + ex anum = it->second.subs(repls, subs_options::no_pattern); + ex aden = it->first.subs(repls, subs_options::no_pattern); + ex result = (*this)*power(anum/aden, nummatches); + return (ex_to(result)).subs_one_level(m, options); + } } return subs_one_level(m, options); @@ -666,12 +677,23 @@ ex power::eval_ncmul(const exvector & v) const ex power::conjugate() const { - ex newbasis = basis.conjugate(); - ex newexponent = exponent.conjugate(); - if (are_ex_trivially_equal(basis, newbasis) && are_ex_trivially_equal(exponent, newexponent)) { - return *this; + // conjugate(pow(x,y))==pow(conjugate(x),conjugate(y)) unless on the + // branch cut which runs along the negative real axis. + if (basis.info(info_flags::positive)) { + ex newexponent = exponent.conjugate(); + if (are_ex_trivially_equal(exponent, newexponent)) { + return *this; + } + return (new power(basis, newexponent))->setflag(status_flags::dynallocated); } - return (new power(newbasis, newexponent))->setflag(status_flags::dynallocated); + if (exponent.info(info_flags::integer)) { + ex newbasis = basis.conjugate(); + if (are_ex_trivially_equal(basis, newbasis)) { + return *this; + } + return (new power(newbasis, exponent))->setflag(status_flags::dynallocated); + } + return conjugate_function(*this).hold(); } ex power::real_part() const @@ -721,8 +743,7 @@ ex power::imag_part() const ex b=basis.imag_part(); ex c=exponent.real_part(); ex d=exponent.imag_part(); - return - power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis))); + return power(abs(basis),c)*exp(-d*atan2(b,a))*sin(c*atan2(b,a)+d*log(abs(basis))); } // protected @@ -765,7 +786,7 @@ unsigned power::return_type() const return basis.return_type(); } -tinfo_t power::return_type_tinfo() const +return_type_t power::return_type_tinfo() const { return basis.return_type_tinfo(); } @@ -864,7 +885,6 @@ ex power::expand_add(const add & a, int n, unsigned options) const intvector k(m-1); intvector k_cum(m-1); // k_cum[l]:=sum(i=0,l,k[l]); intvector upper_limit(m-1); - int l; for (size_t l=0; l(b)); GINAC_ASSERT(!is_exactly_a(b) || @@ -890,7 +910,7 @@ ex power::expand_add(const add & a, int n, unsigned options) const term.push_back(power(b,k[l])); } - const ex & b = a.op(l); + const ex & b = a.op(m - 1); GINAC_ASSERT(!is_exactly_a(b)); GINAC_ASSERT(!is_exactly_a(b) || !is_exactly_a(ex_to(b).exponent) || @@ -904,7 +924,7 @@ ex power::expand_add(const add & a, int n, unsigned options) const term.push_back(power(b,n-k_cum[m-2])); numeric f = binomial(numeric(n),numeric(k[0])); - for (l=1; lsetflag(status_flags::dynallocated)).expand(options)); // increment k[] - l = m-2; - while ((l>=0) && ((++k[l])>upper_limit[l])) { + bool done = false; + std::size_t l = m - 2; + while ((++k[l]) > upper_limit[l]) { k[l] = 0; - --l; + if (l != 0) + --l; + else { + done = true; + break; + } } - if (l<0) break; + if (done) + break; // recalc k_cum[] and upper_limit[] k_cum[l] = (l==0 ? k[0] : k_cum[l-1]+k[l]); @@ -1051,4 +1078,6 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr return result; } +GINAC_BIND_UNARCHIVER(power); + } // namespace GiNaC