ginac/power.cpp

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