X-Git-Url: https://www.ginac.de/ginac.git//ginac.git?p=ginac.git;a=blobdiff_plain;f=ginac%2Fncmul.cpp;h=c4cfe69cb2a93fd3b8afe03f2e571c9a1ba0e34e;hp=cf4f62463218ecb13710ec9070af4f25c9967169;hb=b6e3c62f240698c7e9ed464c57bb6d92741765ba;hpb=9e2d60e206395fc9908e1f9025e50c76b3d7c182 diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp index cf4f6246..c4cfe69c 100644 --- a/ginac/ncmul.cpp +++ b/ginac/ncmul.cpp @@ -3,7 +3,7 @@ * Implementation of GiNaC's non-commutative products of expressions. */ /* - * GiNaC Copyright (C) 1999-2000 Johannes Gutenberg University Mainz, Germany + * GiNaC Copyright (C) 1999-2001 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 @@ -28,13 +28,13 @@ #include "ex.h" #include "add.h" #include "mul.h" +#include "matrix.h" +#include "print.h" #include "archive.h" #include "debugmsg.h" #include "utils.h" -#ifndef NO_NAMESPACE_GINAC namespace GiNaC { -#endif // ndef NO_NAMESPACE_GINAC GINAC_IMPLEMENT_REGISTERED_CLASS(ncmul, exprseq) @@ -42,47 +42,14 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(ncmul, exprseq) // default constructor, destructor, copy constructor assignment operator and helpers ////////// -// public - ncmul::ncmul() { debugmsg("ncmul default constructor",LOGLEVEL_CONSTRUCT); tinfo_key = TINFO_ncmul; } -ncmul::~ncmul() -{ - debugmsg("ncmul destructor",LOGLEVEL_DESTRUCT); - destroy(false); -} - -ncmul::ncmul(const ncmul & other) -{ - debugmsg("ncmul copy constructor",LOGLEVEL_CONSTRUCT); - copy(other); -} - -const ncmul & ncmul::operator=(const ncmul & other) -{ - debugmsg("ncmul operator=",LOGLEVEL_ASSIGNMENT); - if (this != &other) { - destroy(true); - copy(other); - } - return *this; -} - -// protected - -void ncmul::copy(const ncmul & other) -{ - inherited::copy(other); -} - -void ncmul::destroy(bool call_parent) -{ - if (call_parent) inherited::destroy(call_parent); -} +DEFAULT_COPY(ncmul) +DEFAULT_DESTROY(ncmul) ////////// // other constructors @@ -139,24 +106,7 @@ ncmul::ncmul(exvector * vp) : inherited(vp) // archiving ////////// -/** Construct object from archive_node. */ -ncmul::ncmul(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst) -{ - debugmsg("ncmul constructor from archive_node", LOGLEVEL_CONSTRUCT); -} - -/** Unarchive the object. */ -ex ncmul::unarchive(const archive_node &n, const lst &sym_lst) -{ - return (new ncmul(n, sym_lst))->setflag(status_flags::dynallocated); -} - -/** Archive the object. */ -void ncmul::archive(archive_node &n) const -{ - inherited::archive(n); -} - +DEFAULT_ARCHIVING(ncmul) ////////// // functions overriding virtual functions from bases classes @@ -164,42 +114,28 @@ void ncmul::archive(archive_node &n) const // public -basic * ncmul::duplicate() const +void ncmul::print(const print_context & c, unsigned level) const { - debugmsg("ncmul duplicate",LOGLEVEL_ASSIGNMENT); - return new ncmul(*this); -} + debugmsg("ncmul print", LOGLEVEL_PRINT); -void ncmul::print(std::ostream & os, unsigned upper_precedence) const -{ - debugmsg("ncmul print",LOGLEVEL_PRINT); - printseq(os,'(','%',')',precedence,upper_precedence); -} + if (is_of_type(c, print_tree)) { -void ncmul::printraw(std::ostream & os) const -{ - debugmsg("ncmul printraw",LOGLEVEL_PRINT); - os << "%("; - for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) { - (*it).bp->printraw(os); - os << ","; - } - os << ",hash=" << hashvalue << ",flags=" << flags; - os << ")"; -} + inherited::print(c, level); -void ncmul::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const -{ - debugmsg("ncmul print csrc",LOGLEVEL_PRINT); - exvector::const_iterator it; - exvector::const_iterator itend = seq.end()-1; - os << "ncmul("; - for (it=seq.begin(); it!=itend; ++it) { - (*it).bp->printcsrc(os,precedence); - os << ","; - } - (*it).bp->printcsrc(os,precedence); - os << ")"; + } else if (is_of_type(c, print_csrc)) { + + c.s << "ncmul("; + exvector::const_iterator it = seq.begin(), itend = seq.end()-1; + while (it != itend) { + it->print(c, precedence()); + c.s << ","; + it++; + } + it->print(c, precedence()); + c.s << ")"; + + } else + printseq(c, '(', '*', ')', precedence(), level); } bool ncmul::info(unsigned inf) const @@ -211,95 +147,101 @@ typedef std::vector intvector; ex ncmul::expand(unsigned options) const { - exvector sub_expanded_seq; - intvector positions_of_adds; - intvector number_of_add_operands; - - exvector expanded_seq=expandchildren(options); + // First, expand the children + exvector expanded_seq = expandchildren(options); - positions_of_adds.resize(expanded_seq.size()); - number_of_add_operands.resize(expanded_seq.size()); + // Now, look for all the factors that are sums and remember their + // position and number of terms. One remark is in order here: we do not + // take into account the overall_coeff of the add objects. This is + // because in GiNaC, all terms of a sum must be of the same type, so + // a non-zero overall_coeff (which can only be numeric) would imply that + // the sum only has commutative terms. But then it would never appear + // as a factor of an ncmul. + intvector positions_of_adds(expanded_seq.size()); + intvector number_of_add_operands(expanded_seq.size()); - int number_of_adds=0; - int number_of_expanded_terms=1; + int number_of_adds = 0; + int number_of_expanded_terms = 1; - unsigned current_position=0; - exvector::const_iterator last=expanded_seq.end(); + unsigned current_position = 0; + exvector::const_iterator last = expanded_seq.end(); for (exvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) { - if (is_ex_exactly_of_type((*cit),add)) { - positions_of_adds[number_of_adds]=current_position; - const add & expanded_addref=ex_to_add(*cit); - number_of_add_operands[number_of_adds]=expanded_addref.seq.size(); + if (is_ex_exactly_of_type(*cit, add)) { + positions_of_adds[number_of_adds] = current_position; + const add & expanded_addref = ex_to(*cit); + number_of_add_operands[number_of_adds] = expanded_addref.seq.size(); number_of_expanded_terms *= expanded_addref.seq.size(); number_of_adds++; } current_position++; } - if (number_of_adds==0) { - return (new ncmul(expanded_seq,1))->setflag(status_flags::dynallocated || - status_flags::expanded); - } + // If there are no sums, we are done + if (number_of_adds == 0) + return (new ncmul(expanded_seq, true))-> + setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)); + // Now, form all possible products of the terms of the sums with the + // remaining factors, and add them together exvector distrseq; distrseq.reserve(number_of_expanded_terms); - intvector k; - k.resize(number_of_adds); - - int l; - for (l=0; l(expanded_seq[positions_of_adds[i]]); + term[positions_of_adds[i]] = addref.recombine_pair_to_ex(addref.seq[k[i]]); } - distrseq.push_back((new ncmul(term,1))->setflag(status_flags::dynallocated | - status_flags::expanded)); + distrseq.push_back((new ncmul(term, true))-> + setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0))); // increment k[] - l=number_of_adds-1; - while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) { - k[l]=0; + int l = number_of_adds-1; + while ((l>=0) && ((++k[l]) >= number_of_add_operands[l])) { + k[l] = 0; l--; } - if (l<0) break; + if (l<0) + break; } - return (new add(distrseq))->setflag(status_flags::dynallocated | - status_flags::expanded); + return (new add(distrseq))-> + setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0)); } -int ncmul::degree(const symbol & s) const +int ncmul::degree(const ex & s) const { - int deg_sum=0; - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - deg_sum+=(*cit).degree(s); + // Sum up degrees of factors + int deg_sum = 0; + exvector::const_iterator i = seq.begin(), end = seq.end(); + while (i != end) { + deg_sum += i->degree(s); + ++i; } return deg_sum; } -int ncmul::ldegree(const symbol & s) const +int ncmul::ldegree(const ex & s) const { - int deg_sum=0; - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - deg_sum+=(*cit).ldegree(s); + // Sum up degrees of factors + int deg_sum = 0; + exvector::const_iterator i = seq.begin(), end = seq.end(); + while (i != end) { + deg_sum += i->degree(s); + ++i; } return deg_sum; } -ex ncmul::coeff(const symbol & s, int n) const +ex ncmul::coeff(const ex & s, int n) const { exvector coeffseq; coeffseq.reserve(seq.size()); - if (n==0) { + if (n == 0) { // product of individual coeffs // if a non-zero power of s is found, the resulting product will be 0 exvector::const_iterator it=seq.begin(); @@ -310,17 +252,17 @@ ex ncmul::coeff(const symbol & s, int n) const return (new ncmul(coeffseq,1))->setflag(status_flags::dynallocated); } - exvector::const_iterator it=seq.begin(); - bool coeff_found=0; - while (it!=seq.end()) { - ex c=(*it).coeff(s,n); - if (!c.is_zero()) { - coeffseq.push_back(c); - coeff_found=1; + exvector::const_iterator i = seq.begin(), end = seq.end(); + bool coeff_found = false; + while (i != end) { + ex c = i->coeff(s,n); + if (c.is_zero()) { + coeffseq.push_back(*i); } else { - coeffseq.push_back(*it); + coeffseq.push_back(c); + coeff_found = true; } - ++it; + ++i; } if (coeff_found) return (new ncmul(coeffseq,1))->setflag(status_flags::dynallocated); @@ -347,10 +289,8 @@ void ncmul::append_factors(exvector & v, const ex & e) const (is_ex_exactly_of_type(e,ncmul))) { for (unsigned i=0; i unsignedvector; @@ -366,7 +306,7 @@ ex ncmul::eval(int level) const // *(c1,c2,ncmul(...)) (pull out commutative elements) // ncmul(x1,y1,x2,y2) -> *(ncmul(x1,x2),ncmul(y1,y2)) // (collect elements of same type) - // ncmul(x1,x2,x3,...) -> x::eval_ncmul(x1,x2,x3,...) + // ncmul(x1,x2,x3,...) -> x::simplify_ncmul(x1,x2,x3,...) // the following rule would be nice, but produces a recursion, // which must be trapped by introducing a flag that the sub-ncmuls() // are already evaluated (maybe later...) @@ -374,7 +314,7 @@ ex ncmul::eval(int level) const // ncmul(ncmul(x1,x2,...),X,ncmul(y1,y2,...) // (X noncommutative_composite) - if ((level==1)&&(flags & status_flags::evaluated)) { + if ((level==1) && (flags & status_flags::evaluated)) { return *this; } @@ -382,32 +322,33 @@ ex ncmul::eval(int level) const // ncmul(...,*(x1,x2),...,ncmul(x3,x4),...) -> // ncmul(...,x1,x2,...,x3,x4,...) (associativity) - unsigned factors=0; - for (exvector::const_iterator cit=evaledseq.begin(); cit!=evaledseq.end(); ++cit) { - factors += count_factors(*cit); - } - + unsigned factors = 0; + exvector::const_iterator cit = evaledseq.begin(), citend = evaledseq.end(); + while (cit != citend) + factors += count_factors(*cit++); + exvector assocseq; assocseq.reserve(factors); - for (exvector::const_iterator cit=evaledseq.begin(); cit!=evaledseq.end(); ++cit) { - append_factors(assocseq,*cit); - } - + cit = evaledseq.begin(); + while (cit != citend) + append_factors(assocseq, *cit++); + // ncmul(x) -> x if (assocseq.size()==1) return *(seq.begin()); // ncmul() -> 1 - if (assocseq.size()==0) return _ex1(); + if (assocseq.empty()) return _ex1(); // determine return types unsignedvector rettypes; rettypes.reserve(assocseq.size()); - unsigned i=0; + unsigned i = 0; unsigned count_commutative=0; unsigned count_noncommutative=0; unsigned count_noncommutative_composite=0; - for (exvector::const_iterator cit=assocseq.begin(); cit!=assocseq.end(); ++cit) { - switch (rettypes[i]=(*cit).return_type()) { + cit = assocseq.begin(); citend = assocseq.end(); + while (cit != citend) { + switch (rettypes[i] = cit->return_type()) { case return_types::commutative: count_commutative++; break; @@ -420,7 +361,7 @@ ex ncmul::eval(int level) const default: throw(std::logic_error("ncmul::eval(): invalid return type")); } - ++i; + ++i; ++cit; } GINAC_ASSERT(count_commutative+count_noncommutative+count_noncommutative_composite==assocseq.size()); @@ -431,12 +372,12 @@ ex ncmul::eval(int level) const commutativeseq.reserve(count_commutative+1); exvector noncommutativeseq; noncommutativeseq.reserve(assocseq.size()-count_commutative); - for (i=0; isetflag(status_flags::dynallocated)); return (new mul(commutativeseq))->setflag(status_flags::dynallocated); @@ -450,50 +391,52 @@ ex ncmul::eval(int level) const // elements in assocseq GINAC_ASSERT(count_commutative==0); + unsigned assoc_num = assocseq.size(); exvectorvector evv; unsignedvector rttinfos; - evv.reserve(assocseq.size()); - rttinfos.reserve(assocseq.size()); + evv.reserve(assoc_num); + rttinfos.reserve(assoc_num); - for (exvector::const_iterator cit=assocseq.begin(); cit!=assocseq.end(); ++cit) { - unsigned ti=(*cit).return_type_tinfo(); + cit = assocseq.begin(), citend = assocseq.end(); + while (cit != citend) { + unsigned ti = cit->return_type_tinfo(); + unsigned rtt_num = rttinfos.size(); // search type in vector of known types - for (i=0; i=rttinfos.size()) { + if (i >= rtt_num) { // new type rttinfos.push_back(ti); evv.push_back(exvector()); - (*(evv.end()-1)).reserve(assocseq.size()); - (*(evv.end()-1)).push_back(*cit); + (evv.end()-1)->reserve(assoc_num); + (evv.end()-1)->push_back(*cit); } + ++cit; } + unsigned evv_num = evv.size(); #ifdef DO_GINAC_ASSERT - GINAC_ASSERT(evv.size()==rttinfos.size()); - GINAC_ASSERT(evv.size()>0); + GINAC_ASSERT(evv_num == rttinfos.size()); + GINAC_ASSERT(evv_num > 0); unsigned s=0; - for (i=0; isetflag(status_flags::dynallocated)); - } - + return (new mul(splitseq))->setflag(status_flags::dynallocated); } @@ -501,23 +444,34 @@ ex ncmul::eval(int level) const status_flags::evaluated); } -exvector ncmul::get_indices(void) const +ex ncmul::evalm(void) const { - // return union of indices of factors - exvector iv; - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - exvector subiv=(*cit).get_indices(); - iv.reserve(iv.size()+subiv.size()); - for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) { - iv.push_back(*cit2); + // Evaluate children first + exvector *s = new exvector; + s->reserve(seq.size()); + exvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { + s->push_back(it->evalm()); + it++; + } + + // If there are only matrices, simply multiply them + it = s->begin(); itend = s->end(); + if (is_ex_of_type(*it, matrix)) { + matrix prod(ex_to(*it)); + it++; + while (it != itend) { + if (!is_ex_of_type(*it, matrix)) + goto no_matrix; + prod = prod.mul(ex_to(*it)); + it++; } + delete s; + return prod; } - return iv; -} -ex ncmul::subs(const lst & ls, const lst & lr) const -{ - return ncmul(subschildren(ls, lr)); +no_matrix: + return (new ncmul(s))->setflag(status_flags::dynallocated); } ex ncmul::thisexprseq(const exvector & v) const @@ -532,11 +486,24 @@ ex ncmul::thisexprseq(exvector * vp) const // protected -/** Implementation of ex::diff() for a non-commutative product. It always returns 0. +/** Implementation of ex::diff() for a non-commutative product. It applies + * the product rule. * @see ex::diff */ ex ncmul::derivative(const symbol & s) const { - return _ex0(); + unsigned num = seq.size(); + exvector addseq; + addseq.reserve(num); + + // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c) + exvector ncmulseq = seq; + for (unsigned i=0; isetflag(status_flags::dynallocated)); + e.swap(ncmulseq[i]); + } + return (new add(addseq))->setflag(status_flags::dynallocated); } int ncmul::compare_same_type(const basic & other) const @@ -546,30 +513,30 @@ int ncmul::compare_same_type(const basic & other) const unsigned ncmul::return_type(void) const { - if (seq.size()==0) { - // ncmul without factors: should not happen, but commutes + if (seq.empty()) return return_types::commutative; - } - bool all_commutative=1; - unsigned rt; - exvector::const_iterator cit_noncommutative_element; // point to first found nc element + bool all_commutative = true; + exvector::const_iterator noncommutative_element; // point to first found nc element - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - rt=(*cit).return_type(); - if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc - if ((rt==return_types::noncommutative)&&(all_commutative)) { + exvector::const_iterator i = seq.begin(), end = seq.end(); + while (i != end) { + unsigned rt = i->return_type(); + if (rt == return_types::noncommutative_composite) + return rt; // one ncc -> mul also ncc + if ((rt == return_types::noncommutative) && (all_commutative)) { // first nc element found, remember position - cit_noncommutative_element=cit; - all_commutative=0; + noncommutative_element = i; + all_commutative = false; } - if ((rt==return_types::noncommutative)&&(!all_commutative)) { + if ((rt == return_types::noncommutative) && (!all_commutative)) { // another nc element found, compare type_infos - if ((*cit_noncommutative_element).return_type_tinfo()!=(*cit).return_type_tinfo()) { + if (noncommutative_element->return_type_tinfo() != i->return_type_tinfo()) { // diffent types -> mul is ncc return return_types::noncommutative_composite; } } + ++i; } // all factors checked GINAC_ASSERT(!all_commutative); // not all factors should commute, because this is a ncmul(); @@ -578,16 +545,17 @@ unsigned ncmul::return_type(void) const unsigned ncmul::return_type_tinfo(void) const { - if (seq.size()==0) { - // mul without factors: should not happen + if (seq.empty()) return tinfo_key; - } + // return type_info of first noncommutative element - for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) { - if ((*cit).return_type()==return_types::noncommutative) { - return (*cit).return_type_tinfo(); - } + exvector::const_iterator i = seq.begin(), end = seq.end(); + while (i != end) { + if (i->return_type() == return_types::noncommutative) + return i->return_type_tinfo(); + ++i; } + // no noncommutative element found, should not happen return tinfo_key; } @@ -606,9 +574,10 @@ exvector ncmul::expandchildren(unsigned options) const { exvector s; s.reserve(seq.size()); - - for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) { - s.push_back((*it).expand(options)); + exvector::const_iterator it = seq.begin(), itend = seq.end(); + while (it != itend) { + s.push_back(it->expand(options)); + it++; } return s; } @@ -618,22 +587,6 @@ const exvector & ncmul::get_factors(void) const return seq; } -////////// -// static member variables -////////// - -// protected - -unsigned ncmul::precedence=50; - - -////////// -// global constants -////////// - -const ncmul some_ncmul; -const std::type_info & typeid_ncmul = typeid(some_ncmul); - ////////// // friend functions ////////// @@ -645,15 +598,13 @@ ex nonsimplified_ncmul(const exvector & v) ex simplified_ncmul(const exvector & v) { - if (v.size()==0) { + if (v.empty()) return _ex1(); - } else if (v.size()==1) { + else if (v.size() == 1) return v[0]; - } - return (new ncmul(v))->setflag(status_flags::dynallocated | - status_flags::evaluated); + else + return (new ncmul(v))->setflag(status_flags::dynallocated | + status_flags::evaluated); } -#ifndef NO_NAMESPACE_GINAC } // namespace GiNaC -#endif // ndef NO_NAMESPACE_GINAC