ginac/power.cpp

   1 /** @file power.cpp
   2  *
   3  *  Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
   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 <iostream>
  25 #include <stdexcept>
  26
  27 #include "power.h"
  28 #include "expairseq.h"
  29 #include "add.h"
  30 #include "mul.h"
  31 #include "numeric.h"
  32 #include "relational.h"
  33 #include "symbol.h"
  34 #include "debugmsg.h"
  35
  36 #ifndef NO_GINAC_NAMESPACE
  37 namespace GiNaC {
  38 #endif // ndef NO_GINAC_NAMESPACE
  39
  40 typedef vector<int> intvector;
  41
  42 //////////
  43 // default constructor, destructor, copy constructor assignment operator and helpers
  44 //////////
  45
  46 // public
  47
  48 power::power() : basic(TINFO_power)
  49 {
  50     debugmsg("power default constructor",LOGLEVEL_CONSTRUCT);
  51 }
  52
  53 power::~power()
  54 {
  55     debugmsg("power destructor",LOGLEVEL_DESTRUCT);
  56     destroy(0);
  57 }
  58
  59 power::power(power const & other)
  60 {
  61     debugmsg("power copy constructor",LOGLEVEL_CONSTRUCT);
  62     copy(other);
  63 }
  64
  65 power const & power::operator=(power const & other)
  66 {
  67     debugmsg("power operator=",LOGLEVEL_ASSIGNMENT);
  68     if (this != &other) {
  69         destroy(1);
  70         copy(other);
  71     }
  72     return *this;
  73 }
  74
  75 // protected
  76
  77 void power::copy(power const & other)
  78 {
  79     basic::copy(other);
  80     basis=other.basis;
  81     exponent=other.exponent;
  82 }
  83
  84 void power::destroy(bool call_parent)
  85 {
  86     if (call_parent) basic::destroy(call_parent);
  87 }
  88
  89 //////////
  90 // other constructors
  91 //////////
  92
  93 // public
  94
  95 power::power(ex const & lh, ex const & rh) : basic(TINFO_power), basis(lh), exponent(rh)
  96 {
  97     debugmsg("power constructor from ex,ex",LOGLEVEL_CONSTRUCT);
  98     GINAC_ASSERT(basis.return_type()==return_types::commutative);
  99 }
 100
 101 power::power(ex const & lh, numeric const & rh) : basic(TINFO_power), basis(lh), exponent(rh)
 102 {
 103     debugmsg("power constructor from ex,numeric",LOGLEVEL_CONSTRUCT);
 104     GINAC_ASSERT(basis.return_type()==return_types::commutative);
 105 }
 106
 107 //////////
 108 // functions overriding virtual functions from bases classes
 109 //////////
 110
 111 // public
 112
 113 basic * power::duplicate() const
 114 {
 115     debugmsg("power duplicate",LOGLEVEL_DUPLICATE);
 116     return new power(*this);
 117 }
 118
 119 bool power::info(unsigned inf) const
 120 {
 121     if (inf==info_flags::polynomial || inf==info_flags::integer_polynomial || inf==info_flags::rational_polynomial) {
 122         return exponent.info(info_flags::nonnegint);
 123     } else if (inf==info_flags::rational_function) {
 124         return exponent.info(info_flags::integer);
 125     } else {
 126         return basic::info(inf);
 127     }
 128 }
 129
 130 int power::nops() const
 131 {
 132     return 2;
 133 }
 134
 135 ex & power::let_op(int const i)
 136 {
 137     GINAC_ASSERT(i>=0);
 138     GINAC_ASSERT(i<2);
 139
 140     return i==0 ? basis : exponent;
 141 }
 142
 143 int power::degree(symbol const & s) const
 144 {
 145     if (is_exactly_of_type(*exponent.bp,numeric)) {
 146         if ((*basis.bp).compare(s)==0)
 147             return ex_to_numeric(exponent).to_int();
 148         else
 149             return basis.degree(s) * ex_to_numeric(exponent).to_int();
 150     }
 151     return 0;
 152 }
 153
 154 int power::ldegree(symbol const & s) const
 155 {
 156     if (is_exactly_of_type(*exponent.bp,numeric)) {
 157         if ((*basis.bp).compare(s)==0)
 158             return ex_to_numeric(exponent).to_int();
 159         else
 160             return basis.ldegree(s) * ex_to_numeric(exponent).to_int();
 161     }
 162     return 0;
 163 }
 164
 165 ex power::coeff(symbol const & s, int const n) const
 166 {
 167     if ((*basis.bp).compare(s)!=0) {
 168         // basis not equal to s
 169         if (n==0) {
 170             return *this;
 171         } else {
 172             return exZERO();
 173         }
 174     } else if (is_exactly_of_type(*exponent.bp,numeric)&&
 175                (static_cast<numeric const &>(*exponent.bp).compare(numeric(n))==0)) {
 176         return exONE();
 177     }
 178
 179     return exZERO();
 180 }
 181
 182 ex power::eval(int level) const
 183 {
 184     // simplifications: ^(x,0) -> 1 (0^0 handled here)
 185     //                  ^(x,1) -> x
 186     //                  ^(0,x) -> 0 (except if x is real and negative, in which case an exception is thrown)
 187     //                  ^(1,x) -> 1
 188     //                  ^(c1,c2) -> *(c1^n,c1^(c2-n)) (c1, c2 numeric(), 0<(c2-n)<1 except if c1,c2 are rational, but c1^c2 is not)
 189     //                  ^(^(x,c1),c2) -> ^(x,c1*c2) (c1, c2 numeric(), c2 integer or -1 < c1 <= 1, case c1=1 should not happen, see below!)
 190     //                  ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
 191     //                  ^(*(x,c1),c2) -> ^(x,c2)*c1^c2 (c1, c2 numeric(), c1>0)
 192     //                  ^(*(x,c1),c2) -> ^(-x,c2)*c1^c2 (c1, c2 numeric(), c1<0)
 193
 194     debugmsg("power eval",LOGLEVEL_MEMBER_FUNCTION);
 195
 196     if ((level==1)&&(flags & status_flags::evaluated)) {
 197         return *this;
 198     } else if (level == -max_recursion_level) {
 199         throw(std::runtime_error("max recursion level reached"));
 200     }
 201
 202     ex const & ebasis    = level==1 ? basis    : basis.eval(level-1);
 203     ex const & eexponent = level==1 ? exponent : exponent.eval(level-1);
 204
 205     bool basis_is_numerical=0;
 206     bool exponent_is_numerical=0;
 207     numeric * num_basis;
 208     numeric * num_exponent;
 209
 210     if (is_exactly_of_type(*ebasis.bp,numeric)) {
 211         basis_is_numerical=1;
 212         num_basis=static_cast<numeric *>(ebasis.bp);
 213     }
 214     if (is_exactly_of_type(*eexponent.bp,numeric)) {
 215         exponent_is_numerical=1;
 216         num_exponent=static_cast<numeric *>(eexponent.bp);
 217     }
 218
 219     // ^(x,0) -> 1 (0^0 also handled here)
 220     if (eexponent.is_zero())
 221         return exONE();
 222
 223     // ^(x,1) -> x
 224     if (eexponent.is_equal(exONE()))
 225         return ebasis;
 226
 227     // ^(0,x) -> 0 (except if x is real and negative)
 228     if (ebasis.is_zero()) {
 229         if (exponent_is_numerical && num_exponent->is_negative()) {
 230             throw(std::overflow_error("power::eval(): division by zero"));
 231         } else
 232             return exZERO();
 233     }
 234
 235     // ^(1,x) -> 1
 236     if (ebasis.is_equal(exONE()))
 237         return exONE();
 238
 239     if (basis_is_numerical && exponent_is_numerical) {
 240         // ^(c1,c2) -> c1^c2 (c1, c2 numeric(),
 241         // except if c1,c2 are rational, but c1^c2 is not)
 242         bool basis_is_rational = num_basis->is_rational();
 243         bool exponent_is_rational = num_exponent->is_rational();
 244         numeric res = (*num_basis).power(*num_exponent);
 245
 246         if ((!basis_is_rational || !exponent_is_rational)
 247             || res.is_rational()) {
 248             return res;
 249         }
 250         GINAC_ASSERT(!num_exponent->is_integer());  // has been handled by now
 251         // ^(c1,n/m) -> *(c1^q,c1^(n/m-q)), 0<(n/m-h)<1, q integer
 252         if (basis_is_rational && exponent_is_rational
 253             && num_exponent->is_real()
 254             && !num_exponent->is_integer()) {
 255             numeric r, q, n, m;
 256             n = num_exponent->numer();
 257             m = num_exponent->denom();
 258             q = iquo(n, m, r);
 259             if (r.is_negative()) {
 260                 r = r.add(m);
 261                 q = q.sub(numONE());
 262             }
 263             if (q.is_zero())  // the exponent was in the allowed range 0<(n/m)<1
 264                 return this->hold();
 265             else {
 266                 epvector res(2);
 267                 res.push_back(expair(ebasis,r.div(m)));
 268                 res.push_back(expair(ex(num_basis->power(q)),exONE()));
 269                 return (new mul(res))->setflag(status_flags::dynallocated | status_flags::evaluated);
 270                 /*return mul(num_basis->power(q),
 271                            power(ex(*num_basis),ex(r.div(m)))).hold();
 272                 */
 273                 /* return (new mul(num_basis->power(q),
 274                    power(*num_basis,r.div(m)).hold()))->setflag(status_flags::dynallocated | status_flags::evaluated);
 275                 */
 276             }
 277         }
 278     }
 279
 280     // ^(^(x,c1),c2) -> ^(x,c1*c2)
 281     // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1,
 282     // case c1=1 should not happen, see below!)
 283     if (exponent_is_numerical && is_ex_exactly_of_type(ebasis,power)) {
 284         power const & sub_power=ex_to_power(ebasis);
 285         ex const & sub_basis=sub_power.basis;
 286         ex const & sub_exponent=sub_power.exponent;
 287         if (is_ex_exactly_of_type(sub_exponent,numeric)) {
 288             numeric const & num_sub_exponent=ex_to_numeric(sub_exponent);
 289             GINAC_ASSERT(num_sub_exponent!=numeric(1));
 290             if (num_exponent->is_integer() || abs(num_sub_exponent)<1) {
 291                 return power(sub_basis,num_sub_exponent.mul(*num_exponent));
 292             }
 293         }
 294     }
 295
 296     // ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
 297     if (exponent_is_numerical && num_exponent->is_integer() &&
 298         is_ex_exactly_of_type(ebasis,mul)) {
 299         return expand_mul(ex_to_mul(ebasis), *num_exponent);
 300     }
 301
 302     // ^(*(...,x;c1),c2) -> ^(*(...,x;1),c2)*c1^c2 (c1, c2 numeric(), c1>0)
 303     // ^(*(...,x,c1),c2) -> ^(*(...,x;-1),c2)*(-c1)^c2 (c1, c2 numeric(), c1<0)
 304     if (exponent_is_numerical && is_ex_exactly_of_type(ebasis,mul)) {
 305         GINAC_ASSERT(!num_exponent->is_integer()); // should have been handled above
 306         mul const & mulref=ex_to_mul(ebasis);
 307         if (!mulref.overall_coeff.is_equal(exONE())) {
 308             numeric const & num_coeff=ex_to_numeric(mulref.overall_coeff);
 309             if (num_coeff.is_real()) {
 310                 if (num_coeff.is_positive()>0) {
 311                     mul * mulp=new mul(mulref);
 312                     mulp->overall_coeff=exONE();
 313                     mulp->clearflag(status_flags::evaluated);
 314                     mulp->clearflag(status_flags::hash_calculated);
 315                     return (new mul(power(*mulp,exponent),
 316                                     power(num_coeff,*num_exponent)))->
 317                         setflag(status_flags::dynallocated);
 318                 } else {
 319                     GINAC_ASSERT(num_coeff.compare(numZERO())<0);
 320                     if (num_coeff.compare(numMINUSONE())!=0) {
 321                         mul * mulp=new mul(mulref);
 322                         mulp->overall_coeff=exMINUSONE();
 323                         mulp->clearflag(status_flags::evaluated);
 324                         mulp->clearflag(status_flags::hash_calculated);
 325                         return (new mul(power(*mulp,exponent),
 326                                         power(abs(num_coeff),*num_exponent)))->
 327                             setflag(status_flags::dynallocated);
 328                     }
 329                 }
 330             }
 331         }
 332     }
 333
 334     if (are_ex_trivially_equal(ebasis,basis) &&
 335         are_ex_trivially_equal(eexponent,exponent)) {
 336         return this->hold();
 337     }
 338     return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated |
 339                                                    status_flags::evaluated);
 340 }
 341
 342 ex power::evalf(int level) const
 343 {
 344     debugmsg("power evalf",LOGLEVEL_MEMBER_FUNCTION);
 345
 346     ex ebasis;
 347     ex eexponent;
 348
 349     if (level==1) {
 350         ebasis=basis;
 351         eexponent=exponent;
 352     } else if (level == -max_recursion_level) {
 353         throw(std::runtime_error("max recursion level reached"));
 354     } else {
 355         ebasis=basis.evalf(level-1);
 356         eexponent=exponent.evalf(level-1);
 357     }
 358
 359     return power(ebasis,eexponent);
 360 }
 361
 362 ex power::subs(lst const & ls, lst const & lr) const
 363 {
 364     ex const & subsed_basis=basis.subs(ls,lr);
 365     ex const & subsed_exponent=exponent.subs(ls,lr);
 366
 367     if (are_ex_trivially_equal(basis,subsed_basis)&&
 368         are_ex_trivially_equal(exponent,subsed_exponent)) {
 369         return *this;
 370     }
 371
 372     return power(subsed_basis, subsed_exponent);
 373 }
 374
 375 ex power::simplify_ncmul(exvector const & v) const
 376 {
 377     return basic::simplify_ncmul(v);
 378 }
 379
 380 // protected
 381
 382 int power::compare_same_type(basic const & other) const
 383 {
 384     GINAC_ASSERT(is_exactly_of_type(other, power));
 385     power const & o=static_cast<power const &>(const_cast<basic &>(other));
 386
 387     int cmpval;
 388     cmpval=basis.compare(o.basis);
 389     if (cmpval==0) {
 390         return exponent.compare(o.exponent);
 391     }
 392     return cmpval;
 393 }
 394
 395 unsigned power::return_type(void) const
 396 {
 397     return basis.return_type();
 398 }
 399
 400 unsigned power::return_type_tinfo(void) const
 401 {
 402     return basis.return_type_tinfo();
 403 }
 404
 405 ex power::expand(unsigned options) const
 406 {
 407     ex expanded_basis=basis.expand(options);
 408
 409     if (!is_ex_exactly_of_type(exponent,numeric)||
 410         !ex_to_numeric(exponent).is_integer()) {
 411         if (are_ex_trivially_equal(basis,expanded_basis)) {
 412             return this->hold();
 413         } else {
 414             return (new power(expanded_basis,exponent))->
 415                     setflag(status_flags::dynallocated);
 416         }
 417     }
 418
 419     // integer numeric exponent
 420     numeric const & num_exponent=ex_to_numeric(exponent);
 421     int int_exponent = num_exponent.to_int();
 422
 423     if (int_exponent > 0 && is_ex_exactly_of_type(expanded_basis,add)) {
 424         return expand_add(ex_to_add(expanded_basis), int_exponent);
 425     }
 426
 427     if (is_ex_exactly_of_type(expanded_basis,mul)) {
 428         return expand_mul(ex_to_mul(expanded_basis), num_exponent);
 429     }
 430
 431     // cannot expand further
 432     if (are_ex_trivially_equal(basis,expanded_basis)) {
 433         return this->hold();
 434     } else {
 435         return (new power(expanded_basis,exponent))->
 436                setflag(status_flags::dynallocated);
 437     }
 438 }
 439
 440 //////////
 441 // new virtual functions which can be overridden by derived classes
 442 //////////
 443
 444 // none
 445
 446 //////////
 447 // non-virtual functions in this class
 448 //////////
 449
 450 ex power::expand_add(add const & a, int const n) const
 451 {
 452     // expand a^n where a is an add and n is an integer
 453
 454     if (n==2) {
 455         return expand_add_2(a);
 456     }
 457
 458     int m=a.nops();
 459     exvector sum;
 460     sum.reserve((n+1)*(m-1));
 461     intvector k(m-1);
 462     intvector k_cum(m-1); // k_cum[l]:=sum(i=0,l,k[l]);
 463     intvector upper_limit(m-1);
 464     int l;
 465
 466     for (int l=0; l<m-1; l++) {
 467         k[l]=0;
 468         k_cum[l]=0;
 469         upper_limit[l]=n;
 470     }
 471
 472     while (1) {
 473         exvector term;
 474         term.reserve(m+1);
 475         for (l=0; l<m-1; l++) {
 476             ex const & b=a.op(l);
 477             GINAC_ASSERT(!is_ex_exactly_of_type(b,add));
 478             GINAC_ASSERT(!is_ex_exactly_of_type(b,power)||
 479                    !is_ex_exactly_of_type(ex_to_power(b).exponent,numeric)||
 480                    !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer());
 481             if (is_ex_exactly_of_type(b,mul)) {
 482                 term.push_back(expand_mul(ex_to_mul(b),numeric(k[l])));
 483             } else {
 484                 term.push_back(power(b,k[l]));
 485             }
 486         }
 487
 488         ex const & b=a.op(l);
 489         GINAC_ASSERT(!is_ex_exactly_of_type(b,add));
 490         GINAC_ASSERT(!is_ex_exactly_of_type(b,power)||
 491                !is_ex_exactly_of_type(ex_to_power(b).exponent,numeric)||
 492                !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer());
 493         if (is_ex_exactly_of_type(b,mul)) {
 494             term.push_back(expand_mul(ex_to_mul(b),numeric(n-k_cum[m-2])));
 495         } else {
 496             term.push_back(power(b,n-k_cum[m-2]));
 497         }
 498
 499         numeric f=binomial(numeric(n),numeric(k[0]));
 500         for (l=1; l<m-1; l++) {
 501             f=f*binomial(numeric(n-k_cum[l-1]),numeric(k[l]));
 502         }
 503         term.push_back(f);
 504
 505         /*
 506         cout << "begin term" << endl;
 507         for (int i=0; i<m-1; i++) {
 508             cout << "k[" << i << "]=" << k[i] << endl;
 509             cout << "k_cum[" << i << "]=" << k_cum[i] << endl;
 510             cout << "upper_limit[" << i << "]=" << upper_limit[i] << endl;
 511         }
 512         for (exvector::const_iterator cit=term.begin(); cit!=term.end(); ++cit) {
 513             cout << *cit << endl;
 514         }
 515         cout << "end term" << endl;
 516         */
 517
 518         // TODO: optimize this
 519         sum.push_back((new mul(term))->setflag(status_flags::dynallocated));
 520
 521         // increment k[]
 522         l=m-2;
 523         while ((l>=0)&&((++k[l])>upper_limit[l])) {
 524             k[l]=0;
 525             l--;
 526         }
 527         if (l<0) break;
 528
 529         // recalc k_cum[] and upper_limit[]
 530         if (l==0) {
 531             k_cum[0]=k[0];
 532         } else {
 533             k_cum[l]=k_cum[l-1]+k[l];
 534         }
 535         for (int i=l+1; i<m-1; i++) {
 536             k_cum[i]=k_cum[i-1]+k[i];
 537         }
 538
 539         for (int i=l+1; i<m-1; i++) {
 540             upper_limit[i]=n-k_cum[i-1];
 541         }
 542     }
 543     return (new add(sum))->setflag(status_flags::dynallocated);
 544 }
 545
 546 /*
 547 ex power::expand_add_2(add const & a) const
 548 {
 549     // special case: expand a^2 where a is an add
 550
 551     epvector sum;
 552     sum.reserve((a.seq.size()*(a.seq.size()+1))/2);
 553     epvector::const_iterator last=a.seq.end();
 554
 555     for (epvector::const_iterator cit0=a.seq.begin(); cit0!=last; ++cit0) {
 556         ex const & b=a.recombine_pair_to_ex(*cit0);
 557         GINAC_ASSERT(!is_ex_exactly_of_type(b,add));
 558         GINAC_ASSERT(!is_ex_exactly_of_type(b,power)||
 559                !is_ex_exactly_of_type(ex_to_power(b).exponent,numeric)||
 560                !ex_to_numeric(ex_to_power(b).exponent).is_pos_integer());
 561         if (is_ex_exactly_of_type(b,mul)) {
 562             sum.push_back(a.split_ex_to_pair(expand_mul(ex_to_mul(b),numTWO())));
 563         } else {
 564             sum.push_back(a.split_ex_to_pair((new power(b,exTWO()))->
 565                                               setflag(status_flags::dynallocated)));
 566         }
 567         for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) {
 568             sum.push_back(a.split_ex_to_pair((new mul(a.recombine_pair_to_ex(*cit0),
 569                                                       a.recombine_pair_to_ex(*cit1)))->
 570                                               setflag(status_flags::dynallocated),
 571                                              exTWO()));
 572         }
 573     }
 574
 575     GINAC_ASSERT(sum.size()==(a.seq.size()*(a.seq.size()+1))/2);
 576
 577     return (new add(sum))->setflag(status_flags::dynallocated);
 578 }
 579 */
 580
 581 ex power::expand_add_2(add const & a) const
 582 {
 583     // special case: expand a^2 where a is an add
 584
 585     epvector sum;
 586     unsigned a_nops=a.nops();
 587     sum.reserve((a_nops*(a_nops+1))/2);
 588     epvector::const_iterator last=a.seq.end();
 589
 590     // power(+(x,...,z;c),2)=power(+(x,...,z;0),2)+2*c*+(x,...,z;0)+c*c
 591     // first part: ignore overall_coeff and expand other terms
 592     for (epvector::const_iterator cit0=a.seq.begin(); cit0!=last; ++cit0) {
 593         ex const & r=(*cit0).rest;
 594         ex const & c=(*cit0).coeff;
 595
 596         GINAC_ASSERT(!is_ex_exactly_of_type(r,add));
 597         GINAC_ASSERT(!is_ex_exactly_of_type(r,power)||
 598                !is_ex_exactly_of_type(ex_to_power(r).exponent,numeric)||
 599                !ex_to_numeric(ex_to_power(r).exponent).is_pos_integer()||
 600                !is_ex_exactly_of_type(ex_to_power(r).basis,add)||
 601                !is_ex_exactly_of_type(ex_to_power(r).basis,mul)||
 602                !is_ex_exactly_of_type(ex_to_power(r).basis,power));
 603
 604         if (are_ex_trivially_equal(c,exONE())) {
 605             if (is_ex_exactly_of_type(r,mul)) {
 606                 sum.push_back(expair(expand_mul(ex_to_mul(r),numTWO()),exONE()));
 607             } else {
 608                 sum.push_back(expair((new power(r,exTWO()))->setflag(status_flags::dynallocated),
 609                                      exONE()));
 610             }
 611         } else {
 612             if (is_ex_exactly_of_type(r,mul)) {
 613                 sum.push_back(expair(expand_mul(ex_to_mul(r),numTWO()),
 614                                      ex_to_numeric(c).power_dyn(numTWO())));
 615             } else {
 616                 sum.push_back(expair((new power(r,exTWO()))->setflag(status_flags::dynallocated),
 617                                      ex_to_numeric(c).power_dyn(numTWO())));
 618             }
 619         }
 620
 621         for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) {
 622             ex const & r1=(*cit1).rest;
 623             ex const & c1=(*cit1).coeff;
 624             sum.push_back(a.combine_ex_with_coeff_to_pair((new mul(r,r1))->setflag(status_flags::dynallocated),
 625                                                           numTWO().mul(ex_to_numeric(c)).mul_dyn(ex_to_numeric(c1))));
 626         }
 627     }
 628
 629     GINAC_ASSERT(sum.size()==(a.seq.size()*(a.seq.size()+1))/2);
 630
 631     // second part: add terms coming from overall_factor (if != 0)
 632     if (!a.overall_coeff.is_equal(exZERO())) {
 633         for (epvector::const_iterator cit=a.seq.begin(); cit!=a.seq.end(); ++cit) {
 634             sum.push_back(a.combine_pair_with_coeff_to_pair(*cit,ex_to_numeric(a.overall_coeff).mul_dyn(numTWO())));
 635         }
 636         sum.push_back(expair(ex_to_numeric(a.overall_coeff).power_dyn(numTWO()),exONE()));
 637     }
 638
 639     GINAC_ASSERT(sum.size()==(a_nops*(a_nops+1))/2);
 640
 641     return (new add(sum))->setflag(status_flags::dynallocated);
 642 }
 643
 644 ex power::expand_mul(mul const & m, numeric const & n) const
 645 {
 646     // expand m^n where m is a mul and n is and integer
 647
 648     if (n.is_equal(numZERO())) {
 649         return exONE();
 650     }
 651
 652     epvector distrseq;
 653     distrseq.reserve(m.seq.size());
 654     epvector::const_iterator last=m.seq.end();
 655     epvector::const_iterator cit=m.seq.begin();
 656     while (cit!=last) {
 657         if (is_ex_exactly_of_type((*cit).rest,numeric)) {
 658             distrseq.push_back(m.combine_pair_with_coeff_to_pair(*cit,n));
 659         } else {
 660             // it is safe not to call mul::combine_pair_with_coeff_to_pair()
 661             // since n is an integer
 662             distrseq.push_back(expair((*cit).rest,
 663                                       ex_to_numeric((*cit).coeff).mul(n)));
 664         }
 665         ++cit;
 666     }
 667     return (new mul(distrseq,ex_to_numeric(m.overall_coeff).power_dyn(n)))
 668                  ->setflag(status_flags::dynallocated);
 669 }
 670
 671 /*
 672 ex power::expand_commutative_3(ex const & basis, numeric const & exponent,
 673                              unsigned options) const
 674 {
 675     // obsolete
 676
 677     exvector distrseq;
 678     epvector splitseq;
 679
 680     add const & addref=static_cast<add const &>(*basis.bp);
 681
 682     splitseq=addref.seq;
 683     splitseq.pop_back();
 684     ex first_operands=add(splitseq);
 685     ex last_operand=addref.recombine_pair_to_ex(*(addref.seq.end()-1));
 686
 687     int n=exponent.to_int();
 688     for (int k=0; k<=n; k++) {
 689         distrseq.push_back(binomial(n,k)*power(first_operands,numeric(k))*
 690                            power(last_operand,numeric(n-k)));
 691     }
 692     return ex((new add(distrseq))->setflag(status_flags::sub_expanded |
 693                                            status_flags::expanded |
 694                                            status_flags::dynallocated  )).
 695            expand(options);
 696 }
 697 */
 698
 699 /*
 700 ex power::expand_noncommutative(ex const & basis, numeric const & exponent,
 701                                 unsigned options) const
 702 {
 703     ex rest_power=ex(power(basis,exponent.add(numMINUSONE()))).
 704                   expand(options | expand_options::internal_do_not_expand_power_operands);
 705
 706     return ex(mul(rest_power,basis),0).
 707            expand(options | expand_options::internal_do_not_expand_mul_operands);
 708 }
 709 */
 710
 711 //////////
 712 // static member variables
 713 //////////
 714
 715 // protected
 716
 717 unsigned power::precedence=60;
 718
 719 //////////
 720 // global constants
 721 //////////
 722
 723 const power some_power;
 724 type_info const & typeid_power=typeid(some_power);
 725
 726 #ifndef NO_GINAC_NAMESPACE
 727 } // namespace GiNaC
 728 #endif // ndef NO_GINAC_NAMESPACE