X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpower.cpp;h=c6bf7b9770d2fa939eb242fef780f6aa84286dc5;hp=75fd2ff51a9d12bd09d867d30fdf50b860ec6240;hb=1d73fed17debbc88819dc7f84a4a2e608d8eb978;hpb=690cd58cc13ad5052eb5851c573984965d0c40c1 diff --git a/ginac/power.cpp b/ginac/power.cpp index 75fd2ff5..c6bf7b97 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-2006 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2010 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" @@ -41,6 +36,12 @@ #include "archive.h" #include "utils.h" #include "relational.h" +#include "compiler.h" + +#include +#include +#include +#include namespace GiNaC { @@ -49,7 +50,8 @@ GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic, print_func(&power::do_print_latex). print_func(&power::do_print_csrc). print_func(&power::do_print_python). - print_func(&power::do_print_python_repr)) + print_func(&power::do_print_python_repr). + print_func(&power::do_print_csrc_cl_N)) typedef std::vector intvector; @@ -57,7 +59,7 @@ typedef std::vector intvector; // default constructor ////////// -power::power() : inherited(&power::tinfo_static) { } +power::power() { } ////////// // other constructors @@ -69,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); } @@ -82,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 ////////// @@ -161,6 +162,21 @@ static void print_sym_pow(const print_context & c, const symbol &x, int exp) } } +void power::do_print_csrc_cl_N(const print_csrc_cl_N& c, unsigned level) const +{ + if (exponent.is_equal(_ex_1)) { + c.s << "recip("; + basis.print(c); + c.s << ')'; + return; + } + c.s << "expt("; + basis.print(c); + c.s << ", "; + exponent.print(c); + c.s << ')'; +} + void power::do_print_csrc(const print_csrc & c, unsigned level) const { // Integer powers of symbols are printed in a special, optimized way @@ -171,29 +187,20 @@ void power::do_print_csrc(const print_csrc & c, unsigned level) const c.s << '('; else { exp = -exp; - if (is_a(c)) - c.s << "recip("; - else - c.s << "1.0/("; + c.s << "1.0/("; } print_sym_pow(c, ex_to(basis), exp); c.s << ')'; // ^-1 is printed as "1.0/" or with the recip() function of CLN } else if (exponent.is_equal(_ex_1)) { - if (is_a(c)) - c.s << "recip("; - else - c.s << "1.0/("; + c.s << "1.0/("; basis.print(c); c.s << ')'; - // Otherwise, use the pow() or expt() (CLN) functions + // Otherwise, use the pow() function } else { - if (is_a(c)) - c.s << "expt("; - else - c.s << "pow("; + c.s << "pow("; basis.print(c); c.s << ','; exponent.print(c); @@ -231,6 +238,25 @@ bool power::info(unsigned inf) const case info_flags::algebraic: return !exponent.info(info_flags::integer) || 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::has_indices: { + if (flags & status_flags::has_indices) + return true; + else if (flags & status_flags::has_no_indices) + return false; + else if (basis.info(info_flags::has_indices)) { + setflag(status_flags::has_indices); + clearflag(status_flags::has_no_indices); + return true; + } else { + clearflag(status_flags::has_indices); + setflag(status_flags::has_no_indices); + return false; + } + } } return inherited::info(inf); } @@ -335,6 +361,7 @@ ex power::coeff(const ex & s, int n) const * - ^(0,c) -> 0 or exception (depending on the real part of c) * - ^(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,y,z),c) -> *(x^c,y^c,z^c) (if c integer) * - ^(*(x,c1),c2) -> ^(x,c2)*c1^c2 (c1>0) @@ -351,17 +378,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); } @@ -378,7 +401,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()) @@ -395,11 +418,15 @@ ex power::eval(int level) const if (is_exactly_a(ebasis)) return ex_to(ebasis).power(eexponent); - if (exponent_is_numerical) { + // Turn (x^c)^d into x^(c*d) in the case that x is positive and c is real. + 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 ( 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) { @@ -462,8 +489,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)); + } } } @@ -471,7 +499,35 @@ 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)) { + numeric icont = ebasis.integer_content(); + const numeric lead_coeff = + ex_to(ex_to(ebasis).seq.begin()->coeff).div(icont); + + const bool canonicalizable = lead_coeff.is_integer(); + const bool unit_normal = lead_coeff.is_pos_integer(); + if (canonicalizable && (! unit_normal)) + icont = icont.mul(*_num_1_p); + + if (canonicalizable && (icont != *_num1_p)) { + const add& addref = ex_to(ebasis); + add* addp = new add(addref); + addp->setflag(status_flags::dynallocated); + addp->clearflag(status_flags::hash_calculated); + addp->overall_coeff = ex_to(addp->overall_coeff).div_dyn(icont); + for (epvector::iterator i = addp->seq.begin(); i != addp->seq.end(); ++i) + i->coeff = ex_to(i->coeff).div_dyn(icont); + + const numeric c = icont.power(*num_exponent); + if (likely(c != *_num1_p)) + return (new mul(power(*addp, *num_exponent), c))->setflag(status_flags::dynallocated); + else + return power(*addp, *num_exponent); + } + } + // ^(*(...,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)) { @@ -483,6 +539,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), @@ -492,6 +549,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), @@ -576,7 +634,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 { @@ -592,9 +650,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); @@ -607,12 +669,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); + } + 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 (new power(newbasis, newexponent))->setflag(status_flags::dynallocated); + return conjugate_function(*this).hold(); } ex power::real_part() const @@ -706,16 +779,19 @@ 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(); } ex power::expand(unsigned options) const { - if (options == 0 && (flags & status_flags::expanded)) + if (is_a(basis) && exponent.info(info_flags::integer)) { + // A special case worth optimizing. + setflag(status_flags::expanded); return *this; - + } + const ex expanded_basis = basis.expand(options); const ex expanded_exponent = exponent.expand(options); @@ -802,7 +878,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) || @@ -828,7 +903,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) || @@ -842,7 +917,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]); @@ -938,7 +1020,7 @@ ex power::expand_add_2(const add & a, unsigned options) const return (new add(sum))->setflag(status_flags::dynallocated | status_flags::expanded); } -/** Expand factors of m in m^n where m is a mul and n is and integer. +/** Expand factors of m in m^n where m is a mul and n is an integer. * @see power::expand */ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool from_expand) const { @@ -948,8 +1030,13 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr return _ex1; } + // do not bother to rename indices if there are no any. + if ((!(options & expand_options::expand_rename_idx)) + && m.info(info_flags::has_indices)) + options |= expand_options::expand_rename_idx; // Leave it to multiplication since dummy indices have to be renamed - if (get_all_dummy_indices(m).size() > 0 && n.is_positive()) { + if ((options & expand_options::expand_rename_idx) && + (get_all_dummy_indices(m).size() > 0) && n.is_positive()) { ex result = m; exvector va = get_all_dummy_indices(m); sort(va.begin(), va.end(), ex_is_less()); @@ -966,19 +1053,13 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr epvector::const_iterator last = m.seq.end(); epvector::const_iterator cit = m.seq.begin(); while (cit!=last) { - if (is_exactly_a(cit->rest)) { - distrseq.push_back(m.combine_pair_with_coeff_to_pair(*cit, n)); - } else { - // it is safe not to call mul::combine_pair_with_coeff_to_pair() - // since n is an integer - numeric new_coeff = ex_to(cit->coeff).mul(n); - if (from_expand && is_exactly_a(cit->rest) && new_coeff.is_pos_integer()) { - // this happens when e.g. (a+b)^(1/2) gets squared and - // the resulting product needs to be reexpanded - need_reexpand = true; - } - distrseq.push_back(expair(cit->rest, new_coeff)); + expair p = m.combine_pair_with_coeff_to_pair(*cit, n); + if (from_expand && is_exactly_a(cit->rest) && ex_to(p.coeff).is_pos_integer()) { + // this happens when e.g. (a+b)^(1/2) gets squared and + // the resulting product needs to be reexpanded + need_reexpand = true; } + distrseq.push_back(p); ++cit; } @@ -990,4 +1071,6 @@ ex power::expand_mul(const mul & m, const numeric & n, unsigned options, bool fr return result; } +GINAC_BIND_UNARCHIVER(power); + } // namespace GiNaC