This file records noteworthy changes.
0.6.0 (...)
-* Important interface changes:
- gamma() -> Gamma()
- EulerGamma -> gamma
- beta() -> Beta()
+* IMPORTANT: Several interface changes make programs written with GiNaC
+ much clearer but break compatibility with older versions:
+ - f(x).series(x,p[,o]) -> f(x).series(x==p,o)
+ - series(f(x),x,p[,o]) -> series(f(x),x==p,o)
+ - gamma() -> Gamma()
+ - EulerGamma -> gamma
+ - beta() -> Beta()
0.5.4 (15 March 2000)
* Some algorithms in class matrix (notably determinant) were replaced by
symbol x("x");
ex e, d, ed;
- e = sin(x).series(x, 0, 8);
- d = cos(x).series(x, 0, 7);
+ e = sin(x).series(x==0, 8);
+ d = cos(x).series(x==0, 7);
ed = e.diff(x);
ed = series_to_poly(ed);
d = series_to_poly(d);
static unsigned check_series(const ex &e, const ex &point, const ex &d, int order = 8)
{
- ex es = e.series(x, point, order);
+ ex es = e.series(x==point, order);
ex ep = ex_to_pseries(es).convert_to_poly();
if (!(ep - d).is_zero()) {
clog << "series expansion of " << e << " at " << point
unsigned result = 0;
ex e, d;
- e = pow(sin(x), -1).series(x, 0, 8) + pow(sin(-x), -1).series(x, 0, 12);
+ e = pow(sin(x), -1).series(x==0, 8) + pow(sin(-x), -1).series(x==0, 12);
d = Order(pow(x, 6));
result += check_series(e, 0, d);
unsigned result = 0;
ex e, d;
- e = sin(x).series(x, 0, 8) * pow(sin(x), -1).series(x, 0, 12);
+ e = sin(x).series(x==0, 8) * pow(sin(x), -1).series(x==0, 12);
d = 1 + Order(pow(x, 7));
result += check_series(e, 0, d);
unsigned result = 0;
symbol x;
- ex myseries = series(Gamma(x),x,0,order);
+ ex myseries = series(Gamma(x),x==0,order);
// compute the last coefficient numerically:
ex last_coeff = myseries.coeff(x,order-1).evalf();
// compute a bound for that coefficient using a variation of the leading
inline ex diff(basic const & x, symbol const & y, int z=1) {
return diff(ex(x),(y),z);
}
-inline ex series(basic const & x, symbol const & y, ex const & z, int zz=6) {
- return series(ex(x),(y),(z),zz);
+inline ex series(const basic & x, const relational & y, int z) {
+ return series(ex(x),ex(y),(z));
}
-inline ex series(ex const & x, symbol const & y, basic const & z, int zz=6) {
- return series(ex(x),(y),ex(z),zz);
-}
-inline ex series(basic const & x, symbol const & y, basic const & z, int zz=6) {
- return series(ex(x),(y),ex(z),zz);
+inline ex series(const basic & x, const symbol & y, int z) {
+ return series(ex(x),ex(y),(z));
}
// fixes for ex subs(x,y)
inline ex subs(ex const & x, basic const & y) {
inline_single_function_2p_with_defarg('ex','eval','int','0');
inline_single_function_2p_with_defarg('ex','evalf','int','0');
inline_single_function_3p_with_defarg('ex','diff','basic const &','ex','symbol const &','','int','1');
-inline_single_function_4p_with_defarg('ex','series','basic const &','ex','symbol const &','','ex const &','','int','6');
-inline_single_function_4p_with_defarg('ex','series','ex const &','ex','symbol const &','','basic const &','ex','int','6');
-inline_single_function_4p_with_defarg('ex','series','basic const &','ex','symbol const &','','basic const &','ex','int','6');
+inline_single_function_3p('ex','series','const basic &','ex','const relational &','ex','int','');
+inline_single_function_3p('ex','series','const basic &','ex','const symbol &','ex','int','');
inline_function_2p('ex','subs');
inline_single_function_3p('ex','subs','basic const &','ex','lst const &','','lst const &','');
inline_single_function_2p('ex','op','basic const &','ex','int','');
@end example
You can differentiate functions and expand them as Taylor or Laurent
-series (the third argument of @code{series} is the evaluation point, the
-fourth defines the order):
+series in a very natural syntax (the second argument of @code{series} is
+a relation defining the evaluation point, the third specifies the
+order):
@cindex Zeta function
@example
> diff(tan(x),x);
tan(x)^2+1
-> series(sin(x),x,0,4);
+> series(sin(x),x==0,4);
x-1/6*x^3+Order(x^4)
-> series(1/tan(x),x,0,4);
+> series(1/tan(x),x==0,4);
x^(-1)-1/3*x+Order(x^2)
-> series(Gamma(x),x,0,3);
+> series(Gamma(x),x==0,3);
x^(-1)-gamma+(1/12*Pi^2+1/2*gamma^2)*x+
(-1/3*zeta(3)-1/12*Pi^2*gamma-1/6*gamma^3)*x^2+Order(x^3)
> evalf(");
x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
-(0.90747907608088628905)*x^2+Order(x^3)
-> series(Gamma(2*sin(x)-2),x,Pi/2,6);
+> series(Gamma(2*sin(x)-2),x==Pi/2,6);
-(x-1/2*Pi)^(-2)+(-1/12*Pi^2-1/2*gamma^2-1/240)*(x-1/2*Pi)^2
-gamma-1/12+Order((x-1/2*Pi)^3)
@end example
the @code{.subs()} method show how objects of class relational are used
as arguments. There they provide an intuitive syntax for substitutions.
They can also used for creating systems of equations that are to be
-solved for unknown variables.
+solved for unknown variables. More applications of this class will
+appear throughout the next chapters.
@node Archiving, Important Algorithms, Relations, Basic Concepts
symbol v("v"), c("c");
ex gamma = 1/sqrt(1 - pow(v/c,2));
- ex mass_nonrel = gamma.series(v, 0, 10);
+ ex mass_nonrel = gamma.series(v==0, 10);
cout << "the relativistic mass increase with v is " << endl
<< mass_nonrel << endl;
cout << "the inverse square of this series is " << endl
- << pow(mass_nonrel,-2).series(v, 0, 10) << endl;
+ << pow(mass_nonrel,-2).series(v==0, 10) << endl;
// ...
@}
ex mechain_pi(int degr)
@{
symbol x;
- ex pi_expansion = series_to_poly(atan(x).series(x,0,degr));
+ ex pi_expansion = series_to_poly(atan(x).series(x,degr));
ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
-4*pi_expansion.subs(x==numeric(1,239));
return pi_approx;
@}
@end example
-When you run this program, it will type out:
+Note how we just called @code{.series(x,degr)} instead of
+@code{.series(x==0,degr)}. This is a simple shortcut for @code{ex}'s
+method @code{series()}: if the first argument is a symbol the expression
+is expanded in that symbol around point @code{0}. When you run this
+program, it will type out:
@example
2: 3804/1195
# Files which are generated by perl scripts
$(srcdir)/function.h $(srcdir)/function.cpp: $(srcdir)/function.pl
- cd $(srcdir) && perl function.pl
+ cd $(srcdir) && perl -w function.pl
$(srcdir)/lst.h $(srcdir)/lst.cpp: $(srcdir)/container.pl
- cd $(srcdir) && perl container.pl lst
+ cd $(srcdir) && perl -w container.pl lst
$(srcdir)/exprseq.h $(srcdir)/exprseq.cpp: $(srcdir)/container.pl
- cd $(srcdir) && perl container.pl exprseq
+ cd $(srcdir) && perl -w container.pl exprseq
# Force build of headers before compilation
$(srcdir)/add.cpp: $(srcdir)/function.h $(srcdir)/lst.h $(srcdir)/exprseq.h
int ldegree(const symbol & s) const;
ex coeff(const symbol & s, int n=1) const;
ex eval(int level=0) const;
- ex series(const symbol & s, const ex & point, int order) const;
+ ex series(const relational & r, int order) const;
ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
numeric integer_content(void) const;
ex smod(const numeric &xi) const;
throw(std::logic_error("differentiation not supported by this type"));
}
+/** Returns order relation between two objects of same type. Needs to be
+ * implemented by each class. */
int basic::compare_same_type(const basic & other) const
{
return compare_pointers(this, &other);
}
+/** Returns true if two objects of same type are equal. Normally needs
+ * not be reimplemented as long as it wasn't overwritten by some parent
+ * class, since it just calls complare_same_type(). */
bool basic::is_equal_same_type(const basic & other) const
{
return compare_same_type(other)==0;
class symbol;
class lst;
class numeric;
+class relational;
class archive_node;
//typedef vector<ex> exvector;
virtual ex collect(const symbol & s) const;
virtual ex eval(int level=0) const;
virtual ex evalf(int level=0) const;
- virtual ex series(const symbol & s, const ex & point, int order) const;
+ virtual ex series(const relational & r, int order) const;
virtual ex subs(const lst & ls, const lst & lr) const;
virtual ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
virtual numeric integer_content(void) const;
-#!/usr/bin/perl -w
-
if (($#ARGV!=0) and ($#ARGV!=1)) {
die 'usage: container.pl type [maxargs] (type=lst or exprseq)';
}
ex eval(int level = 0) const;
ex evalf(int level = 0) const;
ex diff(const symbol & s, unsigned nth = 1) const;
- ex series(const symbol & s, const ex & point, int order = 6) const;
+ ex series(const ex & r, int order) const;
#ifdef CINT_CONVERSION_WORKAROUND
- ex series(const symbol & s, const basic & point, int order = 6) const { return series(s,ex(point),order); }
+ ex series(const basic & r, int order) const { return series(ex(r),order); }
#endif // def CINT_CONVERSION_WORKAROUND
ex subs(const lst & ls, const lst & lr) const;
ex subs(const ex & e) const;
inline ex diff(const ex & thisex, const symbol & s, unsigned nth = 1)
{ return thisex.diff(s, nth); }
-inline ex series(const ex & thisex, const symbol & s, const ex & point, int order = 6)
-{ return thisex.series(s, point, order); }
+inline ex series(const ex & thisex, const ex & r, int order)
+{ return thisex.series(r, order); }
inline ex subs(const ex & thisex, const ex & e)
{ return thisex.subs(e); }
-#!/usr/bin/perl -w
-
$maxargs=13;
sub generate_seq {
'const ex &','');
$typedef_series_funcp=generate(
-'typedef ex (* series_funcp_${N})(${SEQ1}, const symbol &, const ex &, int);'."\n",
+'typedef ex (* series_funcp_${N})(${SEQ1}, const relational &, int);'."\n",
'const ex &','');
$eval_func_interface=generate(' function_options & eval_func(eval_funcp_${N} e);'."\n",'','');
<<'END_OF_SERIES_SWITCH_STATEMENT','seq[${N}-1]','');
case ${N}:
try {
- res = ((series_funcp_${N})(registered_functions()[serial].series_f))(${SEQ1},s,point,order);
+ res = ((series_funcp_${N})(registered_functions()[serial].series_f))(${SEQ1},r,order);
} catch (do_taylor) {
- res = basic::series(s, point, order);
+ res = basic::series(r, order);
}
return res;
break;
ex expand(unsigned options=0) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const symbol & s, const ex & point, int order) const;
+ ex series(const relational & r, int order) const;
ex thisexprseq(const exvector & v) const;
ex thisexprseq(exvector * vp) const;
protected:
/** Implementation of ex::series for functions.
* \@see ex::series */
-ex function::series(const symbol & s, const ex & point, int order) const
+ex function::series(const relational & r, int order) const
{
GINAC_ASSERT(serial<registered_functions().size());
if (registered_functions()[serial].series_f==0) {
- return basic::series(s, point, order);
+ return basic::series(r, order);
}
ex res;
switch (registered_functions()[serial].nparams) {
return Order(x).hold();
}
-static ex Order_series(const ex & x, const symbol & s, const ex & point, int order)
+static ex Order_series(const ex & x, const relational & r, int order)
{
// Just wrap the function into a pseries object
epvector new_seq;
- new_seq.push_back(expair(Order(_ex1()), numeric(min(x.ldegree(s), order))));
- return pseries(s, point, new_seq);
+ GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
+ const symbol *s = static_cast<symbol *>(r.lhs().bp);
+ new_seq.push_back(expair(Order(_ex1()), numeric(min(x.ldegree(*s), order))));
+ return pseries(r, new_seq);
}
// Differentiation is handled in function::derivative because of its special requirements
}
-static ex Gamma_series(const ex & x, const symbol & s, const ex & pt, int order)
+static ex Gamma_series(const ex & x, const relational & r, int order)
{
// method:
// Taylor series where there is no pole falls back to psi function
// from which follows
// series(Gamma(x),x,-m,order) ==
// series(Gamma(x+m+1)/(x*(x+1)*...*(x+m)),x,-m,order+1);
- const ex x_pt = x.subs(s==pt);
+ const ex x_pt = x.subs(r);
if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive))
throw do_taylor(); // caught by function::series()
// if we got here we have to care for a simple pole at -m:
ex ser_denom = _ex1();
for (numeric p; p<=m; ++p)
ser_denom *= x+p;
- return (Gamma(x+m+_ex1())/ser_denom).series(s, pt, order+1);
+ return (Gamma(x+m+_ex1())/ser_denom).series(r, order+1);
}
}
-static ex Beta_series(const ex & x, const ex & y, const symbol & s, const ex & pt, int order)
+static ex Beta_series(const ex & x, const ex & y, const relational & r, int order)
{
// method:
// Taylor series where there is no pole of one of the Gamma functions
// falls back to Beta function evaluation. Otherwise, fall back to
// Gamma series directly.
// FIXME: this could need some testing, maybe it's wrong in some cases?
- const ex x_pt = x.subs(s==pt);
- const ex y_pt = y.subs(s==pt);
+ const ex x_pt = x.subs(r);
+ const ex y_pt = y.subs(r);
+ GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
+ const symbol *s = static_cast<symbol *>(r.lhs().bp);
ex x_ser, y_ser, xy_ser;
if ((!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive)) &&
(!y_pt.info(info_flags::integer) || y_pt.info(info_flags::positive)))
throw do_taylor(); // caught by function::series()
// trap the case where x is on a pole directly:
if (x.info(info_flags::integer) && !x.info(info_flags::positive))
- x_ser = Gamma(x+s).series(s,pt,order);
+ x_ser = Gamma(x+*s).series(r,order);
else
- x_ser = Gamma(x).series(s,pt,order);
+ x_ser = Gamma(x).series(r,order);
// trap the case where y is on a pole directly:
if (y.info(info_flags::integer) && !y.info(info_flags::positive))
- y_ser = Gamma(y+s).series(s,pt,order);
+ y_ser = Gamma(y+*s).series(r,order);
else
- y_ser = Gamma(y).series(s,pt,order);
+ y_ser = Gamma(y).series(r,order);
// trap the case where y is on a pole directly:
if ((x+y).info(info_flags::integer) && !(x+y).info(info_flags::positive))
- xy_ser = Gamma(y+x+s).series(s,pt,order);
+ xy_ser = Gamma(y+x+*s).series(r,order);
else
- xy_ser = Gamma(y+x).series(s,pt,order);
+ xy_ser = Gamma(y+x).series(r,order);
// compose the result:
- return (x_ser*y_ser/xy_ser).series(s,pt,order);
+ return (x_ser*y_ser/xy_ser).series(r,order);
}
return psi(_ex1(), x);
}
-static ex psi1_series(const ex & x, const symbol & s, const ex & pt, int order)
+static ex psi1_series(const ex & x, const relational & r, int order)
{
// method:
// Taylor series where there is no pole falls back to polygamma function
// from which follows
// series(psi(x),x,-m,order) ==
// series(psi(x+m+1) - 1/x - 1/(x+1) - 1/(x+m)),x,-m,order);
- const ex x_pt = x.subs(s==pt);
+ const ex x_pt = x.subs(r);
if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive))
throw do_taylor(); // caught by function::series()
// if we got here we have to care for a simple pole at -m:
ex recur;
for (numeric p; p<=m; ++p)
recur += power(x+p,_ex_1());
- return (psi(x+m+_ex1())-recur).series(s, pt, order);
+ return (psi(x+m+_ex1())-recur).series(r, order);
}
const unsigned function_index_psi1 =
return psi(n+_ex1(), x);
}
-static ex psi2_series(const ex & n, const ex & x, const symbol & s, const ex & pt, int order)
+static ex psi2_series(const ex & n, const ex & x, const relational & r, int order)
{
// method:
// Taylor series where there is no pole falls back to polygamma function
// series(psi(x),x,-m,order) ==
// series(psi(x+m+1) - (-1)^n * n! * ((x)^(-n-1) + (x+1)^(-n-1) + ...
// ... + (x+m)^(-n-1))),x,-m,order);
- const ex x_pt = x.subs(s==pt);
+ const ex x_pt = x.subs(r);
if (!x_pt.info(info_flags::integer) || x_pt.info(info_flags::positive))
throw do_taylor(); // caught by function::series()
// if we got here we have to care for a pole of order n+1 at -m:
for (numeric p; p<=m; ++p)
recur += power(x+p,-n+_ex_1());
recur *= factorial(n)*power(_ex_1(),n);
- return (psi(n, x+m+_ex1())-recur).series(s, pt, order);
+ return (psi(n, x+m+_ex1())-recur).series(r, order);
}
const unsigned function_index_psi2 =
return power(x, _ex_1());
}
-static ex log_series(const ex & x, const symbol & s, const ex & pt, int order)
+static ex log_series(const ex &x, const relational &r, int order)
{
- if (x.subs(s == pt).is_zero()) {
+ if (x.subs(r).is_zero()) {
epvector seq;
seq.push_back(expair(log(x), _ex0()));
- return pseries(s, pt, seq);
+ return pseries(r, seq);
} else
throw do_taylor();
}
return (_ex1()+power(tan(x),_ex2()));
}
-static ex tan_series(const ex & x, const symbol & s, const ex & pt, int order)
+static ex tan_series(const ex &x, const relational &r, int order)
{
// method:
// Taylor series where there is no pole falls back to tan_deriv.
// On a pole simply expand sin(x)/cos(x).
- const ex x_pt = x.subs(s==pt);
+ const ex x_pt = x.subs(r);
if (!(2*x_pt/Pi).info(info_flags::odd))
throw do_taylor(); // caught by function::series()
// if we got here we have to care for a simple pole
- return (sin(x)/cos(x)).series(s, pt, order+2);
+ return (sin(x)/cos(x)).series(r, order+2);
}
REGISTER_FUNCTION(tan, eval_func(tan_eval).
return _ex1()-power(tanh(x),_ex2());
}
-static ex tanh_series(const ex & x, const symbol & s, const ex & pt, int order)
+static ex tanh_series(const ex &x, const relational &r, int order)
{
// method:
// Taylor series where there is no pole falls back to tanh_deriv.
// On a pole simply expand sinh(x)/cosh(x).
- const ex x_pt = x.subs(s==pt);
+ const ex x_pt = x.subs(r);
if (!(2*I*x_pt/Pi).info(info_flags::odd))
throw do_taylor(); // caught by function::series()
// if we got here we have to care for a simple pole
- return (sinh(x)/cosh(x)).series(s, pt, order+2);
+ return (sinh(x)/cosh(x)).series(r, order+2);
}
REGISTER_FUNCTION(tanh, eval_func(tanh_eval).
ex coeff(const symbol & s, int n=1) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const symbol & s, const ex & point, int order) const;
+ ex series(const relational & s, int order) const;
ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
numeric integer_content(void) const;
ex smod(const numeric &xi) const;
new_seq.push_back(expair(it->rest.normal(), it->coeff));
it++;
}
- ex n = pseries(var, point, new_seq);
+ ex n = pseries(relational(var,point), new_seq);
return (new lst(replace_with_symbol(n, sym_lst, repl_lst), _ex1()))->setflag(status_flags::dynallocated);
}
ex coeff(const symbol & s, int n=1) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const symbol & s, const ex & point, int order) const;
+ ex series(const relational & s, int order) const;
ex subs(const lst & ls, const lst & lr) const;
ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
ex simplify_ncmul(const exvector & v) const;
* the last coefficient can be Order(_ex1()) to represent a truncated,
* non-terminating series.
*
- * @param var_ series variable (must hold a symbol)
- * @param point_ expansion point
+ * @param rel__ expansion variable and point (must hold a relational)
* @param ops_ vector of {coefficient, power} pairs (coefficient must not be zero)
* @return newly constructed pseries */
-pseries::pseries(const ex &var_, const ex &point_, const epvector &ops_)
- : basic(TINFO_pseries), seq(ops_), var(var_), point(point_)
+pseries::pseries(const ex &rel_, const epvector &ops_)
+ : basic(TINFO_pseries), seq(ops_)
{
- debugmsg("pseries constructor from ex,ex,epvector", LOGLEVEL_CONSTRUCT);
- GINAC_ASSERT(is_ex_exactly_of_type(var_, symbol));
+ debugmsg("pseries constructor from rel,epvector", LOGLEVEL_CONSTRUCT);
+ GINAC_ASSERT(is_ex_exactly_of_type(rel_, relational));
+ GINAC_ASSERT(is_ex_exactly_of_type(rel_.lhs(),symbol));
+ point = rel_.rhs();
+ var = *static_cast<symbol *>(rel_.lhs().bp);
}
new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
it++;
}
- return (new pseries(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
+ return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
}
/** Evaluate coefficients numerically. */
new_seq.push_back(expair(it->rest.evalf(level-1), it->coeff));
it++;
}
- return (new pseries(var, point, new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
+ return (new pseries(relational(var,point), new_seq))->setflag(status_flags::dynallocated | status_flags::evaluated);
}
ex pseries::subs(const lst & ls, const lst & lr) const
new_seq.push_back(expair(it->rest.subs(ls, lr), it->coeff));
it++;
}
- return (new pseries(var, point.subs(ls, lr), new_seq))->setflag(status_flags::dynallocated);
+ return (new pseries(relational(var,point.subs(ls, lr)), new_seq))->setflag(status_flags::dynallocated);
}
/** Implementation of ex::diff() for a power series. It treats the series as a
}
it++;
}
- return pseries(var, point, new_seq);
+ return pseries(relational(var,point), new_seq);
} else {
return *this;
}
/** Default implementation of ex::series(). This performs Taylor expansion.
* @see ex::series */
-ex basic::series(const symbol & s, const ex & point, int order) const
+ex basic::series(const relational & r, int order) const
{
epvector seq;
numeric fac(1);
ex deriv = *this;
- ex coeff = deriv.subs(s == point);
+ ex coeff = deriv.subs(r);
+ const symbol *s = static_cast<symbol *>(r.lhs().bp);
+
if (!coeff.is_zero())
seq.push_back(expair(coeff, numeric(0)));
int n;
for (n=1; n<order; n++) {
fac = fac.mul(numeric(n));
- deriv = deriv.diff(s).expand();
+ deriv = deriv.diff(*s).expand();
if (deriv.is_zero()) {
// Series terminates
- return pseries(s, point, seq);
+ return pseries(r, seq);
}
- coeff = fac.inverse() * deriv.subs(s == point);
+ coeff = fac.inverse() * deriv.subs(r);
if (!coeff.is_zero())
seq.push_back(expair(coeff, numeric(n)));
}
// Higher-order terms, if present
- deriv = deriv.diff(s);
+ deriv = deriv.diff(*s);
if (!deriv.is_zero())
seq.push_back(expair(Order(_ex1()), numeric(n)));
- return pseries(s, point, seq);
+ return pseries(r, seq);
}
/** Implementation of ex::series() for symbols.
* @see ex::series */
-ex symbol::series(const symbol & s, const ex & point, int order) const
+ex symbol::series(const relational & r, int order) const
{
epvector seq;
- if (is_equal(s)) {
+ const ex point = r.rhs();
+ GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
+ const symbol *s = static_cast<symbol *>(r.lhs().bp);
+
+ if (this->is_equal(*s)) {
if (order > 0 && !point.is_zero())
seq.push_back(expair(point, _ex0()));
if (order > 1)
seq.push_back(expair(Order(_ex1()), numeric(order)));
} else
seq.push_back(expair(*this, _ex0()));
- return pseries(s, point, seq);
+ return pseries(r, seq);
}
if (!is_compatible_to(other)) {
epvector nul;
nul.push_back(expair(Order(_ex1()), _ex0()));
- return pseries(var, point, nul);
+ return pseries(relational(var,point), nul);
}
// Series addition
}
}
}
- return pseries(var, point, new_seq);
+ return pseries(relational(var,point), new_seq);
}
/** Implementation of ex::series() for sums. This performs series addition when
* adding pseries objects.
* @see ex::series */
-ex add::series(const symbol & s, const ex & point, int order) const
+ex add::series(const relational & r, int order) const
{
ex acc; // Series accumulator
// Get first term from overall_coeff
- acc = overall_coeff.series(s, point, order);
+ acc = overall_coeff.series(r, order);
// Add remaining terms
epvector::const_iterator it = seq.begin();
if (is_ex_exactly_of_type(it->rest, pseries))
op = it->rest;
else
- op = it->rest.series(s, point, order);
+ op = it->rest.series(r, order);
if (!it->coeff.is_equal(_ex1()))
op = ex_to_pseries(op).mul_const(ex_to_numeric(it->coeff));
new_seq.push_back(*it);
it++;
}
- return pseries(var, point, new_seq);
+ return pseries(relational(var,point), new_seq);
}
if (!is_compatible_to(other)) {
epvector nul;
nul.push_back(expair(Order(_ex1()), _ex0()));
- return pseries(var, point, nul);
+ return pseries(relational(var,point), nul);
}
// Series multiplication
}
if (higher_order_c < INT_MAX)
new_seq.push_back(expair(Order(_ex1()), numeric(higher_order_c)));
- return pseries(var, point, new_seq);
+ return pseries(relational(var,point), new_seq);
}
/** Implementation of ex::series() for product. This performs series
* multiplication when multiplying series.
* @see ex::series */
-ex mul::series(const symbol & s, const ex & point, int order) const
+ex mul::series(const relational & r, int order) const
{
ex acc; // Series accumulator
// Get first term from overall_coeff
- acc = overall_coeff.series(s, point, order);
+ acc = overall_coeff.series(r, order);
// Multiply with remaining terms
epvector::const_iterator it = seq.begin();
acc = ex_to_pseries(acc).mul_const(ex_to_numeric(f));
continue;
} else if (!is_ex_exactly_of_type(op, pseries))
- op = op.series(s, point, order);
+ op = op.series(r, order);
if (!it->coeff.is_equal(_ex1()))
op = ex_to_pseries(op).power_const(ex_to_numeric(it->coeff), order);
}
if (!higher_order && !all_sums_zero)
new_seq.push_back(expair(Order(_ex1()), numeric(deg) + p * ldeg));
- return pseries(var, point, new_seq);
+ return pseries(relational(var,point), new_seq);
}
/** Implementation of ex::series() for powers. This performs Laurent expansion
* of reciprocals of series at singularities.
* @see ex::series */
-ex power::series(const symbol & s, const ex & point, int order) const
+ex power::series(const relational & r, int order) const
{
ex e;
if (!is_ex_exactly_of_type(basis, pseries)) {
// Basis is not a series, may there be a singulary?
if (!exponent.info(info_flags::negint))
- return basic::series(s, point, order);
+ return basic::series(r, order);
// Expression is of type something^(-int), check for singularity
- if (!basis.subs(s == point).is_zero())
- return basic::series(s, point, order);
+ if (!basis.subs(r).is_zero())
+ return basic::series(r, order);
// Singularity encountered, expand basis into series
- e = basis.series(s, point, order);
+ e = basis.series(r, order);
} else {
// Basis is a series
e = basis;
/** Re-expansion of a pseries object. */
-ex pseries::series(const symbol & s, const ex & p, int order) const
+ex pseries::series(const relational & r, int order) const
{
- if (var.is_equal(s) && point.is_equal(p)) {
- if (order > degree(s))
+ const ex p = r.rhs();
+ GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
+ const symbol *s = static_cast<symbol *>(r.lhs().bp);
+
+ if (var.is_equal(*s) && point.is_equal(p)) {
+ if (order > degree(*s))
return *this;
else {
epvector new_seq;
new_seq.push_back(*it);
it++;
}
- return pseries(var, point, new_seq);
+ return pseries(r, new_seq);
}
} else
- return convert_to_poly().series(s, p, order);
+ return convert_to_poly().series(r, order);
}
/** Compute the truncated series expansion of an expression.
- * This function returns an expression containing an object of class pseries to
- * represent the series. If the series does not terminate within the given
+ * This function returns an expression containing an object of class pseries
+ * to represent the series. If the series does not terminate within the given
* truncation order, the last term of the series will be an order term.
*
- * @param s expansion variable
- * @param point expansion point
+ * @param r expansion relation, lhs holds variable and rhs holds point
* @param order truncation order of series calculations
* @return an expression holding a pseries object */
-ex ex::series(const symbol &s, const ex &point, int order) const
+ex ex::series(const ex & r, int order) const
{
GINAC_ASSERT(bp!=0);
- ex e;
- try {
- e = bp->series(s, point, order);
- } catch (exception &x) {
- throw (std::logic_error(string("unable to compute series (") + x.what() + ")"));
- }
- return e;
+ ex e;
+ relational rel_;
+
+ if (is_ex_exactly_of_type(r,relational))
+ rel_ = ex_to_relational(r);
+ else if (is_ex_exactly_of_type(r,symbol))
+ rel_ = relational(r,_ex0());
+ else
+ throw (std::logic_error("ex::series(): expansion point has unknown type"));
+
+ try {
+ e = bp->series(rel_, order);
+ } catch (exception &x) {
+ throw (std::logic_error(string("unable to compute series (") + x.what() + ")"));
+ }
+ return e;
}
// other constructors
public:
- pseries(const ex &var_, const ex &point_, const epvector &ops_);
+ pseries(const ex &rel_, const epvector &ops_);
// functions overriding virtual functions from base classes
public:
ex collect(const symbol &s) const;
ex eval(int level=0) const;
ex evalf(int level=0) const;
- ex series(const symbol & s, const ex & p, int order) const;
+ ex series(const relational & r, int order) const;
ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
ex subs(const lst & ls, const lst & lr) const;
protected:
int ldegree(const symbol & s) const;
ex coeff(const symbol & s, int n = 1) const;
ex eval(int level = 0) const;
- ex series(const symbol & s, const ex & point, int order) const;
+ ex series(const relational & s, int order) const;
ex normal(lst &sym_lst, lst &repl_lst, int level=0) const;
ex subs(const lst & ls, const lst & lr) const;
protected:
.BI rem( expression ", " expression ", " symbol )
\- remainder of polynomials
.br
-.BI series( expression ", " "symbol [" ", " "point [" ", " order]] )
+.BI series( expression ", " relation-or-symbol ", " order )
\- series expansion
.br
.BI sqrfree( expression ", " symbol )
return rem(e[0], e[1], ex_to_symbol(e[2]));
}
-static ex f_series2(const exprseq &e)
+static ex f_series(const exprseq &e)
{
- CHECK_ARG(1, symbol, series);
- return e[0].series(ex_to_symbol(e[1]), ex(0));
-}
-
-static ex f_series3(const exprseq &e)
-{
- CHECK_ARG(1, symbol, series);
- return e[0].series(ex_to_symbol(e[1]), e[2]);
-}
-
-static ex f_series4(const exprseq &e)
-{
- CHECK_ARG(1, symbol, series);
- CHECK_ARG(3, numeric, series);
- return e[0].series(ex_to_symbol(e[1]), e[2], ex_to_numeric(e[3]).to_int());
+ CHECK_ARG(2, numeric, series);
+ return e[0].series(e[1], ex_to_numeric(e[2]).to_int());
}
static ex f_sqrfree(const exprseq &e)
{"primpart", fcn_desc(f_primpart, 2)},
{"quo", fcn_desc(f_quo, 3)},
{"rem", fcn_desc(f_rem, 3)},
- {"series", fcn_desc(f_series2, 2)},
- {"series", fcn_desc(f_series3, 3)},
- {"series", fcn_desc(f_series4, 4)},
+ {"series", fcn_desc(f_series, 3)},
{"sqrfree", fcn_desc(f_sqrfree, 2)},
{"sqrt", fcn_desc(f_sqrt, 1)},
{"subs", fcn_desc(f_subs2, 2)},