* Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
/*
- * GiNaC Copyright (C) 1999-2008 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2009 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
//////////
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) {
return basis.return_type();
}
-tinfo_t power::return_type_tinfo() const
+return_type_t power::return_type_tinfo() const
{
return basis.return_type_tinfo();
}
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]);
return result;
}
+GINAC_BIND_UNARCHIVER(power);
+
} // namespace GiNaC