ginac/mul.cpp

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