From f293ecba8b6026a7754795256b2f23910bf70507 Mon Sep 17 00:00:00 2001 From: Richard Kreckel Date: Thu, 23 Mar 2000 16:50:17 +0000 Subject: [PATCH] - We now write f(x).series(x==3,5) instead of f(x).series(x,3,5) and f(x).series(x,4) instead of f(x).series(x,0,4). We also don't allow default arguments any more. --- NEWS | 11 +-- check/exam_differentiation.cpp | 4 +- check/exam_pseries.cpp | 6 +- check/time_gammaseries.cpp | 2 +- cint/dummies.h | 11 ++- cint/dummies.pl | 5 +- doc/tutorial/ginac.texi | 28 +++++--- ginac/Makefile.am | 6 +- ginac/add.h | 2 +- ginac/basic.cpp | 5 ++ ginac/basic.h | 3 +- ginac/container.pl | 2 - ginac/ex.h | 8 +-- ginac/function.pl | 14 ++-- ginac/inifcns.cpp | 8 ++- ginac/inifcns_gamma.cpp | 40 ++++++----- ginac/inifcns_trans.cpp | 18 ++--- ginac/mul.h | 2 +- ginac/normal.cpp | 2 +- ginac/power.h | 2 +- ginac/pseries.cpp | 128 +++++++++++++++++++-------------- ginac/pseries.h | 4 +- ginac/symbol.h | 2 +- ginsh/ginsh.1 | 2 +- ginsh/ginsh_parser.yy | 23 ++---- 25 files changed, 177 insertions(+), 161 deletions(-) diff --git a/NEWS b/NEWS index 0824232a..47dfa00e 100644 --- a/NEWS +++ b/NEWS @@ -1,10 +1,13 @@ 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 diff --git a/check/exam_differentiation.cpp b/check/exam_differentiation.cpp index 0b1e7f19..49eddd66 100644 --- a/check/exam_differentiation.cpp +++ b/check/exam_differentiation.cpp @@ -254,8 +254,8 @@ static unsigned exam_differentiation6(void) 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); diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp index a5625efa..ba7b3193 100644 --- a/check/exam_pseries.cpp +++ b/check/exam_pseries.cpp @@ -26,7 +26,7 @@ static symbol x("x"); 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 @@ -107,7 +107,7 @@ static unsigned exam_series2(void) 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); @@ -120,7 +120,7 @@ static unsigned exam_series3(void) 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); diff --git a/check/time_gammaseries.cpp b/check/time_gammaseries.cpp index 5fb3ddc9..c92f29f7 100644 --- a/check/time_gammaseries.cpp +++ b/check/time_gammaseries.cpp @@ -27,7 +27,7 @@ unsigned Gammaseries(unsigned order) 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 diff --git a/cint/dummies.h b/cint/dummies.h index 08095215..e594a0b4 100644 --- a/cint/dummies.h +++ b/cint/dummies.h @@ -787,14 +787,11 @@ inline ex evalf(basic const & x, int y=0) { 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) { diff --git a/cint/dummies.pl b/cint/dummies.pl index 69828728..890bafe8 100644 --- a/cint/dummies.pl +++ b/cint/dummies.pl @@ -188,9 +188,8 @@ inline_single_function_2p('ex','collect','basic const &','ex','symbol const &',' 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',''); diff --git a/doc/tutorial/ginac.texi b/doc/tutorial/ginac.texi index 89a7b29d..ae19fbaa 100644 --- a/doc/tutorial/ginac.texi +++ b/doc/tutorial/ginac.texi @@ -376,24 +376,25 @@ polynomials): @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 @@ -1192,7 +1193,8 @@ They are created by simply using the C++ operators @code{==}, @code{!=}, 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 @@ -1604,13 +1606,13 @@ int main() 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; // ... @} @@ -1648,7 +1650,7 @@ using namespace GiNaC; 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; @@ -1666,7 +1668,11 @@ int main() @} @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 diff --git a/ginac/Makefile.am b/ginac/Makefile.am index 0c67fd2d..c17aea20 100644 --- a/ginac/Makefile.am +++ b/ginac/Makefile.am @@ -22,13 +22,13 @@ EXTRA_DIST = container.pl function.pl structure.pl # 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 diff --git a/ginac/add.h b/ginac/add.h index 4e5ec2bb..a3c34e9c 100644 --- a/ginac/add.h +++ b/ginac/add.h @@ -70,7 +70,7 @@ public: 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; diff --git a/ginac/basic.cpp b/ginac/basic.cpp index 5c3c147f..5da77cc6 100644 --- a/ginac/basic.cpp +++ b/ginac/basic.cpp @@ -343,11 +343,16 @@ ex basic::derivative(const symbol & s) 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; diff --git a/ginac/basic.h b/ginac/basic.h index c54cc381..574fa37c 100644 --- a/ginac/basic.h +++ b/ginac/basic.h @@ -44,6 +44,7 @@ class ex; class symbol; class lst; class numeric; +class relational; class archive_node; //typedef vector exvector; @@ -137,7 +138,7 @@ public: // only const functions please (may break reference counting) 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; diff --git a/ginac/container.pl b/ginac/container.pl index bab6a88c..86fda5b4 100755 --- a/ginac/container.pl +++ b/ginac/container.pl @@ -1,5 +1,3 @@ -#!/usr/bin/perl -w - if (($#ARGV!=0) and ($#ARGV!=1)) { die 'usage: container.pl type [maxargs] (type=lst or exprseq)'; } diff --git a/ginac/ex.h b/ginac/ex.h index 40dd07e0..94512e08 100644 --- a/ginac/ex.h +++ b/ginac/ex.h @@ -250,9 +250,9 @@ public: 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; @@ -413,8 +413,8 @@ inline ex evalf(const ex & thisex, int level = 0) 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); } diff --git a/ginac/function.pl b/ginac/function.pl index caa4175c..bb6c974f 100755 --- a/ginac/function.pl +++ b/ginac/function.pl @@ -1,5 +1,3 @@ -#!/usr/bin/perl -w - $maxargs=13; sub generate_seq { @@ -155,7 +153,7 @@ $typedef_derivative_funcp=generate( '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",'',''); @@ -205,9 +203,9 @@ $series_switch_statement=generate( <<'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; @@ -479,7 +477,7 @@ public: 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: @@ -947,12 +945,12 @@ ex function::thisexprseq(exvector * vp) const /** 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(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 diff --git a/ginac/inifcns_gamma.cpp b/ginac/inifcns_gamma.cpp index fe059a44..5b35963f 100644 --- a/ginac/inifcns_gamma.cpp +++ b/ginac/inifcns_gamma.cpp @@ -104,7 +104,7 @@ static ex Gamma_deriv(const ex & x, unsigned deriv_param) } -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 @@ -114,7 +114,7 @@ static ex Gamma_series(const ex & x, const symbol & s, const ex & pt, int order) // 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: @@ -122,7 +122,7 @@ static ex Gamma_series(const ex & x, const symbol & s, const ex & pt, int order) 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); } @@ -199,36 +199,38 @@ static ex Beta_deriv(const ex & x, const ex & y, unsigned deriv_param) } -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(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); } @@ -304,7 +306,7 @@ static ex psi1_deriv(const ex & x, unsigned deriv_param) 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 @@ -314,7 +316,7 @@ static ex psi1_series(const ex & x, const symbol & s, const ex & pt, int order) // 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: @@ -322,7 +324,7 @@ static ex psi1_series(const ex & x, const symbol & s, const ex & pt, int order) 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 = @@ -424,7 +426,7 @@ static ex psi2_deriv(const ex & n, const ex & x, unsigned deriv_param) 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 @@ -435,7 +437,7 @@ static ex psi2_series(const ex & n, const ex & x, const symbol & s, const ex & p // 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: @@ -444,7 +446,7 @@ static ex psi2_series(const ex & n, const ex & x, const symbol & s, const ex & p 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 = diff --git a/ginac/inifcns_trans.cpp b/ginac/inifcns_trans.cpp index 87a77cb2..12fabc0e 100644 --- a/ginac/inifcns_trans.cpp +++ b/ginac/inifcns_trans.cpp @@ -144,12 +144,12 @@ static ex log_deriv(const ex & x, unsigned deriv_param) 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(); } @@ -395,16 +395,16 @@ static ex tan_deriv(const ex & x, unsigned deriv_param) 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). @@ -752,16 +752,16 @@ static ex tanh_deriv(const ex & x, unsigned deriv_param) 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). diff --git a/ginac/mul.h b/ginac/mul.h index 01b0b349..a0a46002 100644 --- a/ginac/mul.h +++ b/ginac/mul.h @@ -72,7 +72,7 @@ public: 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; diff --git a/ginac/normal.cpp b/ginac/normal.cpp index d23b2289..85240060 100644 --- a/ginac/normal.cpp +++ b/ginac/normal.cpp @@ -1622,7 +1622,7 @@ ex pseries::normal(lst &sym_lst, lst &repl_lst, int level) 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); } diff --git a/ginac/power.h b/ginac/power.h index 7a76217d..afb5b094 100644 --- a/ginac/power.h +++ b/ginac/power.h @@ -73,7 +73,7 @@ public: 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; diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp index 4faa7200..035551d8 100644 --- a/ginac/pseries.cpp +++ b/ginac/pseries.cpp @@ -97,15 +97,17 @@ void pseries::destroy(bool call_parent) * 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(rel_.lhs().bp); } @@ -292,7 +294,7 @@ ex pseries::eval(int level) const 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. */ @@ -312,7 +314,7 @@ ex pseries::evalf(int level) const 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 @@ -332,7 +334,7 @@ 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 @@ -355,7 +357,7 @@ ex pseries::derivative(const symbol & s) const } it++; } - return pseries(var, point, new_seq); + return pseries(relational(var,point), new_seq); } else { return *this; } @@ -392,42 +394,48 @@ ex pseries::convert_to_poly(bool no_order) const /** 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(r.lhs().bp); + if (!coeff.is_zero()) seq.push_back(expair(coeff, numeric(0))); int n; for (n=1; n(r.lhs().bp); + + if (this->is_equal(*s)) { if (order > 0 && !point.is_zero()) seq.push_back(expair(point, _ex0())); if (order > 1) @@ -436,7 +444,7 @@ ex symbol::series(const symbol & s, const ex & point, int order) const seq.push_back(expair(Order(_ex1()), numeric(order))); } else seq.push_back(expair(*this, _ex0())); - return pseries(s, point, seq); + return pseries(r, seq); } @@ -452,7 +460,7 @@ ex pseries::add_series(const pseries &other) const 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 @@ -510,19 +518,19 @@ ex pseries::add_series(const pseries &other) const } } } - 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(); @@ -532,7 +540,7 @@ ex add::series(const symbol & s, const ex & point, int order) const 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)); @@ -561,7 +569,7 @@ ex pseries::mul_const(const numeric &other) const new_seq.push_back(*it); it++; } - return pseries(var, point, new_seq); + return pseries(relational(var,point), new_seq); } @@ -577,7 +585,7 @@ ex pseries::mul_series(const pseries &other) const 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 @@ -615,19 +623,19 @@ ex pseries::mul_series(const pseries &other) const } 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(); @@ -640,7 +648,7 @@ ex mul::series(const symbol & s, const ex & point, int order) const 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); @@ -695,27 +703,27 @@ ex pseries::power_const(const numeric &p, int deg) const } 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; @@ -727,10 +735,14 @@ ex power::series(const symbol & s, const ex & point, int order) const /** 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(r.lhs().bp); + + if (var.is_equal(*s) && point.is_equal(p)) { + if (order > degree(*s)) return *this; else { epvector new_seq; @@ -744,32 +756,40 @@ ex pseries::series(const symbol & s, const ex & p, int order) const 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; } diff --git a/ginac/pseries.h b/ginac/pseries.h index bf05af9b..db79e141 100644 --- a/ginac/pseries.h +++ b/ginac/pseries.h @@ -50,7 +50,7 @@ protected: // 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: @@ -66,7 +66,7 @@ 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: diff --git a/ginac/symbol.h b/ginac/symbol.h index 6ab13d12..3a81f55d 100644 --- a/ginac/symbol.h +++ b/ginac/symbol.h @@ -81,7 +81,7 @@ public: 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: diff --git a/ginsh/ginsh.1 b/ginsh/ginsh.1 index 414f95b4..53e5bcdb 100644 --- a/ginsh/ginsh.1 +++ b/ginsh/ginsh.1 @@ -306,7 +306,7 @@ detail here. Please refer to the GiNaC documentation. .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 ) diff --git a/ginsh/ginsh_parser.yy b/ginsh/ginsh_parser.yy index 77d21e72..6931d189 100644 --- a/ginsh/ginsh_parser.yy +++ b/ginsh/ginsh_parser.yy @@ -420,23 +420,10 @@ static ex f_rem(const exprseq &e) 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) @@ -529,9 +516,7 @@ static const fcn_init builtin_fcns[] = { {"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)}, -- 2.47.0