#include "relational.h"
#include "operators.h"
#include "symbol.h"
+#include "integral.h"
#include "archive.h"
#include "utils.h"
if (i >= seq.size())
throw (std::out_of_range("op() out of range"));
+ if (is_order_function(seq[i].rest))
+ return Order(power(var-point, seq[i].coeff));
return seq[i].rest * power(var - point, seq[i].coeff);
}
return result;
}
+ex pseries::eval_integ() const
+{
+ epvector *newseq = NULL;
+ for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+ if (newseq) {
+ newseq->push_back(expair(i->rest.eval_integ(), i->coeff));
+ continue;
+ }
+ ex newterm = i->rest.eval_integ();
+ if (!are_ex_trivially_equal(newterm, i->rest)) {
+ newseq = new epvector;
+ newseq->reserve(seq.size());
+ for (epvector::const_iterator j=seq.begin(); j!=i; ++j)
+ newseq->push_back(*j);
+ newseq->push_back(expair(newterm, i->coeff));
+ }
+ }
+
+ ex newpoint = point.eval_integ();
+ if (newseq || !are_ex_trivially_equal(newpoint, point))
+ return (new pseries(var==newpoint, *newseq))
+ ->setflag(status_flags::dynallocated);
+ return *this;
+}
+
ex pseries::subs(const exmap & m, unsigned options) const
{
// If expansion variable is being substituted, convert the series to a
return seq.empty() || !is_order_function((seq.end()-1)->rest);
}
+ex pseries::coeffop(size_t i) const
+{
+ if (i >=nops())
+ throw (std::out_of_range("coeffop() out of range"));
+ return seq[i].rest;
+}
+
+ex pseries::exponop(size_t i) const
+{
+ if (i >= nops())
+ throw (std::out_of_range("exponop() out of range"));
+ return seq[i].coeff;
+}
+
/*
* Implementations of series expansion
nul.push_back(expair(Order(_ex1), _ex0));
return pseries(relational(var,point), nul);
}
+
+ if (seq.empty() || other.seq.empty()) {
+ return (new pseries(var==point, epvector()))
+ ->setflag(status_flags::dynallocated);
+ }
// Series multiplication
epvector new_seq;
const epvector::const_iterator itend = seq.end();
for (epvector::const_iterator it=itbeg; it!=itend; ++it) {
- ex buf = recombine_pair_to_ex(*it);
+ ex expon = it->coeff;
+ int factor = 1;
+ ex buf;
+ if (expon.info(info_flags::integer)) {
+ buf = it->rest;
+ factor = ex_to<numeric>(expon).to_int();
+ } else {
+ buf = recombine_pair_to_ex(*it);
+ }
int real_ldegree = 0;
try {
real_ldegree = buf.expand().ldegree(sym-r.rhs());
- }
- catch (std::runtime_error) {}
+ } catch (std::runtime_error) {}
if (real_ldegree == 0) {
int orderloop = 0;
} while (real_ldegree == orderloop);
}
- ldegrees.push_back(real_ldegree);
+ ldegrees.push_back(factor * real_ldegree);
}
int degsum = std::accumulate(ldegrees.begin(), ldegrees.end(), 0);
- if (degsum>order) {
+ if (degsum >= order) {
epvector epv;
epv.push_back(expair(Order(_ex1), order));
return (new pseries(r, epv))->setflag(status_flags::dynallocated);
ex op = recombine_pair_to_ex(*it).series(r, order-degsum+(*itd), options);
// Series multiplication
- if (it==itbeg)
+ if (it == itbeg)
acc = ex_to<pseries>(op);
else
acc = ex_to<pseries>(acc.mul_series(ex_to<pseries>(op)));
throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
// adjust number of coefficients
- deg = deg - p.to_int()*ldeg;
+ int numcoeff = deg - (p*ldeg).to_int();
+ if (numcoeff <= 0) {
+ epvector epv;
+ epv.reserve(1);
+ epv.push_back(expair(Order(_ex1), deg));
+ return (new pseries(relational(var,point), epv))
+ ->setflag(status_flags::dynallocated);
+ }
// O(x^n)^(-m) is undefined
if (seq.size() == 1 && is_order_function(seq[0].rest) && p.real().is_negative())
// Compute coefficients of the powered series
exvector co;
- co.reserve(deg);
+ co.reserve(numcoeff);
co.push_back(power(coeff(var, ldeg), p));
- for (int i=1; i<deg; ++i) {
+ for (int i=1; i<numcoeff; ++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<numcoeff; ++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)
- new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
+ new_seq.push_back(expair(Order(_ex1), p * ldeg + numcoeff));
return pseries(relational(var,point), new_seq);
}
}
// Is the expression of type something^(-int)?
- if (!must_expand_basis && !exponent.info(info_flags::negint))
+ if (!must_expand_basis && !exponent.info(info_flags::negint)
+ && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
return basic::series(r, order, options);
// Is the expression of type 0^something?
- if (!must_expand_basis && !basis.subs(r, subs_options::no_pattern).is_zero())
+ if (!must_expand_basis && !basis.subs(r, subs_options::no_pattern).is_zero()
+ && (!is_a<add>(basis) || !is_a<numeric>(exponent)))
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();
+ numeric numexp;
+ if (is_a<numeric>(exponent)) {
+ numexp = ex_to<numeric>(exponent);
+ } else {
+ numexp = 0;
+ }
const ex& sym = r.lhs();
// find existing minimal degree
int real_ldegree = basis.expand().ldegree(sym-r.rhs());
} while (real_ldegree == orderloop);
}
- ex e = basis.series(r, order + real_ldegree*(1-intexp), options);
+ if (!(real_ldegree*numexp).is_integer())
+ throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
+ ex e = basis.series(r, (order + real_ldegree*(1-numexp)).to_int(), options);
ex result;
try {
- result = ex_to<pseries>(e).power_const(intexp, order);
- }
- catch (pole_error) {
+ result = ex_to<pseries>(e).power_const(numexp, order);
+ } catch (pole_error) {
epvector ser;
ser.push_back(expair(Order(_ex1), order));
result = pseries(r, ser);
return convert_to_poly().series(r, order, options);
}
+ex integral::series(const relational & r, int order, unsigned options) const
+{
+ if (x.subs(r) != x)
+ throw std::logic_error("Cannot series expand wrt dummy variable");
+
+ // Expanding integrant with r substituted taken in boundaries.
+ ex fseries = f.series(r, order, options);
+ epvector fexpansion;
+ fexpansion.reserve(fseries.nops());
+ for (size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ fexpansion.push_back(expair(
+ currcoeff == Order(_ex1)
+ ? currcoeff
+ : integral(x, a.subs(r), b.subs(r), currcoeff),
+ ex_to<pseries>(fseries).exponop(i)
+ ));
+ }
+
+ // Expanding lower boundary
+ ex result = (new pseries(r, fexpansion))->setflag(status_flags::dynallocated);
+ ex aseries = (a-a.subs(r)).series(r, order, options);
+ fseries = f.series(x == (a.subs(r)), order, options);
+ for (size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ if (is_order_function(currcoeff))
+ break;
+ ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ currcoeff=currcoeff.series(r, orderforf);
+ ex term = ex_to<pseries>(aseries).power_const(ex_to<numeric>(currexpon+1),order);
+ term = ex_to<pseries>(term).mul_const(ex_to<numeric>(-1/(currexpon+1)));
+ term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ }
+
+ // Expanding upper boundary
+ ex bseries = (b-b.subs(r)).series(r, order, options);
+ fseries = f.series(x == (b.subs(r)), order, options);
+ for (size_t i=0; i<fseries.nops(); ++i) {
+ ex currcoeff = ex_to<pseries>(fseries).coeffop(i);
+ if (is_order_function(currcoeff))
+ break;
+ ex currexpon = ex_to<pseries>(fseries).exponop(i);
+ int orderforf = order-ex_to<numeric>(currexpon).to_int()-1;
+ currcoeff=currcoeff.series(r, orderforf);
+ ex term = ex_to<pseries>(bseries).power_const(ex_to<numeric>(currexpon+1),order);
+ term = ex_to<pseries>(term).mul_const(ex_to<numeric>(1/(currexpon+1)));
+ term = ex_to<pseries>(term).mul_series(ex_to<pseries>(currcoeff));
+ result = ex_to<pseries>(result).add_series(ex_to<pseries>(term));
+ }
+
+ return result;
+}
+
/** Compute the truncated series expansion of an expression.
* This function returns an expression containing an object of class pseries
git://www.ginac.de / ginac.git / blobdiff