ginac/power.cpp

   1 /** @file power.cpp
   2  *
   3  *  Implementation of GiNaC's symbolic exponentiation (basis^exponent). */
   4
   5 /*
   6  *  GiNaC Copyright (C) 1999-2002 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 "ncmul.h"
  32 #include "numeric.h"
  33 #include "constant.h"
  34 #include "inifcns.h" // for log() in power::derivative()
  35 #include "matrix.h"
  36 #include "indexed.h"
  37 #include "symbol.h"
  38 #include "print.h"
  39 #include "archive.h"
  40 #include "utils.h"
  41
  42 namespace GiNaC {
  43
  44 GINAC_IMPLEMENT_REGISTERED_CLASS(power, basic)
  45
  46 typedef std::vector<int> intvector;
  47
  48 //////////
  49 // default ctor, dtor, copy ctor, assignment operator and helpers
  50 //////////
  51
  52 power::power() : inherited(TINFO_power) { }
  53
  54 void power::copy(const power & other)
  55 {
  56         inherited::copy(other);
  57         basis = other.basis;
  58         exponent = other.exponent;
  59 }
  60
  61 DEFAULT_DESTROY(power)
  62
  63 //////////
  64 // other ctors
  65 //////////
  66
  67 // all inlined
  68
  69 //////////
  70 // archiving
  71 //////////
  72
  73 power::power(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
  74 {
  75         n.find_ex("basis", basis, sym_lst);
  76         n.find_ex("exponent", exponent, sym_lst);
  77 }
  78
  79 void power::archive(archive_node &n) const
  80 {
  81         inherited::archive(n);
  82         n.add_ex("basis", basis);
  83         n.add_ex("exponent", exponent);
  84 }
  85
  86 DEFAULT_UNARCHIVE(power)
  87
  88 //////////
  89 // functions overriding virtual functions from base classes
  90 //////////
  91
  92 // public
  93
  94 static void print_sym_pow(const print_context & c, const symbol &x, int exp)
  95 {
  96         // Optimal output of integer powers of symbols to aid compiler CSE.
  97         // C.f. ISO/IEC 14882:1998, section 1.9 [intro execution], paragraph 15
  98         // to learn why such a parenthisation is really necessary.
  99         if (exp == 1) {
 100                 x.print(c);
 101         } else if (exp == 2) {
 102                 x.print(c);
 103                 c.s << "*";
 104                 x.print(c);
 105         } else if (exp & 1) {
 106                 x.print(c);
 107                 c.s << "*";
 108                 print_sym_pow(c, x, exp-1);
 109         } else {
 110                 c.s << "(";
 111                 print_sym_pow(c, x, exp >> 1);
 112                 c.s << ")*(";
 113                 print_sym_pow(c, x, exp >> 1);
 114                 c.s << ")";
 115         }
 116 }
 117
 118 void power::print(const print_context & c, unsigned level) const
 119 {
 120         if (is_a<print_tree>(c)) {
 121
 122                 inherited::print(c, level);
 123
 124         } else if (is_a<print_csrc>(c)) {
 125
 126                 // Integer powers of symbols are printed in a special, optimized way
 127                 if (exponent.info(info_flags::integer)
 128                  && (is_exactly_a<symbol>(basis) || is_exactly_a<constant>(basis))) {
 129                         int exp = ex_to<numeric>(exponent).to_int();
 130                         if (exp > 0)
 131                                 c.s << '(';
 132                         else {
 133                                 exp = -exp;
 134                                 if (is_a<print_csrc_cl_N>(c))
 135                                         c.s << "recip(";
 136                                 else
 137                                         c.s << "1.0/(";
 138                         }
 139                         print_sym_pow(c, ex_to<symbol>(basis), exp);
 140                         c.s << ')';
 141
 142                 // <expr>^-1 is printed as "1.0/<expr>" or with the recip() function of CLN
 143                 } else if (exponent.is_equal(_ex_1)) {
 144                         if (is_a<print_csrc_cl_N>(c))
 145                                 c.s << "recip(";
 146                         else
 147                                 c.s << "1.0/(";
 148                         basis.print(c);
 149                         c.s << ')';
 150
 151                 // Otherwise, use the pow() or expt() (CLN) functions
 152                 } else {
 153                         if (is_a<print_csrc_cl_N>(c))
 154                                 c.s << "expt(";
 155                         else
 156                                 c.s << "pow(";
 157                         basis.print(c);
 158                         c.s << ',';
 159                         exponent.print(c);
 160                         c.s << ')';
 161                 }
 162
 163         } else if (is_a<print_python_repr>(c)) {
 164
 165                 c.s << class_name() << '(';
 166                 basis.print(c);
 167                 c.s << ',';
 168                 exponent.print(c);
 169                 c.s << ')';
 170
 171         } else {
 172
 173                 bool is_tex = is_a<print_latex>(c);
 174
 175                 if (is_tex && is_a<numeric>(exponent) && ex_to<numeric>(exponent).is_negative()) {
 176
 177                         // Powers with negative numeric exponents are printed as fractions in TeX
 178                         c.s << "\\frac{1}{";
 179                         power(basis, -exponent).eval().print(c);
 180                         c.s << "}";
 181
 182                 } else if (exponent.is_equal(_ex1_2)) {
 183
 184                         // Square roots are printed in a special way
 185                         c.s << (is_tex ? "\\sqrt{" : "sqrt(");
 186                         basis.print(c);
 187                         c.s << (is_tex ? '}' : ')');
 188
 189                 } else {
 190
 191                         // Ordinary output of powers using '^' or '**'
 192                         if (precedence() <= level)
 193                                 c.s << (is_tex ? "{(" : "(");
 194                         basis.print(c, precedence());
 195                         if (is_a<print_python>(c))
 196                                 c.s << "**";
 197                         else
 198                                 c.s << '^';
 199                         if (is_tex)
 200                                 c.s << '{';
 201                         exponent.print(c, precedence());
 202                         if (is_tex)
 203                                 c.s << '}';
 204                         if (precedence() <= level)
 205                                 c.s << (is_tex ? ")}" : ")");
 206                 }
 207         }
 208 }
 209
 210 bool power::info(unsigned inf) const
 211 {
 212         switch (inf) {
 213                 case info_flags::polynomial:
 214                 case info_flags::integer_polynomial:
 215                 case info_flags::cinteger_polynomial:
 216                 case info_flags::rational_polynomial:
 217                 case info_flags::crational_polynomial:
 218                         return exponent.info(info_flags::nonnegint);
 219                 case info_flags::rational_function:
 220                         return exponent.info(info_flags::integer);
 221                 case info_flags::algebraic:
 222                         return (!exponent.info(info_flags::integer) ||
 223                                         basis.info(inf));
 224         }
 225         return inherited::info(inf);
 226 }
 227
 228 unsigned power::nops() const
 229 {
 230         return 2;
 231 }
 232
 233 ex & power::let_op(int i)
 234 {
 235         GINAC_ASSERT(i>=0);
 236         GINAC_ASSERT(i<2);
 237
 238         return i==0 ? basis : exponent;
 239 }
 240
 241 ex power::map(map_function & f) const
 242 {
 243         return (new power(f(basis), f(exponent)))->setflag(status_flags::dynallocated);
 244 }
 245
 246 int power::degree(const ex & s) const
 247 {
 248         if (is_equal(ex_to<basic>(s)))
 249                 return 1;
 250         else if (is_ex_exactly_of_type(exponent, numeric) && ex_to<numeric>(exponent).is_integer()) {
 251                 if (basis.is_equal(s))
 252                         return ex_to<numeric>(exponent).to_int();
 253                 else
 254                         return basis.degree(s) * ex_to<numeric>(exponent).to_int();
 255         } else if (basis.has(s))
 256                 throw(std::runtime_error("power::degree(): undefined degree because of non-integer exponent"));
 257         else
 258                 return 0;
 259 }
 260
 261 int power::ldegree(const ex & s) const
 262 {
 263         if (is_equal(ex_to<basic>(s)))
 264                 return 1;
 265         else if (is_ex_exactly_of_type(exponent, numeric) && ex_to<numeric>(exponent).is_integer()) {
 266                 if (basis.is_equal(s))
 267                         return ex_to<numeric>(exponent).to_int();
 268                 else
 269                         return basis.ldegree(s) * ex_to<numeric>(exponent).to_int();
 270         } else if (basis.has(s))
 271                 throw(std::runtime_error("power::ldegree(): undefined degree because of non-integer exponent"));
 272         else
 273                 return 0;
 274 }
 275
 276 ex power::coeff(const ex & s, int n) const
 277 {
 278         if (is_equal(ex_to<basic>(s)))
 279                 return n==1 ? _ex1 : _ex0;
 280         else if (!basis.is_equal(s)) {
 281                 // basis not equal to s
 282                 if (n == 0)
 283                         return *this;
 284                 else
 285                         return _ex0;
 286         } else {
 287                 // basis equal to s
 288                 if (is_ex_exactly_of_type(exponent, numeric) && ex_to<numeric>(exponent).is_integer()) {
 289                         // integer exponent
 290                         int int_exp = ex_to<numeric>(exponent).to_int();
 291                         if (n == int_exp)
 292                                 return _ex1;
 293                         else
 294                                 return _ex0;
 295                 } else {
 296                         // non-integer exponents are treated as zero
 297                         if (n == 0)
 298                                 return *this;
 299                         else
 300                                 return _ex0;
 301                 }
 302         }
 303 }
 304
 305 /** Perform automatic term rewriting rules in this class.  In the following
 306  *  x, x1, x2,... stand for a symbolic variables of type ex and c, c1, c2...
 307  *  stand for such expressions that contain a plain number.
 308  *  - ^(x,0) -> 1  (also handles ^(0,0))
 309  *  - ^(x,1) -> x
 310  *  - ^(0,c) -> 0 or exception  (depending on the real part of c)
 311  *  - ^(1,x) -> 1
 312  *  - ^(c1,c2) -> *(c1^n,c1^(c2-n))  (so that 0<(c2-n)<1, try to evaluate roots, possibly in numerator and denominator of c1)
 313  *  - ^(^(x,c1),c2) -> ^(x,c1*c2)  (c2 integer or -1 < c1 <= 1, case c1=1 should not happen, see below!)
 314  *  - ^(*(x,y,z),c) -> *(x^c,y^c,z^c)  (if c integer)
 315  *  - ^(*(x,c1),c2) -> ^(x,c2)*c1^c2  (c1>0)
 316  *  - ^(*(x,c1),c2) -> ^(-x,c2)*c1^c2  (c1<0)
 317  *
 318  *  @param level cut-off in recursive evaluation */
 319 ex power::eval(int level) const
 320 {
 321         if ((level==1) && (flags & status_flags::evaluated))
 322                 return *this;
 323         else if (level == -max_recursion_level)
 324                 throw(std::runtime_error("max recursion level reached"));
 325
 326         const ex & ebasis    = level==1 ? basis    : basis.eval(level-1);
 327         const ex & eexponent = level==1 ? exponent : exponent.eval(level-1);
 328
 329         bool basis_is_numerical = false;
 330         bool exponent_is_numerical = false;
 331         const numeric *num_basis;
 332         const numeric *num_exponent;
 333
 334         if (is_ex_exactly_of_type(ebasis, numeric)) {
 335                 basis_is_numerical = true;
 336                 num_basis = &ex_to<numeric>(ebasis);
 337         }
 338         if (is_ex_exactly_of_type(eexponent, numeric)) {
 339                 exponent_is_numerical = true;
 340                 num_exponent = &ex_to<numeric>(eexponent);
 341         }
 342
 343         // ^(x,0) -> 1  (0^0 also handled here)
 344         if (eexponent.is_zero()) {
 345                 if (ebasis.is_zero())
 346                         throw (std::domain_error("power::eval(): pow(0,0) is undefined"));
 347                 else
 348                         return _ex1;
 349         }
 350
 351         // ^(x,1) -> x
 352         if (eexponent.is_equal(_ex1))
 353                 return ebasis;
 354
 355         // ^(0,c1) -> 0 or exception  (depending on real value of c1)
 356         if (ebasis.is_zero() && exponent_is_numerical) {
 357                 if ((num_exponent->real()).is_zero())
 358                         throw (std::domain_error("power::eval(): pow(0,I) is undefined"));
 359                 else if ((num_exponent->real()).is_negative())
 360                         throw (pole_error("power::eval(): division by zero",1));
 361                 else
 362                         return _ex0;
 363         }
 364
 365         // ^(1,x) -> 1
 366         if (ebasis.is_equal(_ex1))
 367                 return _ex1;
 368
 369         if (exponent_is_numerical) {
 370
 371                 // ^(c1,c2) -> c1^c2  (c1, c2 numeric(),
 372                 // except if c1,c2 are rational, but c1^c2 is not)
 373                 if (basis_is_numerical) {
 374                         const bool basis_is_crational = num_basis->is_crational();
 375                         const bool exponent_is_crational = num_exponent->is_crational();
 376                         if (!basis_is_crational || !exponent_is_crational) {
 377                                 // return a plain float
 378                                 return (new numeric(num_basis->power(*num_exponent)))->setflag(status_flags::dynallocated |
 379                                                                                                status_flags::evaluated |
 380                                                                                                status_flags::expanded);
 381                         }
 382
 383                         const numeric res = num_basis->power(*num_exponent);
 384                         if (res.is_crational()) {
 385                                 return res;
 386                         }
 387                         GINAC_ASSERT(!num_exponent->is_integer());  // has been handled by now
 388
 389                         // ^(c1,n/m) -> *(c1^q,c1^(n/m-q)), 0<(n/m-q)<1, q integer
 390                         if (basis_is_crational && exponent_is_crational
 391                             && num_exponent->is_real()
 392                             && !num_exponent->is_integer()) {
 393                                 const numeric n = num_exponent->numer();
 394                                 const numeric m = num_exponent->denom();
 395                                 numeric r;
 396                                 numeric q = iquo(n, m, r);
 397                                 if (r.is_negative()) {
 398                                         r += m;
 399                                         --q;
 400                                 }
 401                                 if (q.is_zero()) {  // the exponent was in the allowed range 0<(n/m)<1
 402                                         if (num_basis->is_rational() && !num_basis->is_integer()) {
 403                                                 // try it for numerator and denominator separately, in order to
 404                                                 // partially simplify things like (5/8)^(1/3) -> 1/2*5^(1/3)
 405                                                 const numeric bnum = num_basis->numer();
 406                                                 const numeric bden = num_basis->denom();
 407                                                 const numeric res_bnum = bnum.power(*num_exponent);
 408                                                 const numeric res_bden = bden.power(*num_exponent);
 409                                                 if (res_bnum.is_integer())
 410                                                         return (new mul(power(bden,-*num_exponent),res_bnum))->setflag(status_flags::dynallocated | status_flags::evaluated);
 411                                                 if (res_bden.is_integer())
 412                                                         return (new mul(power(bnum,*num_exponent),res_bden.inverse()))->setflag(status_flags::dynallocated | status_flags::evaluated);
 413                                         }
 414                                         return this->hold();
 415                                 } else {
 416                                         // assemble resulting product, but allowing for a re-evaluation,
 417                                         // because otherwise we'll end up with something like
 418                                         //    (7/8)^(4/3)  ->  7/8*(1/2*7^(1/3))
 419                                         // instead of 7/16*7^(1/3).
 420                                         ex prod = power(*num_basis,r.div(m));
 421                                         return prod*power(*num_basis,q);
 422                                 }
 423                         }
 424                 }
 425
 426                 // ^(^(x,c1),c2) -> ^(x,c1*c2)
 427                 // (c1, c2 numeric(), c2 integer or -1 < c1 <= 1,
 428                 // case c1==1 should not happen, see below!)
 429                 if (is_ex_exactly_of_type(ebasis,power)) {
 430                         const power & sub_power = ex_to<power>(ebasis);
 431                         const ex & sub_basis = sub_power.basis;
 432                         const ex & sub_exponent = sub_power.exponent;
 433                         if (is_ex_exactly_of_type(sub_exponent,numeric)) {
 434                                 const numeric & num_sub_exponent = ex_to<numeric>(sub_exponent);
 435                                 GINAC_ASSERT(num_sub_exponent!=numeric(1));
 436                                 if (num_exponent->is_integer() || (abs(num_sub_exponent) - _num1).is_negative())
 437                                         return power(sub_basis,num_sub_exponent.mul(*num_exponent));
 438                         }
 439                 }
 440
 441                 // ^(*(x,y,z),c1) -> *(x^c1,y^c1,z^c1) (c1 integer)
 442                 if (num_exponent->is_integer() && is_ex_exactly_of_type(ebasis,mul)) {
 443                         return expand_mul(ex_to<mul>(ebasis), *num_exponent);
 444                 }
 445
 446                 // ^(*(...,x;c1),c2) -> *(^(*(...,x;1),c2),c1^c2)  (c1, c2 numeric(), c1>0)
 447                 // ^(*(...,x;c1),c2) -> *(^(*(...,x;-1),c2),(-c1)^c2)  (c1, c2 numeric(), c1<0)
 448                 if (is_ex_exactly_of_type(ebasis,mul)) {
 449                         GINAC_ASSERT(!num_exponent->is_integer()); // should have been handled above
 450                         const mul & mulref = ex_to<mul>(ebasis);
 451                         if (!mulref.overall_coeff.is_equal(_ex1)) {
 452                                 const numeric & num_coeff = ex_to<numeric>(mulref.overall_coeff);
 453                                 if (num_coeff.is_real()) {
 454                                         if (num_coeff.is_positive()) {
 455                                                 mul *mulp = new mul(mulref);
 456                                                 mulp->overall_coeff = _ex1;
 457                                                 mulp->clearflag(status_flags::evaluated);
 458                                                 mulp->clearflag(status_flags::hash_calculated);
 459                                                 return (new mul(power(*mulp,exponent),
 460                                                                 power(num_coeff,*num_exponent)))->setflag(status_flags::dynallocated);
 461                                         } else {
 462                                                 GINAC_ASSERT(num_coeff.compare(_num0)<0);
 463                                                 if (!num_coeff.is_equal(_num_1)) {
 464                                                         mul *mulp = new mul(mulref);
 465                                                         mulp->overall_coeff = _ex_1;
 466                                                         mulp->clearflag(status_flags::evaluated);
 467                                                         mulp->clearflag(status_flags::hash_calculated);
 468                                                         return (new mul(power(*mulp,exponent),
 469                                                                         power(abs(num_coeff),*num_exponent)))->setflag(status_flags::dynallocated);
 470                                                 }
 471                                         }
 472                                 }
 473                         }
 474                 }
 475
 476                 // ^(nc,c1) -> ncmul(nc,nc,...) (c1 positive integer, unless nc is a matrix)
 477                 if (num_exponent->is_pos_integer() &&
 478                     ebasis.return_type() != return_types::commutative &&
 479                     !is_ex_of_type(ebasis,matrix)) {
 480                         return ncmul(exvector(num_exponent->to_int(), ebasis), true);
 481                 }
 482         }
 483
 484         if (are_ex_trivially_equal(ebasis,basis) &&
 485             are_ex_trivially_equal(eexponent,exponent)) {
 486                 return this->hold();
 487         }
 488         return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated |
 489                                                        status_flags::evaluated);
 490 }
 491
 492 ex power::evalf(int level) const
 493 {
 494         ex ebasis;
 495         ex eexponent;
 496
 497         if (level==1) {
 498                 ebasis = basis;
 499                 eexponent = exponent;
 500         } else if (level == -max_recursion_level) {
 501                 throw(std::runtime_error("max recursion level reached"));
 502         } else {
 503                 ebasis = basis.evalf(level-1);
 504                 if (!is_exactly_a<numeric>(exponent))
 505                         eexponent = exponent.evalf(level-1);
 506                 else
 507                         eexponent = exponent;
 508         }
 509
 510         return power(ebasis,eexponent);
 511 }
 512
 513 ex power::evalm(void) const
 514 {
 515         const ex ebasis = basis.evalm();
 516         const ex eexponent = exponent.evalm();
 517         if (is_ex_of_type(ebasis,matrix)) {
 518                 if (is_ex_of_type(eexponent,numeric)) {
 519                         return (new matrix(ex_to<matrix>(ebasis).pow(eexponent)))->setflag(status_flags::dynallocated);
 520                 }
 521         }
 522         return (new power(ebasis, eexponent))->setflag(status_flags::dynallocated);
 523 }
 524
 525 ex power::subs(const lst & ls, const lst & lr, bool no_pattern) const
 526 {
 527         const ex &subsed_basis = basis.subs(ls, lr, no_pattern);
 528         const ex &subsed_exponent = exponent.subs(ls, lr, no_pattern);
 529
 530         if (are_ex_trivially_equal(basis, subsed_basis)
 531          && are_ex_trivially_equal(exponent, subsed_exponent))
 532                 return basic::subs(ls, lr, no_pattern);
 533         else
 534                 return power(subsed_basis, subsed_exponent).basic::subs(ls, lr, no_pattern);
 535 }
 536
 537 ex power::simplify_ncmul(const exvector & v) const
 538 {
 539         return inherited::simplify_ncmul(v);
 540 }
 541
 542 // protected
 543
 544 /** Implementation of ex::diff() for a power.
 545  *  @see ex::diff */
 546 ex power::derivative(const symbol & s) const
 547 {
 548         if (exponent.info(info_flags::real)) {
 549                 // D(b^r) = r * b^(r-1) * D(b) (faster than the formula below)
 550                 epvector newseq;
 551                 newseq.reserve(2);
 552                 newseq.push_back(expair(basis, exponent - _ex1));
 553                 newseq.push_back(expair(basis.diff(s), _ex1));
 554                 return mul(newseq, exponent);
 555         } else {
 556                 // D(b^e) = b^e * (D(e)*ln(b) + e*D(b)/b)
 557                 return mul(*this,
 558                            add(mul(exponent.diff(s), log(basis)),
 559                            mul(mul(exponent, basis.diff(s)), power(basis, _ex_1))));
 560         }
 561 }
 562
 563 int power::compare_same_type(const basic & other) const
 564 {
 565         GINAC_ASSERT(is_exactly_a<power>(other));
 566         const power &o = static_cast<const power &>(other);
 567
 568         int cmpval = basis.compare(o.basis);
 569         if (cmpval)
 570                 return cmpval;
 571         else
 572                 return exponent.compare(o.exponent);
 573 }
 574
 575 unsigned power::return_type(void) const
 576 {
 577         return basis.return_type();
 578 }
 579
 580 unsigned power::return_type_tinfo(void) const
 581 {
 582         return basis.return_type_tinfo();
 583 }
 584
 585 ex power::expand(unsigned options) const
 586 {
 587         if (options == 0 && (flags & status_flags::expanded))
 588                 return *this;
 589
 590         const ex expanded_basis = basis.expand(options);
 591         const ex expanded_exponent = exponent.expand(options);
 592
 593         // x^(a+b) -> x^a * x^b
 594         if (is_ex_exactly_of_type(expanded_exponent, add)) {
 595                 const add &a = ex_to<add>(expanded_exponent);
 596                 exvector distrseq;
 597                 distrseq.reserve(a.seq.size() + 1);
 598                 epvector::const_iterator last = a.seq.end();
 599                 epvector::const_iterator cit = a.seq.begin();
 600                 while (cit!=last) {
 601                         distrseq.push_back(power(expanded_basis, a.recombine_pair_to_ex(*cit)));
 602                         ++cit;
 603                 }
 604
 605                 // Make sure that e.g. (x+y)^(2+a) expands the (x+y)^2 factor
 606                 if (ex_to<numeric>(a.overall_coeff).is_integer()) {
 607                         const numeric &num_exponent = ex_to<numeric>(a.overall_coeff);
 608                         int int_exponent = num_exponent.to_int();
 609                         if (int_exponent > 0 && is_ex_exactly_of_type(expanded_basis, add))
 610                                 distrseq.push_back(expand_add(ex_to<add>(expanded_basis), int_exponent));
 611                         else
 612                                 distrseq.push_back(power(expanded_basis, a.overall_coeff));
 613                 } else
 614                         distrseq.push_back(power(expanded_basis, a.overall_coeff));
 615
 616                 // Make sure that e.g. (x+y)^(1+a) -> x*(x+y)^a + y*(x+y)^a
 617                 ex r = (new mul(distrseq))->setflag(status_flags::dynallocated);
 618                 return r.expand();
 619         }
 620
 621         if (!is_ex_exactly_of_type(expanded_exponent, numeric) ||
 622                 !ex_to<numeric>(expanded_exponent).is_integer()) {
 623                 if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent)) {
 624                         return this->hold();
 625                 } else {
 626                         return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
 627                 }
 628         }
 629
 630         // integer numeric exponent
 631         const numeric & num_exponent = ex_to<numeric>(expanded_exponent);
 632         int int_exponent = num_exponent.to_int();
 633
 634         // (x+y)^n, n>0
 635         if (int_exponent > 0 && is_ex_exactly_of_type(expanded_basis,add))
 636                 return expand_add(ex_to<add>(expanded_basis), int_exponent);
 637
 638         // (x*y)^n -> x^n * y^n
 639         if (is_ex_exactly_of_type(expanded_basis,mul))
 640                 return expand_mul(ex_to<mul>(expanded_basis), num_exponent);
 641
 642         // cannot expand further
 643         if (are_ex_trivially_equal(basis,expanded_basis) && are_ex_trivially_equal(exponent,expanded_exponent))
 644                 return this->hold();
 645         else
 646                 return (new power(expanded_basis,expanded_exponent))->setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
 647 }
 648
 649 //////////
 650 // new virtual functions which can be overridden by derived classes
 651 //////////
 652
 653 // none
 654
 655 //////////
 656 // non-virtual functions in this class
 657 //////////
 658
 659 /** expand a^n where a is an add and n is a positive integer.
 660  *  @see power::expand */
 661 ex power::expand_add(const add & a, int n) const
 662 {
 663         if (n==2)
 664                 return expand_add_2(a);
 665
 666         const int m = a.nops();
 667         exvector result;
 668         // The number of terms will be the number of combinatorial compositions,
 669         // i.e. the number of unordered arrangement of m nonnegative integers
 670         // which sum up to n.  It is frequently written as C_n(m) and directly
 671         // related with binomial coefficients:
 672         result.reserve(binomial(numeric(n+m-1), numeric(m-1)).to_int());
 673         intvector k(m-1);
 674         intvector k_cum(m-1); // k_cum[l]:=sum(i=0,l,k[l]);
 675         intvector upper_limit(m-1);
 676         int l;
 677
 678         for (int l=0; l<m-1; ++l) {
 679                 k[l] = 0;
 680                 k_cum[l] = 0;
 681                 upper_limit[l] = n;
 682         }
 683
 684         while (true) {
 685                 exvector term;
 686                 term.reserve(m+1);
 687                 for (l=0; l<m-1; ++l) {
 688                         const ex & b = a.op(l);
 689                         GINAC_ASSERT(!is_exactly_a<add>(b));
 690                         GINAC_ASSERT(!is_exactly_a<power>(b) ||
 691                                      !is_exactly_a<numeric>(ex_to<power>(b).exponent) ||
 692                                      !ex_to<numeric>(ex_to<power>(b).exponent).is_pos_integer() ||
 693                                      !is_exactly_a<add>(ex_to<power>(b).basis) ||
 694                                      !is_exactly_a<mul>(ex_to<power>(b).basis) ||
 695                                      !is_exactly_a<power>(ex_to<power>(b).basis));
 696                         if (is_ex_exactly_of_type(b,mul))
 697                                 term.push_back(expand_mul(ex_to<mul>(b),numeric(k[l])));
 698                         else
 699                                 term.push_back(power(b,k[l]));
 700                 }
 701
 702                 const ex & b = a.op(l);
 703                 GINAC_ASSERT(!is_exactly_a<add>(b));
 704                 GINAC_ASSERT(!is_exactly_a<power>(b) ||
 705                              !is_exactly_a<numeric>(ex_to<power>(b).exponent) ||
 706                              !ex_to<numeric>(ex_to<power>(b).exponent).is_pos_integer() ||
 707                              !is_exactly_a<add>(ex_to<power>(b).basis) ||
 708                              !is_exactly_a<mul>(ex_to<power>(b).basis) ||
 709                              !is_exactly_a<power>(ex_to<power>(b).basis));
 710                 if (is_ex_exactly_of_type(b,mul))
 711                         term.push_back(expand_mul(ex_to<mul>(b),numeric(n-k_cum[m-2])));
 712                 else
 713                         term.push_back(power(b,n-k_cum[m-2]));
 714
 715                 numeric f = binomial(numeric(n),numeric(k[0]));
 716                 for (l=1; l<m-1; ++l)
 717                         f *= binomial(numeric(n-k_cum[l-1]),numeric(k[l]));
 718
 719                 term.push_back(f);
 720
 721                 result.push_back((new mul(term))->setflag(status_flags::dynallocated));
 722
 723                 // increment k[]
 724                 l = m-2;
 725                 while ((l>=0) && ((++k[l])>upper_limit[l])) {
 726                         k[l] = 0;
 727                         --l;
 728                 }
 729                 if (l<0) break;
 730
 731                 // recalc k_cum[] and upper_limit[]
 732                 k_cum[l] = (l==0 ? k[0] : k_cum[l-1]+k[l]);
 733
 734                 for (int i=l+1; i<m-1; ++i)
 735                         k_cum[i] = k_cum[i-1]+k[i];
 736
 737                 for (int i=l+1; i<m-1; ++i)
 738                         upper_limit[i] = n-k_cum[i-1];
 739         }
 740
 741         return (new add(result))->setflag(status_flags::dynallocated |
 742                                           status_flags::expanded);
 743 }
 744
 745
 746 /** Special case of power::expand_add. Expands a^2 where a is an add.
 747  *  @see power::expand_add */
 748 ex power::expand_add_2(const add & a) const
 749 {
 750         epvector sum;
 751         unsigned a_nops = a.nops();
 752         sum.reserve((a_nops*(a_nops+1))/2);
 753         epvector::const_iterator last = a.seq.end();
 754
 755         // power(+(x,...,z;c),2)=power(+(x,...,z;0),2)+2*c*+(x,...,z;0)+c*c
 756         // first part: ignore overall_coeff and expand other terms
 757         for (epvector::const_iterator cit0=a.seq.begin(); cit0!=last; ++cit0) {
 758                 const ex & r = cit0->rest;
 759                 const ex & c = cit0->coeff;
 760
 761                 GINAC_ASSERT(!is_exactly_a<add>(r));
 762                 GINAC_ASSERT(!is_exactly_a<power>(r) ||
 763                              !is_exactly_a<numeric>(ex_to<power>(r).exponent) ||
 764                              !ex_to<numeric>(ex_to<power>(r).exponent).is_pos_integer() ||
 765                              !is_exactly_a<add>(ex_to<power>(r).basis) ||
 766                              !is_exactly_a<mul>(ex_to<power>(r).basis) ||
 767                              !is_exactly_a<power>(ex_to<power>(r).basis));
 768
 769                 if (are_ex_trivially_equal(c,_ex1)) {
 770                         if (is_ex_exactly_of_type(r,mul)) {
 771                                 sum.push_back(expair(expand_mul(ex_to<mul>(r),_num2),
 772                                                      _ex1));
 773                         } else {
 774                                 sum.push_back(expair((new power(r,_ex2))->setflag(status_flags::dynallocated),
 775                                                      _ex1));
 776                         }
 777                 } else {
 778                         if (is_ex_exactly_of_type(r,mul)) {
 779                                 sum.push_back(a.combine_ex_with_coeff_to_pair(expand_mul(ex_to<mul>(r),_num2),
 780                                                      ex_to<numeric>(c).power_dyn(_num2)));
 781                         } else {
 782                                 sum.push_back(a.combine_ex_with_coeff_to_pair((new power(r,_ex2))->setflag(status_flags::dynallocated),
 783                                                      ex_to<numeric>(c).power_dyn(_num2)));
 784                         }
 785                 }
 786
 787                 for (epvector::const_iterator cit1=cit0+1; cit1!=last; ++cit1) {
 788                         const ex & r1 = cit1->rest;
 789                         const ex & c1 = cit1->coeff;
 790                         sum.push_back(a.combine_ex_with_coeff_to_pair((new mul(r,r1))->setflag(status_flags::dynallocated),
 791                                                                       _num2.mul(ex_to<numeric>(c)).mul_dyn(ex_to<numeric>(c1))));
 792                 }
 793         }
 794
 795         GINAC_ASSERT(sum.size()==(a.seq.size()*(a.seq.size()+1))/2);
 796
 797         // second part: add terms coming from overall_factor (if != 0)
 798         if (!a.overall_coeff.is_zero()) {
 799                 epvector::const_iterator i = a.seq.begin(), end = a.seq.end();
 800                 while (i != end) {
 801                         sum.push_back(a.combine_pair_with_coeff_to_pair(*i, ex_to<numeric>(a.overall_coeff).mul_dyn(_num2)));
 802                         ++i;
 803                 }
 804                 sum.push_back(expair(ex_to<numeric>(a.overall_coeff).power_dyn(_num2),_ex1));
 805         }
 806
 807         GINAC_ASSERT(sum.size()==(a_nops*(a_nops+1))/2);
 808
 809         return (new add(sum))->setflag(status_flags::dynallocated | status_flags::expanded);
 810 }
 811
 812 /** Expand factors of m in m^n where m is a mul and n is and integer.
 813  *  @see power::expand */
 814 ex power::expand_mul(const mul & m, const numeric & n) const
 815 {
 816         GINAC_ASSERT(n.is_integer());
 817
 818         if (n.is_zero())
 819                 return _ex1;
 820
 821         epvector distrseq;
 822         distrseq.reserve(m.seq.size());
 823         epvector::const_iterator last = m.seq.end();
 824         epvector::const_iterator cit = m.seq.begin();
 825         while (cit!=last) {
 826                 if (is_ex_exactly_of_type(cit->rest,numeric)) {
 827                         distrseq.push_back(m.combine_pair_with_coeff_to_pair(*cit, n));
 828                 } else {
 829                         // it is safe not to call mul::combine_pair_with_coeff_to_pair()
 830                         // since n is an integer
 831                         distrseq.push_back(expair(cit->rest, ex_to<numeric>(cit->coeff).mul(n)));
 832                 }
 833                 ++cit;
 834         }
 835         return (new mul(distrseq, ex_to<numeric>(m.overall_coeff).power_dyn(n)))->setflag(status_flags::dynallocated);
 836 }
 837
 838 } // namespace GiNaC