* Implementation of GiNaC's products of expressions. */
/*
- * GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2002 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., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
+#include <iostream>
#include <vector>
#include <stdexcept>
#include "power.h"
#include "matrix.h"
#include "archive.h"
-#include "debugmsg.h"
#include "utils.h"
namespace GiNaC {
GINAC_IMPLEMENT_REGISTERED_CLASS(mul, expairseq)
//////////
-// default ctor, dctor, copy ctor assignment operator and helpers
+// default ctor, dtor, copy ctor, assignment operator and helpers
//////////
mul::mul()
{
- debugmsg("mul default ctor",LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_mul;
}
mul::mul(const ex & lh, const ex & rh)
{
- debugmsg("mul ctor from ex,ex",LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_mul;
- overall_coeff = _ex1();
+ overall_coeff = _ex1;
construct_from_2_ex(lh,rh);
GINAC_ASSERT(is_canonical());
}
mul::mul(const exvector & v)
{
- debugmsg("mul ctor from exvector",LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_mul;
- overall_coeff = _ex1();
+ overall_coeff = _ex1;
construct_from_exvector(v);
GINAC_ASSERT(is_canonical());
}
mul::mul(const epvector & v)
{
- debugmsg("mul ctor from epvector",LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_mul;
- overall_coeff = _ex1();
+ overall_coeff = _ex1;
construct_from_epvector(v);
GINAC_ASSERT(is_canonical());
}
mul::mul(const epvector & v, const ex & oc)
{
- debugmsg("mul ctor from epvector,ex",LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_mul;
overall_coeff = oc;
construct_from_epvector(v);
mul::mul(epvector * vp, const ex & oc)
{
- debugmsg("mul ctor from epvector *,ex",LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_mul;
GINAC_ASSERT(vp!=0);
overall_coeff = oc;
mul::mul(const ex & lh, const ex & mh, const ex & rh)
{
- debugmsg("mul ctor from ex,ex,ex",LOGLEVEL_CONSTRUCT);
tinfo_key = TINFO_mul;
exvector factors;
factors.reserve(3);
factors.push_back(lh);
factors.push_back(mh);
factors.push_back(rh);
- overall_coeff = _ex1();
+ overall_coeff = _ex1;
construct_from_exvector(factors);
GINAC_ASSERT(is_canonical());
}
void mul::print(const print_context & c, unsigned level) const
{
- debugmsg("mul print", LOGLEVEL_PRINT);
-
if (is_a<print_tree>(c)) {
inherited::print(c, level);
if (precedence() <= level)
c.s << "(";
- if (!overall_coeff.is_equal(_ex1())) {
+ if (!overall_coeff.is_equal(_ex1)) {
overall_coeff.print(c, precedence());
c.s << "*";
}
while (it != itend) {
// If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
- if (it == seq.begin() && ex_to<numeric>(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0) {
+ if (it == seq.begin() && ex_to<numeric>(it->coeff).is_integer() && it->coeff.info(info_flags::negative)) {
if (is_a<print_csrc_cl_N>(c))
c.s << "recip(";
else
}
// If the exponent is 1 or -1, it is left out
- if (it->coeff.compare(_ex1()) == 0 || it->coeff.compare(_num_1()) == 0)
+ if (it->coeff.is_equal(_ex1) || it->coeff.is_equal(_ex_1))
it->rest.print(c, precedence());
else {
// Outer parens around ex needed for broken gcc-2.95 parser:
// Separator is "/" for negative integer powers, "*" otherwise
++it;
if (it != itend) {
- if (ex_to<numeric>(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0)
+ if (ex_to<numeric>(it->coeff).is_integer() && it->coeff.info(info_flags::negative))
c.s << "/";
else
c.s << "*";
if (precedence() <= level)
c.s << ")";
+ } else if (is_a<print_python_repr>(c)) {
+ c.s << class_name() << '(';
+ op(0).print(c);
+ for (unsigned i=1; i<nops(); ++i) {
+ c.s << ',';
+ op(i).print(c);
+ }
+ c.s << ')';
} else {
if (precedence() <= level) {
numeric coeff = ex_to<numeric>(overall_coeff);
if (coeff.csgn() == -1)
c.s << '-';
- if (!coeff.is_equal(_num1()) &&
- !coeff.is_equal(_num_1())) {
+ if (!coeff.is_equal(_num1) &&
+ !coeff.is_equal(_num_1)) {
if (coeff.is_rational()) {
if (coeff.is_negative())
(-coeff).print(c);
return (new mul(coeffseq))->setflag(status_flags::dynallocated);
}
- return _ex0();
+ return _ex0;
}
/** Perform automatic term rewriting rules in this class. In the following
* @param level cut-off in recursive evaluation */
ex mul::eval(int level) const
{
- debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
-
epvector *evaled_seqp = evalchildren(level);
if (evaled_seqp) {
// do more evaluation later
#ifdef DO_GINAC_ASSERT
epvector::const_iterator i = seq.begin(), end = seq.end();
while (i != end) {
- GINAC_ASSERT((!is_ex_exactly_of_type(i->rest, mul)) ||
+ GINAC_ASSERT((!is_exactly_a<mul>(i->rest)) ||
(!(ex_to<numeric>(i->coeff).is_integer())));
GINAC_ASSERT(!(i->is_canonical_numeric()));
if (is_ex_exactly_of_type(recombine_pair_to_ex(*i), numeric))
print(print_tree(std::cerr));
- GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*i), numeric));
+ GINAC_ASSERT(!is_exactly_a<numeric>(recombine_pair_to_ex(*i)));
/* for paranoia */
expair p = split_ex_to_pair(recombine_pair_to_ex(*i));
GINAC_ASSERT(p.rest.is_equal(i->rest));
if (flags & status_flags::evaluated) {
GINAC_ASSERT(seq.size()>0);
- GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1()));
+ GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1));
return *this;
}
int seq_size = seq.size();
if (overall_coeff.is_zero()) {
// *(...,x;0) -> 0
- return _ex0();
+ return _ex0;
} else if (seq_size==0) {
// *(;c) -> c
return overall_coeff;
- } else if (seq_size==1 && overall_coeff.is_equal(_ex1())) {
+ } else if (seq_size==1 && overall_coeff.is_equal(_ex1)) {
// *(x;1) -> x
return recombine_pair_to_ex(*(seq.begin()));
} else if ((seq_size==1) &&
is_ex_exactly_of_type((*seq.begin()).rest,add) &&
- ex_to<numeric>((*seq.begin()).coeff).is_equal(_num1())) {
+ ex_to<numeric>((*seq.begin()).coeff).is_equal(_num1)) {
// *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
const add & addref = ex_to<add>((*seq.begin()).rest);
epvector *distrseq = new epvector();
ex mul::evalm(void) const
{
// numeric*matrix
- if (seq.size() == 1 && seq[0].coeff.is_equal(_ex1())
+ if (seq.size() == 1 && seq[0].coeff.is_equal(_ex1)
&& is_ex_of_type(seq[0].rest, matrix))
return ex_to<matrix>(seq[0].rest).mul(ex_to<numeric>(overall_coeff));
epvector::const_iterator i = seq.begin(), end = seq.end();
epvector::iterator i2 = mulseq.begin();
while (i != end) {
- expair ep = split_ex_to_pair(power(i->rest, i->coeff - _ex1()) *
+ expair ep = split_ex_to_pair(power(i->rest, i->coeff - _ex1) *
i->rest.diff(s));
ep.swap(*i2);
addseq.push_back((new mul(mulseq, overall_coeff * i->coeff))->setflag(status_flags::dynallocated));
if (is_ex_exactly_of_type(powerref.exponent,numeric))
return expair(powerref.basis,powerref.exponent);
}
- return expair(e,_ex1());
+ return expair(e,_ex1);
}
expair mul::combine_ex_with_coeff_to_pair(const ex & e,
// we create a temporary power object
// otherwise it would be hard to correctly simplify
// expression like (4^(1/3))^(3/2)
- if (are_ex_trivially_equal(c,_ex1()))
+ if (are_ex_trivially_equal(c,_ex1))
return split_ex_to_pair(e);
return split_ex_to_pair(power(e,c));
// we create a temporary power object
// otherwise it would be hard to correctly simplify
// expression like (4^(1/3))^(3/2)
- if (are_ex_trivially_equal(c,_ex1()))
+ if (are_ex_trivially_equal(c,_ex1))
return p;
return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
ex mul::recombine_pair_to_ex(const expair & p) const
{
- if (ex_to<numeric>(p.coeff).is_equal(_num1()))
+ if (ex_to<numeric>(p.coeff).is_equal(_num1))
return p.rest;
else
return power(p.rest,p.coeff);
*it = ep;
return true;
}
- if (ex_to<numeric>((*it).coeff).is_equal(_num1())) {
+ if (ex_to<numeric>((*it).coeff).is_equal(_num1)) {
// combined pair has coeff 1 and must be moved to the end
return true;
}
ex mul::default_overall_coeff(void) const
{
- return _ex1();
+ return _ex1;
}
void mul::combine_overall_coeff(const ex & c)
{
- GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
- GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
+ GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
+ GINAC_ASSERT(is_exactly_a<numeric>(c));
overall_coeff = ex_to<numeric>(overall_coeff).mul_dyn(ex_to<numeric>(c));
}
void mul::combine_overall_coeff(const ex & c1, const ex & c2)
{
- GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
- GINAC_ASSERT(is_ex_exactly_of_type(c1,numeric));
- GINAC_ASSERT(is_ex_exactly_of_type(c2,numeric));
+ GINAC_ASSERT(is_exactly_a<numeric>(overall_coeff));
+ GINAC_ASSERT(is_exactly_a<numeric>(c1));
+ GINAC_ASSERT(is_exactly_a<numeric>(c2));
overall_coeff = ex_to<numeric>(overall_coeff).mul_dyn(ex_to<numeric>(c1).power(ex_to<numeric>(c2)));
}
bool mul::can_make_flat(const expair & p) const
{
- GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
+ GINAC_ASSERT(is_exactly_a<numeric>(p.coeff));
// this assertion will probably fail somewhere
// it would require a more careful make_flat, obeying the power laws
// probably should return true only if p.coeff is integer
- return ex_to<numeric>(p.coeff).is_equal(_num1());
+ return ex_to<numeric>(p.coeff).is_equal(_num1);
}
ex mul::expand(unsigned options) const
// with the next one that is found while collecting the factors which are
// not sums
int number_of_adds = 0;
- ex last_expanded = _ex1();
+ ex last_expanded = _ex1;
epvector non_adds;
non_adds.reserve(expanded_seq.size());
epvector::const_iterator cit = expanded_seq.begin(), last = expanded_seq.end();
while (cit != last) {
if (is_ex_exactly_of_type(cit->rest, add) &&
- (cit->coeff.is_equal(_ex1()))) {
+ (cit->coeff.is_equal(_ex1))) {
++number_of_adds;
if (is_ex_exactly_of_type(last_expanded, add)) {
+#if 0
+ // Expand a product of two sums, simple and robust version.
const add & add1 = ex_to<add>(last_expanded);
const add & add2 = ex_to<add>(cit->rest);
- int n1 = add1.nops();
- int n2 = add2.nops();
+ const int n1 = add1.nops();
+ const int n2 = add2.nops();
+ ex tmp_accu;
exvector distrseq;
- distrseq.reserve(n1*n2);
+ distrseq.reserve(n2);
for (int i1=0; i1<n1; ++i1) {
- for (int i2=0; i2<n2; ++i2) {
- distrseq.push_back(add1.op(i1) * add2.op(i2));
+ distrseq.clear();
+ // cache the first operand (for efficiency):
+ const ex op1 = add1.op(i1);
+ for (int i2=0; i2<n2; ++i2)
+ distrseq.push_back(op1 * add2.op(i2));
+ tmp_accu += (new add(distrseq))->
+ setflag(status_flags::dynallocated);
+ }
+ last_expanded = tmp_accu;
+#else
+ // Expand a product of two sums, aggressive version.
+ // Caring for the overall coefficients in separate loops can
+ // sometimes give a performance gain of up to 15%!
+
+ const int sizedifference = ex_to<add>(last_expanded).seq.size()-ex_to<add>(cit->rest).seq.size();
+ // add2 is for the inner loop and should be the bigger of the two sums
+ // in the presence of asymptotically good sorting:
+ const add& add1 = (sizedifference<0 ? ex_to<add>(last_expanded) : ex_to<add>(cit->rest));
+ const add& add2 = (sizedifference<0 ? ex_to<add>(cit->rest) : ex_to<add>(last_expanded));
+ const epvector::const_iterator add1begin = add1.seq.begin();
+ const epvector::const_iterator add1end = add1.seq.end();
+ const epvector::const_iterator add2begin = add2.seq.begin();
+ const epvector::const_iterator add2end = add2.seq.end();
+ epvector distrseq;
+ distrseq.reserve(add1.seq.size()+add2.seq.size());
+ // Multiply add2 with the overall coefficient of add1 and append it to distrseq:
+ if (!add1.overall_coeff.is_zero()) {
+ if (add1.overall_coeff.is_equal(_ex1))
+ distrseq.insert(distrseq.end(),add2begin,add2end);
+ else
+ for (epvector::const_iterator i=add2begin; i!=add2end; ++i)
+ distrseq.push_back(expair(i->rest, ex_to<numeric>(i->coeff).mul_dyn(ex_to<numeric>(add1.overall_coeff))));
+ }
+ // Multiply add1 with the overall coefficient of add2 and append it to distrseq:
+ if (!add2.overall_coeff.is_zero()) {
+ if (add2.overall_coeff.is_equal(_ex1))
+ distrseq.insert(distrseq.end(),add1begin,add1end);
+ else
+ for (epvector::const_iterator i=add1begin; i!=add1end; ++i)
+ distrseq.push_back(expair(i->rest, ex_to<numeric>(i->coeff).mul_dyn(ex_to<numeric>(add2.overall_coeff))));
+ }
+ // Compute the new overall coefficient and put it together:
+ ex tmp_accu = (new add(distrseq, add1.overall_coeff*add2.overall_coeff))->setflag(status_flags::dynallocated);
+ // Multiply explicitly all non-numeric terms of add1 and add2:
+ for (epvector::const_iterator i1=add1begin; i1!=add1end; ++i1) {
+ // We really have to combine terms here in order to compactify
+ // the result. Otherwise it would become waayy tooo bigg.
+ numeric oc;
+ distrseq.clear();
+ for (epvector::const_iterator i2=add2begin; i2!=add2end; ++i2) {
+ // Don't push_back expairs which might have a rest that evaluates to a numeric,
+ // since that would violate an invariant of expairseq:
+ const ex rest = (new mul(i1->rest, i2->rest))->setflag(status_flags::dynallocated);
+ if (is_ex_exactly_of_type(rest, numeric))
+ oc += ex_to<numeric>(rest).mul(ex_to<numeric>(i1->coeff).mul(ex_to<numeric>(i2->coeff)));
+ else
+ distrseq.push_back(expair(rest, ex_to<numeric>(i1->coeff).mul_dyn(ex_to<numeric>(i2->coeff))));
}
+ tmp_accu += (new add(distrseq, oc))->setflag(status_flags::dynallocated);
}
- last_expanded = (new add(distrseq))->
- setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
+ last_expanded = tmp_accu;
+#endif
} else {
non_adds.push_back(split_ex_to_pair(last_expanded));
last_expanded = cit->rest;
* pointer, if sequence is unchanged. */
epvector * mul::expandchildren(unsigned options) const
{
- epvector::const_iterator last = seq.end();
+ const epvector::const_iterator last = seq.end();
epvector::const_iterator cit = seq.begin();
while (cit!=last) {
const ex & factor = recombine_pair_to_ex(*cit);
git://www.ginac.de / ginac.git / blobdiff