* Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
/*
- * GiNaC Copyright (C) 1999-2007 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2015 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
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
-#include <vector>
-#include <iostream>
-#include <stdexcept>
-#include <limits>
-
#include "power.h"
#include "expairseq.h"
#include "add.h"
#include "relational.h"
#include "compiler.h"
+#include <iostream>
+#include <limits>
+#include <stdexcept>
+#include <vector>
+
namespace GiNaC {
GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(power, basic,
// default constructor
//////////
-power::power() : inherited(&power::tinfo_static) { }
+power::power() { }
//////////
// other constructors
// 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);
}
n.add_ex("exponent", exponent);
}
-DEFAULT_UNARCHIVE(power)
-
//////////
// functions overriding virtual functions from base classes
//////////
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::positive) && exponent.info(info_flags::real)) ||
+ (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;
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
* - ^(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)
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<numeric>(ebasis)) {
- basis_is_numerical = true;
num_basis = &ex_to<numeric>(ebasis);
}
if (is_exactly_a<numeric>(eexponent)) {
- exponent_is_numerical = true;
num_exponent = &ex_to<numeric>(eexponent);
}
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())
if (is_exactly_a<power>(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) {
}
// ^(^(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<power>(ebasis)) {
const power & sub_power = ex_to<power>(ebasis);
if (is_exactly_a<numeric>(sub_exponent)) {
const numeric & num_sub_exponent = ex_to<numeric>(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));
}
}
// (2*x + 6*y)^(-4) -> 1/16*(x + 3*y)^(-4)
if (num_exponent->is_integer() && is_exactly_a<add>(ebasis)) {
numeric icont = ebasis.integer_content();
- const numeric& lead_coeff =
- ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div_dyn(icont);
+ const numeric lead_coeff =
+ ex_to<numeric>(ex_to<add>(ebasis).seq.begin()->coeff).div(icont);
const bool canonicalizable = lead_coeff.is_integer();
const bool unit_normal = lead_coeff.is_pos_integer();
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),
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),
}
// 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
{
for (exmap::const_iterator it = m.begin(); it != m.end(); ++it) {
int nummatches = std::numeric_limits<int>::max();
- lst repls;
- if (tryfactsubs(*this, it->first, nummatches, repls))
- return (ex_to<basic>((*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<basic>(result)).subs_one_level(m, options);
+ }
}
return subs_one_level(m, options);
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
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
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<symbol>(basis) && exponent.info(info_flags::integer)) {
+ // A special case worth optimizing.
+ setflag(status_flags::expanded);
return *this;
-
+ }
+
+ // (x*p)^c -> x^c * p^c, if p>0
+ // makes sense before expanding the basis
+ if (is_exactly_a<mul>(basis) && !basis.info(info_flags::indefinite)) {
+ const mul &m = ex_to<mul>(basis);
+ exvector prodseq;
+ epvector 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;
+
+ // search for positive/negative factors
+ while (cit!=last) {
+ ex e=m.recombine_pair_to_ex(*cit);
+ if (e.info(info_flags::positive))
+ prodseq.push_back(pow(e, exponent).expand(options));
+ else if (e.info(info_flags::negative)) {
+ prodseq.push_back(pow(-e, exponent).expand(options));
+ possign = !possign;
+ } else
+ powseq.push_back(*cit);
+ ++cit;
+ }
+
+ // take care on the numeric coefficient
+ ex coeff=(possign? _ex1 : _ex_1);
+ if (m.overall_coeff.info(info_flags::positive) && m.overall_coeff != _ex1)
+ prodseq.push_back(power(m.overall_coeff, exponent));
+ else if (m.overall_coeff.info(info_flags::negative) && m.overall_coeff != _ex_1)
+ prodseq.push_back(power(-m.overall_coeff, exponent));
+ else
+ coeff *= m.overall_coeff;
+
+ // If positive/negative factors are found, then extract them.
+ // In either case we set a flag to avoid the second run on a part
+ // which does not have positive/negative terms.
+ if (prodseq.size() > 0) {
+ ex newbasis = coeff*mul(powseq);
+ ex_to<basic>(newbasis).setflag(status_flags::purely_indefinite);
+ return ((new mul(prodseq))->setflag(status_flags::dynallocated)*(new power(newbasis, exponent))->setflag(status_flags::dynallocated).expand(options)).expand(options);
+ } else
+ ex_to<basic>(basis).setflag(status_flags::purely_indefinite);
+ }
+
const ex expanded_basis = basis.expand(options);
const ex expanded_exponent = exponent.expand(options);
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<m-1; ++l) {
k[l] = 0;
while (true) {
exvector term;
term.reserve(m+1);
- for (l=0; l<m-1; ++l) {
+ for (std::size_t l = 0; l < m - 1; ++l) {
const ex & b = a.op(l);
GINAC_ASSERT(!is_exactly_a<add>(b));
GINAC_ASSERT(!is_exactly_a<power>(b) ||
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<add>(b));
GINAC_ASSERT(!is_exactly_a<power>(b) ||
!is_exactly_a<numeric>(ex_to<power>(b).exponent) ||
term.push_back(power(b,n-k_cum[m-2]));
numeric f = binomial(numeric(n),numeric(k[0]));
- for (l=1; l<m-1; ++l)
+ for (std::size_t l = 1; l < m - 1; ++l)
f *= binomial(numeric(n-k_cum[l-1]),numeric(k[l]));
term.push_back(f);
result.push_back(ex((new mul(term))->setflag(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]);
if (c.is_equal(_ex1)) {
if (is_exactly_a<mul>(r)) {
- sum.push_back(expair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
- _ex1));
+ sum.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r), *_num2_p, options, true),
+ _ex1));
} else {
- sum.push_back(expair((new power(r,_ex2))->setflag(status_flags::dynallocated),
- _ex1));
+ sum.push_back(a.combine_ex_with_coeff_to_pair((new power(r,_ex2))->setflag(status_flags::dynallocated),
+ _ex1));
}
} else {
if (is_exactly_a<mul>(r)) {
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());
return result;
}
+GINAC_BIND_UNARCHIVER(power);
+
} // namespace GiNaC
git://www.ginac.de / ginac.git / blobdiff