f95bbf6df4b0bc6fab2d8ef0c387f83added262c
[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-2011 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
22  */
23
24 #include "pseries.h"
25 #include "add.h"
26 #include "inifcns.h" // for Order function
27 #include "lst.h"
28 #include "mul.h"
29 #include "power.h"
30 #include "relational.h"
31 #include "operators.h"
32 #include "symbol.h"
33 #include "integral.h"
34 #include "archive.h"
35 #include "utils.h"
36
37 #include <limits>
38 #include <numeric>
39 #include <stdexcept>
40
41 namespace GiNaC {
42
43 GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(pseries, basic,
44   print_func<print_context>(&pseries::do_print).
45   print_func<print_latex>(&pseries::do_print_latex).
46   print_func<print_tree>(&pseries::do_print_tree).
47   print_func<print_python>(&pseries::do_print_python).
48   print_func<print_python_repr>(&pseries::do_print_python_repr))
49
50
51 /*
52  *  Default constructor
53  */
54
55 pseries::pseries() { }
56
57
58 /*
59  *  Other ctors
60  */
61
62 /** Construct pseries from a vector of coefficients and powers.
63  *  expair.rest holds the coefficient, expair.coeff holds the power.
64  *  The powers must be integers (positive or negative) and in ascending order;
65  *  the last coefficient can be Order(_ex1) to represent a truncated,
66  *  non-terminating series.
67  *
68  *  @param rel_  expansion variable and point (must hold a relational)
69  *  @param ops_  vector of {coefficient, power} pairs (coefficient must not be zero)
70  *  @return newly constructed pseries */
71 pseries::pseries(const ex &rel_, const epvector &ops_) : seq(ops_)
72 {
73         GINAC_ASSERT(is_a<relational>(rel_));
74         GINAC_ASSERT(is_a<symbol>(rel_.lhs()));
75         point = rel_.rhs();
76         var = rel_.lhs();
77 }
78
79
80 /*
81  *  Archiving
82  */
83
84 void pseries::read_archive(const archive_node &n, lst &sym_lst) 
85 {
86         inherited::read_archive(n, sym_lst);
87         archive_node::archive_node_cit first = n.find_first("coeff");
88         archive_node::archive_node_cit last = n.find_last("power");
89         ++last;
90         seq.reserve((last-first)/2);
91
92         for (archive_node::archive_node_cit loc = first; loc < last;) {
93                 ex rest;
94                 ex coeff;
95                 n.find_ex_by_loc(loc++, rest, sym_lst);
96                 n.find_ex_by_loc(loc++, coeff, sym_lst);
97                 seq.push_back(expair(rest, coeff));
98         }
99
100         n.find_ex("var", var, sym_lst);
101         n.find_ex("point", point, sym_lst);
102 }
103
104 void pseries::archive(archive_node &n) const
105 {
106         inherited::archive(n);
107         epvector::const_iterator i = seq.begin(), iend = seq.end();
108         while (i != iend) {
109                 n.add_ex("coeff", i->rest);
110                 n.add_ex("power", i->coeff);
111                 ++i;
112         }
113         n.add_ex("var", var);
114         n.add_ex("point", point);
115 }
116
117
118 //////////
119 // functions overriding virtual functions from base classes
120 //////////
121
122 void pseries::print_series(const print_context & c, const char *openbrace, const char *closebrace, const char *mul_sym, const char *pow_sym, unsigned level) const
123 {
124         if (precedence() <= level)
125                 c.s << '(';
126                 
127         // objects of type pseries must not have any zero entries, so the
128         // trivial (zero) pseries needs a special treatment here:
129         if (seq.empty())
130                 c.s << '0';
131
132         epvector::const_iterator i = seq.begin(), end = seq.end();
133         while (i != end) {
134
135                 // print a sign, if needed
136                 if (i != seq.begin())
137                         c.s << '+';
138
139                 if (!is_order_function(i->rest)) {
140
141                         // print 'rest', i.e. the expansion coefficient
142                         if (i->rest.info(info_flags::numeric) &&
143                                 i->rest.info(info_flags::positive)) {
144                                 i->rest.print(c);
145                         } else {
146                                 c.s << openbrace << '(';
147                                 i->rest.print(c);
148                                 c.s << ')' << closebrace;
149                         }
150
151                         // print 'coeff', something like (x-1)^42
152                         if (!i->coeff.is_zero()) {
153                                 c.s << mul_sym;
154                                 if (!point.is_zero()) {
155                                         c.s << openbrace << '(';
156                                         (var-point).print(c);
157                                         c.s << ')' << closebrace;
158                                 } else
159                                         var.print(c);
160                                 if (i->coeff.compare(_ex1)) {
161                                         c.s << pow_sym;
162                                         c.s << openbrace;
163                                         if (i->coeff.info(info_flags::negative)) {
164                                                 c.s << '(';
165                                                 i->coeff.print(c);
166                                                 c.s << ')';
167                                         } else
168                                                 i->coeff.print(c);
169                                         c.s << closebrace;
170                                 }
171                         }
172                 } else
173                         Order(power(var-point,i->coeff)).print(c);
174                 ++i;
175         }
176
177         if (precedence() <= level)
178                 c.s << ')';
179 }
180
181 void pseries::do_print(const print_context & c, unsigned level) const
182 {
183         print_series(c, "", "", "*", "^", level);
184 }
185
186 void pseries::do_print_latex(const print_latex & c, unsigned level) const
187 {
188         print_series(c, "{", "}", " ", "^", level);
189 }
190
191 void pseries::do_print_python(const print_python & c, unsigned level) const
192 {
193         print_series(c, "", "", "*", "**", level);
194 }
195
196 void pseries::do_print_tree(const print_tree & c, unsigned level) const
197 {
198         c.s << std::string(level, ' ') << class_name() << " @" << this
199             << std::hex << ", hash=0x" << hashvalue << ", flags=0x" << flags << std::dec
200             << std::endl;
201         size_t num = seq.size();
202         for (size_t i=0; i<num; ++i) {
203                 seq[i].rest.print(c, level + c.delta_indent);
204                 seq[i].coeff.print(c, level + c.delta_indent);
205                 c.s << std::string(level + c.delta_indent, ' ') << "-----" << std::endl;
206         }
207         var.print(c, level + c.delta_indent);
208         point.print(c, level + c.delta_indent);
209 }
210
211 void pseries::do_print_python_repr(const print_python_repr & c, unsigned level) const
212 {
213         c.s << class_name() << "(relational(";
214         var.print(c);
215         c.s << ',';
216         point.print(c);
217         c.s << "),[";
218         size_t num = seq.size();
219         for (size_t i=0; i<num; ++i) {
220                 if (i)
221                         c.s << ',';
222                 c.s << '(';
223                 seq[i].rest.print(c);
224                 c.s << ',';
225                 seq[i].coeff.print(c);
226                 c.s << ')';
227         }
228         c.s << "])";
229 }
230
231 int pseries::compare_same_type(const basic & other) const
232 {
233         GINAC_ASSERT(is_a<pseries>(other));
234         const pseries &o = static_cast<const pseries &>(other);
235         
236         // first compare the lengths of the series...
237         if (seq.size()>o.seq.size())
238                 return 1;
239         if (seq.size()<o.seq.size())
240                 return -1;
241         
242         // ...then the expansion point...
243         int cmpval = var.compare(o.var);
244         if (cmpval)
245                 return cmpval;
246         cmpval = point.compare(o.point);
247         if (cmpval)
248                 return cmpval;
249         
250         // ...and if that failed the individual elements
251         epvector::const_iterator it = seq.begin(), o_it = o.seq.begin();
252         while (it!=seq.end() && o_it!=o.seq.end()) {
253                 cmpval = it->compare(*o_it);
254                 if (cmpval)
255                         return cmpval;
256                 ++it;
257                 ++o_it;
258         }
259
260         // so they are equal.
261         return 0;
262 }
263
264 /** Return the number of operands including a possible order term. */
265 size_t pseries::nops() const
266 {
267         return seq.size();
268 }
269
270 /** Return the ith term in the series when represented as a sum. */
271 ex pseries::op(size_t i) const
272 {
273         if (i >= seq.size())
274                 throw (std::out_of_range("op() out of range"));
275
276         if (is_order_function(seq[i].rest))
277                 return Order(power(var-point, seq[i].coeff));
278         return seq[i].rest * power(var - point, seq[i].coeff);
279 }
280
281 /** Return degree of highest power of the series.  This is usually the exponent
282  *  of the Order term.  If s is not the expansion variable of the series, the
283  *  series is examined termwise. */
284 int pseries::degree(const ex &s) const
285 {
286         if (var.is_equal(s)) {
287                 // Return last exponent
288                 if (seq.size())
289                         return ex_to<numeric>((seq.end()-1)->coeff).to_int();
290                 else
291                         return 0;
292         } else {
293                 epvector::const_iterator it = seq.begin(), itend = seq.end();
294                 if (it == itend)
295                         return 0;
296                 int max_pow = std::numeric_limits<int>::min();
297                 while (it != itend) {
298                         int pow = it->rest.degree(s);
299                         if (pow > max_pow)
300                                 max_pow = pow;
301                         ++it;
302                 }
303                 return max_pow;
304         }
305 }
306
307 /** Return degree of lowest power of the series.  This is usually the exponent
308  *  of the leading term.  If s is not the expansion variable of the series, the
309  *  series is examined termwise.  If s is the expansion variable but the
310  *  expansion point is not zero the series is not expanded to find the degree.
311  *  I.e.: (1-x) + (1-x)^2 + Order((1-x)^3) has ldegree(x) 1, not 0. */
312 int pseries::ldegree(const ex &s) const
313 {
314         if (var.is_equal(s)) {
315                 // Return first exponent
316                 if (seq.size())
317                         return ex_to<numeric>((seq.begin())->coeff).to_int();
318                 else
319                         return 0;
320         } else {
321                 epvector::const_iterator it = seq.begin(), itend = seq.end();
322                 if (it == itend)
323                         return 0;
324                 int min_pow = std::numeric_limits<int>::max();
325                 while (it != itend) {
326                         int pow = it->rest.ldegree(s);
327                         if (pow < min_pow)
328                                 min_pow = pow;
329                         ++it;
330                 }
331                 return min_pow;
332         }
333 }
334
335 /** Return coefficient of degree n in power series if s is the expansion
336  *  variable.  If the expansion point is nonzero, by definition the n=1
337  *  coefficient in s of a+b*(s-z)+c*(s-z)^2+Order((s-z)^3) is b (assuming
338  *  the expansion took place in the s in the first place).
339  *  If s is not the expansion variable, an attempt is made to convert the
340  *  series to a polynomial and return the corresponding coefficient from
341  *  there. */
342 ex pseries::coeff(const ex &s, int n) const
343 {
344         if (var.is_equal(s)) {
345                 if (seq.empty())
346                         return _ex0;
347                 
348                 // Binary search in sequence for given power
349                 numeric looking_for = numeric(n);
350                 int lo = 0, hi = seq.size() - 1;
351                 while (lo <= hi) {
352                         int mid = (lo + hi) / 2;
353                         GINAC_ASSERT(is_exactly_a<numeric>(seq[mid].coeff));
354                         int cmp = ex_to<numeric>(seq[mid].coeff).compare(looking_for);
355                         switch (cmp) {
356                                 case -1:
357                                         lo = mid + 1;
358                                         break;
359                                 case 0:
360                                         return seq[mid].rest;
361                                 case 1:
362                                         hi = mid - 1;
363                                         break;
364                                 default:
365                                         throw(std::logic_error("pseries::coeff: compare() didn't return -1, 0 or 1"));
366                         }
367                 }
368                 return _ex0;
369         } else
370                 return convert_to_poly().coeff(s, n);
371 }
372
373 /** Does nothing. */
374 ex pseries::collect(const ex &s, bool distributed) const
375 {
376         return *this;
377 }
378
379 /** Perform coefficient-wise automatic term rewriting rules in this class. */
380 ex pseries::eval(int level) const
381 {
382         if (level == 1)
383                 return this->hold();
384         
385         if (level == -max_recursion_level)
386                 throw (std::runtime_error("pseries::eval(): recursion limit exceeded"));
387         
388         // Construct a new series with evaluated coefficients
389         epvector new_seq;
390         new_seq.reserve(seq.size());
391         epvector::const_iterator it = seq.begin(), itend = seq.end();
392         while (it != itend) {
393                 new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
394                 ++it;
395         }
396         return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
397 }
398
399 /** Evaluate coefficients numerically. */
400 ex pseries::evalf(int level) const
401 {
402         if (level == 1)
403                 return *this;
404         
405         if (level == -max_recursion_level)
406                 throw (std::runtime_error("pseries::evalf(): recursion limit exceeded"));
407         
408         // Construct a new series with evaluated coefficients
409         epvector new_seq;
410         new_seq.reserve(seq.size());
411         epvector::const_iterator it = seq.begin(), itend = seq.end();
412         while (it != itend) {
413                 new_seq.push_back(expair(it->rest.evalf(level-1), it->coeff));
414                 ++it;
415         }
416         return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
417 }
418
419 ex pseries::conjugate() const
420 {
421         if(!var.info(info_flags::real))
422                 return conjugate_function(*this).hold();
423
424         epvector * newseq = conjugateepvector(seq);
425         ex newpoint = point.conjugate();
426
427         if (!newseq     && are_ex_trivially_equal(point, newpoint)) {
428                 return *this;
429         }
430
431         ex result = (new pseries(var==newpoint, newseq ? *newseq : seq))->setflag(status_flags::dynallocated);
432         if (newseq) {
433                 delete newseq;
434         }
435         return result;
436 }
437
438 ex pseries::real_part() const
439 {
440         if(!var.info(info_flags::real))
441                 return real_part_function(*this).hold();
442         ex newpoint = point.real_part();
443         if(newpoint != point)
444                 return real_part_function(*this).hold();
445
446         epvector v;
447         v.reserve(seq.size());
448         for(epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i)
449                 v.push_back(expair((i->rest).real_part(), i->coeff));
450         return (new pseries(var==point, v))->setflag(status_flags::dynallocated);
451 }
452
453 ex pseries::imag_part() const
454 {
455         if(!var.info(info_flags::real))
456                 return imag_part_function(*this).hold();
457         ex newpoint = point.real_part();
458         if(newpoint != point)
459                 return imag_part_function(*this).hold();
460
461         epvector v;
462         v.reserve(seq.size());
463         for(epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i)
464                 v.push_back(expair((i->rest).imag_part(), i->coeff));
465         return (new pseries(var==point, v))->setflag(status_flags::dynallocated);
466 }
467
468 ex pseries::eval_integ() const
469 {
470         epvector *newseq = NULL;
471         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
472                 if (newseq) {
473                         newseq->push_back(expair(i->rest.eval_integ(), i->coeff));
474                         continue;
475                 }
476                 ex newterm = i->rest.eval_integ();
477                 if (!are_ex_trivially_equal(newterm, i->rest)) {
478                         newseq = new epvector;
479                         newseq->reserve(seq.size());
480                         for (epvector::const_iterator j=seq.begin(); j!=i; ++j)
481                                 newseq->push_back(*j);
482                         newseq->push_back(expair(newterm, i->coeff));
483                 }
484         }
485
486         ex newpoint = point.eval_integ();
487         if (newseq || !are_ex_trivially_equal(newpoint, point))
488                 return (new pseries(var==newpoint, *newseq))
489                        ->setflag(status_flags::dynallocated);
490         return *this;
491 }
492
493 ex pseries::evalm() const
494 {
495         // evalm each coefficient
496         epvector newseq;
497         bool something_changed = false;
498         for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
499                 if (something_changed) {
500                         ex newcoeff = i->rest.evalm();
501                         if (!newcoeff.is_zero())
502                                 newseq.push_back(expair(newcoeff, i->coeff));
503                 }
504                 else {
505                         ex newcoeff = i->rest.evalm();
506                         if (!are_ex_trivially_equal(newcoeff, i->rest)) {
507                                 something_changed = true;
508                                 newseq.reserve(seq.size());
509                                 std::copy(seq.begin(), i, std::back_inserter<epvector>(newseq));
510                                 if (!newcoeff.is_zero())
511                                         newseq.push_back(expair(newcoeff, i->coeff));
512                         }
513                 }
514         }
515         if (something_changed)
516                 return (new pseries(var==point, newseq))->setflag(status_flags::dynallocated);
517         else
518                 return *this;
519 }
520
521 ex pseries::subs(const exmap & m, unsigned options) const
522 {
523         // If expansion variable is being substituted, convert the series to a
524         // polynomial and do the substitution there because the result might
525         // no longer be a power series
526         if (m.find(var) != m.end())
527                 return convert_to_poly(true).subs(m, options);
528         
529         // Otherwise construct a new series with substituted coefficients and
530         // expansion point
531         epvector newseq;
532         newseq.reserve(seq.size());
533         epvector::const_iterator it = seq.begin(), itend = seq.end();
534         while (it != itend) {
535                 newseq.push_back(expair(it->rest.subs(m, options), it->coeff));
536                 ++it;
537         }
538         return (new pseries(relational(var,point.subs(m, options)), newseq))->setflag(status_flags::dynallocated);
539 }
540
541 /** Implementation of ex::expand() for a power series.  It expands all the
542  *  terms individually and returns the resulting series as a new pseries. */
543 ex pseries::expand(unsigned options) const
544 {
545         epvector newseq;
546         epvector::const_iterator i = seq.begin(), end = seq.end();
547         while (i != end) {
548                 ex restexp = i->rest.expand();
549                 if (!restexp.is_zero())
550                         newseq.push_back(expair(restexp, i->coeff));
551                 ++i;
552         }
553         return (new pseries(relational(var,point), newseq))
554                 ->setflag(status_flags::dynallocated | (options == 0 ? status_flags::expanded : 0));
555 }
556
557 /** Implementation of ex::diff() for a power series.
558  *  @see ex::diff */
559 ex pseries::derivative(const symbol & s) const
560 {
561         epvector new_seq;
562         epvector::const_iterator it = seq.begin(), itend = seq.end();
563
564         if (s == var) {
565                 
566                 // FIXME: coeff might depend on var
567                 while (it != itend) {
568                         if (is_order_function(it->rest)) {
569                                 new_seq.push_back(expair(it->rest, it->coeff - 1));
570                         } else {
571                                 ex c = it->rest * it->coeff;
572                                 if (!c.is_zero())
573                                         new_seq.push_back(expair(c, it->coeff - 1));
574                         }
575                         ++it;
576                 }
577
578         } else {
579
580                 while (it != itend) {
581                         if (is_order_function(it->rest)) {
582                                 new_seq.push_back(*it);
583                         } else {
584                                 ex c = it->rest.diff(s);
585                                 if (!c.is_zero())
586                                         new_seq.push_back(expair(c, it->coeff));
587                         }
588                         ++it;
589                 }
590         }
591
592         return pseries(relational(var,point), new_seq);
593 }
594
595 ex pseries::convert_to_poly(bool no_order) const
596 {
597         ex e;
598         epvector::const_iterator it = seq.begin(), itend = seq.end();
599         
600         while (it != itend) {
601                 if (is_order_function(it->rest)) {
602                         if (!no_order)
603                                 e += Order(power(var - point, it->coeff));
604                 } else
605                         e += it->rest * power(var - point, it->coeff);
606                 ++it;
607         }
608         return e;
609 }
610
611 bool pseries::is_terminating() const
612 {
613         return seq.empty() || !is_order_function((seq.end()-1)->rest);
614 }
615
616 ex pseries::coeffop(size_t i) const
617 {
618         if (i >=nops())
619                 throw (std::out_of_range("coeffop() out of range"));
620         return seq[i].rest;
621 }
622
623 ex pseries::exponop(size_t i) const
624 {
625         if (i >= nops())
626                 throw (std::out_of_range("exponop() out of range"));
627         return seq[i].coeff;
628 }
629
630
631 /*
632  *  Implementations of series expansion
633  */
634
635 /** Default implementation of ex::series(). This performs Taylor expansion.
636  *  @see ex::series */
637 ex basic::series(const relational & r, int order, unsigned options) const
638 {
639         epvector seq;
640         const symbol &s = ex_to<symbol>(r.lhs());
641
642         // default for order-values that make no sense for Taylor expansion
643         if ((order <= 0) && this->has(s)) {
644                 seq.push_back(expair(Order(_ex1), order));
645                 return pseries(r, seq);
646         }
647
648         // do Taylor expansion
649         numeric fac = 1;
650         ex deriv = *this;
651         ex coeff = deriv.subs(r, subs_options::no_pattern);
652
653         if (!coeff.is_zero()) {
654                 seq.push_back(expair(coeff, _ex0));
655         }
656
657         int n;
658         for (n=1; n<order; ++n) {
659                 fac = fac.mul(n);
660                 // We need to test for zero in order to see if the series terminates.
661                 // The problem is that there is no such thing as a perfect test for
662                 // zero.  Expanding the term occasionally helps a little...
663                 deriv = deriv.diff(s).expand();
664                 if (deriv.is_zero())  // Series terminates
665                         return pseries(r, seq);
666
667                 coeff = deriv.subs(r, subs_options::no_pattern);
668                 if (!coeff.is_zero())
669                         seq.push_back(expair(fac.inverse() * coeff, n));
670         }
671         
672         // Higher-order terms, if present
673         deriv = deriv.diff(s);
674         if (!deriv.expand().is_zero())
675                 seq.push_back(expair(Order(_ex1), n));
676         return pseries(r, seq);
677 }
678
679
680 /** Implementation of ex::series() for symbols.
681  *  @see ex::series */
682 ex symbol::series(const relational & r, int order, unsigned options) const
683 {
684         epvector seq;
685         const ex point = r.rhs();
686         GINAC_ASSERT(is_a<symbol>(r.lhs()));
687
688         if (this->is_equal_same_type(ex_to<symbol>(r.lhs()))) {
689                 if (order > 0 && !point.is_zero())
690                         seq.push_back(expair(point, _ex0));
691                 if (order > 1)
692                         seq.push_back(expair(_ex1, _ex1));
693                 else
694                         seq.push_back(expair(Order(_ex1), numeric(order)));
695         } else
696                 seq.push_back(expair(*this, _ex0));
697         return pseries(r, seq);
698 }
699
700
701 /** Add one series object to another, producing a pseries object that
702  *  represents the sum.
703  *
704  *  @param other  pseries object to add with
705  *  @return the sum as a pseries */
706 ex pseries::add_series(const pseries &other) const
707 {
708         // Adding two series with different variables or expansion points
709         // results in an empty (constant) series 
710         if (!is_compatible_to(other)) {
711                 epvector nul;
712                 nul.push_back(expair(Order(_ex1), _ex0));
713                 return pseries(relational(var,point), nul);
714         }
715         
716         // Series addition
717         epvector new_seq;
718         epvector::const_iterator a = seq.begin();
719         epvector::const_iterator b = other.seq.begin();
720         epvector::const_iterator a_end = seq.end();
721         epvector::const_iterator b_end = other.seq.end();
722         int pow_a = std::numeric_limits<int>::max(), pow_b = std::numeric_limits<int>::max();
723         for (;;) {
724                 // If a is empty, fill up with elements from b and stop
725                 if (a == a_end) {
726                         while (b != b_end) {
727                                 new_seq.push_back(*b);
728                                 ++b;
729                         }
730                         break;
731                 } else
732                         pow_a = ex_to<numeric>((*a).coeff).to_int();
733                 
734                 // If b is empty, fill up with elements from a and stop
735                 if (b == b_end) {
736                         while (a != a_end) {
737                                 new_seq.push_back(*a);
738                                 ++a;
739                         }
740                         break;
741                 } else
742                         pow_b = ex_to<numeric>((*b).coeff).to_int();
743                 
744                 // a and b are non-empty, compare powers
745                 if (pow_a < pow_b) {
746                         // a has lesser power, get coefficient from a
747                         new_seq.push_back(*a);
748                         if (is_order_function((*a).rest))
749                                 break;
750                         ++a;
751                 } else if (pow_b < pow_a) {
752                         // b has lesser power, get coefficient from b
753                         new_seq.push_back(*b);
754                         if (is_order_function((*b).rest))
755                                 break;
756                         ++b;
757                 } else {
758                         // Add coefficient of a and b
759                         if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
760                                 new_seq.push_back(expair(Order(_ex1), (*a).coeff));
761                                 break;  // Order term ends the sequence
762                         } else {
763                                 ex sum = (*a).rest + (*b).rest;
764                                 if (!(sum.is_zero()))
765                                         new_seq.push_back(expair(sum, numeric(pow_a)));
766                                 ++a;
767                                 ++b;
768                         }
769                 }
770         }
771         return pseries(relational(var,point), new_seq);
772 }
773
774
775 /** Implementation of ex::series() for sums. This performs series addition when
776  *  adding pseries objects.
777  *  @see ex::series */
778 ex add::series(const relational & r, int order, unsigned options) const
779 {
780         ex acc; // Series accumulator
781         
782         // Get first term from overall_coeff
783         acc = overall_coeff.series(r, order, options);
784         
785         // Add remaining terms
786         epvector::const_iterator it = seq.begin();
787         epvector::const_iterator itend = seq.end();
788         for (; it!=itend; ++it) {
789                 ex op;
790                 if (is_exactly_a<pseries>(it->rest))
791                         op = it->rest;
792                 else
793                         op = it->rest.series(r, order, options);
794                 if (!it->coeff.is_equal(_ex1))
795                         op = ex_to<pseries>(op).mul_const(ex_to<numeric>(it->coeff));
796                 
797                 // Series addition
798                 acc = ex_to<pseries>(acc).add_series(ex_to<pseries>(op));
799         }
800         return acc;
801 }
802
803
804 /** Multiply a pseries object with a numeric constant, producing a pseries
805  *  object that represents the product.
806  *
807  *  @param other  constant to multiply with
808  *  @return the product as a pseries */
809 ex pseries::mul_const(const numeric &other) const
810 {
811         epvector new_seq;
812         new_seq.reserve(seq.size());
813         
814         epvector::const_iterator it = seq.begin(), itend = seq.end();
815         while (it != itend) {
816                 if (!is_order_function(it->rest))
817                         new_seq.push_back(expair(it->rest * other, it->coeff));
818                 else
819                         new_seq.push_back(*it);
820                 ++it;
821         }
822         return pseries(relational(var,point), new_seq);
823 }
824
825
826 /** Multiply one pseries object to another, producing a pseries object that
827  *  represents the product.
828  *
829  *  @param other  pseries object to multiply with
830  *  @return the product as a pseries */
831 ex pseries::mul_series(const pseries &other) const
832 {
833         // Multiplying two series with different variables or expansion points
834         // results in an empty (constant) series 
835         if (!is_compatible_to(other)) {
836                 epvector nul;
837                 nul.push_back(expair(Order(_ex1), _ex0));
838                 return pseries(relational(var,point), nul);
839         }
840
841         if (seq.empty() || other.seq.empty()) {
842                 return (new pseries(var==point, epvector()))
843                        ->setflag(status_flags::dynallocated);
844         }
845         
846         // Series multiplication
847         epvector new_seq;
848         int a_max = degree(var);
849         int b_max = other.degree(var);
850         int a_min = ldegree(var);
851         int b_min = other.ldegree(var);
852         int cdeg_min = a_min + b_min;
853         int cdeg_max = a_max + b_max;
854         
855         int higher_order_a = std::numeric_limits<int>::max();
856         int higher_order_b = std::numeric_limits<int>::max();
857         if (is_order_function(coeff(var, a_max)))
858                 higher_order_a = a_max + b_min;
859         if (is_order_function(other.coeff(var, b_max)))
860                 higher_order_b = b_max + a_min;
861         int higher_order_c = std::min(higher_order_a, higher_order_b);
862         if (cdeg_max >= higher_order_c)
863                 cdeg_max = higher_order_c - 1;
864         
865         for (int cdeg=cdeg_min; cdeg<=cdeg_max; ++cdeg) {
866                 ex co = _ex0;
867                 // c(i)=a(0)b(i)+...+a(i)b(0)
868                 for (int i=a_min; cdeg-i>=b_min; ++i) {
869                         ex a_coeff = coeff(var, i);
870                         ex b_coeff = other.coeff(var, cdeg-i);
871                         if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
872                                 co += a_coeff * b_coeff;
873                 }
874                 if (!co.is_zero())
875                         new_seq.push_back(expair(co, numeric(cdeg)));
876         }
877         if (higher_order_c < std::numeric_limits<int>::max())
878                 new_seq.push_back(expair(Order(_ex1), numeric(higher_order_c)));
879         return pseries(relational(var, point), new_seq);
880 }
881
882
883 /** Implementation of ex::series() for product. This performs series
884  *  multiplication when multiplying series.
885  *  @see ex::series */
886 ex mul::series(const relational & r, int order, unsigned options) const
887 {
888         pseries acc; // Series accumulator
889
890         GINAC_ASSERT(is_a<symbol>(r.lhs()));
891         const ex& sym = r.lhs();
892                 
893         // holds ldegrees of the series of individual factors
894         std::vector<int> ldegrees;
895         std::vector<bool> ldegree_redo;
896
897         // find minimal degrees
898         const epvector::const_iterator itbeg = seq.begin();
899         const epvector::const_iterator itend = seq.end();
900         // first round: obtain a bound up to which minimal degrees have to be
901         // considered
902         for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
903
904                 ex expon = it->coeff;
905                 int factor = 1;
906                 ex buf;
907                 if (expon.info(info_flags::integer)) {
908                         buf = it->rest;
909                         factor = ex_to<numeric>(expon).to_int();
910                 } else {
911                         buf = recombine_pair_to_ex(*it);
912                 }
913
914                 int real_ldegree = 0;
915                 bool flag_redo = false;
916                 try {
917                         real_ldegree = buf.expand().ldegree(sym-r.rhs());
918                 } catch (std::runtime_error) {}
919
920                 if (real_ldegree == 0) {
921                         if ( factor < 0 ) {
922                                 // This case must terminate, otherwise we would have division by
923                                 // zero.
924                                 int orderloop = 0;
925                                 do {
926                                         orderloop++;
927                                         real_ldegree = buf.series(r, orderloop, options).ldegree(sym);
928                                 } while (real_ldegree == orderloop);
929                         } else {
930                                 // Here it is possible that buf does not have a ldegree, therefore
931                                 // check only if ldegree is negative, otherwise reconsider the case
932                                 // in the second round.
933                                 real_ldegree = buf.series(r, 0, options).ldegree(sym);
934                                 if (real_ldegree == 0)
935                                         flag_redo = true;
936                         }
937                 }
938
939                 ldegrees.push_back(factor * real_ldegree);
940                 ldegree_redo.push_back(flag_redo);
941         }
942
943         int degbound = order-std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
944         // Second round: determine the remaining positive ldegrees by the series
945         // method.
946         // here we can ignore ldegrees larger than degbound
947         size_t j = 0;
948         for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
949                 if ( ldegree_redo[j] ) {
950                         ex expon = it->coeff;
951                         int factor = 1;
952                         ex buf;
953                         if (expon.info(info_flags::integer)) {
954                                 buf = it->rest;
955                                 factor = ex_to<numeric>(expon).to_int();
956                         } else {
957                                 buf = recombine_pair_to_ex(*it);
958                         }
959                         int real_ldegree = 0;
960                         int orderloop = 0;
961                         do {
962                                 orderloop++;
963                                 real_ldegree = buf.series(r, orderloop, options).ldegree(sym);
964                         } while ((real_ldegree == orderloop)
965                                         && ( factor*real_ldegree < degbound));
966                         ldegrees[j] = factor * real_ldegree;
967                         degbound -= factor * real_ldegree;
968                 }
969                 j++;
970         }
971
972         int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
973
974         if (degsum >= order) {
975                 epvector epv;
976                 epv.push_back(expair(Order(_ex1), order));
977                 return (new pseries(r, epv))->setflag(status_flags::dynallocated);
978         }
979
980         // Multiply with remaining terms
981         std::vector<int>::const_iterator itd = ldegrees.begin();
982         for (epvector::const_iterator it=itbeg; it!=itend; ++it, ++itd) {
983
984                 // do series expansion with adjusted order
985                 ex op = recombine_pair_to_ex(*it).series(r, order-degsum+(*itd), options);
986
987                 // Series multiplication
988                 if (it == itbeg)
989                         acc = ex_to<pseries>(op);
990                 else
991                         acc = ex_to<pseries>(acc.mul_series(ex_to<pseries>(op)));
992         }
993
994         return acc.mul_const(ex_to<numeric>(overall_coeff));
995 }
996
997
998 /** Compute the p-th power of a series.
999  *
1000  *  @param p  power to compute
1001  *  @param deg  truncation order of series calculation */
1002 ex pseries::power_const(const numeric &p, int deg) const
1003 {
1004         // method:
1005         // (due to Leonhard Euler)
1006         // let A(x) be this series and for the time being let it start with a
1007         // constant (later we'll generalize):
1008         //     A(x) = a_0 + a_1*x + a_2*x^2 + ...
1009         // We want to compute
1010         //     C(x) = A(x)^p
1011         //     C(x) = c_0 + c_1*x + c_2*x^2 + ...
1012         // Taking the derivative on both sides and multiplying with A(x) one
1013         // immediately arrives at
1014         //     C'(x)*A(x) = p*C(x)*A'(x)
1015         // Multiplying this out and comparing coefficients we get the recurrence
1016         // formula
1017         //     c_i = (i*p*a_i*c_0 + ((i-1)*p-1)*a_{i-1}*c_1 + ...
1018         //                    ... + (p-(i-1))*a_1*c_{i-1})/(a_0*i)
1019         // which can easily be solved given the starting value c_0 = (a_0)^p.
1020         // For the more general case where the leading coefficient of A(x) is not
1021         // a constant, just consider A2(x) = A(x)*x^m, with some integer m and
1022         // repeat the above derivation.  The leading power of C2(x) = A2(x)^2 is
1023         // then of course x^(p*m) but the recurrence formula still holds.
1024         
1025         if (seq.empty()) {
1026                 // as a special case, handle the empty (zero) series honoring the
1027                 // usual power laws such as implemented in power::eval()
1028                 if (p.real().is_zero())
1029                         throw std::domain_error("pseries::power_const(): pow(0,I) is undefined");
1030                 else if (p.real().is_negative())
1031                         throw pole_error("pseries::power_const(): division by zero",1);
1032                 else
1033                         return *this;
1034         }
1035         
1036         const int ldeg = ldegree(var);
1037         if (!(p*ldeg).is_integer())
1038                 throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
1039
1040         // adjust number of coefficients
1041         int numcoeff = deg - (p*ldeg).to_int();
1042         if (numcoeff <= 0) {
1043                 epvector epv;
1044                 epv.reserve(1);
1045                 epv.push_back(expair(Order(_ex1), deg));
1046                 return (new pseries(relational(var,point), epv))
1047                        ->setflag(status_flags::dynallocated);
1048         }
1049         
1050         // O(x^n)^(-m) is undefined
1051         if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
1052                 throw pole_error("pseries::power_const(): division by zero",1);
1053         
1054         // Compute coefficients of the powered series
1055         exvector co;
1056         co.reserve(numcoeff);
1057         co.push_back(power(coeff(var, ldeg), p));
1058         for (int i=1; i<numcoeff; ++i) {
1059                 ex sum = _ex0;
1060                 for (int j=1; j<=i; ++j) {
1061                         ex c = coeff(var, j + ldeg);
1062                         if (is_order_function(c)) {
1063                                 co.push_back(Order(_ex1));
1064                                 break;
1065                         } else
1066                                 sum += (p * j - (i - j)) * co[i - j] * c;
1067                 }
1068                 co.push_back(sum / coeff(var, ldeg) / i);
1069         }
1070         
1071         // Construct new series (of non-zero coefficients)
1072         epvector new_seq;
1073         bool higher_order = false;
1074         for (int i=0; i<numcoeff; ++i) {
1075                 if (!co[i].is_zero())
1076                         new_seq.push_back(expair(co[i], p * ldeg + i));
1077                 if (is_order_function(co[i])) {
1078                         higher_order = true;
1079                         break;
1080                 }
1081         }
1082         if (!higher_order)
1083                 new_seq.push_back(expair(Order(_ex1), p * ldeg + numcoeff));
1084
1085         return pseries(relational(var,point), new_seq);
1086 }
1087
1088
1089 /** Return a new pseries object with the powers shifted by deg. */
1090 pseries pseries::shift_exponents(int deg) const
1091 {
1092         epvector newseq = seq;
1093         epvector::iterator i = newseq.begin(), end  = newseq.end();
1094         while (i != end) {
1095                 i->coeff += deg;
1096                 ++i;
1097         }
1098         return pseries(relational(var, point), newseq);
1099 }
1100
1101
1102 /** Implementation of ex::series() for powers. This performs Laurent expansion
1103  *  of reciprocals of series at singularities.
1104  *  @see ex::series */
1105 ex power::series(const relational & r, int order, unsigned options) const
1106 {
1107         // If basis is already a series, just power it
1108         if (is_exactly_a<pseries>(basis))
1109                 return ex_to<pseries>(basis).power_const(ex_to<numeric>(exponent), order);
1110
1111         // Basis is not a series, may there be a singularity?
1112         bool must_expand_basis = false;
1113         try {
1114                 basis.subs(r, subs_options::no_pattern);
1115         } catch (pole_error) {
1116                 must_expand_basis = true;
1117         }
1118
1119         bool exponent_is_regular = true;
1120         try {
1121                 exponent.subs(r, subs_options::no_pattern);
1122         } catch (pole_error) {
1123                 exponent_is_regular = false;
1124         }
1125
1126         if (!exponent_is_regular) {
1127                 ex l = exponent*log(basis);
1128                 // this == exp(l);
1129                 ex le = l.series(r, order, options);
1130                 // Note: expanding exp(l) won't help, since that will attempt
1131                 // Taylor expansion, and fail (because exponent is "singular")
1132                 // Still l itself might be expanded in Taylor series.
1133                 // Examples:
1134                 // sin(x)/x*log(cos(x))
1135                 // 1/x*log(1 + x)
1136                 return exp(le).series(r, order, options);
1137                 // Note: if l happens to have a Laurent expansion (with
1138                 // negative powers of (var - point)), expanding exp(le)
1139                 // will barf (which is The Right Thing).
1140         }
1141
1142         // Is the expression of type something^(-int)?
1143         if (!must_expand_basis && !exponent.info(info_flags::negint)
1144          && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
1145                 return basic::series(r, order, options);
1146
1147         // Is the expression of type 0^something?
1148         if (!must_expand_basis && !basis.subs(r, subs_options::no_pattern).is_zero()
1149          && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
1150                 return basic::series(r, order, options);
1151
1152         // Singularity encountered, is the basis equal to (var - point)?
1153         if (basis.is_equal(r.lhs() - r.rhs())) {
1154                 epvector new_seq;
1155                 if (ex_to<numeric>(exponent).to_int() < order)
1156                         new_seq.push_back(expair(_ex1, exponent));
1157                 else
1158                         new_seq.push_back(expair(Order(_ex1), exponent));
1159                 return pseries(r, new_seq);
1160         }
1161
1162         // No, expand basis into series
1163
1164         numeric numexp;
1165         if (is_a<numeric>(exponent)) {
1166                 numexp = ex_to<numeric>(exponent);
1167         } else {
1168                 numexp = 0;
1169         }
1170         const ex& sym = r.lhs();
1171         // find existing minimal degree
1172         ex eb = basis.expand();
1173         int real_ldegree = 0;
1174         if (eb.info(info_flags::rational_function))
1175                 real_ldegree = eb.ldegree(sym-r.rhs());
1176         if (real_ldegree == 0) {
1177                 int orderloop = 0;
1178                 do {
1179                         orderloop++;
1180                         real_ldegree = basis.series(r, orderloop, options).ldegree(sym);
1181                 } while (real_ldegree == orderloop);
1182         }
1183
1184         if (!(real_ldegree*numexp).is_integer())
1185                 throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
1186         ex e = basis.series(r, (order + real_ldegree*(1-numexp)).to_int(), options);
1187         
1188         ex result;
1189         try {
1190                 result = ex_to<pseries>(e).power_const(numexp, order);
1191         } catch (pole_error) {
1192                 epvector ser;
1193                 ser.push_back(expair(Order(_ex1), order));
1194                 result = pseries(r, ser);
1195         }
1196
1197         return result;
1198 }
1199
1200
1201 /** Re-expansion of a pseries object. */
1202 ex pseries::series(const relational & r, int order, unsigned options) const
1203 {
1204         const ex p = r.rhs();
1205         GINAC_ASSERT(is_a<symbol>(r.lhs()));
1206         const symbol &s = ex_to<symbol>(r.lhs());
1207         
1208         if (var.is_equal(s) && point.is_equal(p)) {
1209                 if (order > degree(s))
1210                         return *this;
1211                 else {
1212                         epvector new_seq;
1213                         epvector::const_iterator it = seq.begin(), itend = seq.end();
1214                         while (it != itend) {
1215                                 int o = ex_to<numeric>(it->coeff).to_int();
1216                                 if (o >= order) {
1217                                         new_seq.push_back(expair(Order(_ex1), o));
1218                                         break;
1219                                 }
1220                                 new_seq.push_back(*it);
1221                                 ++it;
1222                         }
1223                         return pseries(r, new_seq);
1224                 }
1225         } else
1226                 return convert_to_poly().series(r, order, options);
1227 }
1228
1229 ex integral::series(const relational & r, int order, unsigned options) const
1230 {
1231         if (x.subs(r) != x)
1232                 throw std::logic_error("Cannot series expand wrt dummy variable");
1233         
1234         // Expanding integrant with r substituted taken in boundaries.
1235         ex fseries = f.series(r, order, options);
1236         epvector fexpansion;
1237         fexpansion.reserve(fseries.nops());
1238         for (size_t i=0; i<fseries.nops(); ++i) {
1239                 ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
1240                 currcoeff = (currcoeff == Order(_ex1))
1241                         ? currcoeff
1242                         : integral(x, a.subs(r), b.subs(r), currcoeff);
1243                 if (currcoeff != 0)
1244                         fexpansion.push_back(
1245                                 expair(currcoeff, ex_to<pseries>(fseries).exponop(i)));
1246         }
1247
1248         // Expanding lower boundary
1249         ex result = (new pseries(r, fexpansion))->setflag(status_flags::dynallocated);
1250         ex aseries = (a-a.subs(r)).series(r, order, options);
1251         fseries = f.series(x == (a.subs(r)), order, options);
1252         for (size_t i=0; i<fseries.nops(); ++i) {
1253                 ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
1254                 if (is_order_function(currcoeff))
1255                         break;
1256                 ex currexpon = ex_to<pseries>(fseries).exponop(i);
1257                 int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
1258                 currcoeff = currcoeff.series(r, orderforf);
1259                 ex term = ex_to<pseries>(aseries).power_const(ex_to<numeric>(currexpon+1),order);
1260                 term = ex_to<pseries>(term).mul_const(ex_to<numeric>(-1/(currexpon+1)));
1261                 term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
1262                 result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
1263         }
1264
1265         // Expanding upper boundary
1266         ex bseries = (b-b.subs(r)).series(r, order, options);
1267         fseries = f.series(x == (b.subs(r)), order, options);
1268         for (size_t i=0; i<fseries.nops(); ++i) {
1269                 ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
1270                 if (is_order_function(currcoeff))
1271                         break;
1272                 ex currexpon = ex_to<pseries>(fseries).exponop(i);
1273                 int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
1274                 currcoeff = currcoeff.series(r, orderforf);
1275                 ex term = ex_to<pseries>(bseries).power_const(ex_to<numeric>(currexpon+1),order);
1276                 term = ex_to<pseries>(term).mul_const(ex_to<numeric>(1/(currexpon+1)));
1277                 term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
1278                 result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
1279         }
1280
1281         return result;
1282 }
1283
1284
1285 /** Compute the truncated series expansion of an expression.
1286  *  This function returns an expression containing an object of class pseries 
1287  *  to represent the series. If the series does not terminate within the given
1288  *  truncation order, the last term of the series will be an order term.
1289  *
1290  *  @param r  expansion relation, lhs holds variable and rhs holds point
1291  *  @param order  truncation order of series calculations
1292  *  @param options  of class series_options
1293  *  @return an expression holding a pseries object */
1294 ex ex::series(const ex & r, int order, unsigned options) const
1295 {
1296         ex e;
1297         relational rel_;
1298         
1299         if (is_a<relational>(r))
1300                 rel_ = ex_to<relational>(r);
1301         else if (is_a<symbol>(r))
1302                 rel_ = relational(r,_ex0);
1303         else
1304                 throw (std::logic_error("ex::series(): expansion point has unknown type"));
1305         
1306         e = bp->series(rel_, order, options);
1307         return e;
1308 }
1309
1310 GINAC_BIND_UNARCHIVER(pseries);
1311
1312 } // namespace GiNaC