- We now write f(x).series(x==3,5) instead of f(x).series(x,3,5) and
authorRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Thu, 23 Mar 2000 16:50:17 +0000 (16:50 +0000)
committerRichard Kreckel <Richard.Kreckel@uni-mainz.de>
Thu, 23 Mar 2000 16:50:17 +0000 (16:50 +0000)
  f(x).series(x,4) instead of f(x).series(x,0,4).  We also don't allow
  default arguments any more.

25 files changed:
NEWS
check/exam_differentiation.cpp
check/exam_pseries.cpp
check/time_gammaseries.cpp
cint/dummies.h
cint/dummies.pl
doc/tutorial/ginac.texi
ginac/Makefile.am
ginac/add.h
ginac/basic.cpp
ginac/basic.h
ginac/container.pl
ginac/ex.h
ginac/function.pl
ginac/inifcns.cpp
ginac/inifcns_gamma.cpp
ginac/inifcns_trans.cpp
ginac/mul.h
ginac/normal.cpp
ginac/power.h
ginac/pseries.cpp
ginac/pseries.h
ginac/symbol.h
ginsh/ginsh.1
ginsh/ginsh_parser.yy

diff --git a/NEWS b/NEWS
index 0824232ae0651befc14547b75e65a33715549406..47dfa00eb7699768ed27fdd9f0c22736bb45051b 100644 (file)
--- 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
index 0b1e7f19b636686874582f05b98b4c3d8fc7d3e0..49eddd66b5dc57040e86810ffbf8a7551892dd60 100644 (file)
@@ -254,8 +254,8 @@ static unsigned exam_differentiation6(void)
     symbol x("x");
     ex e, d, ed;
     
-    e = sin(x).series(x0, 8);
-    d = cos(x).series(x0, 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);
index a5625efa2c011afbc832c476a22e3069e3fd6e86..ba7b31933b1f17409dea89604bc553fe05157baf 100644 (file)
@@ -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(xpoint, 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);
     
index 5fb3ddc95af6f7b7a5dfd5abb6c54cb03562b2fa..c92f29f7521309e2c5e09febfcf2794615a1a9ff 100644 (file)
@@ -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
index 08095215956720c3bfb3794e9dea6ee8683d4b1e..e594a0b4d15a70162849b399a325438a3d50f990 100644 (file)
@@ -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) {
index 69828728c8003d0e11701170eb84e6bf0222326b..890bafe8dd258c071d8e2946401af6594eeb52e3 100644 (file)
@@ -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','');
index 89a7b29da32d4df2d78f6bf9654f5c15efc38a02..ae19fbaaaf9915903156a7c507cd4e32ee4a68cc 100644 (file)
@@ -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(v0, 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(v0, 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
index 0c67fd2de5d4d7a5a658d72d0381199696060f05..c17aea204d633aebcf6ab085d26c96e451e0c582 100644 (file)
@@ -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
index 4e5ec2bbe014e3c2ba21c9dce799855985125de2..a3c34e9ccd592fa436a13e8144ef7bb4f9177101 100644 (file)
@@ -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;
index 5c3c147faf57d039fe57e94c9be0e2f92112d7bb..5da77cc6e95b3ef523cb69621ee80a8296740e7a 100644 (file)
@@ -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;
index c54cc3816f49b67ab184b07ac9202d7948ea619d..574fa37c006872c99c487186c1305fc3eb852c38 100644 (file)
@@ -44,6 +44,7 @@ class ex;
 class symbol;
 class lst;
 class numeric;
+class relational;
 class archive_node;
 
 //typedef vector<ex> 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;
index bab6a88c34919824a46bd7dd0202ca65eb1ed533..86fda5b48d693e2f109ce896bdafb39cc7f5f844 100755 (executable)
@@ -1,5 +1,3 @@
-#!/usr/bin/perl -w
-
 if (($#ARGV!=0) and ($#ARGV!=1)) {
     die 'usage: container.pl type [maxargs] (type=lst or exprseq)';
 }
index 40dd07e02e77043b3009bcd9d09c783a928af0c0..94512e087e2b0ab85c1b884ff95e1f87cc3f9039 100644 (file)
@@ -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); }
index caa4175c3cc269b77d6619a654c5a8bd783a78f7..bb6c974f6f83ee6b6634cec5b5c14f28d0962a16 100755 (executable)
@@ -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<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) {
index fcf2d05dc088052d60e7d12f9e577b7a103cd4f9..3e74166d601ed46fdf451c229baea45d20e7a699 100644 (file)
@@ -158,12 +158,14 @@ static ex Order_eval(const ex & x)
        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
index fe059a4449003ee08099fb79cbcd274efe19588c..5b35963f078a8ea975b705bc3dc4326a83a2eedc 100644 (file)
@@ -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<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);
 }
 
 
@@ -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 =
index 87a77cb2b97e1c0a42755b09ec2184f7923b920f..12fabc0e98499ea5e4c11ff196f753b6d32b3872 100644 (file)
@@ -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).
index 01b0b3495f4b5481798972d4f1fb627f7a89fb89..a0a460027630db6f48871ec4c7ac1693e09a2809 100644 (file)
@@ -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;
index d23b2289d520aa535bd0dd6fb43b6e43080536a7..852400608fe263832822931758c949725c80fe06 100644 (file)
@@ -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);
 }
 
index 7a76217d2aa94f70f314c9436ae6ef6cfc3cfef3..afb5b0940835aa265d08e06f961a7c9a5dfcfa0b 100644 (file)
@@ -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;
index 4faa7200c8551a161fa700ae80c0a5b4c0b09342..035551d8889ab0f2bef95853d18bd8117c53dd70 100644 (file)
@@ -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<symbol *>(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<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)
@@ -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<symbol *>(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;
 }
 
 
index bf05af9b26b5a14175b33a10e15b6a146743f6bf..db79e141cddf55f2e727761ee7e5a15d50ebdd3e 100644 (file)
@@ -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:
index 6ab13d12ae4a0be94b349de540b2fd30193e2d3f..3a81f55d68f2e2bd52925c0117a26dea1928d97b 100644 (file)
@@ -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:
index 414f95b46de6d813d86d60bd2834bd99aa0be3df..53e5bcdb467e7bf1d4e06d2bf1530f186d691c2a 100644 (file)
@@ -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 )
index 77d21e72d0d873750fceb71e6c06ec39c83d4e80..6931d18944f3d64dbada3089b77723bb17db51fe 100644 (file)
@@ -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)},