- '?time' now works
[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 <stdexcept>
25
26 #include "pseries.h"
27 #include "add.h"
28 #include "inifcns.h"
29 #include "lst.h"
30 #include "mul.h"
31 #include "power.h"
32 #include "relational.h"
33 #include "symbol.h"
34 #include "archive.h"
35 #include "utils.h"
36 #include "debugmsg.h"
37
38 #ifndef NO_NAMESPACE_GINAC
39 namespace GiNaC {
40 #endif // ndef NO_NAMESPACE_GINAC
41
42 GINAC_IMPLEMENT_REGISTERED_CLASS(pseries, basic)
43
44 /*
45  *  Default constructor, destructor, copy constructor, assignment operator and helpers
46  */
47
48 pseries::pseries() : basic(TINFO_pseries)
49 {
50     debugmsg("pseries default constructor", LOGLEVEL_CONSTRUCT);
51 }
52
53 pseries::~pseries()
54 {
55     debugmsg("pseries destructor", LOGLEVEL_DESTRUCT);
56     destroy(false);
57 }
58
59 pseries::pseries(const pseries &other)
60 {
61     debugmsg("pseries copy constructor", LOGLEVEL_CONSTRUCT);
62     copy(other);
63 }
64
65 const pseries &pseries::operator=(const pseries & other)
66 {
67     debugmsg("pseries operator=", LOGLEVEL_ASSIGNMENT);
68     if (this != &other) {
69         destroy(true);
70         copy(other);
71     }
72     return *this;
73 }
74
75 void pseries::copy(const pseries &other)
76 {
77     inherited::copy(other);
78     seq = other.seq;
79     var = other.var;
80     point = other.point;
81 }
82
83 void pseries::destroy(bool call_parent)
84 {
85     if (call_parent)
86         inherited::destroy(call_parent);
87 }
88
89
90 /*
91  *  Other constructors
92  */
93
94 /** Construct pseries from a vector of coefficients and powers.
95  *  expair.rest holds the coefficient, expair.coeff holds the power.
96  *  The powers must be integers (positive or negative) and in ascending order;
97  *  the last coefficient can be Order(_ex1()) to represent a truncated,
98  *  non-terminating series.
99  *
100  *  @param rel__  expansion variable and point (must hold a relational)
101  *  @param ops_  vector of {coefficient, power} pairs (coefficient must not be zero)
102  *  @return newly constructed pseries */
103 pseries::pseries(const ex &rel_, const epvector &ops_)
104     : basic(TINFO_pseries), seq(ops_)
105 {
106     debugmsg("pseries constructor from rel,epvector", LOGLEVEL_CONSTRUCT);
107     GINAC_ASSERT(is_ex_exactly_of_type(rel_, relational));
108     GINAC_ASSERT(is_ex_exactly_of_type(rel_.lhs(),symbol));
109     point = rel_.rhs();
110     var = *static_cast<symbol *>(rel_.lhs().bp);
111 }
112
113
114 /*
115  *  Archiving
116  */
117
118 /** Construct object from archive_node. */
119 pseries::pseries(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
120 {
121     debugmsg("pseries constructor from archive_node", LOGLEVEL_CONSTRUCT);
122     for (unsigned int i=0; true; i++) {
123         ex rest;
124         ex coeff;
125         if (n.find_ex("coeff", rest, sym_lst, i) && n.find_ex("power", coeff, sym_lst, i))
126             seq.push_back(expair(rest, coeff));
127         else
128             break;
129     }
130     n.find_ex("var", var, sym_lst);
131     n.find_ex("point", point, sym_lst);
132 }
133
134 /** Unarchive the object. */
135 ex pseries::unarchive(const archive_node &n, const lst &sym_lst)
136 {
137     return (new pseries(n, sym_lst))->setflag(status_flags::dynallocated);
138 }
139
140 /** Archive the object. */
141 void pseries::archive(archive_node &n) const
142 {
143     inherited::archive(n);
144     epvector::const_iterator i = seq.begin(), iend = seq.end();
145     while (i != iend) {
146         n.add_ex("coeff", i->rest);
147         n.add_ex("power", i->coeff);
148         i++;
149     }
150     n.add_ex("var", var);
151     n.add_ex("point", point);
152 }
153
154
155 /*
156  *  Functions overriding virtual functions from base classes
157  */
158
159 basic *pseries::duplicate() const
160 {
161     debugmsg("pseries duplicate", LOGLEVEL_DUPLICATE);
162     return new pseries(*this);
163 }
164
165 void pseries::print(ostream &os, unsigned upper_precedence) const
166 {
167     debugmsg("pseries print", LOGLEVEL_PRINT);
168     convert_to_poly().print(os, upper_precedence);
169 }
170
171 void pseries::printraw(ostream &os) const
172 {
173         debugmsg("pseries printraw", LOGLEVEL_PRINT);
174         os << "pseries(" << var << ";" << point << ";";
175         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) {
176                 os << "(" << (*i).rest << "," << (*i).coeff << "),";
177         }
178         os << ")";
179 }
180
181 unsigned pseries::nops(void) const
182 {
183     return seq.size();
184 }
185
186 ex pseries::op(int i) const
187 {
188     if (i < 0 || unsigned(i) >= seq.size())
189         throw (std::out_of_range("op() out of range"));
190     return seq[i].rest * power(var - point, seq[i].coeff);
191 }
192
193 ex &pseries::let_op(int i)
194 {
195     throw (std::logic_error("let_op not defined for pseries"));
196 }
197
198 int pseries::degree(const symbol &s) const
199 {
200     if (var.is_equal(s)) {
201         // Return last exponent
202         if (seq.size())
203             return ex_to_numeric((*(seq.end() - 1)).coeff).to_int();
204         else
205             return 0;
206     } else {
207         epvector::const_iterator it = seq.begin(), itend = seq.end();
208         if (it == itend)
209             return 0;
210         int max_pow = INT_MIN;
211         while (it != itend) {
212             int pow = it->rest.degree(s);
213             if (pow > max_pow)
214                 max_pow = pow;
215             it++;
216         }
217         return max_pow;
218     }
219 }
220
221 int pseries::ldegree(const symbol &s) const
222 {
223     if (var.is_equal(s)) {
224         // Return first exponent
225         if (seq.size())
226             return ex_to_numeric((*(seq.begin())).coeff).to_int();
227         else
228             return 0;
229     } else {
230         epvector::const_iterator it = seq.begin(), itend = seq.end();
231         if (it == itend)
232             return 0;
233         int min_pow = INT_MAX;
234         while (it != itend) {
235             int pow = it->rest.ldegree(s);
236             if (pow < min_pow)
237                 min_pow = pow;
238             it++;
239         }
240         return min_pow;
241     }
242 }
243
244 ex pseries::coeff(const symbol &s, int n) const
245 {
246     if (var.is_equal(s)) {
247                 if (seq.size() == 0)
248                         return _ex0();
249
250                 // Binary search in sequence for given power
251                 numeric looking_for = numeric(n);
252                 int lo = 0, hi = seq.size() - 1;
253                 while (lo <= hi) {
254                         int mid = (lo + hi) / 2;
255                         GINAC_ASSERT(is_ex_exactly_of_type(seq[mid].coeff, numeric));
256                         int cmp = ex_to_numeric(seq[mid].coeff).compare(looking_for);
257                         switch (cmp) {
258                                 case -1:
259                                         lo = mid + 1;
260                                         break;
261                                 case 0:
262                                         return seq[mid].rest;
263                                 case 1:
264                                         hi = mid - 1;
265                                         break;
266                                 default:
267                                         throw(std::logic_error("pseries::coeff: compare() didn't return -1, 0 or 1"));
268                         }
269                 }
270                 return _ex0();
271     } else
272         return convert_to_poly().coeff(s, n);
273 }
274
275 ex pseries::collect(const symbol &s) const
276 {
277         return *this;
278 }
279
280 /** Evaluate coefficients. */
281 ex pseries::eval(int level) const
282 {
283     if (level == 1)
284         return this->hold();
285
286         if (level == -max_recursion_level)
287                 throw (std::runtime_error("pseries::eval(): recursion limit exceeded"));
288     
289     // Construct a new series with evaluated coefficients
290     epvector new_seq;
291     new_seq.reserve(seq.size());
292     epvector::const_iterator it = seq.begin(), itend = seq.end();
293     while (it != itend) {
294         new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
295         it++;
296     }
297     return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
298 }
299
300 /** Evaluate coefficients numerically. */
301 ex pseries::evalf(int level) const
302 {
303         if (level == 1)
304                 return *this;
305
306         if (level == -max_recursion_level)
307                 throw (std::runtime_error("pseries::evalf(): recursion limit exceeded"));
308
309     // Construct a new series with evaluated coefficients
310     epvector new_seq;
311     new_seq.reserve(seq.size());
312     epvector::const_iterator it = seq.begin(), itend = seq.end();
313     while (it != itend) {
314         new_seq.push_back(expair(it->rest.evalf(level-1), it->coeff));
315         it++;
316     }
317     return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
318 }
319
320 ex pseries::subs(const lst & ls, const lst & lr) const
321 {
322         // If expansion variable is being substituted, convert the series to a
323         // polynomial and do the substitution there because the result might
324         // no longer be a power series
325         if (ls.has(var))
326                 return convert_to_poly(true).subs(ls, lr);
327
328         // Otherwise construct a new series with substituted coefficients and
329         // expansion point
330         epvector new_seq;
331         new_seq.reserve(seq.size());
332         epvector::const_iterator it = seq.begin(), itend = seq.end();
333         while (it != itend) {
334                 new_seq.push_back(expair(it->rest.subs(ls, lr), it->coeff));
335                 it++;
336         }
337     return (new pseries(relational(var,point.subs(ls, lr)), new_seq))->setflag(status_flags::dynallocated);
338 }
339
340 /** Implementation of ex::diff() for a power series.  It treats the series as a
341  *  polynomial.
342  *  @see ex::diff */
343 ex pseries::derivative(const symbol & s) const
344 {
345     if (s == var) {
346         epvector new_seq;
347         epvector::const_iterator it = seq.begin(), itend = seq.end();
348         
349         // FIXME: coeff might depend on var
350         while (it != itend) {
351             if (is_order_function(it->rest)) {
352                 new_seq.push_back(expair(it->rest, it->coeff - 1));
353             } else {
354                 ex c = it->rest * it->coeff;
355                 if (!c.is_zero())
356                     new_seq.push_back(expair(c, it->coeff - 1));
357             }
358             it++;
359         }
360         return pseries(relational(var,point), new_seq);
361     } else {
362         return *this;
363     }
364 }
365
366
367 /*
368  *  Construct ordinary polynomial out of series
369  */
370
371 /** Convert a pseries object to an ordinary polynomial.
372  *
373  *  @param no_order flag: discard higher order terms */
374 ex pseries::convert_to_poly(bool no_order) const
375 {
376     ex e;
377     epvector::const_iterator it = seq.begin(), itend = seq.end();
378     
379     while (it != itend) {
380         if (is_order_function(it->rest)) {
381             if (!no_order)
382                 e += Order(power(var - point, it->coeff));
383         } else
384             e += it->rest * power(var - point, it->coeff);
385         it++;
386     }
387     return e;
388 }
389
390
391 /*
392  *  Implementation of series expansion
393  */
394
395 /** Default implementation of ex::series(). This performs Taylor expansion.
396  *  @see ex::series */
397 ex basic::series(const relational & r, int order) const
398 {
399     epvector seq;
400     numeric fac(1);
401     ex deriv = *this;
402     ex coeff = deriv.subs(r);
403     const symbol *s = static_cast<symbol *>(r.lhs().bp);
404     
405     if (!coeff.is_zero())
406         seq.push_back(expair(coeff, numeric(0)));
407     
408     int n;
409     for (n=1; n<order; n++) {
410         fac = fac.mul(numeric(n));
411         deriv = deriv.diff(*s).expand();
412         if (deriv.is_zero()) {
413             // Series terminates
414             return pseries(r, seq);
415         }
416         coeff = fac.inverse() * deriv.subs(r);
417         if (!coeff.is_zero())
418             seq.push_back(expair(coeff, numeric(n)));
419     }
420     
421     // Higher-order terms, if present
422     deriv = deriv.diff(*s);
423     if (!deriv.is_zero())
424         seq.push_back(expair(Order(_ex1()), numeric(n)));
425     return pseries(r, seq);
426 }
427
428
429 /** Implementation of ex::series() for symbols.
430  *  @see ex::series */
431 ex symbol::series(const relational & r, int order) const
432 {
433         epvector seq;
434     const ex point = r.rhs();
435     GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
436     const symbol *s = static_cast<symbol *>(r.lhs().bp);
437     
438         if (this->is_equal(*s)) {
439                 if (order > 0 && !point.is_zero())
440                         seq.push_back(expair(point, _ex0()));
441                 if (order > 1)
442                         seq.push_back(expair(_ex1(), _ex1()));
443                 else
444                         seq.push_back(expair(Order(_ex1()), numeric(order)));
445         } else
446                 seq.push_back(expair(*this, _ex0()));
447         return pseries(r, seq);
448 }
449
450
451 /** Add one series object to another, producing a pseries object that
452  *  represents the sum.
453  *
454  *  @param other  pseries object to add with
455  *  @return the sum as a pseries */
456 ex pseries::add_series(const pseries &other) const
457 {
458     // Adding two series with different variables or expansion points
459     // results in an empty (constant) series 
460     if (!is_compatible_to(other)) {
461         epvector nul;
462         nul.push_back(expair(Order(_ex1()), _ex0()));
463         return pseries(relational(var,point), nul);
464     }
465     
466     // Series addition
467     epvector new_seq;
468     epvector::const_iterator a = seq.begin();
469     epvector::const_iterator b = other.seq.begin();
470     epvector::const_iterator a_end = seq.end();
471     epvector::const_iterator b_end = other.seq.end();
472     int pow_a = INT_MAX, pow_b = INT_MAX;
473     for (;;) {
474         // If a is empty, fill up with elements from b and stop
475         if (a == a_end) {
476             while (b != b_end) {
477                 new_seq.push_back(*b);
478                 b++;
479             }
480             break;
481         } else
482             pow_a = ex_to_numeric((*a).coeff).to_int();
483         
484         // If b is empty, fill up with elements from a and stop
485         if (b == b_end) {
486             while (a != a_end) {
487                 new_seq.push_back(*a);
488                 a++;
489             }
490             break;
491         } else
492             pow_b = ex_to_numeric((*b).coeff).to_int();
493         
494         // a and b are non-empty, compare powers
495         if (pow_a < pow_b) {
496             // a has lesser power, get coefficient from a
497             new_seq.push_back(*a);
498             if (is_order_function((*a).rest))
499                 break;
500             a++;
501         } else if (pow_b < pow_a) {
502             // b has lesser power, get coefficient from b
503             new_seq.push_back(*b);
504             if (is_order_function((*b).rest))
505                 break;
506             b++;
507         } else {
508             // Add coefficient of a and b
509             if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
510                 new_seq.push_back(expair(Order(_ex1()), (*a).coeff));
511                 break;  // Order term ends the sequence
512             } else {
513                 ex sum = (*a).rest + (*b).rest;
514                 if (!(sum.is_zero()))
515                     new_seq.push_back(expair(sum, numeric(pow_a)));
516                 a++;
517                 b++;
518             }
519         }
520     }
521     return pseries(relational(var,point), new_seq);
522 }
523
524
525 /** Implementation of ex::series() for sums. This performs series addition when
526  *  adding pseries objects.
527  *  @see ex::series */
528 ex add::series(const relational & r, int order) const
529 {
530     ex acc; // Series accumulator
531     
532     // Get first term from overall_coeff
533     acc = overall_coeff.series(r, order);
534
535     // Add remaining terms
536     epvector::const_iterator it = seq.begin();
537     epvector::const_iterator itend = seq.end();
538     for (; it!=itend; it++) {
539         ex op;
540         if (is_ex_exactly_of_type(it->rest, pseries))
541             op = it->rest;
542         else
543             op = it->rest.series(r, order);
544         if (!it->coeff.is_equal(_ex1()))
545             op = ex_to_pseries(op).mul_const(ex_to_numeric(it->coeff));
546         
547         // Series addition
548         acc = ex_to_pseries(acc).add_series(ex_to_pseries(op));
549     }
550     return acc;
551 }
552
553
554 /** Multiply a pseries object with a numeric constant, producing a pseries
555  *  object that represents the product.
556  *
557  *  @param other  constant to multiply with
558  *  @return the product as a pseries */
559 ex pseries::mul_const(const numeric &other) const
560 {
561     epvector new_seq;
562     new_seq.reserve(seq.size());
563     
564     epvector::const_iterator it = seq.begin(), itend = seq.end();
565     while (it != itend) {
566         if (!is_order_function(it->rest))
567             new_seq.push_back(expair(it->rest * other, it->coeff));
568         else
569             new_seq.push_back(*it);
570         it++;
571     }
572     return pseries(relational(var,point), new_seq);
573 }
574
575
576 /** Multiply one pseries object to another, producing a pseries object that
577  *  represents the product.
578  *
579  *  @param other  pseries object to multiply with
580  *  @return the product as a pseries */
581 ex pseries::mul_series(const pseries &other) const
582 {
583     // Multiplying two series with different variables or expansion points
584     // results in an empty (constant) series 
585     if (!is_compatible_to(other)) {
586         epvector nul;
587         nul.push_back(expair(Order(_ex1()), _ex0()));
588         return pseries(relational(var,point), nul);
589     }
590
591     // Series multiplication
592     epvector new_seq;
593     
594     const symbol *s = static_cast<symbol *>(var.bp);
595     int a_max = degree(*s);
596     int b_max = other.degree(*s);
597     int a_min = ldegree(*s);
598     int b_min = other.ldegree(*s);
599     int cdeg_min = a_min + b_min;
600     int cdeg_max = a_max + b_max;
601     
602     int higher_order_a = INT_MAX;
603     int higher_order_b = INT_MAX;
604     if (is_order_function(coeff(*s, a_max)))
605         higher_order_a = a_max + b_min;
606     if (is_order_function(other.coeff(*s, b_max)))
607         higher_order_b = b_max + a_min;
608     int higher_order_c = min(higher_order_a, higher_order_b);
609     if (cdeg_max >= higher_order_c)
610         cdeg_max = higher_order_c - 1;
611     
612     for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) {
613         ex co = _ex0();
614         // c(i)=a(0)b(i)+...+a(i)b(0)
615         for (int i=a_min; cdeg-i>=b_min; i++) {
616             ex a_coeff = coeff(*s, i);
617             ex b_coeff = other.coeff(*s, cdeg-i);
618             if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
619                 co += coeff(*s, i) * other.coeff(*s, cdeg-i);
620         }
621         if (!co.is_zero())
622             new_seq.push_back(expair(co, numeric(cdeg)));
623     }
624     if (higher_order_c < INT_MAX)
625         new_seq.push_back(expair(Order(_ex1()), numeric(higher_order_c)));
626     return pseries(relational(var,point), new_seq);
627 }
628
629
630 /** Implementation of ex::series() for product. This performs series
631  *  multiplication when multiplying series.
632  *  @see ex::series */
633 ex mul::series(const relational & r, int order) const
634 {
635     ex acc; // Series accumulator
636     
637     // Get first term from overall_coeff
638     acc = overall_coeff.series(r, order);
639     
640     // Multiply with remaining terms
641     epvector::const_iterator it = seq.begin();
642     epvector::const_iterator itend = seq.end();
643     for (; it!=itend; it++) {
644         ex op = it->rest;
645         if (op.info(info_flags::numeric)) {
646             // series * const (special case, faster)
647             ex f = power(op, it->coeff);
648             acc = ex_to_pseries(acc).mul_const(ex_to_numeric(f));
649             continue;
650         } else if (!is_ex_exactly_of_type(op, pseries))
651             op = op.series(r, order);
652         if (!it->coeff.is_equal(_ex1()))
653             op = ex_to_pseries(op).power_const(ex_to_numeric(it->coeff), order);
654
655         // Series multiplication
656         acc = ex_to_pseries(acc).mul_series(ex_to_pseries(op));
657     }
658     return acc;
659 }
660
661
662 /** Compute the p-th power of a series.
663  *
664  *  @param p  power to compute
665  *  @param deg  truncation order of series calculation */
666 ex pseries::power_const(const numeric &p, int deg) const
667 {
668     int i;
669     const symbol *s = static_cast<symbol *>(var.bp);
670     int ldeg = ldegree(*s);
671     
672     // Calculate coefficients of powered series
673     exvector co;
674     co.reserve(deg);
675     ex co0;
676     co.push_back(co0 = power(coeff(*s, ldeg), p));
677     bool all_sums_zero = true;
678     for (i=1; i<deg; i++) {
679         ex sum = _ex0();
680         for (int j=1; j<=i; j++) {
681             ex c = coeff(*s, j + ldeg);
682             if (is_order_function(c)) {
683                 co.push_back(Order(_ex1()));
684                 break;
685             } else
686                 sum += (p * j - (i - j)) * co[i - j] * c;
687         }
688         if (!sum.is_zero())
689             all_sums_zero = false;
690         co.push_back(co0 * sum / numeric(i));
691     }
692     
693     // Construct new series (of non-zero coefficients)
694     epvector new_seq;
695     bool higher_order = false;
696     for (i=0; i<deg; i++) {
697         if (!co[i].is_zero())
698             new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
699         if (is_order_function(co[i])) {
700             higher_order = true;
701             break;
702         }
703     }
704     if (!higher_order && !all_sums_zero)
705         new_seq.push_back(expair(Order(_ex1()), numeric(deg) + p * ldeg));
706     return pseries(relational(var,point), new_seq);
707 }
708
709
710 /** Implementation of ex::series() for powers. This performs Laurent expansion
711  *  of reciprocals of series at singularities.
712  *  @see ex::series */
713 ex power::series(const relational & r, int order) const
714 {
715     ex e;
716     if (!is_ex_exactly_of_type(basis, pseries)) {
717         // Basis is not a series, may there be a singulary?
718         if (!exponent.info(info_flags::negint))
719             return basic::series(r, order);
720         
721         // Expression is of type something^(-int), check for singularity
722         if (!basis.subs(r).is_zero())
723             return basic::series(r, order);
724         
725         // Singularity encountered, expand basis into series
726         e = basis.series(r, order);
727     } else {
728         // Basis is a series
729         e = basis;
730     }
731     
732     // Power e
733     return ex_to_pseries(e).power_const(ex_to_numeric(exponent), order);
734 }
735
736
737 /** Re-expansion of a pseries object. */
738 ex pseries::series(const relational & r, int order) const
739 {
740     const ex p = r.rhs();
741     GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
742     const symbol *s = static_cast<symbol *>(r.lhs().bp);
743     
744         if (var.is_equal(*s) && point.is_equal(p)) {
745                 if (order > degree(*s))
746                         return *this;
747                 else {
748                 epvector new_seq;
749                 epvector::const_iterator it = seq.begin(), itend = seq.end();
750                         while (it != itend) {
751                                 int o = ex_to_numeric(it->coeff).to_int();
752                                 if (o >= order) {
753                                         new_seq.push_back(expair(Order(_ex1()), o));
754                                         break;
755                                 }
756                                 new_seq.push_back(*it);
757                                 it++;
758                         }
759                         return pseries(r, new_seq);
760                 }
761         } else
762                 return convert_to_poly().series(r, order);
763 }
764
765
766 /** Compute the truncated series expansion of an expression.
767  *  This function returns an expression containing an object of class pseries 
768  *  to represent the series. If the series does not terminate within the given
769  *  truncation order, the last term of the series will be an order term.
770  *
771  *  @param r  expansion relation, lhs holds variable and rhs holds point
772  *  @param order  truncation order of series calculations
773  *  @return an expression holding a pseries object */
774 ex ex::series(const ex & r, int order) const
775 {
776     GINAC_ASSERT(bp!=0);
777     ex e;
778     relational rel_;
779     
780     if (is_ex_exactly_of_type(r,relational))
781         rel_ = ex_to_relational(r);
782     else if (is_ex_exactly_of_type(r,symbol))
783         rel_ = relational(r,_ex0());
784     else
785         throw (std::logic_error("ex::series(): expansion point has unknown type"));
786     
787     try {
788         e = bp->series(rel_, order);
789     } catch (exception &x) {
790         throw (std::logic_error(string("unable to compute series (") + x.what() + ")"));
791     }
792     return e;
793 }
794
795
796 // Global constants
797 const pseries some_pseries;
798 const type_info & typeid_pseries = typeid(some_pseries);
799
800 #ifndef NO_NAMESPACE_GINAC
801 } // namespace GiNaC
802 #endif // ndef NO_NAMESPACE_GINAC