- switched to automake build environment
[ginac.git] / 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);