X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fpower.cpp;h=ec339b9a660afc6aac21fa9307b6bec5e26144ce;hp=9e720b70b04280a150de8d6ab99e6932f4d39699;hb=b21923803208edf7ec3b0bfaa51f813cb25a0d86;hpb=619d77d2676f7f1a562fb9fefc0ba6754fe2d750 diff --git a/ginac/power.cpp b/ginac/power.cpp index 9e720b70..ec339b9a 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-2007 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2008 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 @@ -41,6 +41,7 @@ #include "archive.h" #include "utils.h" #include "relational.h" +#include "compiler.h" 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; @@ -161,6 +163,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 +188,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 +239,23 @@ 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::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); } @@ -467,8 +492,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 +502,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)) { @@ -718,9 +772,12 @@ tinfo_t power::return_type_tinfo() const 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); @@ -943,7 +1000,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 { @@ -953,8 +1010,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());