3 * Implementation of class for extended truncated power-series and
4 * methods for series expansion. */
7 * GiNaC Copyright (C) 1999 Johannes Gutenberg University Mainz, Germany
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.
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.
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
29 #include "relational.h"
34 * Default constructor, destructor, copy constructor, assignment operator and helpers
37 series::series() : basic(TINFO_series)
39 debugmsg("series default constructor", LOGLEVEL_CONSTRUCT);
44 debugmsg("series destructor", LOGLEVEL_DESTRUCT);
48 series::series(series const &other)
50 debugmsg("series copy constructor", LOGLEVEL_CONSTRUCT);
54 series const &series::operator=(series const & other)
56 debugmsg("series operator=", LOGLEVEL_ASSIGNMENT);
64 void series::copy(series const &other)
66 inherited::copy(other);
72 void series::destroy(bool call_parent)
75 inherited::destroy(call_parent);
83 /** Construct series from a vector of coefficients and powers.
84 * expair.rest holds the coefficient, expair.coeff holds the power.
85 * The powers must be integers (positive or negative) and in ascending order;
86 * the last coefficient can be Order(exONE()) to represent a truncated,
87 * non-terminating series.
89 * @param var_ series variable (must hold a symbol)
90 * @param point_ expansion point
91 * @param ops_ vector of {coefficient, power} pairs (coefficient must not be zero)
92 * @return newly constructed series */
93 series::series(ex const &var_, ex const &point_, epvector const &ops_)
94 : basic(TINFO_series), seq(ops_), var(var_), point(point_)
96 debugmsg("series constructor from ex,ex,epvector", LOGLEVEL_CONSTRUCT);
97 ASSERT(is_ex_exactly_of_type(var_, symbol));
102 * Functions overriding virtual functions from base classes
105 basic *series::duplicate() const
107 debugmsg("series duplicate", LOGLEVEL_DUPLICATE);
108 return new series(*this);
111 // Highest degree of variable
112 int series::degree(symbol const &s) const
114 if (var.is_equal(s)) {
115 // Return last exponent
117 return ex_to_numeric((*(seq.end() - 1)).coeff).to_int();
121 epvector::const_iterator it = seq.begin(), itend = seq.end();
124 int max_pow = INT_MIN;
125 while (it != itend) {
126 int pow = it->rest.degree(s);
135 // Lowest degree of variable
136 int series::ldegree(symbol const &s) const
138 if (var.is_equal(s)) {
139 // Return first exponent
141 return ex_to_numeric((*(seq.begin())).coeff).to_int();
145 epvector::const_iterator it = seq.begin(), itend = seq.end();
148 int min_pow = INT_MAX;
149 while (it != itend) {
150 int pow = it->rest.ldegree(s);
159 // Coefficient of variable
160 ex series::coeff(symbol const &s, int n) const
162 if (var.is_equal(s)) {
163 epvector::const_iterator it = seq.begin(), itend = seq.end();
164 while (it != itend) {
165 int pow = ex_to_numeric(it->coeff).to_int();
174 return convert_to_poly().coeff(s, n);
177 ex series::eval(int level) const
182 // Construct a new series with evaluated coefficients
184 new_seq.reserve(seq.size());
185 epvector::const_iterator it = seq.begin(), itend = seq.end();
186 while (it != itend) {
187 new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
190 return (new series(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
193 ex series::evalf(int level) const
195 return convert_to_poly().evalf(level);
200 * Construct expression (polynomial) out of series
203 /** Convert a series object to an ordinary polynomial.
205 * @param no_order flag: discard higher order terms */
206 ex series::convert_to_poly(bool no_order) const
209 epvector::const_iterator it = seq.begin(), itend = seq.end();
211 while (it != itend) {
212 if (is_order_function(it->rest)) {
214 e += Order(power(var - point, it->coeff));
216 e += it->rest * power(var - point, it->coeff);
224 * Implementation of series expansion
227 /** Default implementation of ex::series(). This performs Taylor expansion.
229 ex basic::series(symbol const & s, ex const & point, int order) const
234 ex coeff = deriv.subs(s == point);
235 if (!coeff.is_zero())
236 seq.push_back(expair(coeff, numeric(0)));
239 for (n=1; n<order; n++) {
240 fac = fac.mul(numeric(n));
241 deriv = deriv.diff(s).expand();
242 if (deriv.is_zero()) {
244 return series::series(s, point, seq);
246 coeff = power(fac, -1) * deriv.subs(s == point);
247 if (!coeff.is_zero())
248 seq.push_back(expair(coeff, numeric(n)));
251 // Higher-order terms, if present
252 deriv = deriv.diff(s);
253 if (!deriv.is_zero())
254 seq.push_back(expair(Order(exONE()), numeric(n)));
255 return series::series(s, point, seq);
259 /** Add one series object to another, producing a series object that represents
262 * @param other series object to add with
263 * @return the sum as a series */
264 ex series::add_series(const series &other) const
266 // Adding two series with different variables or expansion points
267 // results in an empty (constant) series
268 if (!is_compatible_to(other)) {
270 nul.push_back(expair(Order(exONE()), exZERO()));
271 return series(var, point, nul);
276 epvector::const_iterator a = seq.begin();
277 epvector::const_iterator b = other.seq.begin();
278 epvector::const_iterator a_end = seq.end();
279 epvector::const_iterator b_end = other.seq.end();
280 int pow_a = INT_MAX, pow_b = INT_MAX;
282 // If a is empty, fill up with elements from b and stop
285 new_seq.push_back(*b);
290 pow_a = ex_to_numeric((*a).coeff).to_int();
292 // If b is empty, fill up with elements from a and stop
295 new_seq.push_back(*a);
300 pow_b = ex_to_numeric((*b).coeff).to_int();
302 // a and b are non-empty, compare powers
304 // a has lesser power, get coefficient from a
305 new_seq.push_back(*a);
306 if (is_order_function((*a).rest))
309 } else if (pow_b < pow_a) {
310 // b has lesser power, get coefficient from b
311 new_seq.push_back(*b);
312 if (is_order_function((*b).rest))
316 // Add coefficient of a and b
317 if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
318 new_seq.push_back(expair(Order(exONE()), (*a).coeff));
319 break; // Order term ends the sequence
321 ex sum = (*a).rest + (*b).rest;
322 if (!(sum.is_zero()))
323 new_seq.push_back(expair(sum, numeric(pow_a)));
329 return series(var, point, new_seq);
333 /** Implementation of ex::series() for sums. This performs series addition when
334 * adding series objects.
337 ex add::series(symbol const & s, ex const & point, int order) const
339 ex acc; // Series accumulator
342 epvector::const_iterator it = seq.begin();
343 epvector::const_iterator itend = seq.end();
345 if (is_ex_exactly_of_type(it->rest, series))
348 acc = it->rest.series(s, point, order);
349 if (!it->coeff.is_equal(exONE()))
350 acc = ex_to_series(acc).mul_const(ex_to_numeric(it->coeff));
354 // Add remaining terms
355 for (; it!=itend; it++) {
357 if (is_ex_exactly_of_type(it->rest, series))
360 op = it->rest.series(s, point, order);
361 if (!it->coeff.is_equal(exONE()))
362 op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
365 acc = ex_to_series(acc).add_series(ex_to_series(op));
370 ex add::series(symbol const & s, ex const & point, int order) const
372 ex acc; // Series accumulator
374 // Get first term from overall_coeff
375 acc = overall_coeff.series(s,point,order);
377 // Add remaining terms
378 epvector::const_iterator it = seq.begin();
379 epvector::const_iterator itend = seq.end();
380 for (; it!=itend; it++) {
382 if (is_ex_exactly_of_type(it->rest, series))
385 op = it->rest.series(s, point, order);
386 if (!it->coeff.is_equal(exONE()))
387 op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
390 acc = ex_to_series(acc).add_series(ex_to_series(op));
396 /** Multiply a series object with a numeric constant, producing a series object
397 * that represents the product.
399 * @param other constant to multiply with
400 * @return the product as a series */
401 ex series::mul_const(const numeric &other) const
404 new_seq.reserve(seq.size());
406 epvector::const_iterator it = seq.begin(), itend = seq.end();
407 while (it != itend) {
408 if (!is_order_function(it->rest))
409 new_seq.push_back(expair(it->rest * other, it->coeff));
411 new_seq.push_back(*it);
414 return series(var, point, new_seq);
418 /** Multiply one series object to another, producing a series object that
419 * represents the product.
421 * @param other series object to multiply with
422 * @return the product as a series */
423 ex series::mul_series(const series &other) const
425 // Multiplying two series with different variables or expansion points
426 // results in an empty (constant) series
427 if (!is_compatible_to(other)) {
429 nul.push_back(expair(Order(exONE()), exZERO()));
430 return series(var, point, nul);
433 // Series multiplication
436 const symbol *s = static_cast<symbol *>(var.bp);
437 int a_max = degree(*s);
438 int b_max = other.degree(*s);
439 int a_min = ldegree(*s);
440 int b_min = other.ldegree(*s);
441 int cdeg_min = a_min + b_min;
442 int cdeg_max = a_max + b_max;
444 int higher_order_a = INT_MAX;
445 int higher_order_b = INT_MAX;
446 if (is_order_function(coeff(*s, a_max)))
447 higher_order_a = a_max + b_min;
448 if (is_order_function(other.coeff(*s, b_max)))
449 higher_order_b = b_max + a_min;
450 int higher_order_c = min(higher_order_a, higher_order_b);
451 if (cdeg_max >= higher_order_c)
452 cdeg_max = higher_order_c - 1;
454 for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) {
456 // c(i)=a(0)b(i)+...+a(i)b(0)
457 for (int i=a_min; cdeg-i>=b_min; i++) {
458 ex a_coeff = coeff(*s, i);
459 ex b_coeff = other.coeff(*s, cdeg-i);
460 if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
461 co += coeff(*s, i) * other.coeff(*s, cdeg-i);
464 new_seq.push_back(expair(co, numeric(cdeg)));
466 if (higher_order_c < INT_MAX)
467 new_seq.push_back(expair(Order(exONE()), numeric(higher_order_c)));
468 return series::series(var, point, new_seq);
472 /** Implementation of ex::series() for product. This performs series multiplication when multiplying series.
475 ex mul::series(symbol const & s, ex const & point, int order) const
477 ex acc; // Series accumulator
480 epvector::const_iterator it = seq.begin();
481 epvector::const_iterator itend = seq.end();
483 if (is_ex_exactly_of_type(it->rest, series))
486 acc = it->rest.series(s, point, order);
487 if (!it->coeff.is_equal(exONE()))
488 acc = ex_to_series(acc).power_const(ex_to_numeric(it->coeff), order);
492 // Multiply with remaining terms
493 for (; it!=itend; it++) {
495 if (op.info(info_flags::numeric)) {
496 // series * const (special case, faster)
497 ex f = power(op, it->coeff);
498 acc = ex_to_series(acc).mul_const(ex_to_numeric(f));
500 } else if (!is_ex_exactly_of_type(op, series))
501 op = op.series(s, point, order);
502 if (!it->coeff.is_equal(exONE()))
503 op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
505 // Series multiplication
506 acc = ex_to_series(acc).mul_series(ex_to_series(op));
512 ex mul::series(symbol const & s, ex const & point, int order) const
514 ex acc; // Series accumulator
516 // Get first term from overall_coeff
517 acc = overall_coeff.series(s, point, order);
519 // Multiply with remaining terms
520 epvector::const_iterator it = seq.begin();
521 epvector::const_iterator itend = seq.end();
522 for (; it!=itend; it++) {
524 if (op.info(info_flags::numeric)) {
525 // series * const (special case, faster)
526 ex f = power(op, it->coeff);
527 acc = ex_to_series(acc).mul_const(ex_to_numeric(f));
529 } else if (!is_ex_exactly_of_type(op, series))
530 op = op.series(s, point, order);
531 if (!it->coeff.is_equal(exONE()))
532 op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
534 // Series multiplication
535 acc = ex_to_series(acc).mul_series(ex_to_series(op));
541 /** Compute the p-th power of a series.
543 * @param p power to compute
544 * @param deg truncation order of series calculation */
545 ex series::power_const(const numeric &p, int deg) const
548 const symbol *s = static_cast<symbol *>(var.bp);
549 int ldeg = ldegree(*s);
551 // Calculate coefficients of powered series
555 co.push_back(co0 = power(coeff(*s, ldeg), p));
556 bool all_sums_zero = true;
557 for (i=1; i<deg; i++) {
559 for (int j=1; j<=i; j++) {
560 ex c = coeff(*s, j + ldeg);
561 if (is_order_function(c)) {
562 co.push_back(Order(exONE()));
565 sum += (p * j - (i - j)) * co[i - j] * c;
568 all_sums_zero = false;
569 co.push_back(co0 * sum / numeric(i));
572 // Construct new series (of non-zero coefficients)
574 bool higher_order = false;
575 for (i=0; i<deg; i++) {
576 if (!co[i].is_zero())
577 new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
578 if (is_order_function(co[i])) {
583 if (!higher_order && !all_sums_zero)
584 new_seq.push_back(expair(Order(exONE()), numeric(deg) + p * ldeg));
585 return series::series(var, point, new_seq);
589 /** Implementation of ex::series() for powers. This performs Laurent expansion
590 * of reciprocals of series at singularities.
592 ex power::series(symbol const & s, ex const & point, int order) const
595 if (!is_ex_exactly_of_type(basis, series)) {
596 // Basis is not a series, may there be a singulary?
597 if (!exponent.info(info_flags::negint))
598 return basic::series(s, point, order);
600 // Expression is of type something^(-int), check for singularity
601 if (!basis.subs(s == point).is_zero())
602 return basic::series(s, point, order);
604 // Singularity encountered, expand basis into series
605 e = basis.series(s, point, order);
612 return ex_to_series(e).power_const(ex_to_numeric(exponent), order);
616 /** Compute the truncated series expansion of an expression.
617 * This function returns an expression containing an object of class series to
618 * represent the series. If the series does not terminate within the given
619 * truncation order, the last term of the series will be an order term.
621 * @param s expansion variable
622 * @param point expansion point
623 * @param order truncation order of series calculations
624 * @return an expression holding a series object */
625 ex ex::series(symbol const &s, ex const &point, int order) const
628 return bp->series(s, point, order);
633 const series some_series;
634 type_info const & typeid_series = typeid(some_series);