* methods for series expansion. */
/*
- * GiNaC Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
*/
-#include <iostream>
+#include <numeric>
#include <stdexcept>
#include "pseries.h"
epvector seq;
numeric fac = 1;
ex deriv = *this;
- ex coeff = deriv.subs(r);
+ ex coeff = deriv.subs(r, subs_options::no_pattern);
const symbol &s = ex_to<symbol>(r.lhs());
if (!coeff.is_zero())
seq.push_back(expair(coeff, _ex0));
int n;
- for (n=1; n<order; ++n) {
+ for (n=1; n<=order; ++n) {
fac = fac.mul(n);
// We need to test for zero in order to see if the series terminates.
// The problem is that there is no such thing as a perfect test for
if (deriv.is_zero()) // Series terminates
return pseries(r, seq);
- coeff = deriv.subs(r);
+ coeff = deriv.subs(r, subs_options::no_pattern);
if (!coeff.is_zero())
seq.push_back(expair(fac.inverse() * coeff, n));
}
// Higher-order terms, if present
- deriv = deriv.diff(s);
- if (!deriv.expand().is_zero())
- seq.push_back(expair(Order(_ex1), n));
+ int ldeg;
+ try {
+ ldeg = std::abs(deriv.ldegree(s));
+ }
+ catch (std::runtime_error) {
+ ldeg = 0;
+ }
+ if (ldeg != 0) {
+ // pure polynomial
+ if (!deriv.subs(r, subs_options::no_pattern).is_zero()) {
+ seq.push_back(expair(Order(_ex1), n));
+ } else {
+ seq.push_back(expair(Order(_ex1), n+ldeg-1));
+ }
+ } else {
+ // something more complicated -> loop until next coefficient is found
+ for (;; ++n) {
+ deriv = deriv.diff(s).expand();
+ if (deriv.is_zero()) {
+ break;
+ }
+ if (!deriv.subs(r, subs_options::no_pattern).is_zero()) {
+ seq.push_back(expair(Order(_ex1), n));
+ break;
+ }
+ }
+ }
+
return pseries(r, seq);
}
{
pseries acc; // Series accumulator
+ GINAC_ASSERT(is_a<symbol>(r.lhs()));
+ const ex& sym = r.lhs();
+
+ // holds ldegrees of the series of individual factors
+ std::vector<int> ldegrees;
+ // flag if series have to be re-calculated
+ bool negldegree = false;
+
// Multiply with remaining terms
const epvector::const_iterator itbeg = seq.begin();
const epvector::const_iterator itend = seq.end();
for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
ex op = recombine_pair_to_ex(*it).series(r, order, options);
+ int ldeg = op.ldegree(sym);
+ ldegrees.push_back(ldeg);
+ if (ldeg < 0) {
+ negldegree = true;
+ }
+
// Series multiplication
if (it==itbeg)
acc = ex_to<pseries>(op);
else
acc = ex_to<pseries>(acc.mul_series(ex_to<pseries>(op)));
}
- return acc.mul_const(ex_to<numeric>(overall_coeff));
+
+ if (seq.size() > 1) {
+
+ // re-calculation of series
+
+ pseries newacc;
+
+ const int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
+
+ // Multiply with remaining terms
+ const epvector::const_reverse_iterator itbeg = seq.rbegin();
+ const epvector::const_reverse_iterator itend = seq.rend();
+ std::vector<int>::const_reverse_iterator itd = ldegrees.rbegin();
+ for (epvector::const_reverse_iterator it=itbeg; it!=itend; ++it, ++itd) {
+
+ ex op = recombine_pair_to_ex(*it).series(r, order-degsum+(*itd), options);
+
+ // Series multiplication
+ if (it==itbeg)
+ newacc = ex_to<pseries>(op);
+ else
+ newacc = ex_to<pseries>(newacc.mul_series(ex_to<pseries>(op)));
+ }
+ return newacc.mul_const(ex_to<numeric>(overall_coeff));
+ } else {
+ return acc.mul_const(ex_to<numeric>(overall_coeff));
+ }
}
if (!(p*ldeg).is_integer())
throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
+ // adjust number of coefficients
+ deg = deg - p.to_int()*ldeg;
+
// O(x^n)^(-m) is undefined
if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
throw pole_error("pseries::power_const(): division by zero",1);
co.reserve(deg);
co.push_back(power(coeff(var, ldeg), p));
bool all_sums_zero = true;
- for (int i=1; i<deg; ++i) {
+ for (int i=1; i<=deg; ++i) {
ex sum = _ex0;
for (int j=1; j<=i; ++j) {
ex c = coeff(var, j + ldeg);
// Construct new series (of non-zero coefficients)
epvector new_seq;
bool higher_order = false;
- for (int i=0; i<deg; ++i) {
+ for (int i=0; i<=deg; ++i) {
if (!co[i].is_zero())
new_seq.push_back(expair(co[i], p * ldeg + i));
if (is_order_function(co[i])) {
}
}
if (!higher_order && !all_sums_zero)
- new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
+ new_seq.push_back(expair(Order(_ex1), p * ldeg + deg + 1));
+
return pseries(relational(var,point), new_seq);
}
// Basis is not a series, may there be a singularity?
bool must_expand_basis = false;
try {
- basis.subs(r);
+ basis.subs(r, subs_options::no_pattern);
} catch (pole_error) {
must_expand_basis = true;
}
// Is the expression of type something^(-int)?
if (!must_expand_basis && !exponent.info(info_flags::negint))
return basic::series(r, order, options);
-
+
// Is the expression of type 0^something?
- if (!must_expand_basis && !basis.subs(r).is_zero())
+ if (!must_expand_basis && !basis.subs(r, subs_options::no_pattern).is_zero())
return basic::series(r, order, options);
// Singularity encountered, is the basis equal to (var - point)?
}
// No, expand basis into series
+ int intexp = ex_to<numeric>(exponent).to_int();
ex e = basis.series(r, order, options);
- return ex_to<pseries>(e).power_const(ex_to<numeric>(exponent), order);
+ int ldeg = ex_to<pseries>(e).ldegree(r.lhs());
+ if (intexp * ldeg < 0) {
+ e = basis.series(r, order + ldeg*(1-intexp), options);
+ }
+ return ex_to<pseries>(e).power_const(intexp, order);
}
git://www.ginac.de / ginac.git / blobdiff