* Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
/*
- * GiNaC Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2006 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
#include "lst.h"
#include "archive.h"
#include "utils.h"
+#include "relational.h"
namespace GiNaC {
// default constructor
//////////
-power::power() : inherited(TINFO_power) { }
+power::power() : inherited(&power::tinfo_static) { }
//////////
// other constructors
return *this;
}
+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);
+}
+
int power::degree(const ex & s) const
{
if (is_equal(ex_to<basic>(s)))
if (ebasis.is_equal(_ex1))
return _ex1;
+ // power of a function calculated by separate rules defined for this function
+ if (is_exactly_a<function>(ebasis))
+ return ex_to<function>(ebasis).power(eexponent);
+
if (exponent_is_numerical) {
// ^(c1,c2) -> c1^c2 (c1, c2 numeric(),
return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
}
+bool power::has(const ex & other, unsigned options) const
+{
+ if (!(options & has_options::algebraic))
+ return basic::has(other, options);
+ if (!is_a<power>(other))
+ return basic::has(other, options);
+ if (!exponent.info(info_flags::integer)
+ || !other.op(1).info(info_flags::integer))
+ return basic::has(other, options);
+ if (exponent.info(info_flags::posint)
+ && other.op(1).info(info_flags::posint)
+ && ex_to<numeric>(exponent).to_int()
+ > ex_to<numeric>(other.op(1)).to_int()
+ && basis.match(other.op(0)))
+ return true;
+ if (exponent.info(info_flags::negint)
+ && other.op(1).info(info_flags::negint)
+ && ex_to<numeric>(exponent).to_int()
+ < ex_to<numeric>(other.op(1)).to_int()
+ && basis.match(other.op(0)))
+ return true;
+ return basic::has(other, options);
+}
+
// from mul.cpp
extern bool tryfactsubs(const ex &, const ex &, int &, lst &);
return (new power(newbasis, newexponent))->setflag(status_flags::dynallocated);
}
+ex power::real_part() const
+{
+ if (exponent.info(info_flags::integer)) {
+ ex basis_real = basis.real_part();
+ if (basis_real == basis)
+ return *this;
+ realsymbol a("a"),b("b");
+ ex result;
+ if (exponent.info(info_flags::posint))
+ result = power(a+I*b,exponent);
+ else
+ result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
+ result = result.expand();
+ result = result.real_part();
+ result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+ return result;
+ }
+
+ ex a = basis.real_part();
+ 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))*cos(c*atan2(b,a)+d*log(abs(basis)));
+}
+
+ex power::imag_part() const
+{
+ if (exponent.info(info_flags::integer)) {
+ ex basis_real = basis.real_part();
+ if (basis_real == basis)
+ return 0;
+ realsymbol a("a"),b("b");
+ ex result;
+ if (exponent.info(info_flags::posint))
+ result = power(a+I*b,exponent);
+ else
+ result = power(a/(a*a+b*b)-I*b/(a*a+b*b),-exponent);
+ result = result.expand();
+ result = result.imag_part();
+ result = result.subs(lst( a==basis_real, b==basis.imag_part() ));
+ return result;
+ }
+
+ ex a=basis.real_part();
+ 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)));
+}
+
+// protected
+
// protected
/** Implementation of ex::diff() for a power.
* @see ex::diff */
ex power::derivative(const symbol & s) const
{
- if (exponent.info(info_flags::real)) {
+ if (is_a<numeric>(exponent)) {
// D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
epvector newseq;
newseq.reserve(2);
return basis.return_type();
}
-unsigned power::return_type_tinfo() const
+tinfo_t power::return_type_tinfo() const
{
return basis.return_type_tinfo();
}
const size_t m = a.nops();
exvector result;
// The number of terms will be the number of combinatorial compositions,
- // i.e. the number of unordered arrangement of m nonnegative integers
+ // i.e. the number of unordered arrangements of m nonnegative integers
// which sum up to n. It is frequently written as C_n(m) and directly
// related with binomial coefficients:
result.reserve(binomial(numeric(n+m-1), numeric(m-1)).to_int());
}
// Leave it to multiplication since dummy indices have to be renamed
- if (get_all_dummy_indices(m).size() > 0) {
+ if (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());
+
for (int i=1; i < n.to_int(); i++)
- result *= rename_dummy_indices_uniquely(m,m);
+ result *= rename_dummy_indices_uniquely(va, m);
return result;
}