index e775d82e9e2d85207dd1330bd494d801d195ab9e..c4cfe69cb2a93fd3b8afe03f2e571c9a1ba0e34e 100644 (file)
@@ -28,6 +28,7 @@
#include "ex.h"
#include "mul.h"
+#include "matrix.h"
#include "print.h"
#include "archive.h"
#include "debugmsg.h"
@@ -126,15 +127,15 @@ void ncmul::print(const print_context & c, unsigned level) const
c.s << "ncmul(";
exvector::const_iterator it = seq.begin(), itend = seq.end()-1;
while (it != itend) {
-                       it->print(c, precedence);
+                       it->print(c, precedence());
c.s << ",";
it++;
}
-               it->print(c, precedence);
+               it->print(c, precedence());
c.s << ")";

} else
-               printseq(c, '(', '*', ')', precedence, level);
+               printseq(c, '(', '*', ')', precedence(), level);
}

bool ncmul::info(unsigned inf) const
@@ -146,85 +147,91 @@ typedef std::vector<int> intvector;

ex ncmul::expand(unsigned options) const
{
-       exvector sub_expanded_seq;
-
-       exvector expanded_seq=expandchildren(options);
-
-
-       int number_of_expanded_terms=1;
-
-       unsigned current_position=0;
-       exvector::const_iterator last=expanded_seq.end();
+       // First, expand the children
+       exvector expanded_seq = expandchildren(options);
+
+       // 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.
+
+       int number_of_expanded_terms = 1;
+
+       unsigned current_position = 0;
+       exvector::const_iterator last = expanded_seq.end();
for (exvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) {
}
current_position++;
}

-               return (new ncmul(expanded_seq,1))->setflag(status_flags::dynallocated ||
-                                                                                                       status_flags::expanded);
-       }
+       // If there are no sums, we are done
+               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;
-
-       int l;
-       for (l=0; l<number_of_adds; l++) {
-               k[l]=0;
-       }

-       while (1) {
-               exvector term;
-               term=expanded_seq;
-               for (l=0; l<number_of_adds; l++) {
+       while (true) {
+               exvector term = expanded_seq;
+               for (int i=0; i<number_of_adds; 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[]
-                       k[l]=0;
+               while ((l>=0) && ((++k[l]) >= number_of_add_operands[l])) {
+                       k[l] = 0;
l--;
}
-               if (l<0) break;
+               if (l<0)
+                       break;
}

-                                                                               status_flags::expanded);
+               setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
}

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 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;
}
@@ -234,7 +241,7 @@ 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();
@@ -245,17 +252,17 @@ ex ncmul::coeff(const ex & 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);
@@ -282,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<e.nops(); i++)
append_factors(v,e.op(i));
-
-               return;
-       }
-       v.push_back(e);
+       } else
+               v.push_back(e);
}

typedef std::vector<unsigned> unsignedvector;
@@ -317,30 +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;
@@ -353,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());

@@ -364,7 +372,8 @@ ex ncmul::eval(int level) const
commutativeseq.reserve(count_commutative+1);
exvector noncommutativeseq;
noncommutativeseq.reserve(assocseq.size()-count_commutative);
-               for (i=0; i<assocseq.size(); ++i) {
+               unsigned num = assocseq.size();
+               for (unsigned i=0; i<num; ++i) {
if (rettypes[i]==return_types::commutative)
commutativeseq.push_back(assocseq[i]);
else
@@ -382,48 +391,51 @@ 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(); ++i) {
-                               if (ti==rttinfos[i]) {
+                       for (i=0; i<rtt_num; ++i) {
+                               if (ti == rttinfos[i]) {
evv[i].push_back(*cit);
break;
}
}
-                       if (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; i<evv.size(); ++i) {
+               for (i=0; i<evv_num; ++i)
s += evv[i].size();
-               }
-               GINAC_ASSERT(s==assocseq.size());
+               GINAC_ASSERT(s == assoc_num);
#endif // def DO_GINAC_ASSERT

// if all elements are of same type, simplify the string
-               if (evv.size()==1)
+               if (evv_num == 1)
return evv[0][0].simplify_ncmul(evv[0]);

exvector splitseq;
-               splitseq.reserve(evv.size());
-               for (i=0; i<evv.size(); ++i) {
+               splitseq.reserve(evv_num);
+               for (i=0; i<evv_num; ++i)
splitseq.push_back((new ncmul(evv[i]))->setflag(status_flags::dynallocated));
-               }

return (new mul(splitseq))->setflag(status_flags::dynallocated);
}
@@ -432,9 +444,34 @@ ex ncmul::eval(int level) const
status_flags::evaluated);
}

-ex ncmul::subs(const lst & ls, const lst & lr) const
+ex ncmul::evalm(void) const
{
-       return ncmul(subschildren(ls, lr));
+       // 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<matrix>(*it));
+               it++;
+               while (it != itend) {
+                       if (!is_ex_of_type(*it, matrix))
+                               goto no_matrix;
+                       prod = prod.mul(ex_to<matrix>(*it));
+                       it++;
+               }
+               delete s;
+               return prod;
+       }
+
+no_matrix:
+       return (new ncmul(s))->setflag(status_flags::dynallocated);
}

ex ncmul::thisexprseq(const exvector & v) const
@@ -449,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();
+
+       // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c)
+       exvector ncmulseq = seq;
+       for (unsigned i=0; i<num; ++i) {
+               ex e = seq[i].diff(s);
+               e.swap(ncmulseq[i]);
+               e.swap(ncmulseq[i]);
+       }
}

int ncmul::compare_same_type(const basic & other) const
@@ -463,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();
@@ -495,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;
}
@@ -523,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;
}
@@ -535,14 +587,6 @@ const exvector & ncmul::get_factors(void) const
return seq;
}

-//////////
-// static member variables
-//////////
-
-// protected
-
-unsigned ncmul::precedence = 50;
-
//////////
// friend functions
//////////
@@ -554,13 +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);
}

} // namespace GiNaC