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