#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
/** 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)
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
{
}
// Coefficient of variable
-ex series::coeff(symbol const &s, int n) const
+ex series::coeff(symbol const &s, int const n) const
{
if (var.is_equal(s)) {
epvector::const_iterator it = seq.begin(), itend = seq.end();
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;
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);
}
const series some_series;
type_info const & typeid_series = typeid(some_series);
+#ifndef NO_GINAC_NAMESPACE
} // namespace GiNaC
+#endif // ndef NO_GINAC_NAMESPACE