ginac/ncmul.cpp

   1 /** @file ncmul.cpp
   2  *
   3  *  Implementation of GiNaC's non-commutative products of expressions. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include <algorithm>
  24 #include <iostream>
  25 #include <stdexcept>
  26
  27 #include "ncmul.h"
  28 #include "ex.h"
  29 #include "add.h"
  30 #include "mul.h"
  31 #include "archive.h"
  32 #include "debugmsg.h"
  33 #include "utils.h"
  34
  35 namespace GiNaC {
  36
  37 GINAC_IMPLEMENT_REGISTERED_CLASS(ncmul, exprseq)
  38
  39 //////////
  40 // default constructor, destructor, copy constructor assignment operator and helpers
  41 //////////
  42
  43 ncmul::ncmul()
  44 {
  45         debugmsg("ncmul default constructor",LOGLEVEL_CONSTRUCT);
  46         tinfo_key = TINFO_ncmul;
  47 }
  48
  49 DEFAULT_COPY(ncmul)
  50 DEFAULT_DESTROY(ncmul)
  51
  52 //////////
  53 // other constructors
  54 //////////
  55
  56 // public
  57
  58 ncmul::ncmul(const ex & lh, const ex & rh) : inherited(lh,rh)
  59 {
  60         debugmsg("ncmul constructor from ex,ex",LOGLEVEL_CONSTRUCT);
  61         tinfo_key = TINFO_ncmul;
  62 }
  63
  64 ncmul::ncmul(const ex & f1, const ex & f2, const ex & f3) : inherited(f1,f2,f3)
  65 {
  66         debugmsg("ncmul constructor from 3 ex",LOGLEVEL_CONSTRUCT);
  67         tinfo_key = TINFO_ncmul;
  68 }
  69
  70 ncmul::ncmul(const ex & f1, const ex & f2, const ex & f3,
  71              const ex & f4) : inherited(f1,f2,f3,f4)
  72 {
  73         debugmsg("ncmul constructor from 4 ex",LOGLEVEL_CONSTRUCT);
  74         tinfo_key = TINFO_ncmul;
  75 }
  76
  77 ncmul::ncmul(const ex & f1, const ex & f2, const ex & f3,
  78              const ex & f4, const ex & f5) : inherited(f1,f2,f3,f4,f5)
  79 {
  80         debugmsg("ncmul constructor from 5 ex",LOGLEVEL_CONSTRUCT);
  81         tinfo_key = TINFO_ncmul;
  82 }
  83
  84 ncmul::ncmul(const ex & f1, const ex & f2, const ex & f3,
  85              const ex & f4, const ex & f5, const ex & f6) : inherited(f1,f2,f3,f4,f5,f6)
  86 {
  87         debugmsg("ncmul constructor from 6 ex",LOGLEVEL_CONSTRUCT);
  88         tinfo_key = TINFO_ncmul;
  89 }
  90
  91 ncmul::ncmul(const exvector & v, bool discardable) : inherited(v,discardable)
  92 {
  93         debugmsg("ncmul constructor from exvector,bool",LOGLEVEL_CONSTRUCT);
  94         tinfo_key = TINFO_ncmul;
  95 }
  96
  97 ncmul::ncmul(exvector * vp) : inherited(vp)
  98 {
  99         debugmsg("ncmul constructor from exvector *",LOGLEVEL_CONSTRUCT);
 100         tinfo_key = TINFO_ncmul;
 101 }
 102
 103 //////////
 104 // archiving
 105 //////////
 106
 107 DEFAULT_ARCHIVING(ncmul)
 108
 109 //////////
 110 // functions overriding virtual functions from bases classes
 111 //////////
 112
 113 // public
 114
 115 void ncmul::print(std::ostream & os, unsigned upper_precedence) const
 116 {
 117         debugmsg("ncmul print",LOGLEVEL_PRINT);
 118         printseq(os,'(','*',')',precedence,upper_precedence);
 119 }
 120
 121 void ncmul::printraw(std::ostream & os) const
 122 {
 123         debugmsg("ncmul printraw",LOGLEVEL_PRINT);
 124         os << "ncmul(";
 125         for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
 126                 (*it).bp->printraw(os);
 127                 os << ",";
 128         }
 129         os << ",hash=" << hashvalue << ",flags=" << flags;
 130         os << ")";
 131 }
 132
 133 void ncmul::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
 134 {
 135         debugmsg("ncmul print csrc",LOGLEVEL_PRINT);
 136         exvector::const_iterator it;
 137         exvector::const_iterator itend = seq.end()-1;
 138         os << "ncmul(";
 139         for (it=seq.begin(); it!=itend; ++it) {
 140                 (*it).bp->printcsrc(os,precedence);
 141                 os << ",";
 142         }
 143         (*it).bp->printcsrc(os,precedence);
 144         os << ")";
 145 }
 146
 147 bool ncmul::info(unsigned inf) const
 148 {
 149         throw(std::logic_error("which flags have to be implemented in ncmul::info()?"));
 150 }
 151
 152 typedef std::vector<int> intvector;
 153
 154 ex ncmul::expand(unsigned options) const
 155 {
 156         exvector sub_expanded_seq;
 157         intvector positions_of_adds;
 158         intvector number_of_add_operands;
 159
 160         exvector expanded_seq=expandchildren(options);
 161
 162         positions_of_adds.resize(expanded_seq.size());
 163         number_of_add_operands.resize(expanded_seq.size());
 164
 165         int number_of_adds=0;
 166         int number_of_expanded_terms=1;
 167
 168         unsigned current_position=0;
 169         exvector::const_iterator last=expanded_seq.end();
 170         for (exvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) {
 171                 if (is_ex_exactly_of_type((*cit),add)) {
 172                         positions_of_adds[number_of_adds]=current_position;
 173                         const add & expanded_addref=ex_to_add(*cit);
 174                         number_of_add_operands[number_of_adds]=expanded_addref.seq.size();
 175                         number_of_expanded_terms *= expanded_addref.seq.size();
 176                         number_of_adds++;
 177                 }
 178                 current_position++;
 179         }
 180
 181         if (number_of_adds==0) {
 182                 return (new ncmul(expanded_seq,1))->setflag(status_flags::dynallocated ||
 183                                                                                                         status_flags::expanded);
 184         }
 185
 186         exvector distrseq;
 187         distrseq.reserve(number_of_expanded_terms);
 188
 189         intvector k;
 190         k.resize(number_of_adds);
 191
 192         int l;
 193         for (l=0; l<number_of_adds; l++) {
 194                 k[l]=0;
 195         }
 196
 197         while (1) {
 198                 exvector term;
 199                 term=expanded_seq;
 200                 for (l=0; l<number_of_adds; l++) {
 201                         GINAC_ASSERT(is_ex_exactly_of_type(expanded_seq[positions_of_adds[l]],add));
 202                         const add & addref=ex_to_add(expanded_seq[positions_of_adds[l]]);
 203                         term[positions_of_adds[l]]=addref.recombine_pair_to_ex(addref.seq[k[l]]);
 204                 }
 205                 distrseq.push_back((new ncmul(term,1))->setflag(status_flags::dynallocated |
 206                                                                                                                 status_flags::expanded));
 207
 208                 // increment k[]
 209                 l=number_of_adds-1;
 210                 while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) {
 211                         k[l]=0;
 212                         l--;
 213                 }
 214                 if (l<0) break;
 215         }
 216
 217         return (new add(distrseq))->setflag(status_flags::dynallocated |
 218                                                                                 status_flags::expanded);
 219 }
 220
 221 int ncmul::degree(const symbol & s) const
 222 {
 223         int deg_sum=0;
 224         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 225                 deg_sum+=(*cit).degree(s);
 226         }
 227         return deg_sum;
 228 }
 229
 230 int ncmul::ldegree(const symbol & s) const
 231 {
 232         int deg_sum=0;
 233         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 234                 deg_sum+=(*cit).ldegree(s);
 235         }
 236         return deg_sum;
 237 }
 238
 239 ex ncmul::coeff(const symbol & s, int n) const
 240 {
 241         exvector coeffseq;
 242         coeffseq.reserve(seq.size());
 243
 244         if (n==0) {
 245                 // product of individual coeffs
 246                 // if a non-zero power of s is found, the resulting product will be 0
 247                 exvector::const_iterator it=seq.begin();
 248                 while (it!=seq.end()) {
 249                         coeffseq.push_back((*it).coeff(s,n));
 250                         ++it;
 251                 }
 252                 return (new ncmul(coeffseq,1))->setflag(status_flags::dynallocated);
 253         }
 254
 255         exvector::const_iterator it=seq.begin();
 256         bool coeff_found=0;
 257         while (it!=seq.end()) {
 258                 ex c=(*it).coeff(s,n);
 259                 if (!c.is_zero()) {
 260                         coeffseq.push_back(c);
 261                         coeff_found=1;
 262                 } else {
 263                         coeffseq.push_back(*it);
 264                 }
 265                 ++it;
 266         }
 267
 268         if (coeff_found) return (new ncmul(coeffseq,1))->setflag(status_flags::dynallocated);
 269
 270         return _ex0();
 271 }
 272
 273 unsigned ncmul::count_factors(const ex & e) const
 274 {
 275         if ((is_ex_exactly_of_type(e,mul)&&(e.return_type()!=return_types::commutative))||
 276                 (is_ex_exactly_of_type(e,ncmul))) {
 277                 unsigned factors=0;
 278                 for (unsigned i=0; i<e.nops(); i++)
 279                         factors += count_factors(e.op(i));
 280
 281                 return factors;
 282         }
 283         return 1;
 284 }
 285
 286 void ncmul::append_factors(exvector & v, const ex & e) const
 287 {
 288         if ((is_ex_exactly_of_type(e,mul)&&(e.return_type()!=return_types::commutative))||
 289                 (is_ex_exactly_of_type(e,ncmul))) {
 290                 for (unsigned i=0; i<e.nops(); i++)
 291                         append_factors(v,e.op(i));
 292
 293                 return;
 294         }
 295         v.push_back(e);
 296 }
 297
 298 typedef std::vector<unsigned> unsignedvector;
 299 typedef std::vector<exvector> exvectorvector;
 300
 301 ex ncmul::eval(int level) const
 302 {
 303         // simplifications: ncmul(...,*(x1,x2),...,ncmul(x3,x4),...) ->
 304         //                      ncmul(...,x1,x2,...,x3,x4,...) (associativity)
 305         //                  ncmul(x) -> x
 306         //                  ncmul() -> 1
 307         //                  ncmul(...,c1,...,c2,...)
 308         //                      *(c1,c2,ncmul(...)) (pull out commutative elements)
 309         //                  ncmul(x1,y1,x2,y2) -> *(ncmul(x1,x2),ncmul(y1,y2))
 310         //                      (collect elements of same type)
 311         //                  ncmul(x1,x2,x3,...) -> x::simplify_ncmul(x1,x2,x3,...)
 312         // the following rule would be nice, but produces a recursion,
 313         // which must be trapped by introducing a flag that the sub-ncmuls()
 314         // are already evaluated (maybe later...)
 315         //                  ncmul(x1,x2,...,X,y1,y2,...) ->
 316         //                      ncmul(ncmul(x1,x2,...),X,ncmul(y1,y2,...)
 317         //                      (X noncommutative_composite)
 318
 319         if ((level==1) && (flags & status_flags::evaluated)) {
 320                 return *this;
 321         }
 322
 323         exvector evaledseq=evalchildren(level);
 324
 325         // ncmul(...,*(x1,x2),...,ncmul(x3,x4),...) ->
 326         //     ncmul(...,x1,x2,...,x3,x4,...) (associativity)
 327         unsigned factors=0;
 328         for (exvector::const_iterator cit=evaledseq.begin(); cit!=evaledseq.end(); ++cit)
 329                 factors += count_factors(*cit);
 330
 331         exvector assocseq;
 332         assocseq.reserve(factors);
 333         for (exvector::const_iterator cit=evaledseq.begin(); cit!=evaledseq.end(); ++cit)
 334                 append_factors(assocseq,*cit);
 335
 336         // ncmul(x) -> x
 337         if (assocseq.size()==1) return *(seq.begin());
 338
 339         // ncmul() -> 1
 340         if (assocseq.size()==0) return _ex1();
 341
 342         // determine return types
 343         unsignedvector rettypes;
 344         rettypes.reserve(assocseq.size());
 345         unsigned i=0;
 346         unsigned count_commutative=0;
 347         unsigned count_noncommutative=0;
 348         unsigned count_noncommutative_composite=0;
 349         for (exvector::const_iterator cit=assocseq.begin(); cit!=assocseq.end(); ++cit) {
 350                 switch (rettypes[i]=(*cit).return_type()) {
 351                 case return_types::commutative:
 352                         count_commutative++;
 353                         break;
 354                 case return_types::noncommutative:
 355                         count_noncommutative++;
 356                         break;
 357                 case return_types::noncommutative_composite:
 358                         count_noncommutative_composite++;
 359                         break;
 360                 default:
 361                         throw(std::logic_error("ncmul::eval(): invalid return type"));
 362                 }
 363                 ++i;
 364         }
 365         GINAC_ASSERT(count_commutative+count_noncommutative+count_noncommutative_composite==assocseq.size());
 366
 367         // ncmul(...,c1,...,c2,...) ->
 368         //     *(c1,c2,ncmul(...)) (pull out commutative elements)
 369         if (count_commutative!=0) {
 370                 exvector commutativeseq;
 371                 commutativeseq.reserve(count_commutative+1);
 372                 exvector noncommutativeseq;
 373                 noncommutativeseq.reserve(assocseq.size()-count_commutative);
 374                 for (i=0; i<assocseq.size(); ++i) {
 375                         if (rettypes[i]==return_types::commutative)
 376                                 commutativeseq.push_back(assocseq[i]);
 377                         else
 378                                 noncommutativeseq.push_back(assocseq[i]);
 379                 }
 380                 commutativeseq.push_back((new ncmul(noncommutativeseq,1))->setflag(status_flags::dynallocated));
 381                 return (new mul(commutativeseq))->setflag(status_flags::dynallocated);
 382         }
 383
 384         // ncmul(x1,y1,x2,y2) -> *(ncmul(x1,x2),ncmul(y1,y2))
 385         //     (collect elements of same type)
 386
 387         if (count_noncommutative_composite==0) {
 388                 // there are neither commutative nor noncommutative_composite
 389                 // elements in assocseq
 390                 GINAC_ASSERT(count_commutative==0);
 391
 392                 exvectorvector evv;
 393                 unsignedvector rttinfos;
 394                 evv.reserve(assocseq.size());
 395                 rttinfos.reserve(assocseq.size());
 396
 397                 for (exvector::const_iterator cit=assocseq.begin(); cit!=assocseq.end(); ++cit) {
 398                         unsigned ti=(*cit).return_type_tinfo();
 399                         // search type in vector of known types
 400                         for (i=0; i<rttinfos.size(); ++i) {
 401                                 if (ti==rttinfos[i]) {
 402                                         evv[i].push_back(*cit);
 403                                         break;
 404                                 }
 405                         }
 406                         if (i>=rttinfos.size()) {
 407                                 // new type
 408                                 rttinfos.push_back(ti);
 409                                 evv.push_back(exvector());
 410                                 (*(evv.end()-1)).reserve(assocseq.size());
 411                                 (*(evv.end()-1)).push_back(*cit);
 412                         }
 413                 }
 414
 415 #ifdef DO_GINAC_ASSERT
 416                 GINAC_ASSERT(evv.size()==rttinfos.size());
 417                 GINAC_ASSERT(evv.size()>0);
 418                 unsigned s=0;
 419                 for (i=0; i<evv.size(); ++i) {
 420                         s += evv[i].size();
 421                 }
 422                 GINAC_ASSERT(s==assocseq.size());
 423 #endif // def DO_GINAC_ASSERT
 424
 425                 // if all elements are of same type, simplify the string
 426                 if (evv.size()==1)
 427                         return evv[0][0].simplify_ncmul(evv[0]);
 428
 429                 exvector splitseq;
 430                 splitseq.reserve(evv.size());
 431                 for (i=0; i<evv.size(); ++i) {
 432                         splitseq.push_back((new ncmul(evv[i]))->setflag(status_flags::dynallocated));
 433                 }
 434
 435                 return (new mul(splitseq))->setflag(status_flags::dynallocated);
 436         }
 437
 438         return (new ncmul(assocseq))->setflag(status_flags::dynallocated |
 439                                                                                   status_flags::evaluated);
 440 }
 441
 442 ex ncmul::subs(const lst & ls, const lst & lr) const
 443 {
 444         return ncmul(subschildren(ls, lr));
 445 }
 446
 447 ex ncmul::thisexprseq(const exvector & v) const
 448 {
 449         return (new ncmul(v))->setflag(status_flags::dynallocated);
 450 }
 451
 452 ex ncmul::thisexprseq(exvector * vp) const
 453 {
 454         return (new ncmul(vp))->setflag(status_flags::dynallocated);
 455 }
 456
 457 // protected
 458
 459 /** Implementation of ex::diff() for a non-commutative product. It always returns 0.
 460  *  @see ex::diff */
 461 ex ncmul::derivative(const symbol & s) const
 462 {
 463         return _ex0();
 464 }
 465
 466 int ncmul::compare_same_type(const basic & other) const
 467 {
 468         return inherited::compare_same_type(other);
 469 }
 470
 471 unsigned ncmul::return_type(void) const
 472 {
 473         if (seq.size()==0) {
 474                 // ncmul without factors: should not happen, but commutes
 475                 return return_types::commutative;
 476         }
 477
 478         bool all_commutative=1;
 479         unsigned rt;
 480         exvector::const_iterator cit_noncommutative_element; // point to first found nc element
 481
 482         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 483                 rt=(*cit).return_type();
 484                 if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc
 485                 if ((rt==return_types::noncommutative)&&(all_commutative)) {
 486                         // first nc element found, remember position
 487                         cit_noncommutative_element=cit;
 488                         all_commutative=0;
 489                 }
 490                 if ((rt==return_types::noncommutative)&&(!all_commutative)) {
 491                         // another nc element found, compare type_infos
 492                         if ((*cit_noncommutative_element).return_type_tinfo()!=(*cit).return_type_tinfo()) {
 493                                 // diffent types -> mul is ncc
 494                                 return return_types::noncommutative_composite;
 495                         }
 496                 }
 497         }
 498         // all factors checked
 499         GINAC_ASSERT(!all_commutative); // not all factors should commute, because this is a ncmul();
 500         return all_commutative ? return_types::commutative : return_types::noncommutative;
 501 }
 502
 503 unsigned ncmul::return_type_tinfo(void) const
 504 {
 505         if (seq.size()==0) {
 506                 // mul without factors: should not happen
 507                 return tinfo_key;
 508         }
 509         // return type_info of first noncommutative element
 510         for (exvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 511                 if ((*cit).return_type()==return_types::noncommutative) {
 512                         return (*cit).return_type_tinfo();
 513                 }
 514         }
 515         // no noncommutative element found, should not happen
 516         return tinfo_key;
 517 }
 518
 519 //////////
 520 // new virtual functions which can be overridden by derived classes
 521 //////////
 522
 523 // none
 524
 525 //////////
 526 // non-virtual functions in this class
 527 //////////
 528
 529 exvector ncmul::expandchildren(unsigned options) const
 530 {
 531         exvector s;
 532         s.reserve(seq.size());
 533
 534         for (exvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
 535                 s.push_back((*it).expand(options));
 536         }
 537         return s;
 538 }
 539
 540 const exvector & ncmul::get_factors(void) const
 541 {
 542         return seq;
 543 }
 544
 545 //////////
 546 // static member variables
 547 //////////
 548
 549 // protected
 550
 551 unsigned ncmul::precedence = 50;
 552
 553 //////////
 554 // friend functions
 555 //////////
 556
 557 ex nonsimplified_ncmul(const exvector & v)
 558 {
 559         return (new ncmul(v))->setflag(status_flags::dynallocated);
 560 }
 561
 562 ex simplified_ncmul(const exvector & v)
 563 {
 564         if (v.size()==0) {
 565                 return _ex1();
 566         } else if (v.size()==1) {
 567                 return v[0];
 568         }
 569         return (new ncmul(v))->setflag(status_flags::dynallocated |
 570                                        status_flags::evaluated);
 571 }
 572
 573 } // namespace GiNaC