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