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