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