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