- enforced GiNaC coding standards :-)
[ginac.git] / ginac / series.cpp
1 /** @file series.cpp
2  *
3  *  Implementation of class for extended truncated power-series and
4  *  methods for series expansion.
5  *
6  *  GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
7  *
8  *  This program is free software; you can redistribute it and/or modify
9  *  it under the terms of the GNU General Public License as published by
10  *  the Free Software Foundation; either version 2 of the License, or
11  *  (at your option) any later version.
12  *
13  *  This program is distributed in the hope that it will be useful,
14  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16  *  GNU General Public License for more details.
17  *
18  *  You should have received a copy of the GNU General Public License
19  *  along with this program; if not, write to the Free Software
20  *  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
21  */
22
23 #include "ginac.h"
24
25
26 /*
27  *  Default constructor, destructor, copy constructor, assignment operator and helpers
28  */
29
30 series::series() : basic(TINFO_SERIES)
31 {
32     debugmsg("series default constructor", LOGLEVEL_CONSTRUCT);
33 }
34
35 series::~series()
36 {
37     debugmsg("series destructor", LOGLEVEL_DESTRUCT);
38     destroy(false);
39 }
40
41 series::series(series const &other)
42 {
43     debugmsg("series copy constructor", LOGLEVEL_CONSTRUCT);
44     copy(other);
45 }
46
47 series const &series::operator=(series const & other)
48 {
49     debugmsg("series operator=", LOGLEVEL_ASSIGNMENT);
50     if (this != &other) {
51         destroy(true);
52         copy(other);
53     }
54     return *this;
55 }
56
57 void series::copy(series const &other)
58 {
59     inherited::copy(other);
60     seq = other.seq;
61     var = other.var;
62     point = other.point;
63 }
64
65 void series::destroy(bool call_parent)
66 {
67     if (call_parent)
68         inherited::destroy(call_parent);
69 }
70
71
72 /*
73  *  Other constructors
74  */
75
76 /** Construct series from a vector of coefficients and powers.
77  *  expair.rest holds the coefficient, expair.coeff holds the power.
78  *  The powers must be integers (positive or negative) and in ascending order;
79  *  the last coefficient can be Order(exONE()) to represent a truncated,
80  *  non-terminating series.
81  *
82  *  @param var_  series variable (must hold a symbol)
83  *  @param point_  expansion point
84  *  @param ops_  vector of {coefficient, power} pairs (coefficient must not be zero)
85  *  @return newly constructed series */
86 series::series(ex const &var_, ex const &point_, epvector const &ops_)
87     : basic(TINFO_SERIES), seq(ops_), var(var_), point(point_)
88 {
89     debugmsg("series constructor from ex,ex,epvector", LOGLEVEL_CONSTRUCT);
90     ASSERT(is_ex_exactly_of_type(var_, symbol));
91 }
92
93
94 /*
95  *  Functions overriding virtual functions from base classes
96  */
97
98 basic *series::duplicate() const
99 {
100     debugmsg("series duplicate", LOGLEVEL_DUPLICATE);
101     return new series(*this);
102 }
103
104 // Highest degree of variable
105 int series::degree(symbol const &s) const
106 {
107     if (var == s) {
108         // Return last exponent
109         if (seq.size())
110             return ex_to_numeric((*(seq.end() - 1)).coeff).to_int();
111         else
112             return 0;
113     } else {
114         epvector::const_iterator it = seq.begin(), itend = seq.end();
115         if (it == itend)
116             return 0;
117         int max_pow = INT_MIN;
118         while (it != itend) {
119             int pow = it->rest.degree(s);
120             if (pow > max_pow)
121                 max_pow = pow;
122             it++;
123         }
124         return max_pow;
125     }
126 }
127
128 // Lowest degree of variable
129 int series::ldegree(symbol const &s) const
130 {
131     if (var == s) {
132         // Return first exponent
133         if (seq.size())
134             return ex_to_numeric((*(seq.begin())).coeff).to_int();
135         else
136             return 0;
137     } else {
138         epvector::const_iterator it = seq.begin(), itend = seq.end();
139         if (it == itend)
140             return 0;
141         int min_pow = INT_MAX;
142         while (it != itend) {
143             int pow = it->rest.ldegree(s);
144             if (pow < min_pow)
145                 min_pow = pow;
146             it++;
147         }
148         return min_pow;
149     }
150 }
151
152 // Coefficient of variable
153 ex series::coeff(symbol const &s, int n) const
154 {
155     if (var == s) {
156         epvector::const_iterator it = seq.begin(), itend = seq.end();
157         while (it != itend) {
158             int pow = ex_to_numeric(it->coeff).to_int();
159             if (pow == n)
160                 return it->rest;
161             if (pow > n)
162                 return exZERO();
163             it++;
164         }
165         return exZERO();
166     } else
167         return convert_to_poly().coeff(s, n);
168 }
169
170 ex series::eval(int level) const
171 {
172     if (level == 1)
173         return this->hold();
174
175     // Construct a new series with evaluated coefficients
176     epvector new_seq;
177     new_seq.reserve(seq.size());
178     epvector::const_iterator it = seq.begin(), itend = seq.end();
179     while (it != itend) {
180         new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
181         it++;
182     }
183     return (new series(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
184 }
185
186 ex series::evalf(int level) const
187 {
188     return convert_to_poly().evalf(level);
189 }
190
191
192 /*
193  *  Construct expression (polynomial) out of series
194  */
195
196 /** Convert a series object to an ordinary polynomial.
197  *
198  *  @param no_order flag: discard higher order terms */
199 ex series::convert_to_poly(bool no_order) const
200 {
201     ex e;
202     epvector::const_iterator it = seq.begin(), itend = seq.end();
203
204     while (it != itend) {
205         if (is_order_function(it->rest)) {
206             if (!no_order)
207                 e += Order(power(var - point, it->coeff));
208         } else
209             e += it->rest * power(var - point, it->coeff);
210         it++;
211     }
212     return e;
213 }
214
215
216 /*
217  *  Implementation of series expansion
218  */
219
220 /** Default implementation of ex::series(). This performs Taylor expansion.
221  *  @see ex::series */
222 ex basic::series(symbol const & s, ex const & point, int order) const
223 {
224     epvector seq;
225     numeric fac(1);
226     ex deriv = *this;
227     ex coeff = deriv.subs(s == point);
228     if (!coeff.is_zero())
229         seq.push_back(expair(coeff, numeric(0)));
230
231     int n;
232     for (n=1; n<order; n++) {
233         fac = fac.mul(numeric(n));
234         deriv = deriv.diff(s).expand();
235         if (deriv.is_zero()) {
236             // Series terminates
237             return series::series(s, point, seq);
238         }
239         coeff = power(fac, -1) * deriv.subs(s == point);
240         if (!coeff.is_zero())
241             seq.push_back(expair(coeff, numeric(n)));
242     }
243
244     // Higher-order terms, if present
245     deriv = deriv.diff(s);
246     if (!deriv.is_zero())
247         seq.push_back(expair(Order(exONE()), numeric(n)));
248     return series::series(s, point, seq);
249 }
250
251
252 /** Add one series object to another, producing a series object that represents
253  *  the sum.
254  *
255  *  @param other  series object to add with
256  *  @return the sum as a series */
257 ex series::add_series(const series &other) const
258 {
259     // Adding two series with different variables or expansion points
260     // results in an empty (constant) series 
261     if (!is_compatible_to(other)) {
262         epvector nul;
263         nul.push_back(expair(Order(exONE()), exZERO()));
264         return series(var, point, nul);
265     }
266
267     // Series addition
268     epvector new_seq;
269     epvector::const_iterator a = seq.begin();
270     epvector::const_iterator b = other.seq.begin();
271     epvector::const_iterator a_end = seq.end();
272     epvector::const_iterator b_end = other.seq.end();
273     int pow_a = INT_MAX, pow_b = INT_MAX;
274     for (;;) {
275         // If a is empty, fill up with elements from b and stop
276         if (a == a_end) {
277             while (b != b_end) {
278                 new_seq.push_back(*b);
279                 b++;
280             }
281             break;
282         } else
283             pow_a = ex_to_numeric((*a).coeff).to_int();
284
285         // If b is empty, fill up with elements from a and stop
286         if (b == b_end) {
287             while (a != a_end) {
288                 new_seq.push_back(*a);
289                 a++;
290             }
291             break;
292         } else
293             pow_b = ex_to_numeric((*b).coeff).to_int();
294
295         // a and b are non-empty, compare powers
296         if (pow_a < pow_b) {
297             // a has lesser power, get coefficient from a
298             new_seq.push_back(*a);
299             if (is_order_function((*a).rest))
300                 break;
301             a++;
302         } else if (pow_b < pow_a) {
303             // b has lesser power, get coefficient from b
304             new_seq.push_back(*b);
305             if (is_order_function((*b).rest))
306                 break;
307             b++;
308         } else {
309             // Add coefficient of a and b
310             if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
311                 new_seq.push_back(expair(Order(exONE()), (*a).coeff));
312                 break;  // Order term ends the sequence
313             } else {
314                 ex sum = (*a).rest + (*b).rest;
315                 if (!(sum == exZERO()))
316                     new_seq.push_back(expair(sum, numeric(pow_a)));
317                 a++;
318                 b++;
319             }
320         }
321     }
322     return series(var, point, new_seq);
323 }
324
325
326 /** Implementation of ex::series() for sums. This performs series addition when
327  *  adding series objects.
328  *  @see ex::series */
329 /*
330 ex add::series(symbol const & s, ex const & point, int order) const
331 {
332     ex acc; // Series accumulator
333
334     // Get first term
335     epvector::const_iterator it = seq.begin();
336     epvector::const_iterator itend = seq.end();
337     if (it != itend) {
338         if (is_ex_exactly_of_type(it->rest, series))
339             acc = it->rest;
340         else
341             acc = it->rest.series(s, point, order);
342         if (!it->coeff.is_equal(exONE()))
343             acc = ex_to_series(acc).mul_const(ex_to_numeric(it->coeff));
344         it++;
345     }
346
347     // Add remaining terms
348     for (; it!=itend; it++) {
349         ex op;
350         if (is_ex_exactly_of_type(it->rest, series))
351             op = it->rest;
352         else
353             op = it->rest.series(s, point, order);
354         if (!it->coeff.is_equal(exONE()))
355             op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
356
357         // Series addition
358         acc = ex_to_series(acc).add_series(ex_to_series(op));
359     }
360     return acc;
361 }
362 */
363 ex add::series(symbol const & s, ex const & point, int order) const
364 {
365     ex acc; // Series accumulator
366
367     // Get first term from overall_coeff
368     acc = overall_coeff.series(s,point,order);
369
370     // Add remaining terms
371     epvector::const_iterator it = seq.begin();
372     epvector::const_iterator itend = seq.end();
373     for (; it!=itend; it++) {
374         ex op;
375         if (is_ex_exactly_of_type(it->rest, series))
376             op = it->rest;
377         else
378             op = it->rest.series(s, point, order);
379         if (!it->coeff.is_equal(exONE()))
380             op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
381
382         // Series addition
383         acc = ex_to_series(acc).add_series(ex_to_series(op));
384     }
385     return acc;
386 }
387
388
389 /** Multiply a series object with a numeric constant, producing a series object
390  *  that represents the product.
391  *
392  *  @param other  constant to multiply with
393  *  @return the product as a series */
394 ex series::mul_const(const numeric &other) const
395 {
396     epvector new_seq;
397     new_seq.reserve(seq.size());
398
399     epvector::const_iterator it = seq.begin(), itend = seq.end();
400     while (it != itend) {
401         if (!is_order_function(it->rest))
402             new_seq.push_back(expair(it->rest * other, it->coeff));
403         else
404             new_seq.push_back(*it);
405         it++;
406     }
407     return series(var, point, new_seq);
408 }
409
410
411 /** Multiply one series object to another, producing a series object that
412  *  represents the product.
413  *
414  *  @param other  series object to multiply with
415  *  @return the product as a series */
416 ex series::mul_series(const series &other) const
417 {
418     // Multiplying two series with different variables or expansion points
419     // results in an empty (constant) series 
420     if (!is_compatible_to(other)) {
421         epvector nul;
422         nul.push_back(expair(Order(exONE()), exZERO()));
423         return series(var, point, nul);
424     }
425
426     // Series multiplication
427     epvector new_seq;
428
429     const symbol *s = static_cast<symbol *>(var.bp);
430     int a_max = degree(*s);
431     int b_max = other.degree(*s);
432     int a_min = ldegree(*s);
433     int b_min = other.ldegree(*s);
434     int cdeg_min = a_min + b_min;
435     int cdeg_max = a_max + b_max;
436
437     int higher_order_a = INT_MAX;
438     int higher_order_b = INT_MAX;
439     if (is_order_function(coeff(*s, a_max)))
440         higher_order_a = a_max + b_min;
441     if (is_order_function(other.coeff(*s, b_max)))
442         higher_order_b = b_max + a_min;
443     int higher_order_c = min(higher_order_a, higher_order_b);
444     if (cdeg_max >= higher_order_c)
445         cdeg_max = higher_order_c - 1;
446
447     for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) {
448         ex co = exZERO();
449         // c(i)=a(0)b(i)+...+a(i)b(0)
450         for (int i=a_min; cdeg-i>=b_min; i++) {
451             ex a_coeff = coeff(*s, i);
452             ex b_coeff = other.coeff(*s, cdeg-i);
453             if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
454                 co += coeff(*s, i) * other.coeff(*s, cdeg-i);
455         }
456         if (!co.is_zero())
457             new_seq.push_back(expair(co, numeric(cdeg)));
458     }
459     if (higher_order_c < INT_MAX)
460         new_seq.push_back(expair(Order(exONE()), numeric(higher_order_c)));
461     return series::series(var, point, new_seq);
462 }
463
464
465 /** Implementation of ex::series() for product. This performs series multiplication when multiplying series.
466  *  @see ex::series */
467 /*
468 ex mul::series(symbol const & s, ex const & point, int order) const
469 {
470     ex acc; // Series accumulator
471
472     // Get first term
473     epvector::const_iterator it = seq.begin();
474     epvector::const_iterator itend = seq.end();
475     if (it != itend) {
476         if (is_ex_exactly_of_type(it->rest, series))
477             acc = it->rest;
478         else
479             acc = it->rest.series(s, point, order);
480         if (!it->coeff.is_equal(exONE()))
481             acc = ex_to_series(acc).power_const(ex_to_numeric(it->coeff), order);
482         it++;
483     }
484
485     // Multiply with remaining terms
486     for (; it!=itend; it++) {
487         ex op = it->rest;
488         if (op.info(info_flags::numeric)) {
489             // series * const (special case, faster)
490             ex f = power(op, it->coeff);
491             acc = ex_to_series(acc).mul_const(ex_to_numeric(f));
492             continue;
493         } else if (!is_ex_exactly_of_type(op, series))
494             op = op.series(s, point, order);
495         if (!it->coeff.is_equal(exONE()))
496             op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
497
498         // Series multiplication
499         acc = ex_to_series(acc).mul_series(ex_to_series(op));
500     }
501     return acc;
502 }
503 */
504
505 ex mul::series(symbol const & s, ex const & point, int order) const
506 {
507     ex acc; // Series accumulator
508
509     // Get first term from overall_coeff
510     acc = overall_coeff.series(s, point, order);
511
512     // Multiply with remaining terms
513     epvector::const_iterator it = seq.begin();
514     epvector::const_iterator itend = seq.end();
515     for (; it!=itend; it++) {
516         ex op = it->rest;
517         if (op.info(info_flags::numeric)) {
518             // series * const (special case, faster)
519             ex f = power(op, it->coeff);
520             acc = ex_to_series(acc).mul_const(ex_to_numeric(f));
521             continue;
522         } else if (!is_ex_exactly_of_type(op, series))
523             op = op.series(s, point, order);
524         if (!it->coeff.is_equal(exONE()))
525             op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
526
527         // Series multiplication
528         acc = ex_to_series(acc).mul_series(ex_to_series(op));
529     }
530     return acc;
531 }
532
533
534 /** Compute the p-th power of a series.
535  *
536  *  @param p  power to compute
537  *  @param deg  truncation order of series calculation */
538 ex series::power_const(const numeric &p, int deg) const
539 {
540     int i;
541     const symbol *s = static_cast<symbol *>(var.bp);
542     int ldeg = ldegree(*s);
543
544     // Calculate coefficients of powered series
545     exvector co;
546     co.reserve(deg);
547     ex co0;
548     co.push_back(co0 = power(coeff(*s, ldeg), p));
549     bool all_sums_zero = true;
550     for (i=1; i<deg; i++) {
551         ex sum = exZERO();
552         for (int j=1; j<=i; j++) {
553             ex c = coeff(*s, j + ldeg);
554             if (is_order_function(c)) {
555                 co.push_back(Order(exONE()));
556                 break;
557             } else
558                 sum += (p * j - (i - j)) * co[i - j] * c;
559         }
560         if (!sum.is_zero())
561             all_sums_zero = false;
562         co.push_back(co0 * sum / numeric(i));
563     }
564
565     // Construct new series (of non-zero coefficients)
566     epvector new_seq;
567     bool higher_order = false;
568     for (i=0; i<deg; i++) {
569         if (!co[i].is_zero())
570             new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
571         if (is_order_function(co[i])) {
572             higher_order = true;
573             break;
574         }
575     }
576     if (!higher_order && !all_sums_zero)
577         new_seq.push_back(expair(Order(exONE()), numeric(deg) + p * ldeg));
578     return series::series(var, point, new_seq);
579 }
580
581
582 /** Implementation of ex::series() for powers. This performs Laurent expansion
583  *  of reciprocals of series at singularities.
584  *  @see ex::series */
585 ex power::series(symbol const & s, ex const & point, int order) const
586 {
587     ex e;
588     if (!is_ex_exactly_of_type(basis, series)) {
589         // Basis is not a series, may there be a singulary?
590         if (!exponent.info(info_flags::negint))
591             return basic::series(s, point, order);
592
593         // Expression is of type something^(-int), check for singularity
594         if (!basis.subs(s == point).is_zero())
595             return basic::series(s, point, order);
596
597         // Singularity encountered, expand basis into series
598         e = basis.series(s, point, order);
599     } else {
600         // Basis is a series
601         e = basis;
602     }
603
604     // Power e
605     return ex_to_series(e).power_const(ex_to_numeric(exponent), order);
606 }
607
608
609 /** Compute the truncated series expansion of an expression.
610  *  This function returns an expression containing an object of class series to
611  *  represent the series. If the series does not terminate within the given
612  *  truncation order, the last term of the series will be an order term.
613  *
614  *  @param s  expansion variable
615  *  @param point  expansion point
616  *  @param order  truncation order of series calculations
617  *  @return an expression holding a series object */
618 ex ex::series(symbol const &s, ex const &point, int order) const
619 {
620     ASSERT(bp!=0);
621     return bp->series(s, point, order);
622 }
623
624
625 // Global constants
626 const series some_series;
627 type_info const & typeid_series = typeid(some_series);