* methods for series expansion. */
/*
- * GiNaC Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+ * GiNaC Copyright (C) 1999-2005 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
#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;
throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
// adjust number of coefficients
- deg = deg - (p*ldeg).to_int();
+ 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) && !is_a<add>(basis))
+ 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() && !is_a<add>(basis))
+ 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)?
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);
+ currcoeff = (currcoeff == Order(_ex1))
+ ? currcoeff
+ : integral(x, a.subs(r), b.subs(r), currcoeff);
+ if (currcoeff != 0)
+ fexpansion.push_back(
+ expair(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