* Implementation of class for extended truncated power-series and
* methods for series expansion. */
-#include "ginac.h"
+/*
+ * GiNaC Copyright (C) 1999 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
+ * the Free Software Foundation; either version 2 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program; if not, write to the Free Software
+ * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+ */
+#include "series.h"
+#include "add.h"
+#include "inifcns.h"
+#include "mul.h"
+#include "power.h"
+#include "relational.h"
+#include "symbol.h"
+#include "debugmsg.h"
+#include "utils.h"
+
+#ifndef NO_GINAC_NAMESPACE
+namespace GiNaC {
+#endif // ndef NO_GINAC_NAMESPACE
/*
* Default constructor, destructor, copy constructor, assignment operator and helpers
*/
-series::series() : basic(TINFO_SERIES)
+series::series() : basic(TINFO_series)
{
debugmsg("series default constructor", LOGLEVEL_CONSTRUCT);
}
/** Construct series from a vector of coefficients and powers.
* expair.rest holds the coefficient, expair.coeff holds the power.
* The powers must be integers (positive or negative) and in ascending order;
- * the last coefficient can be Order(exONE()) to represent a truncated,
+ * the last coefficient can be Order(_ex1()) to represent a truncated,
* non-terminating series.
*
* @param var_ series variable (must hold a symbol)
* @param ops_ vector of {coefficient, power} pairs (coefficient must not be zero)
* @return newly constructed series */
series::series(ex const &var_, ex const &point_, epvector const &ops_)
- : basic(TINFO_SERIES), seq(ops_), var(var_), point(point_)
+ : basic(TINFO_series), seq(ops_), var(var_), point(point_)
{
debugmsg("series constructor from ex,ex,epvector", LOGLEVEL_CONSTRUCT);
- ASSERT(is_ex_exactly_of_type(var_, symbol));
+ GINAC_ASSERT(is_ex_exactly_of_type(var_, symbol));
}
return new series(*this);
}
+void series::print(ostream &os, unsigned upper_precedence) const
+{
+ debugmsg("symbol print", LOGLEVEL_PRINT);
+ convert_to_poly().print(os, upper_precedence);
+}
+
+void series::printraw(ostream &os) const
+{
+ debugmsg("symbol printraw", LOGLEVEL_PRINT);
+ os << "series(" << var << ";" << point << ";";
+ for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) {
+ os << "(" << (*i).rest << "," << (*i).coeff << "),";
+ }
+ os << ")";
+}
+
// Highest degree of variable
int series::degree(symbol const &s) const
{
- if (var == s) {
+ if (var.is_equal(s)) {
// Return last exponent
if (seq.size())
return ex_to_numeric((*(seq.end() - 1)).coeff).to_int();
// Lowest degree of variable
int series::ldegree(symbol const &s) const
{
- if (var == s) {
+ if (var.is_equal(s)) {
// Return first exponent
if (seq.size())
return ex_to_numeric((*(seq.begin())).coeff).to_int();
}
// Coefficient of variable
-ex series::coeff(symbol const &s, int n) const
+ex series::coeff(symbol const &s, int const n) const
{
- if (var == s) {
+ if (var.is_equal(s)) {
epvector::const_iterator it = seq.begin(), itend = seq.end();
while (it != itend) {
int pow = ex_to_numeric(it->coeff).to_int();
if (pow == n)
return it->rest;
if (pow > n)
- return exZERO();
+ return _ex0();
it++;
}
- return exZERO();
+ return _ex0();
} else
return convert_to_poly().coeff(s, n);
}
{
if (level == 1)
return this->hold();
-
+
// Construct a new series with evaluated coefficients
epvector new_seq;
new_seq.reserve(seq.size());
return (new series(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
}
+/** Evaluate numerically. The order term is dropped. */
ex series::evalf(int level) const
{
return convert_to_poly().evalf(level);
}
-
/*
* Construct expression (polynomial) out of series
*/
{
ex e;
epvector::const_iterator it = seq.begin(), itend = seq.end();
-
+
while (it != itend) {
if (is_order_function(it->rest)) {
if (!no_order)
ex coeff = deriv.subs(s == point);
if (!coeff.is_zero())
seq.push_back(expair(coeff, numeric(0)));
-
+
int n;
for (n=1; n<order; n++) {
fac = fac.mul(numeric(n));
if (!coeff.is_zero())
seq.push_back(expair(coeff, numeric(n)));
}
-
+
// Higher-order terms, if present
deriv = deriv.diff(s);
if (!deriv.is_zero())
- seq.push_back(expair(Order(exONE()), numeric(n)));
+ seq.push_back(expair(Order(_ex1()), numeric(n)));
return series::series(s, point, seq);
}
// results in an empty (constant) series
if (!is_compatible_to(other)) {
epvector nul;
- nul.push_back(expair(Order(exONE()), exZERO()));
+ nul.push_back(expair(Order(_ex1()), _ex0()));
return series(var, point, nul);
}
-
+
// Series addition
epvector new_seq;
epvector::const_iterator a = seq.begin();
break;
} else
pow_a = ex_to_numeric((*a).coeff).to_int();
-
+
// If b is empty, fill up with elements from a and stop
if (b == b_end) {
while (a != a_end) {
break;
} else
pow_b = ex_to_numeric((*b).coeff).to_int();
-
+
// a and b are non-empty, compare powers
if (pow_a < pow_b) {
// a has lesser power, get coefficient from a
} else {
// Add coefficient of a and b
if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
- new_seq.push_back(expair(Order(exONE()), (*a).coeff));
+ new_seq.push_back(expair(Order(_ex1()), (*a).coeff));
break; // Order term ends the sequence
} else {
ex sum = (*a).rest + (*b).rest;
- if (!(sum == exZERO()))
+ if (!(sum.is_zero()))
new_seq.push_back(expair(sum, numeric(pow_a)));
a++;
b++;
acc = it->rest;
else
acc = it->rest.series(s, point, order);
- if (!it->coeff.is_equal(exONE()))
+ if (!it->coeff.is_equal(_ex1()))
acc = ex_to_series(acc).mul_const(ex_to_numeric(it->coeff));
it++;
}
op = it->rest;
else
op = it->rest.series(s, point, order);
- if (!it->coeff.is_equal(exONE()))
+ if (!it->coeff.is_equal(_ex1()))
op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
// Series addition
ex add::series(symbol const & s, ex const & point, int order) const
{
ex acc; // Series accumulator
-
+
// Get first term from overall_coeff
acc = overall_coeff.series(s,point,order);
-
+
// Add remaining terms
epvector::const_iterator it = seq.begin();
epvector::const_iterator itend = seq.end();
op = it->rest;
else
op = it->rest.series(s, point, order);
- if (!it->coeff.is_equal(exONE()))
+ if (!it->coeff.is_equal(_ex1()))
op = ex_to_series(op).mul_const(ex_to_numeric(it->coeff));
-
+
// Series addition
acc = ex_to_series(acc).add_series(ex_to_series(op));
}
{
epvector new_seq;
new_seq.reserve(seq.size());
-
+
epvector::const_iterator it = seq.begin(), itend = seq.end();
while (it != itend) {
if (!is_order_function(it->rest))
// results in an empty (constant) series
if (!is_compatible_to(other)) {
epvector nul;
- nul.push_back(expair(Order(exONE()), exZERO()));
+ nul.push_back(expair(Order(_ex1()), _ex0()));
return series(var, point, nul);
}
// Series multiplication
epvector new_seq;
-
+
const symbol *s = static_cast<symbol *>(var.bp);
int a_max = degree(*s);
int b_max = other.degree(*s);
int b_min = other.ldegree(*s);
int cdeg_min = a_min + b_min;
int cdeg_max = a_max + b_max;
-
+
int higher_order_a = INT_MAX;
int higher_order_b = INT_MAX;
if (is_order_function(coeff(*s, a_max)))
int higher_order_c = min(higher_order_a, higher_order_b);
if (cdeg_max >= higher_order_c)
cdeg_max = higher_order_c - 1;
-
+
for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) {
- ex co = exZERO();
+ ex co = _ex0();
// c(i)=a(0)b(i)+...+a(i)b(0)
for (int i=a_min; cdeg-i>=b_min; i++) {
ex a_coeff = coeff(*s, i);
new_seq.push_back(expair(co, numeric(cdeg)));
}
if (higher_order_c < INT_MAX)
- new_seq.push_back(expair(Order(exONE()), numeric(higher_order_c)));
+ new_seq.push_back(expair(Order(_ex1()), numeric(higher_order_c)));
return series::series(var, point, new_seq);
}
-/** Implementation of ex::series() for product. This performs series multiplication when multiplying series.
- * @see ex::series */
/*
ex mul::series(symbol const & s, ex const & point, int order) const
{
acc = it->rest;
else
acc = it->rest.series(s, point, order);
- if (!it->coeff.is_equal(exONE()))
+ if (!it->coeff.is_equal(_ex1()))
acc = ex_to_series(acc).power_const(ex_to_numeric(it->coeff), order);
it++;
}
continue;
} else if (!is_ex_exactly_of_type(op, series))
op = op.series(s, point, order);
- if (!it->coeff.is_equal(exONE()))
+ if (!it->coeff.is_equal(_ex1()))
op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
// Series multiplication
}
*/
+/** Implementation of ex::series() for product. This performs series
+ * multiplication when multiplying series.
+ * @see ex::series */
ex mul::series(symbol const & s, ex const & point, int order) const
{
ex acc; // Series accumulator
-
+
// Get first term from overall_coeff
acc = overall_coeff.series(s, point, order);
-
+
// Multiply with remaining terms
epvector::const_iterator it = seq.begin();
epvector::const_iterator itend = seq.end();
continue;
} else if (!is_ex_exactly_of_type(op, series))
op = op.series(s, point, order);
- if (!it->coeff.is_equal(exONE()))
+ if (!it->coeff.is_equal(_ex1()))
op = ex_to_series(op).power_const(ex_to_numeric(it->coeff), order);
// Series multiplication
int i;
const symbol *s = static_cast<symbol *>(var.bp);
int ldeg = ldegree(*s);
-
+
// Calculate coefficients of powered series
exvector co;
co.reserve(deg);
co.push_back(co0 = power(coeff(*s, ldeg), p));
bool all_sums_zero = true;
for (i=1; i<deg; i++) {
- ex sum = exZERO();
+ ex sum = _ex0();
for (int j=1; j<=i; j++) {
ex c = coeff(*s, j + ldeg);
if (is_order_function(c)) {
- co.push_back(Order(exONE()));
+ co.push_back(Order(_ex1()));
break;
} else
sum += (p * j - (i - j)) * co[i - j] * c;
all_sums_zero = false;
co.push_back(co0 * sum / numeric(i));
}
-
+
// Construct new series (of non-zero coefficients)
epvector new_seq;
bool higher_order = false;
}
}
if (!higher_order && !all_sums_zero)
- new_seq.push_back(expair(Order(exONE()), numeric(deg) + p * ldeg));
+ new_seq.push_back(expair(Order(_ex1()), numeric(deg) + p * ldeg));
return series::series(var, point, new_seq);
}
// Basis is not a series, may there be a singulary?
if (!exponent.info(info_flags::negint))
return basic::series(s, point, order);
-
+
// Expression is of type something^(-int), check for singularity
if (!basis.subs(s == point).is_zero())
return basic::series(s, point, order);
-
+
// Singularity encountered, expand basis into series
e = basis.series(s, point, order);
} else {
// Basis is a series
e = basis;
}
-
+
// Power e
return ex_to_series(e).power_const(ex_to_numeric(exponent), order);
}
* @return an expression holding a series object */
ex ex::series(symbol const &s, ex const &point, int order) const
{
- ASSERT(bp!=0);
+ GINAC_ASSERT(bp!=0);
return bp->series(s, point, order);
}
// Global constants
const series some_series;
type_info const & typeid_series = typeid(some_series);
+
+#ifndef NO_GINAC_NAMESPACE
+} // namespace GiNaC
+#endif // ndef NO_GINAC_NAMESPACE
git://www.ginac.de / ginac.git / blobdiff