ginac/mul.cpp

   1 /** @file mul.cpp
   2  *
   3  *  Implementation of GiNaC's products of expressions. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999 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 <vector>
  24 #include <stdexcept>
  25
  26 #include "mul.h"
  27 #include "add.h"
  28 #include "power.h"
  29 #include "debugmsg.h"
  30
  31 #ifndef NO_GINAC_NAMESPACE
  32 namespace GiNaC {
  33 #endif // ndef NO_GINAC_NAMESPACE
  34
  35 //////////
  36 // default constructor, destructor, copy constructor assignment operator and helpers
  37 //////////
  38
  39 // public
  40
  41 mul::mul()
  42 {
  43     debugmsg("mul default constructor",LOGLEVEL_CONSTRUCT);
  44     tinfo_key = TINFO_mul;
  45 }
  46
  47 mul::~mul()
  48 {
  49     debugmsg("mul destructor",LOGLEVEL_DESTRUCT);
  50     destroy(0);
  51 }
  52
  53 mul::mul(mul const & other)
  54 {
  55     debugmsg("mul copy constructor",LOGLEVEL_CONSTRUCT);
  56     copy(other);
  57 }
  58
  59 mul const & mul::operator=(mul const & other)
  60 {
  61     debugmsg("mul operator=",LOGLEVEL_ASSIGNMENT);
  62     if (this != &other) {
  63         destroy(1);
  64         copy(other);
  65     }
  66     return *this;
  67 }
  68
  69 // protected
  70
  71 void mul::copy(mul const & other)
  72 {
  73     expairseq::copy(other);
  74 }
  75
  76 void mul::destroy(bool call_parent)
  77 {
  78     if (call_parent) expairseq::destroy(call_parent);
  79 }
  80
  81 //////////
  82 // other constructors
  83 //////////
  84
  85 // public
  86
  87 mul::mul(ex const & lh, ex const & rh)
  88 {
  89     debugmsg("mul constructor from ex,ex",LOGLEVEL_CONSTRUCT);
  90     tinfo_key = TINFO_mul;
  91     overall_coeff=exONE();
  92     construct_from_2_ex(lh,rh);
  93     GINAC_ASSERT(is_canonical());
  94 }
  95
  96 mul::mul(exvector const & v)
  97 {
  98     debugmsg("mul constructor from exvector",LOGLEVEL_CONSTRUCT);
  99     tinfo_key = TINFO_mul;
 100     overall_coeff=exONE();
 101     construct_from_exvector(v);
 102     GINAC_ASSERT(is_canonical());
 103 }
 104
 105 /*
 106 mul::mul(epvector const & v, bool do_not_canonicalize)
 107 {
 108     debugmsg("mul constructor from epvector,bool",LOGLEVEL_CONSTRUCT);
 109     tinfo_key = TINFO_mul;
 110     if (do_not_canonicalize) {
 111         seq=v;
 112 #ifdef EXPAIRSEQ_USE_HASHTAB
 113         combine_same_terms(); // to build hashtab
 114 #endif // def EXPAIRSEQ_USE_HASHTAB
 115     } else {
 116         construct_from_epvector(v);
 117     }
 118     GINAC_ASSERT(is_canonical());
 119 }
 120 */
 121
 122 mul::mul(epvector const & v)
 123 {
 124     debugmsg("mul constructor from epvector",LOGLEVEL_CONSTRUCT);
 125     tinfo_key = TINFO_mul;
 126     overall_coeff=exONE();
 127     construct_from_epvector(v);
 128     GINAC_ASSERT(is_canonical());
 129 }
 130
 131 mul::mul(epvector const & v, ex const & oc)
 132 {
 133     debugmsg("mul constructor from epvector,ex",LOGLEVEL_CONSTRUCT);
 134     tinfo_key = TINFO_mul;
 135     overall_coeff=oc;
 136     construct_from_epvector(v);
 137     GINAC_ASSERT(is_canonical());
 138 }
 139
 140 mul::mul(epvector * vp, ex const & oc)
 141 {
 142     debugmsg("mul constructor from epvector *,ex",LOGLEVEL_CONSTRUCT);
 143     tinfo_key = TINFO_mul;
 144     GINAC_ASSERT(vp!=0);
 145     overall_coeff=oc;
 146     construct_from_epvector(*vp);
 147     delete vp;
 148     GINAC_ASSERT(is_canonical());
 149 }
 150
 151 mul::mul(ex const & lh, ex const & mh, ex const & rh)
 152 {
 153     debugmsg("mul constructor from ex,ex,ex",LOGLEVEL_CONSTRUCT);
 154     tinfo_key = TINFO_mul;
 155     exvector factors;
 156     factors.reserve(3);
 157     factors.push_back(lh);
 158     factors.push_back(mh);
 159     factors.push_back(rh);
 160     overall_coeff=exONE();
 161     construct_from_exvector(factors);
 162     GINAC_ASSERT(is_canonical());
 163 }
 164
 165 //////////
 166 // functions overriding virtual functions from bases classes
 167 //////////
 168
 169 // public
 170
 171 basic * mul::duplicate() const
 172 {
 173     debugmsg("mul duplicate",LOGLEVEL_ASSIGNMENT);
 174     return new mul(*this);
 175 }
 176
 177 bool mul::info(unsigned inf) const
 178 {
 179     // TODO: optimize
 180     if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial || inf==info_flags::rational_function) {
 181         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
 182             if (!(recombine_pair_to_ex(*it).info(inf)))
 183                 return false;
 184         }
 185         return true;
 186     } else {
 187         return expairseq::info(inf);
 188     }
 189 }
 190
 191 typedef vector<int> intvector;
 192
 193 int mul::degree(symbol const & s) const
 194 {
 195     int deg_sum=0;
 196     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 197         deg_sum+=(*cit).rest.degree(s) * ex_to_numeric((*cit).coeff).to_int();
 198     }
 199     return deg_sum;
 200 }
 201
 202 int mul::ldegree(symbol const & s) const
 203 {
 204     int deg_sum=0;
 205     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 206         deg_sum+=(*cit).rest.ldegree(s) * ex_to_numeric((*cit).coeff).to_int();
 207     }
 208     return deg_sum;
 209 }
 210
 211 ex mul::coeff(symbol const & s, int const n) const
 212 {
 213     exvector coeffseq;
 214     coeffseq.reserve(seq.size()+1);
 215
 216     if (n==0) {
 217         // product of individual coeffs
 218         // if a non-zero power of s is found, the resulting product will be 0
 219         epvector::const_iterator it=seq.begin();
 220         while (it!=seq.end()) {
 221             coeffseq.push_back(recombine_pair_to_ex(*it).coeff(s,n));
 222             ++it;
 223         }
 224         coeffseq.push_back(overall_coeff);
 225         return (new mul(coeffseq))->setflag(status_flags::dynallocated);
 226     }
 227
 228     epvector::const_iterator it=seq.begin();
 229     bool coeff_found=0;
 230     while (it!=seq.end()) {
 231         ex t=recombine_pair_to_ex(*it);
 232         ex c=t.coeff(s,n);
 233         if (!c.is_zero()) {
 234             coeffseq.push_back(c);
 235             coeff_found=1;
 236         } else {
 237             coeffseq.push_back(t);
 238         }
 239         ++it;
 240     }
 241     if (coeff_found) {
 242         coeffseq.push_back(overall_coeff);
 243         return (new mul(coeffseq))->setflag(status_flags::dynallocated);
 244     }
 245
 246     return exZERO();
 247 }
 248
 249 /*
 250 ex mul::eval(int level) const
 251 {
 252     // simplifications: *(...,x,(c1,1),(c2,1)) -> *(...,x,(c1*c2,1)) (c1, c2 numeric(), move pairs to end first)
 253     //                  *(...,x,1) -> *(...,x)
 254     //                  *(...,x,0) -> 0
 255     //                  *(+(x,y,...),(c,1)) -> *(+(*(x,c),*(y,c),...)) (c numeric())
 256     //                  *(x) -> x
 257     //                  *() -> 1
 258
 259     debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
 260
 261     if ((level==1)&&(flags & status_flags::evaluated)) {
 262 #ifdef DO_GINAC_ASSERT
 263         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 264             GINAC_ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
 265                    (!(ex_to_numeric((*cit).coeff).is_integer())));
 266         }
 267
 268         // test if all numerics were moved to the end and
 269         // all numerics with coeff 1 to the very end
 270         if (seq.size()!=0) {
 271             epvector::const_iterator cit=seq.end();
 272             bool all_coeff_1=true;
 273             bool all_numeric=true;
 274             do {
 275                 cit--;
 276                 if (is_ex_exactly_of_type((*cit).rest,numeric)) {
 277                     GINAC_ASSERT(all_numeric);
 278                     if ((*cit).coeff.is_equal(exONE())) {
 279                         GINAC_ASSERT(all_coeff_1);
 280                     } else {
 281                         all_coeff_1=false;
 282                     }
 283                 } else {
 284                     all_numeric=false;
 285                 }
 286             } while (cit!=seq.begin());
 287         }
 288 #endif // def DO_GINAC_ASSERT
 289         return *this;
 290     }
 291
 292     epvector newseq;
 293     epvector::iterator it1,it2;
 294     bool seq_copied=false;
 295
 296     epvector * evaled_seqp=evalchildren(level);
 297     if (evaled_seqp!=0) {
 298         // do more evaluation later
 299         return (new mul(evaled_seqp))->setflag(status_flags::dynallocated);
 300     }
 301
 302     // combine pairs with coeff 1 (all numerics should be at end, assert below)
 303     if (seq.size()>1) {
 304         // count number of pairs with coeff 1
 305         unsigned num_coeff_1=0;
 306         bool still_numeric=true;
 307         epvector::const_iterator cit=seq.end();
 308         unsigned first_pos;
 309         unsigned second_pos;
 310         do {
 311             cit--;
 312             if (is_ex_exactly_of_type((*cit).rest,numeric)) {
 313                 if ((*cit).coeff.is_equal(exONE())) {
 314                     num_coeff_1++;
 315                 }
 316             } else {
 317                 still_numeric=false;
 318             }
 319         } while ((cit!=seq.begin())&&still_numeric);
 320         if (num_coeff_1>1) {
 321             newseq=seq;
 322
 323     }
 324
 325
 326 #ifdef DO_GINAC_ASSERT
 327     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 328         GINAC_ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
 329                (!(ex_to_numeric((*cit).coeff).is_integer())));
 330     }
 331
 332     // test if all numerics were moved to the end and
 333     // all numerics with coeff 1 to the very end
 334     if (seq.size()!=0) {
 335         epvector::const_iterator cit=seq.end();
 336         bool all_coeff_1=true;
 337         bool all_numeric=true;
 338         do {
 339             cit--;
 340             if (is_ex_exactly_of_type((*cit).rest,numeric)) {
 341                 GINAC_ASSERT(all_numeric);
 342                 if ((*cit).coeff.is_equal(exONE())) {
 343                     GINAC_ASSERT(all_coeff_1);
 344                 } else {
 345                     all_coeff_1=false;
 346                 }
 347             } else {
 348                 all_numeric=false;
 349             }
 350         } while (cit!=seq.begin());
 351     }
 352 #endif // def DO_GINAC_ASSERT
 353
 354     if (flags & status_flags::evaluated) {
 355         return *this;
 356     }
 357
 358     expair const & last_expair=*(seq.end()-1);
 359     expair const & next_to_last_expair=*(seq.end()-2);
 360     int seq_size = seq.size();
 361
 362     // *(...,x,(c1,1),(c2,1)) -> *(...,x,(c1*c2,1)) (c1, c2 numeric())
 363     if ((!seq_copied) && (seq_size>=2) &&
 364         is_ex_exactly_of_type(last_expair.rest,numeric) &&
 365         ex_to_numeric(last_expair.coeff).is_equal(numONE()) &&
 366         is_ex_exactly_of_type(next_to_last_expair.rest,numeric) &&
 367         ex_to_numeric(next_to_last_expair.coeff).is_equal(numONE()) ) {
 368         newseq=seq;
 369         seq_copied=true;
 370         it2=newseq.end()-1;
 371         it1=it2-1;
 372     }
 373     while (seq_copied && (newseq.size()>=2) &&
 374            is_ex_exactly_of_type((*it1).rest,numeric) &&
 375            ex_to_numeric((*it1).coeff).is_equal(numONE()) &&
 376            is_ex_exactly_of_type((*it2).rest,numeric) &&
 377            ex_to_numeric((*it2).coeff).is_equal(numONE()) ) {
 378         *it1=expair(ex_to_numeric((*it1).rest).mul_dyn(ex_to_numeric((*it2).rest)),exONE());
 379         newseq.pop_back();
 380         it2=newseq.end()-1;
 381         it1=it2-1;
 382     }
 383
 384     // *(...,x,1) -> *(...,x)
 385     if ((!seq_copied) && (seq_size>=1) &&
 386         (is_ex_exactly_of_type(last_expair.rest,numeric)) &&
 387         (ex_to_numeric(last_expair.rest).compare(numONE())==0)) {
 388         newseq=seq;
 389         seq_copied=true;
 390         it2=newseq.end()-1;
 391     }
 392     if (seq_copied && (newseq.size()>=1) &&
 393         (is_ex_exactly_of_type((*it2).rest,numeric)) &&
 394         (ex_to_numeric((*it2).rest).compare(numONE())==0)) {
 395         newseq.pop_back();
 396         it2=newseq.end()-1;
 397     }
 398
 399     // *(...,x,0) -> 0
 400     if ((!seq_copied) && (seq_size>=1) &&
 401         (is_ex_exactly_of_type(last_expair.rest,numeric)) &&
 402         (ex_to_numeric(last_expair.rest).is_zero())) {
 403         return exZERO();
 404     }
 405     if (seq_copied && (newseq.size()>=1) &&
 406         (is_ex_exactly_of_type((*it2).rest,numeric)) &&
 407         (ex_to_numeric((*it2).rest).is_zero())) {
 408         return exZERO();
 409     }
 410
 411     // *(+(x,y,...),c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
 412     if ((!seq_copied) && (seq_size==2) &&
 413         is_ex_exactly_of_type(next_to_last_expair.rest,add) &&
 414         is_ex_exactly_of_type(last_expair.rest,numeric) &&
 415         ex_to_numeric(last_expair.coeff).is_equal(numONE()) &&
 416         (ex_to_numeric(next_to_last_expair.coeff).compare(numONE())==0)) {
 417         add const & addref=ex_to_add(next_to_last_expair.rest);
 418         epvector distrseq;
 419         distrseq.reserve(addref.seq.size());
 420         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
 421             distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit,
 422                                    last_expair.rest));
 423         }
 424         // special treatment for the last element if it is numeric (to
 425         // avoid terms like (2/3)*(3/2)) is no longer necessary, this
 426         // is handled in add::combine_pair_with_coeff_to_pair()
 427         return (new add(distrseq,1))->setflag(status_flags::dynallocated  |
 428                                               status_flags::evaluated );
 429     }
 430     if (seq_copied && (newseq.size()==2) &&
 431         is_ex_exactly_of_type(newseq[0].rest,add) &&
 432         is_ex_exactly_of_type(newseq[1].rest,numeric) &&
 433         ex_to_numeric(newseq[1].coeff).is_equal(numONE()) &&
 434         (ex_to_numeric(newseq[0].coeff).compare(numONE())==0)) {
 435         add const & addref=ex_to_add(newseq[0].rest);
 436         epvector distrseq;
 437         distrseq.reserve(addref.seq.size());
 438         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
 439             distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit,
 440                                    newseq[1].rest));
 441         }
 442         // special treatment for the last element if it is numeric (to
 443         // avoid terms like (2/3)*(3/2)) is no longer necessary, this
 444         // is handled in add::combine_pair_with_coeff_to_pair()
 445         return (new add(distrseq,1))->setflag(status_flags::dynallocated  |
 446                                               status_flags::evaluated );
 447     }
 448
 449     // *() -> 1
 450     if ((!seq_copied) && (seq_size==0)) {
 451         return exONE();
 452     } else if (seq_copied && (newseq.size()==0)) {
 453         return exONE();
 454     }
 455
 456     // *(x) -> x
 457     if ((!seq_copied) && (seq_size==1)) {
 458         return recombine_pair_to_ex(*(seq.begin()));
 459     } else if (seq_copied && (newseq.size()==1)) {
 460         return recombine_pair_to_ex(*(newseq.begin()));
 461     }
 462
 463     if (!seq_copied) return this->hold();
 464
 465     return (new mul(newseq,1))->setflag(status_flags::dynallocated  |
 466                                         status_flags::evaluated );
 467 }
 468 */
 469
 470 ex mul::eval(int level) const
 471 {
 472     // simplifications  *(...,x;0) -> 0
 473     //                  *(+(x,y,...);c) -> *(+(*(x,c),*(y,c),...)) (c numeric())
 474     //                  *(x;1) -> x
 475     //                  *(;c) -> c
 476
 477     debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
 478
 479     epvector * evaled_seqp=evalchildren(level);
 480     if (evaled_seqp!=0) {
 481         // do more evaluation later
 482         return (new mul(evaled_seqp,overall_coeff))->
 483                    setflag(status_flags::dynallocated);
 484     }
 485
 486 #ifdef DO_GINAC_ASSERT
 487     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 488         GINAC_ASSERT((!is_ex_exactly_of_type((*cit).rest,mul))||
 489                (!(ex_to_numeric((*cit).coeff).is_integer())));
 490         GINAC_ASSERT(!((*cit).is_numeric_with_coeff_1()));
 491         if (is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric)) {
 492             printtree(cerr,0);
 493         }
 494         GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric));
 495         /* for paranoia */
 496         expair p=split_ex_to_pair(recombine_pair_to_ex(*cit));
 497         GINAC_ASSERT(p.rest.is_equal((*cit).rest));
 498         GINAC_ASSERT(p.coeff.is_equal((*cit).coeff));
 499         /* end paranoia */
 500     }
 501 #endif // def DO_GINAC_ASSERT
 502
 503     if (flags & status_flags::evaluated) {
 504         GINAC_ASSERT(seq.size()>0);
 505         GINAC_ASSERT((seq.size()>1)||!overall_coeff.is_equal(exONE()));
 506         return *this;
 507     }
 508
 509     int seq_size=seq.size();
 510     if (overall_coeff.is_equal(exZERO())) {
 511         // *(...,x;0) -> 0
 512         return exZERO();
 513     } else if (seq_size==0) {
 514         // *(;c) -> c
 515         return overall_coeff;
 516     } else if ((seq_size==1)&&overall_coeff.is_equal(exONE())) {
 517         // *(x;1) -> x
 518         return recombine_pair_to_ex(*(seq.begin()));
 519     } else if ((seq_size==1) &&
 520                is_ex_exactly_of_type((*seq.begin()).rest,add) &&
 521                ex_to_numeric((*seq.begin()).coeff).is_equal(numONE())) {
 522         // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
 523         add const & addref=ex_to_add((*seq.begin()).rest);
 524         epvector distrseq;
 525         distrseq.reserve(addref.seq.size());
 526         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
 527             distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit,
 528                                    overall_coeff));
 529         }
 530         return (new add(distrseq,
 531                         ex_to_numeric(addref.overall_coeff).
 532                         mul_dyn(ex_to_numeric(overall_coeff))))
 533             ->setflag(status_flags::dynallocated  |
 534                       status_flags::evaluated );
 535     }
 536     return this->hold();
 537 }
 538
 539 /*
 540 ex mul::eval(int level) const
 541 {
 542     // simplifications: *(...,x,c1,c2) -> *(...,x,c1*c2) (c1, c2 numeric())
 543     //                  *(...,(c1,c2)) -> (...,(c1^c2,1)) (normalize)
 544     //                  *(...,x,1) -> +(...,x)
 545     //                  *(...,x,0) -> 0
 546     //                  *(+(x,y,...),c) -> *(+(*(x,c),*(y,c),...)) (c numeric())
 547     //                  *(x) -> x
 548     //                  *() -> 1
 549
 550     debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
 551
 552     epvector newseq=seq;
 553     epvector::iterator it1,it2;
 554
 555     // *(...,x,c1,c2) -> *(...,x,c1*c2) (c1, c2 numeric())
 556     it2=newseq.end()-1;
 557     it1=it2-1;
 558     while ((newseq.size()>=2)&&is_exactly_of_type(*(*it1).rest.bp,numeric)&&
 559                                is_exactly_of_type(*(*it2).rest.bp,numeric)) {
 560         *it1=expair(ex_to_numeric((*it1).rest).power(ex_to_numeric((*it1).coeff))
 561                     .mul(ex_to_numeric((*it2).rest).power(ex_to_numeric((*it2).coeff))),exONE());
 562         newseq.pop_back();
 563         it2=newseq.end()-1;
 564         it1=it2-1;
 565     }
 566
 567     if ((newseq.size()>=1)&&is_exactly_of_type(*(*it2).rest.bp,numeric)) {
 568         // *(...,(c1,c2)) -> (...,(c1^c2,1)) (normalize)
 569         *it2=expair(ex_to_numeric((*it2).rest).power(ex_to_numeric((*it2).coeff)),exONE());
 570         // *(...,x,1) -> *(...,x)
 571         if (static_cast<numeric &>(*(*it2).rest.bp).compare(numONE())==0) {
 572             newseq.pop_back();
 573             it2=newseq.end()-1;
 574         }
 575     }
 576
 577     // *(...,x,0) -> 0
 578     if ((newseq.size()>=1)&&is_exactly_of_type(*(*it2).rest.bp,numeric)) {
 579         if (static_cast<numeric &>(*(*it2).rest.bp).is_zero()==0) {
 580             return exZERO();
 581         }
 582     }
 583
 584     // *(+(x,y,...),c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
 585     if ((newseq.size()==2)&&is_ex_exactly_of_type(newseq[0].rest,add)&&
 586         is_ex_exactly_of_type(newseq[1].rest,numeric)&&
 587         (ex_to_numeric(newseq[0].coeff).compare(numONE())==0)) {
 588         add const & addref=ex_to_add(newseq[0].rest);
 589         numeric const & numref=ex_to_numeric(newseq[1].rest);
 590         epvector distrseq;
 591         distrseq.reserve(addref.seq.size());
 592         for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
 593             distrseq.push_back(expair((*cit).rest,ex_to_numeric((*cit).coeff).mul(numref)));
 594         }
 595         return (new add(distrseq,1))->setflag(status_flags::dynallocated  |
 596                                               status_flags::evaluated );
 597     }
 598
 599     if (newseq.size()==0) {
 600         // *() -> 1
 601         return exONE();
 602     } else if (newseq.size()==1) {
 603         // *(x) -> x
 604         return recombine_pair_to_ex(*newseq.begin());
 605     }
 606
 607     return (new mul(newseq,1))->setflag(status_flags::dynallocated  |
 608                                         status_flags::evaluated );
 609 }
 610 */
 611
 612 exvector mul::get_indices(void) const
 613 {
 614     // return union of indices of factors
 615     exvector iv;
 616     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 617         exvector subiv=(*cit).rest.get_indices();
 618         iv.reserve(iv.size()+subiv.size());
 619         for (exvector::const_iterator cit2=subiv.begin(); cit2!=subiv.end(); ++cit2) {
 620             iv.push_back(*cit2);
 621         }
 622     }
 623     return iv;
 624 }
 625
 626 ex mul::simplify_ncmul(exvector const & v) const
 627 {
 628     throw(std::logic_error("mul::simplify_ncmul() should never have been called!"));
 629 }
 630
 631 // protected
 632
 633 int mul::compare_same_type(basic const & other) const
 634 {
 635     return expairseq::compare_same_type(other);
 636 }
 637
 638 bool mul::is_equal_same_type(basic const & other) const
 639 {
 640     return expairseq::is_equal_same_type(other);
 641 }
 642
 643 unsigned mul::return_type(void) const
 644 {
 645     if (seq.size()==0) {
 646         // mul without factors: should not happen, but commutes
 647         return return_types::commutative;
 648     }
 649
 650     bool all_commutative=1;
 651     unsigned rt;
 652     epvector::const_iterator cit_noncommutative_element; // point to first found nc element
 653
 654     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 655         rt=(*cit).rest.return_type();
 656         if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc
 657         if ((rt==return_types::noncommutative)&&(all_commutative)) {
 658             // first nc element found, remember position
 659             cit_noncommutative_element=cit;
 660             all_commutative=0;
 661         }
 662         if ((rt==return_types::noncommutative)&&(!all_commutative)) {
 663                 // another nc element found, compare type_infos
 664             if ((*cit_noncommutative_element).rest.return_type_tinfo()!=(*cit).rest.return_type_tinfo()) {
 665                 // diffent types -> mul is ncc
 666                 return return_types::noncommutative_composite;
 667             }
 668         }
 669     }
 670     // all factors checked
 671     return all_commutative ? return_types::commutative : return_types::noncommutative;
 672 }
 673
 674 unsigned mul::return_type_tinfo(void) const
 675 {
 676     if (seq.size()==0) {
 677         // mul without factors: should not happen
 678         return tinfo_key;
 679     }
 680     // return type_info of first noncommutative element
 681     for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 682         if ((*cit).rest.return_type()==return_types::noncommutative) {
 683             return (*cit).rest.return_type_tinfo();
 684         }
 685     }
 686     // no noncommutative element found, should not happen
 687     return tinfo_key;
 688 }
 689
 690 ex mul::thisexpairseq(epvector const & v, ex const & oc) const
 691 {
 692     return (new mul(v,oc))->setflag(status_flags::dynallocated);
 693 }
 694
 695 ex mul::thisexpairseq(epvector * vp, ex const & oc) const
 696 {
 697     return (new mul(vp,oc))->setflag(status_flags::dynallocated);
 698 }
 699
 700 expair mul::split_ex_to_pair(ex const & e) const
 701 {
 702     if (is_ex_exactly_of_type(e,power)) {
 703         power const & powerref=ex_to_power(e);
 704         if (is_ex_exactly_of_type(powerref.exponent,numeric)) {
 705             return expair(powerref.basis,powerref.exponent);
 706         }
 707     }
 708     return expair(e,exONE());
 709 }
 710
 711 expair mul::combine_ex_with_coeff_to_pair(ex const & e,
 712                                           ex const & c) const
 713 {
 714     // to avoid duplication of power simplification rules,
 715     // we create a temporary power object
 716     // otherwise it would be hard to correctly simplify
 717     // expression like (4^(1/3))^(3/2)
 718     if (are_ex_trivially_equal(c,exONE())) {
 719         return split_ex_to_pair(e);
 720     }
 721     return split_ex_to_pair(power(e,c));
 722 }
 723
 724 expair mul::combine_pair_with_coeff_to_pair(expair const & p,
 725                                             ex const & c) const
 726 {
 727     // to avoid duplication of power simplification rules,
 728     // we create a temporary power object
 729     // otherwise it would be hard to correctly simplify
 730     // expression like (4^(1/3))^(3/2)
 731     if (are_ex_trivially_equal(c,exONE())) {
 732         return p;
 733     }
 734     return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
 735 }
 736
 737 ex mul::recombine_pair_to_ex(expair const & p) const
 738 {
 739     // if (p.coeff.compare(exONE())==0) {
 740     // if (are_ex_trivially_equal(p.coeff,exONE())) {
 741     if (ex_to_numeric(p.coeff).is_equal(numONE())) {
 742         return p.rest;
 743     } else {
 744         return power(p.rest,p.coeff);
 745     }
 746 }
 747
 748 bool mul::expair_needs_further_processing(epp it)
 749 {
 750     if (is_ex_exactly_of_type((*it).rest,mul) &&
 751         ex_to_numeric((*it).coeff).is_integer()) {
 752         // combined pair is product with integer power -> expand it
 753         *it=split_ex_to_pair(recombine_pair_to_ex(*it));
 754         return true;
 755     }
 756     if (is_ex_exactly_of_type((*it).rest,numeric)) {
 757         expair ep=split_ex_to_pair(recombine_pair_to_ex(*it));
 758         if (!ep.is_equal(*it)) {
 759             // combined pair is a numeric power which can be simplified
 760             *it=ep;
 761             return true;
 762         }
 763         if (ex_to_numeric((*it).coeff).is_equal(numONE())) {
 764             // combined pair has coeff 1 and must be moved to the end
 765             return true;
 766         }
 767     }
 768     return false;
 769 }
 770
 771 ex mul::default_overall_coeff(void) const
 772 {
 773     return exONE();
 774 }
 775
 776 void mul::combine_overall_coeff(ex const & c)
 777 {
 778     GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
 779     GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
 780     overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c));
 781 }
 782
 783 void mul::combine_overall_coeff(ex const & c1, ex const & c2)
 784 {
 785     GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
 786     GINAC_ASSERT(is_ex_exactly_of_type(c1,numeric));
 787     GINAC_ASSERT(is_ex_exactly_of_type(c2,numeric));
 788     overall_coeff = ex_to_numeric(overall_coeff).
 789                         mul_dyn(ex_to_numeric(c1).power(ex_to_numeric(c2)));
 790 }
 791
 792 bool mul::can_make_flat(expair const & p) const
 793 {
 794     GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
 795     // this assertion will probably fail somewhere
 796     // it would require a more careful make_flat, obeying the power laws
 797     // probably should return true only if p.coeff is integer
 798     return ex_to_numeric(p.coeff).is_equal(numONE());
 799 }
 800
 801 ex mul::expand(unsigned options) const
 802 {
 803     exvector sub_expanded_seq;
 804     intvector positions_of_adds;
 805     intvector number_of_add_operands;
 806
 807     epvector * expanded_seqp=expandchildren(options);
 808
 809     epvector const & expanded_seq = expanded_seqp==0 ? seq : *expanded_seqp;
 810
 811     positions_of_adds.resize(expanded_seq.size());
 812     number_of_add_operands.resize(expanded_seq.size());
 813
 814     int number_of_adds=0;
 815     int number_of_expanded_terms=1;
 816
 817     unsigned current_position=0;
 818     epvector::const_iterator last=expanded_seq.end();
 819     for (epvector::const_iterator cit=expanded_seq.begin(); cit!=last; ++cit) {
 820         if (is_ex_exactly_of_type((*cit).rest,add)&&
 821             (ex_to_numeric((*cit).coeff).is_equal(numONE()))) {
 822             positions_of_adds[number_of_adds]=current_position;
 823             add const & expanded_addref=ex_to_add((*cit).rest);
 824             int addref_nops=expanded_addref.nops();
 825             number_of_add_operands[number_of_adds]=addref_nops;
 826             number_of_expanded_terms *= addref_nops;
 827             number_of_adds++;
 828         }
 829         current_position++;
 830     }
 831
 832     if (number_of_adds==0) {
 833         if (expanded_seqp==0) {
 834             return this->setflag(status_flags::expanded);
 835         }
 836         return (new mul(expanded_seqp,overall_coeff))->
 837                      setflag(status_flags::dynallocated ||
 838                              status_flags::expanded);
 839     }
 840
 841     exvector distrseq;
 842     distrseq.reserve(number_of_expanded_terms);
 843
 844     intvector k;
 845     k.resize(number_of_adds);
 846
 847     int l;
 848     for (l=0; l<number_of_adds; l++) {
 849         k[l]=0;
 850     }
 851
 852     while (1) {
 853         epvector term;
 854         term=expanded_seq;
 855         for (l=0; l<number_of_adds; l++) {
 856             add const & addref=ex_to_add(expanded_seq[positions_of_adds[l]].rest);
 857             GINAC_ASSERT(term[positions_of_adds[l]].coeff.compare(exONE())==0);
 858             term[positions_of_adds[l]]=split_ex_to_pair(addref.op(k[l]));
 859         }
 860         /*
 861         cout << "mul::expand() term begin" << endl;
 862         for (epvector::const_iterator cit=term.begin(); cit!=term.end(); ++cit) {
 863             cout << "rest" << endl;
 864             (*cit).rest.printtree(cout);
 865             cout << "coeff" << endl;
 866             (*cit).coeff.printtree(cout);
 867         }
 868         cout << "mul::expand() term end" << endl;
 869         */
 870         distrseq.push_back((new mul(term,overall_coeff))->
 871                                 setflag(status_flags::dynallocated |
 872                                         status_flags::expanded));
 873
 874         // increment k[]
 875         l=number_of_adds-1;
 876         while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) {
 877             k[l]=0;
 878             l--;
 879         }
 880         if (l<0) break;
 881     }
 882
 883     if (expanded_seqp!=0) {
 884         delete expanded_seqp;
 885     }
 886     /*
 887     cout << "mul::expand() distrseq begin" << endl;
 888     for (exvector::const_iterator cit=distrseq.begin(); cit!=distrseq.end(); ++cit) {
 889         (*cit).printtree(cout);
 890     }
 891     cout << "mul::expand() distrseq end" << endl;
 892     */
 893
 894     return (new add(distrseq))->setflag(status_flags::dynallocated |
 895                                         status_flags::expanded);
 896 }
 897
 898 /*
 899 ex mul::expand(unsigned options) const
 900 {
 901     exvector sub_expanded_seq;
 902     intvector positions_of_adds;
 903     intvector number_of_add_operands;
 904
 905     sub_expanded_seq.resize(seq.size());
 906     positions_of_adds.resize(seq.size());
 907     number_of_add_operands.reserve(seq.size());
 908
 909     int number_of_adds=0;
 910     int number_of_expanded_terms=1;
 911     for (unsigned current_position=0; current_position<seq.size(); current_position++) {
 912         ex const & expanded_ex=recombine_pair_to_ex(seq[current_position]).expand(options);
 913         if (is_ex_exactly_of_type(expanded_ex,add)) {
 914             positions_of_adds[number_of_adds]=current_position;
 915             add const & expanded_addref=ex_to_add(expanded_ex);
 916             number_of_add_operands[number_of_adds]=expanded_addref.seq.size();
 917             number_of_expanded_terms *= expanded_addref.seq.size();
 918             number_of_adds++;
 919         }
 920         sub_expanded_seq.push_back(expanded_ex);
 921     }
 922
 923     exvector distrseq;
 924     distrseq.reserve(number_of_expanded_terms);
 925
 926     intvector k;
 927     k.resize(number_of_adds);
 928
 929     int l;
 930     for (l=0; l<number_of_adds; l++) {
 931         k[l]=0;
 932     }
 933
 934     while (1) {
 935         exvector term;
 936         term=sub_expanded_seq;
 937         for (l=0; l<number_of_adds; l++) {
 938             add const & addref=ex_to_add(sub_expanded_seq[positions_of_adds[l]]);
 939             term[positions_of_adds[l]]=addref.recombine_pair_to_ex(addref.seq[k[l]]);
 940         }
 941         distrseq.push_back((new mul(term))->setflag(status_flags::dynallocated |
 942                                                     status_flags::expanded));
 943
 944         // increment k[]
 945         l=number_of_adds-1;
 946         while ((l>=0)&&((++k[l])>=number_of_add_operands[l])) {
 947             k[l]=0;
 948             l--;
 949         }
 950         if (l<0) break;
 951     }
 952
 953     return (new add(distrseq))->setflag(status_flags::dynallocated |
 954                                         status_flags::expanded);
 955 }
 956 */
 957
 958 //////////
 959 // new virtual functions which can be overridden by derived classes
 960 //////////
 961
 962 // none
 963
 964 //////////
 965 // non-virtual functions in this class
 966 //////////
 967
 968 epvector * mul::expandchildren(unsigned options) const
 969 {
 970     epvector::const_iterator last=seq.end();
 971     epvector::const_iterator cit=seq.begin();
 972     while (cit!=last) {
 973         ex const & factor=recombine_pair_to_ex(*cit);
 974         ex const & expanded_factor=factor.expand(options);
 975         if (!are_ex_trivially_equal(factor,expanded_factor)) {
 976
 977             // something changed, copy seq, eval and return it
 978             epvector *s=new epvector;
 979             s->reserve(seq.size());
 980
 981             // copy parts of seq which are known not to have changed
 982             epvector::const_iterator cit2=seq.begin();
 983             while (cit2!=cit) {
 984                 s->push_back(*cit2);
 985                 ++cit2;
 986             }
 987             // copy first changed element
 988             s->push_back(split_ex_to_pair(expanded_factor));
 989             ++cit2;
 990             // copy rest
 991             while (cit2!=last) {
 992                 s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options)));
 993                 ++cit2;
 994             }
 995             return s;
 996         }
 997         ++cit;
 998     }
 999
1000     return 0; // nothing has changed
1001 }
1002
1003 //////////
1004 // static member variables
1005 //////////
1006
1007 // protected
1008
1009 unsigned mul::precedence=50;
1010
1011
1012 //////////
1013 // global constants
1014 //////////
1015
1016 const mul some_mul;
1017 type_info const & typeid_mul=typeid(some_mul);
1018
1019 #ifndef NO_GINAC_NAMESPACE
1020 } // namespace GiNaC
1021 #endif // ndef NO_GINAC_NAMESPACE