ginac/mul.cpp

   1 /** @file mul.cpp
   2  *
   3  *  Implementation of GiNaC's products of expressions. */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2001 Johannes Gutenberg University Mainz, Germany
   7  *
   8  *  This program is free software; you can redistribute it and/or modify
   9  *  it under the terms of the GNU General Public License as published by
  10  *  the Free Software Foundation; either version 2 of the License, or
  11  *  (at your option) any later version.
  12  *
  13  *  This program is distributed in the hope that it will be useful,
  14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  16  *  GNU General Public License for more details.
  17  *
  18  *  You should have received a copy of the GNU General Public License
  19  *  along with this program; if not, write to the Free Software
  20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
  21  */
  22
  23 #include <vector>
  24 #include <stdexcept>
  25
  26 #include "mul.h"
  27 #include "add.h"
  28 #include "power.h"
  29 #include "archive.h"
  30 #include "debugmsg.h"
  31 #include "utils.h"
  32
  33 namespace GiNaC {
  34
  35 GINAC_IMPLEMENT_REGISTERED_CLASS(mul, expairseq)
  36
  37 //////////
  38 // default ctor, dctor, copy ctor assignment operator and helpers
  39 //////////
  40
  41 // public
  42
  43 mul::mul()
  44 {
  45         debugmsg("mul default ctor",LOGLEVEL_CONSTRUCT);
  46         tinfo_key = TINFO_mul;
  47 }
  48
  49 // protected
  50
  51 /** For use by copy ctor and assignment operator. */
  52 void mul::copy(const mul & other)
  53 {
  54         inherited::copy(other);
  55 }
  56
  57 void mul::destroy(bool call_parent)
  58 {
  59         if (call_parent) inherited::destroy(call_parent);
  60 }
  61
  62 //////////
  63 // other ctors
  64 //////////
  65
  66 // public
  67
  68 mul::mul(const ex & lh, const ex & rh)
  69 {
  70         debugmsg("mul ctor from ex,ex",LOGLEVEL_CONSTRUCT);
  71         tinfo_key = TINFO_mul;
  72         overall_coeff = _ex1();
  73         construct_from_2_ex(lh,rh);
  74         GINAC_ASSERT(is_canonical());
  75 }
  76
  77 mul::mul(const exvector & v)
  78 {
  79         debugmsg("mul ctor from exvector",LOGLEVEL_CONSTRUCT);
  80         tinfo_key = TINFO_mul;
  81         overall_coeff = _ex1();
  82         construct_from_exvector(v);
  83         GINAC_ASSERT(is_canonical());
  84 }
  85
  86 mul::mul(const epvector & v)
  87 {
  88         debugmsg("mul ctor from epvector",LOGLEVEL_CONSTRUCT);
  89         tinfo_key = TINFO_mul;
  90         overall_coeff = _ex1();
  91         construct_from_epvector(v);
  92         GINAC_ASSERT(is_canonical());
  93 }
  94
  95 mul::mul(const epvector & v, const ex & oc)
  96 {
  97         debugmsg("mul ctor from epvector,ex",LOGLEVEL_CONSTRUCT);
  98         tinfo_key = TINFO_mul;
  99         overall_coeff = oc;
 100         construct_from_epvector(v);
 101         GINAC_ASSERT(is_canonical());
 102 }
 103
 104 mul::mul(epvector * vp, const ex & oc)
 105 {
 106         debugmsg("mul ctor from epvector *,ex",LOGLEVEL_CONSTRUCT);
 107         tinfo_key = TINFO_mul;
 108         GINAC_ASSERT(vp!=0);
 109         overall_coeff = oc;
 110         construct_from_epvector(*vp);
 111         delete vp;
 112         GINAC_ASSERT(is_canonical());
 113 }
 114
 115 mul::mul(const ex & lh, const ex & mh, const ex & rh)
 116 {
 117         debugmsg("mul ctor from ex,ex,ex",LOGLEVEL_CONSTRUCT);
 118         tinfo_key = TINFO_mul;
 119         exvector factors;
 120         factors.reserve(3);
 121         factors.push_back(lh);
 122         factors.push_back(mh);
 123         factors.push_back(rh);
 124         overall_coeff = _ex1();
 125         construct_from_exvector(factors);
 126         GINAC_ASSERT(is_canonical());
 127 }
 128
 129 //////////
 130 // archiving
 131 //////////
 132
 133 /** Construct object from archive_node. */
 134 mul::mul(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 135 {
 136         debugmsg("mul ctor from archive_node", LOGLEVEL_CONSTRUCT);
 137 }
 138
 139 /** Unarchive the object. */
 140 ex mul::unarchive(const archive_node &n, const lst &sym_lst)
 141 {
 142         return (new mul(n, sym_lst))->setflag(status_flags::dynallocated);
 143 }
 144
 145 /** Archive the object. */
 146 void mul::archive(archive_node &n) const
 147 {
 148         inherited::archive(n);
 149 }
 150
 151 //////////
 152 // functions overriding virtual functions from bases classes
 153 //////////
 154
 155 // public
 156
 157 void mul::print(std::ostream & os, unsigned upper_precedence) const
 158 {
 159         debugmsg("mul print",LOGLEVEL_PRINT);
 160         if (precedence<=upper_precedence) os << "(";
 161         bool first = true;
 162         // first print the overall numeric coefficient:
 163         numeric coeff = ex_to_numeric(overall_coeff);
 164         if (coeff.csgn()==-1) os << '-';
 165         if (!coeff.is_equal(_num1()) &&
 166                 !coeff.is_equal(_num_1())) {
 167                 if (coeff.is_rational()) {
 168                         if (coeff.is_negative())
 169                                 os << -coeff;
 170                         else
 171                                 os << coeff;
 172                 } else {
 173                         if (coeff.csgn()==-1)
 174                                 (-coeff).print(os, precedence);
 175                         else
 176                                 coeff.print(os, precedence);
 177                 }
 178                 os << '*';
 179         }
 180         // then proceed with the remaining factors:
 181         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 182                 if (!first) {
 183                         os << '*';
 184                 } else {
 185                         first=false;
 186                 }
 187                 recombine_pair_to_ex(*cit).print(os,precedence);
 188         }
 189         if (precedence<=upper_precedence) os << ")";
 190 }
 191
 192 void mul::printraw(std::ostream & os) const
 193 {
 194         debugmsg("mul printraw",LOGLEVEL_PRINT);
 195
 196         os << "*(";
 197         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
 198                 os << "(";
 199                 (*it).rest.bp->printraw(os);
 200                 os << ",";
 201                 (*it).coeff.bp->printraw(os);
 202                 os << "),";
 203         }
 204         os << ",hash=" << hashvalue << ",flags=" << flags;
 205         os << ")";
 206 }
 207
 208 void mul::printcsrc(std::ostream & os, unsigned type, unsigned upper_precedence) const
 209 {
 210         debugmsg("mul print csrc", LOGLEVEL_PRINT);
 211         if (precedence <= upper_precedence)
 212                 os << "(";
 213
 214         if (!overall_coeff.is_equal(_ex1())) {
 215                 overall_coeff.bp->printcsrc(os,type,precedence);
 216                 os << "*";
 217         }
 218
 219         // Print arguments, separated by "*" or "/"
 220         epvector::const_iterator it = seq.begin();
 221         epvector::const_iterator itend = seq.end();
 222         while (it != itend) {
 223
 224                 // If the first argument is a negative integer power, it gets printed as "1.0/<expr>"
 225                 if (it == seq.begin() && ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0) {
 226                         if (type == csrc_types::ctype_cl_N)
 227                                 os << "recip(";
 228                         else
 229                                 os << "1.0/";
 230                 }
 231
 232                 // If the exponent is 1 or -1, it is left out
 233                 if (it->coeff.compare(_ex1()) == 0 || it->coeff.compare(_num_1()) == 0)
 234                         it->rest.bp->printcsrc(os, type, precedence);
 235                 else
 236                         // outer parens around ex needed for broken gcc-2.95 parser:
 237                         (ex(power(it->rest, abs(ex_to_numeric(it->coeff))))).bp->printcsrc(os, type, upper_precedence);
 238
 239                 // Separator is "/" for negative integer powers, "*" otherwise
 240                 ++it;
 241                 if (it != itend) {
 242                         if (ex_to_numeric(it->coeff).is_integer() && it->coeff.compare(_num0()) < 0)
 243                                 os << "/";
 244                         else
 245                                 os << "*";
 246                 }
 247         }
 248         if (precedence <= upper_precedence)
 249                 os << ")";
 250 }
 251
 252 bool mul::info(unsigned inf) const
 253 {
 254         switch (inf) {
 255                 case info_flags::polynomial:
 256                 case info_flags::integer_polynomial:
 257                 case info_flags::cinteger_polynomial:
 258                 case info_flags::rational_polynomial:
 259                 case info_flags::crational_polynomial:
 260                 case info_flags::rational_function: {
 261                         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
 262                                 if (!(recombine_pair_to_ex(*i).info(inf)))
 263                                         return false;
 264                         }
 265                         return overall_coeff.info(inf);
 266                 }
 267                 case info_flags::algebraic: {
 268                         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
 269                                 if ((recombine_pair_to_ex(*i).info(inf)))
 270                                         return true;
 271                         }
 272                         return false;
 273                 }
 274         }
 275         return inherited::info(inf);
 276 }
 277
 278 int mul::degree(const symbol & s) const
 279 {
 280         int deg_sum = 0;
 281         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 282                 if (ex_to_numeric(cit->coeff).is_integer())
 283                         deg_sum+=cit->rest.degree(s) * ex_to_numeric(cit->coeff).to_int();
 284         }
 285         return deg_sum;
 286 }
 287
 288 int mul::ldegree(const symbol & s) const
 289 {
 290         int deg_sum = 0;
 291         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 292                 if (ex_to_numeric(cit->coeff).is_integer())
 293                         deg_sum+=cit->rest.ldegree(s) * ex_to_numeric(cit->coeff).to_int();
 294         }
 295         return deg_sum;
 296 }
 297
 298 ex mul::coeff(const symbol & s, int n) const
 299 {
 300         exvector coeffseq;
 301         coeffseq.reserve(seq.size()+1);
 302
 303         if (n==0) {
 304                 // product of individual coeffs
 305                 // if a non-zero power of s is found, the resulting product will be 0
 306                 epvector::const_iterator it = seq.begin();
 307                 while (it!=seq.end()) {
 308                         coeffseq.push_back(recombine_pair_to_ex(*it).coeff(s,n));
 309                         ++it;
 310                 }
 311                 coeffseq.push_back(overall_coeff);
 312                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
 313         }
 314
 315         epvector::const_iterator it=seq.begin();
 316         bool coeff_found = 0;
 317         while (it!=seq.end()) {
 318                 ex t = recombine_pair_to_ex(*it);
 319                 ex c = t.coeff(s,n);
 320                 if (!c.is_zero()) {
 321                         coeffseq.push_back(c);
 322                         coeff_found = 1;
 323                 } else {
 324                         coeffseq.push_back(t);
 325                 }
 326                 ++it;
 327         }
 328         if (coeff_found) {
 329                 coeffseq.push_back(overall_coeff);
 330                 return (new mul(coeffseq))->setflag(status_flags::dynallocated);
 331         }
 332
 333         return _ex0();
 334 }
 335
 336 ex mul::eval(int level) const
 337 {
 338         // simplifications  *(...,x;0) -> 0
 339         //                  *(+(x,y,...);c) -> *(+(*(x,c),*(y,c),...)) (c numeric())
 340         //                  *(x;1) -> x
 341         //                  *(;c) -> c
 342
 343         debugmsg("mul eval",LOGLEVEL_MEMBER_FUNCTION);
 344
 345         epvector * evaled_seqp = evalchildren(level);
 346         if (evaled_seqp!=0) {
 347                 // do more evaluation later
 348                 return (new mul(evaled_seqp,overall_coeff))->
 349                            setflag(status_flags::dynallocated);
 350         }
 351
 352 #ifdef DO_GINAC_ASSERT
 353         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 354                 GINAC_ASSERT((!is_ex_exactly_of_type((*cit).rest,mul)) ||
 355                              (!(ex_to_numeric((*cit).coeff).is_integer())));
 356                 GINAC_ASSERT(!(cit->is_canonical_numeric()));
 357                 if (is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric))
 358                     printtree(std::cerr,0);
 359                 GINAC_ASSERT(!is_ex_exactly_of_type(recombine_pair_to_ex(*cit),numeric));
 360                 /* for paranoia */
 361                 expair p = split_ex_to_pair(recombine_pair_to_ex(*cit));
 362                 GINAC_ASSERT(p.rest.is_equal((*cit).rest));
 363                 GINAC_ASSERT(p.coeff.is_equal((*cit).coeff));
 364                 /* end paranoia */
 365         }
 366 #endif // def DO_GINAC_ASSERT
 367
 368         if (flags & status_flags::evaluated) {
 369                 GINAC_ASSERT(seq.size()>0);
 370                 GINAC_ASSERT(seq.size()>1 || !overall_coeff.is_equal(_ex1()));
 371                 return *this;
 372         }
 373
 374         int seq_size = seq.size();
 375         if (overall_coeff.is_equal(_ex0())) {
 376                 // *(...,x;0) -> 0
 377                 return _ex0();
 378         } else if (seq_size==0) {
 379                 // *(;c) -> c
 380                 return overall_coeff;
 381         } else if (seq_size==1 && overall_coeff.is_equal(_ex1())) {
 382                 // *(x;1) -> x
 383                 return recombine_pair_to_ex(*(seq.begin()));
 384         } else if ((seq_size==1) &&
 385                    is_ex_exactly_of_type((*seq.begin()).rest,add) &&
 386                    ex_to_numeric((*seq.begin()).coeff).is_equal(_num1())) {
 387                 // *(+(x,y,...);c) -> +(*(x,c),*(y,c),...) (c numeric(), no powers of +())
 388                 const add & addref = ex_to_add((*seq.begin()).rest);
 389                 epvector distrseq;
 390                 distrseq.reserve(addref.seq.size());
 391                 for (epvector::const_iterator cit=addref.seq.begin(); cit!=addref.seq.end(); ++cit) {
 392                         distrseq.push_back(addref.combine_pair_with_coeff_to_pair(*cit, overall_coeff));
 393                 }
 394                 return (new add(distrseq,
 395                                 ex_to_numeric(addref.overall_coeff).
 396                                 mul_dyn(ex_to_numeric(overall_coeff))))
 397                       ->setflag(status_flags::dynallocated | status_flags::evaluated);
 398         }
 399         return this->hold();
 400 }
 401
 402 ex mul::evalf(int level) const
 403 {
 404         if (level==1)
 405                 return mul(seq,overall_coeff);
 406
 407         if (level==-max_recursion_level)
 408                 throw(std::runtime_error("max recursion level reached"));
 409
 410         epvector s;
 411         s.reserve(seq.size());
 412
 413         --level;
 414         for (epvector::const_iterator it=seq.begin(); it!=seq.end(); ++it) {
 415                 s.push_back(combine_ex_with_coeff_to_pair((*it).rest.evalf(level),
 416                                                           (*it).coeff));
 417         }
 418         return mul(s,overall_coeff.evalf(level));
 419 }
 420
 421 ex mul::simplify_ncmul(const exvector & v) const
 422 {
 423         throw(std::logic_error("mul::simplify_ncmul() should never have been called!"));
 424 }
 425
 426 // protected
 427
 428 /** Implementation of ex::diff() for a product.  It applies the product rule.
 429  *  @see ex::diff */
 430 ex mul::derivative(const symbol & s) const
 431 {
 432         exvector addseq;
 433         addseq.reserve(seq.size());
 434
 435         // D(a*b*c) = D(a)*b*c + a*D(b)*c + a*b*D(c)
 436         for (unsigned i=0; i!=seq.size(); ++i) {
 437                 epvector mulseq = seq;
 438                 mulseq[i] = split_ex_to_pair(power(seq[i].rest,seq[i].coeff - _ex1()) *
 439                                              seq[i].rest.diff(s));
 440                 addseq.push_back((new mul(mulseq,overall_coeff*seq[i].coeff))->setflag(status_flags::dynallocated));
 441         }
 442         return (new add(addseq))->setflag(status_flags::dynallocated);
 443 }
 444
 445 int mul::compare_same_type(const basic & other) const
 446 {
 447         return inherited::compare_same_type(other);
 448 }
 449
 450 bool mul::is_equal_same_type(const basic & other) const
 451 {
 452         return inherited::is_equal_same_type(other);
 453 }
 454
 455 unsigned mul::return_type(void) const
 456 {
 457         if (seq.size()==0) {
 458                 // mul without factors: should not happen, but commutes
 459                 return return_types::commutative;
 460         }
 461
 462         bool all_commutative = 1;
 463         unsigned rt;
 464         epvector::const_iterator cit_noncommutative_element; // point to first found nc element
 465
 466         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 467                 rt=(*cit).rest.return_type();
 468                 if (rt==return_types::noncommutative_composite) return rt; // one ncc -> mul also ncc
 469                 if ((rt==return_types::noncommutative)&&(all_commutative)) {
 470                         // first nc element found, remember position
 471                         cit_noncommutative_element = cit;
 472                         all_commutative = 0;
 473                 }
 474                 if ((rt==return_types::noncommutative)&&(!all_commutative)) {
 475                         // another nc element found, compare type_infos
 476                         if ((*cit_noncommutative_element).rest.return_type_tinfo()!=(*cit).rest.return_type_tinfo()) {
 477                                 // diffent types -> mul is ncc
 478                                 return return_types::noncommutative_composite;
 479                         }
 480                 }
 481         }
 482         // all factors checked
 483         return all_commutative ? return_types::commutative : return_types::noncommutative;
 484 }
 485
 486 unsigned mul::return_type_tinfo(void) const
 487 {
 488         if (seq.size()==0)
 489                 return tinfo_key;  // mul without factors: should not happen
 490
 491         // return type_info of first noncommutative element
 492         for (epvector::const_iterator cit=seq.begin(); cit!=seq.end(); ++cit) {
 493                 if ((*cit).rest.return_type()==return_types::noncommutative)
 494                         return (*cit).rest.return_type_tinfo();
 495         }
 496         // no noncommutative element found, should not happen
 497         return tinfo_key;
 498 }
 499
 500 ex mul::thisexpairseq(const epvector & v, const ex & oc) const
 501 {
 502         return (new mul(v,oc))->setflag(status_flags::dynallocated);
 503 }
 504
 505 ex mul::thisexpairseq(epvector * vp, const ex & oc) const
 506 {
 507         return (new mul(vp,oc))->setflag(status_flags::dynallocated);
 508 }
 509
 510 expair mul::split_ex_to_pair(const ex & e) const
 511 {
 512         if (is_ex_exactly_of_type(e,power)) {
 513                 const power & powerref = ex_to_power(e);
 514                 if (is_ex_exactly_of_type(powerref.exponent,numeric))
 515                         return expair(powerref.basis,powerref.exponent);
 516         }
 517         return expair(e,_ex1());
 518 }
 519
 520 expair mul::combine_ex_with_coeff_to_pair(const ex & e,
 521                                           const ex & c) const
 522 {
 523         // to avoid duplication of power simplification rules,
 524         // we create a temporary power object
 525         // otherwise it would be hard to correctly simplify
 526         // expression like (4^(1/3))^(3/2)
 527         if (are_ex_trivially_equal(c,_ex1()))
 528                 return split_ex_to_pair(e);
 529
 530         return split_ex_to_pair(power(e,c));
 531 }
 532
 533 expair mul::combine_pair_with_coeff_to_pair(const expair & p,
 534                                             const ex & c) const
 535 {
 536         // to avoid duplication of power simplification rules,
 537         // we create a temporary power object
 538         // otherwise it would be hard to correctly simplify
 539         // expression like (4^(1/3))^(3/2)
 540         if (are_ex_trivially_equal(c,_ex1()))
 541                 return p;
 542
 543         return split_ex_to_pair(power(recombine_pair_to_ex(p),c));
 544 }
 545
 546 ex mul::recombine_pair_to_ex(const expair & p) const
 547 {
 548         if (ex_to_numeric(p.coeff).is_equal(_num1()))
 549                 return p.rest;
 550         else
 551                 return power(p.rest,p.coeff);
 552 }
 553
 554 bool mul::expair_needs_further_processing(epp it)
 555 {
 556         if (is_ex_exactly_of_type((*it).rest,mul) &&
 557                 ex_to_numeric((*it).coeff).is_integer()) {
 558                 // combined pair is product with integer power -> expand it
 559                 *it = split_ex_to_pair(recombine_pair_to_ex(*it));
 560                 return true;
 561         }
 562         if (is_ex_exactly_of_type((*it).rest,numeric)) {
 563                 expair ep=split_ex_to_pair(recombine_pair_to_ex(*it));
 564                 if (!ep.is_equal(*it)) {
 565                         // combined pair is a numeric power which can be simplified
 566                         *it = ep;
 567                         return true;
 568                 }
 569                 if (ex_to_numeric((*it).coeff).is_equal(_num1())) {
 570                         // combined pair has coeff 1 and must be moved to the end
 571                         return true;
 572                 }
 573         }
 574         return false;
 575 }
 576
 577 ex mul::default_overall_coeff(void) const
 578 {
 579         return _ex1();
 580 }
 581
 582 void mul::combine_overall_coeff(const ex & c)
 583 {
 584         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
 585         GINAC_ASSERT(is_ex_exactly_of_type(c,numeric));
 586         overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c));
 587 }
 588
 589 void mul::combine_overall_coeff(const ex & c1, const ex & c2)
 590 {
 591         GINAC_ASSERT(is_ex_exactly_of_type(overall_coeff,numeric));
 592         GINAC_ASSERT(is_ex_exactly_of_type(c1,numeric));
 593         GINAC_ASSERT(is_ex_exactly_of_type(c2,numeric));
 594         overall_coeff = ex_to_numeric(overall_coeff).mul_dyn(ex_to_numeric(c1).power(ex_to_numeric(c2)));
 595 }
 596
 597 bool mul::can_make_flat(const expair & p) const
 598 {
 599         GINAC_ASSERT(is_ex_exactly_of_type(p.coeff,numeric));
 600         // this assertion will probably fail somewhere
 601         // it would require a more careful make_flat, obeying the power laws
 602         // probably should return true only if p.coeff is integer
 603         return ex_to_numeric(p.coeff).is_equal(_num1());
 604 }
 605
 606 ex mul::expand(unsigned options) const
 607 {
 608         if (flags & status_flags::expanded)
 609                 return *this;
 610
 611         exvector sub_expanded_seq;
 612
 613         epvector * expanded_seqp = expandchildren(options);
 614
 615         const epvector & expanded_seq = expanded_seqp==0 ? seq : *expanded_seqp;
 616
 617         int number_of_adds = 0;
 618         epvector non_adds;
 619         non_adds.reserve(expanded_seq.size());
 620         epvector::const_iterator cit = expanded_seq.begin();
 621         epvector::const_iterator last = expanded_seq.end();
 622         ex last_expanded = _ex1();
 623         while (cit!=last) {
 624                 if (is_ex_exactly_of_type((*cit).rest,add) &&
 625                         ((*cit).coeff.is_equal(_ex1()))) {
 626                         ++number_of_adds;
 627                         if (is_ex_exactly_of_type(last_expanded,add)) {
 628                                 // expand adds
 629                                 const add & add1 = ex_to_add(last_expanded);
 630                                 const add & add2 = ex_to_add((*cit).rest);
 631                                 int n1 = add1.nops();
 632                                 int n2 = add2.nops();
 633                                 exvector distrseq;
 634                                 distrseq.reserve(n1*n2);
 635                                 for (int i1=0; i1<n1; ++i1) {
 636                                         for (int i2=0; i2<n2; ++i2) {
 637                                                 distrseq.push_back(add1.op(i1)*add2.op(i2));
 638                                         }
 639                                 }
 640                                 last_expanded = (new add(distrseq))->setflag(status_flags::dynallocated | status_flags::expanded);
 641                         } else {
 642                                 non_adds.push_back(split_ex_to_pair(last_expanded));
 643                                 last_expanded = (*cit).rest;
 644                         }
 645                 } else {
 646                         non_adds.push_back(*cit);
 647                 }
 648                 ++cit;
 649         }
 650         if (expanded_seqp)
 651                 delete expanded_seqp;
 652
 653         if (is_ex_exactly_of_type(last_expanded,add)) {
 654                 add const & finaladd = ex_to_add(last_expanded);
 655                 exvector distrseq;
 656                 int n = finaladd.nops();
 657                 distrseq.reserve(n);
 658                 for (int i=0; i<n; ++i) {
 659                         epvector factors = non_adds;
 660                         factors.push_back(split_ex_to_pair(finaladd.op(i)));
 661                         distrseq.push_back((new mul(factors,overall_coeff))->setflag(status_flags::dynallocated | status_flags::expanded));
 662                 }
 663                 return ((new add(distrseq))->
 664                         setflag(status_flags::dynallocated | status_flags::expanded));
 665         }
 666         non_adds.push_back(split_ex_to_pair(last_expanded));
 667         return (new mul(non_adds,overall_coeff))->
 668                 setflag(status_flags::dynallocated | status_flags::expanded);
 669 }
 670
 671
 672 //////////
 673 // new virtual functions which can be overridden by derived classes
 674 //////////
 675
 676 // none
 677
 678 //////////
 679 // non-virtual functions in this class
 680 //////////
 681
 682
 683 /** Member-wise expand the expairs representing this sequence.  This must be
 684  *  overridden from expairseq::expandchildren() and done iteratively in order
 685  *  to allow for early cancallations and thus safe memory.
 686  *
 687  *  @see mul::expand()
 688  *  @return pointer to epvector containing expanded representation or zero
 689  *  pointer, if sequence is unchanged. */
 690 epvector * mul::expandchildren(unsigned options) const
 691 {
 692         epvector::const_iterator last = seq.end();
 693         epvector::const_iterator cit = seq.begin();
 694         while (cit!=last) {
 695                 const ex & factor = recombine_pair_to_ex(*cit);
 696                 const ex & expanded_factor = factor.expand(options);
 697                 if (!are_ex_trivially_equal(factor,expanded_factor)) {
 698
 699                         // something changed, copy seq, eval and return it
 700                         epvector *s = new epvector;
 701                         s->reserve(seq.size());
 702
 703                         // copy parts of seq which are known not to have changed
 704                         epvector::const_iterator cit2 = seq.begin();
 705                         while (cit2!=cit) {
 706                                 s->push_back(*cit2);
 707                                 ++cit2;
 708                         }
 709                         // copy first changed element
 710                         s->push_back(split_ex_to_pair(expanded_factor));
 711                         ++cit2;
 712                         // copy rest
 713                         while (cit2!=last) {
 714                                 s->push_back(split_ex_to_pair(recombine_pair_to_ex(*cit2).expand(options)));
 715                                 ++cit2;
 716                         }
 717                         return s;
 718                 }
 719                 ++cit;
 720         }
 721
 722         return 0; // nothing has changed
 723 }
 724
 725 //////////
 726 // static member variables
 727 //////////
 728
 729 // protected
 730
 731 unsigned mul::precedence = 50;
 732
 733 } // namespace GiNaC