+ // simplifications: ncmul(...,*(x1,x2),...,ncmul(x3,x4),...) ->
+ // ncmul(...,x1,x2,...,x3,x4,...) (associativity)
+ // ncmul(x) -> x
+ // ncmul() -> 1
+ // ncmul(...,c1,...,c2,...)
+ // *(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::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...)
+ // ncmul(x1,x2,...,X,y1,y2,...) ->
+ // ncmul(ncmul(x1,x2,...),X,ncmul(y1,y2,...)
+ // (X noncommutative_composite)
+
+ if ((level==1) && (flags & status_flags::evaluated)) {
+ return *this;
+ }
+
+ exvector evaledseq=evalchildren(level);
+
+ // ncmul(...,*(x1,x2),...,ncmul(x3,x4),...) ->
+ // ncmul(...,x1,x2,...,x3,x4,...) (associativity)
+ unsigned factors = 0;
+ exvector::const_iterator cit = evaledseq.begin(), citend = evaledseq.end();
+ while (cit != citend)
+ factors += count_factors(*cit++);
+
+ exvector assocseq;
+ assocseq.reserve(factors);
+ cit = evaledseq.begin();
+ while (cit != citend)
+ append_factors(assocseq, *cit++);
+
+ // ncmul(x) -> x
+ if (assocseq.size()==1) return *(seq.begin());
+
+ // ncmul() -> 1
+ if (assocseq.empty()) return _ex1();
+
+ // determine return types
+ unsignedvector rettypes;
+ rettypes.reserve(assocseq.size());
+ unsigned i = 0;
+ unsigned count_commutative=0;
+ unsigned count_noncommutative=0;
+ unsigned count_noncommutative_composite=0;
+ cit = assocseq.begin(); citend = assocseq.end();
+ while (cit != citend) {
+ switch (rettypes[i] = cit->return_type()) {
+ case return_types::commutative:
+ count_commutative++;
+ break;
+ case return_types::noncommutative:
+ count_noncommutative++;
+ break;
+ case return_types::noncommutative_composite:
+ count_noncommutative_composite++;
+ break;
+ default:
+ throw(std::logic_error("ncmul::eval(): invalid return type"));
+ }
+ ++i; ++cit;
+ }
+ GINAC_ASSERT(count_commutative+count_noncommutative+count_noncommutative_composite==assocseq.size());
+
+ // ncmul(...,c1,...,c2,...) ->
+ // *(c1,c2,ncmul(...)) (pull out commutative elements)
+ if (count_commutative!=0) {
+ exvector commutativeseq;
+ commutativeseq.reserve(count_commutative+1);
+ exvector noncommutativeseq;
+ noncommutativeseq.reserve(assocseq.size()-count_commutative);
+ unsigned num = assocseq.size();
+ for (unsigned i=0; i<num; ++i) {
+ if (rettypes[i]==return_types::commutative)
+ commutativeseq.push_back(assocseq[i]);
+ else
+ noncommutativeseq.push_back(assocseq[i]);
+ }
+ commutativeseq.push_back((new ncmul(noncommutativeseq,1))->setflag(status_flags::dynallocated));
+ return (new mul(commutativeseq))->setflag(status_flags::dynallocated);
+ }
+
+ // ncmul(x1,y1,x2,y2) -> *(ncmul(x1,x2),ncmul(y1,y2))
+ // (collect elements of same type)
+
+ if (count_noncommutative_composite==0) {
+ // there are neither commutative nor noncommutative_composite
+ // elements in assocseq
+ GINAC_ASSERT(count_commutative==0);
+
+ unsigned assoc_num = assocseq.size();
+ exvectorvector evv;
+ unsignedvector rttinfos;
+ evv.reserve(assoc_num);
+ rttinfos.reserve(assoc_num);
+
+ 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<rtt_num; ++i) {
+ if (ti == rttinfos[i]) {
+ evv[i].push_back(*cit);
+ break;
+ }
+ }
+ if (i >= rtt_num) {
+ // new type
+ rttinfos.push_back(ti);
+ evv.push_back(exvector());
+ (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_num == rttinfos.size());
+ GINAC_ASSERT(evv_num > 0);
+ unsigned s=0;
+ for (i=0; i<evv_num; ++i)
+ s += evv[i].size();
+ GINAC_ASSERT(s == assoc_num);
+#endif // def DO_GINAC_ASSERT
+
+ // if all elements are of same type, simplify the string
+ if (evv_num == 1)
+ return evv[0][0].simplify_ncmul(evv[0]);
+
+ exvector splitseq;
+ 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);
+ }
+
+ return (new ncmul(assocseq))->setflag(status_flags::dynallocated |
+ 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);
git://www.ginac.de / ginac.git / blobdiff