- The default implementations of evalf(), diff(), normal() and expand() use

[ginac.git] / ginac / ncmul.cpp

diff --git a/ginac/ncmul.cpp b/ginac/ncmul.cpp
index abe72ef1f726e64e560ebaca45061d3f4ddc04e9..dbd5732840ccde39e9704e09893f4dc44a236b1f 100644 (file)
--- a/ginac/ncmul.cpp
+++ b/ginac/ncmul.cpp
@@ -28,6 +28,7 @@
 #include "ex.h"
 #include "add.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
@@ -163,7 +164,7 @@ ex ncmul::expand(unsigned options) const
        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);
+                       const add & expanded_addref=ex_to<add>(*cit);
                        number_of_add_operands[number_of_adds]=expanded_addref.seq.size();
                        number_of_expanded_terms *= expanded_addref.seq.size();
                        number_of_adds++;
@@ -192,7 +193,7 @@ ex ncmul::expand(unsigned options) const
                term=expanded_seq;
                for (l=0; l<number_of_adds; l++) {
                        GINAC_ASSERT(is_ex_exactly_of_type(expanded_seq[positions_of_adds[l]],add));
-                       const add & addref=ex_to_add(expanded_seq[positions_of_adds[l]]);
+                       const add & addref=ex_to<add>(expanded_seq[positions_of_adds[l]]);
                        term[positions_of_adds[l]]=addref.recombine_pair_to_ex(addref.seq[k[l]]);
                }
                distrseq.push_back((new ncmul(term,1))->setflag(status_flags::dynallocated |
@@ -432,6 +433,36 @@ ex ncmul::eval(int level) const
                                                                                  status_flags::evaluated);
 }
 
+ex ncmul::evalm(void) const
+{
+       // 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
 {
        return (new ncmul(v))->setflag(status_flags::dynallocated);
@@ -444,11 +475,21 @@ 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();
+       exvector addseq;
+       addseq.reserve(seq.size());
+       
+       // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c)
+       for (unsigned i=0; i!=seq.size(); ++i) {
+               exvector ncmulseq = seq;
+               ncmulseq[i] = seq[i].diff(s);
+               addseq.push_back((new ncmul(ncmulseq))->setflag(status_flags::dynallocated));
+       }
+       return (new add(addseq))->setflag(status_flags::dynallocated);
 }
 
 int ncmul::compare_same_type(const basic & other) const
@@ -518,9 +559,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;
 }
@@ -530,14 +572,6 @@ const exvector & ncmul::get_factors(void) const
        return seq;
 }
 
-//////////
-// static member variables
-//////////
-
-// protected
-
-unsigned ncmul::precedence = 50;
-
 //////////
 // friend functions
 //////////