* 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
* 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::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));
}
}
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),
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
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 result;
}
+GINAC_BIND_UNARCHIVER(power);
+
} // namespace GiNaC
git://www.ginac.de / ginac.git / blobdiff