- 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 0824232..47dfa00 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 0b1e7f1..49eddd6 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 a5625ef..ba7b319 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 5fb3ddc..c92f29f 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 0809521..e594a0b 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 6982872..890bafe 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 89a7b29..ae19fba 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 0c67fd2..c17aea2 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 4e5ec2b..a3c34e9 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 5c3c147..5da77cc 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 c54cc38..574fa37 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 bab6a88..86fda5b 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 40dd07e..94512e0 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 caa4175..bb6c974 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 fcf2d05..3e74166 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 fe059a4..5b35963 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 87a77cb..12fabc0 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 01b0b34..a0a4600 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 d23b228..8524006 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 7a76217..afb5b09 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 4faa720..035551d 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 bf05af9..db79e14 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 6ab13d1..3a81f55 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 414f95b..53e5bcd 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 77d21e7..6931d18 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)},