]> www.ginac.de Git - ginac.git/blobdiff - ginac/pseries.cpp
- Fixed a thinko in atan(const numeric &).
[ginac.git] / ginac / pseries.cpp
index 0f29e5f547698e08323e31525f6dc451a4ec49cc..f55c0fede6079ad75800c87fed184edf4194bf78 100644 (file)
@@ -35,9 +35,9 @@
 #include "utils.h"
 #include "debugmsg.h"
 
-#ifndef NO_GINAC_NAMESPACE
+#ifndef NO_NAMESPACE_GINAC
 namespace GiNaC {
-#endif // ndef NO_GINAC_NAMESPACE
+#endif // ndef NO_NAMESPACE_GINAC
 
 GINAC_IMPLEMENT_REGISTERED_CLASS(pseries, basic)
 
@@ -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 ex,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);
 }
 
 
@@ -117,7 +119,7 @@ pseries::pseries(const ex &var_, const ex &point_, const epvector &ops_)
 pseries::pseries(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
 {
     debugmsg("pseries constructor from archive_node", LOGLEVEL_CONSTRUCT);
-    for (unsigned int i=0; true; i++) {
+    for (unsigned int i=0; true; ++i) {
         ex rest;
         ex coeff;
         if (n.find_ex("coeff", rest, sym_lst, i) && n.find_ex("power", coeff, sym_lst, i))
@@ -143,16 +145,15 @@ void pseries::archive(archive_node &n) const
     while (i != iend) {
         n.add_ex("coeff", i->rest);
         n.add_ex("power", i->coeff);
-        i++;
+        ++i;
     }
     n.add_ex("var", var);
     n.add_ex("point", point);
 }
 
-
-/*
- *  Functions overriding virtual functions from base classes
- */
+//////////
+// functions overriding virtual functions from bases classes
+//////////
 
 basic *pseries::duplicate() const
 {
@@ -160,27 +161,81 @@ basic *pseries::duplicate() const
     return new pseries(*this);
 }
 
-void pseries::print(ostream &os, unsigned upper_precedence) const
+void pseries::print(std::ostream &os, unsigned upper_precedence) const
 {
     debugmsg("pseries print", LOGLEVEL_PRINT);
-    convert_to_poly().print(os, upper_precedence);
+    for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+        // omit zero terms
+        if (i->rest.is_zero())
+            continue;
+        // print a sign, if needed
+        if (i!=seq.begin())
+            os << '+';
+        if (!is_order_function(i->rest)) {
+            // print 'rest', i.e. the expansion coefficient
+            if (i->rest.info(info_flags::numeric) &&
+                i->rest.info(info_flags::positive)) {
+                os << i->rest;
+            } else
+                os << "(" << i->rest << ')';
+            // print 'coeff', something like (x-1)^42
+            if (!i->coeff.is_zero()) {
+                os << '*';
+                if (!point.is_zero())
+                    os << '(' << var-point << ')';
+                else
+                    os << var;
+                if (i->coeff.compare(_ex1())) {
+                    os << '^';
+                    if (i->coeff.info(info_flags::negative))
+                        os << '(' << i->coeff << ')';
+                    else
+                        os << i->coeff;
+                }
+            }
+        } else {
+            os << Order(power(var-point,i->coeff));
+        }
+    }
 }
 
-void pseries::printraw(ostream &os) const
+
+void pseries::printraw(std::ostream &os) const
 {
-       debugmsg("pseries printraw", LOGLEVEL_PRINT);
-       os << "pseries(" << var << ";" << point << ";";
-       for (epvector::const_iterator i=seq.begin(); i!=seq.end(); i++) {
-               os << "(" << (*i).rest << "," << (*i).coeff << "),";
-       }
-       os << ")";
+    debugmsg("pseries printraw", LOGLEVEL_PRINT);
+    os << "pseries(" << var << ";" << point << ";";
+    for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i) {
+        os << "(" << (*i).rest << "," << (*i).coeff << "),";
+    }
+    os << ")";
 }
 
+
+void pseries::printtree(std::ostream & os, unsigned indent) const
+{
+    debugmsg("pseries printtree",LOGLEVEL_PRINT);
+    os << std::string(indent,' ') << "pseries " 
+       << ", hash=" << hashvalue
+       << " (0x" << std::hex << hashvalue << std::dec << ")"
+       << ", flags=" << flags << std::endl;
+    for (unsigned i=0; i<seq.size(); ++i) {
+        seq[i].rest.printtree(os,indent+delta_indent);
+        seq[i].coeff.printtree(os,indent+delta_indent);
+        if (i!=seq.size()-1)
+            os << std::string(indent+delta_indent,' ') << "-----" << std::endl;
+    }
+    var.printtree(os, indent+delta_indent);
+    point.printtree(os, indent+delta_indent);
+}
+
+/** Return the number of operands including a possible order term. */
 unsigned pseries::nops(void) const
 {
     return seq.size();
 }
 
+
+/** Return the ith term in the series when represented as a sum. */
 ex pseries::op(int i) const
 {
     if (i < 0 || unsigned(i) >= seq.size())
@@ -188,11 +243,16 @@ ex pseries::op(int i) const
     return seq[i].rest * power(var - point, seq[i].coeff);
 }
 
+
 ex &pseries::let_op(int i)
 {
     throw (std::logic_error("let_op not defined for pseries"));
 }
 
+
+/** Return degree of highest power of the series.  This is usually the exponent
+ *  of the Order term.  If s is not the expansion variable of the series, the
+ *  series is examined termwise. */
 int pseries::degree(const symbol &s) const
 {
     if (var.is_equal(s)) {
@@ -210,12 +270,17 @@ int pseries::degree(const symbol &s) const
             int pow = it->rest.degree(s);
             if (pow > max_pow)
                 max_pow = pow;
-            it++;
+            ++it;
         }
         return max_pow;
     }
 }
 
+/** Return degree of lowest power of the series.  This is usually the exponent
+ *  of the leading term.  If s is not the expansion variable of the series, the
+ *  series is examined termwise.  If s is the expansion variable but the
+ *  expansion point is not zero the series is not expanded to find the degree.
+ *  I.e.: (1-x) + (1-x)^2 + Order((1-x)^3) has ldegree(x) 1, not 0. */
 int pseries::ldegree(const symbol &s) const
 {
     if (var.is_equal(s)) {
@@ -233,7 +298,7 @@ int pseries::ldegree(const symbol &s) const
             int pow = it->rest.ldegree(s);
             if (pow < min_pow)
                 min_pow = pow;
-            it++;
+            ++it;
         }
         return min_pow;
     }
@@ -242,60 +307,143 @@ int pseries::ldegree(const symbol &s) const
 ex pseries::coeff(const symbol &s, int n) const
 {
     if (var.is_equal(s)) {
-        epvector::const_iterator it = seq.begin(), itend = seq.end();
-        while (it != itend) {
-            int pow = ex_to_numeric(it->coeff).to_int();
-            if (pow == n)
-                return it->rest;
-            if (pow > n)
-                return _ex0();
-            it++;
+        if (seq.size() == 0)
+            return _ex0();
+        
+        // Binary search in sequence for given power
+        numeric looking_for = numeric(n);
+        int lo = 0, hi = seq.size() - 1;
+        while (lo <= hi) {
+            int mid = (lo + hi) / 2;
+            GINAC_ASSERT(is_ex_exactly_of_type(seq[mid].coeff, numeric));
+            int cmp = ex_to_numeric(seq[mid].coeff).compare(looking_for);
+            switch (cmp) {
+                case -1:
+                    lo = mid + 1;
+                    break;
+                case 0:
+                    return seq[mid].rest;
+                case 1:
+                    hi = mid - 1;
+                    break;
+                default:
+                    throw(std::logic_error("pseries::coeff: compare() didn't return -1, 0 or 1"));
+            }
         }
         return _ex0();
     } else
         return convert_to_poly().coeff(s, n);
 }
 
+
+ex pseries::collect(const symbol &s) const
+{
+    return *this;
+}
+
+
+/** Evaluate coefficients. */
 ex pseries::eval(int level) const
 {
     if (level == 1)
         return this->hold();
     
+    if (level == -max_recursion_level)
+        throw (std::runtime_error("pseries::eval(): recursion limit exceeded"));
+    
     // Construct a new series with evaluated coefficients
     epvector new_seq;
     new_seq.reserve(seq.size());
     epvector::const_iterator it = seq.begin(), itend = seq.end();
     while (it != itend) {
         new_seq.push_back(expair(it->rest.eval(level-1), it->coeff));
-        it++;
+        ++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 numerically.  The order term is dropped. */
+
+/** Evaluate coefficients numerically. */
 ex pseries::evalf(int level) const
 {
-    return convert_to_poly().evalf(level);
+    if (level == 1)
+        return *this;
+    
+    if (level == -max_recursion_level)
+        throw (std::runtime_error("pseries::evalf(): recursion limit exceeded"));
+    
+    // Construct a new series with evaluated coefficients
+    epvector new_seq;
+    new_seq.reserve(seq.size());
+    epvector::const_iterator it = seq.begin(), itend = seq.end();
+    while (it != itend) {
+        new_seq.push_back(expair(it->rest.evalf(level-1), it->coeff));
+        ++it;
+    }
+    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
 {
-       // If expansion variable is being substituted, convert the series to a
-       // polynomial and do the substitution there because the result might
-       // no longer be a power series
-       if (ls.has(var))
-               return convert_to_poly(true).subs(ls, lr);
+    // If expansion variable is being substituted, convert the series to a
+    // polynomial and do the substitution there because the result might
+    // no longer be a power series
+    if (ls.has(var))
+        return convert_to_poly(true).subs(ls, lr);
+    
+    // Otherwise construct a new series with substituted coefficients and
+    // expansion point
+    epvector newseq;
+    newseq.reserve(seq.size());
+    epvector::const_iterator it = seq.begin(), itend = seq.end();
+    while (it != itend) {
+        newseq.push_back(expair(it->rest.subs(ls, lr), it->coeff));
+        ++it;
+    }
+    return (new pseries(relational(var,point.subs(ls, lr)), newseq))->setflag(status_flags::dynallocated);
+}
+
+
+/** Implementation of ex::expand() for a power series.  It expands all the
+ *  terms individually and returns the resulting series as a new pseries.
+ *  @see ex::diff */
+ex pseries::expand(unsigned options) const
+{
+    epvector newseq;
+    newseq.reserve(seq.size());
+    for (epvector::const_iterator i=seq.begin(); i!=seq.end(); ++i)
+        newseq.push_back(expair(i->rest.expand(), i->coeff));
+    return (new pseries(relational(var,point), newseq))
+        ->setflag(status_flags::dynallocated |
+                  status_flags::expanded);
+}
+
 
-       // Otherwise construct a new series with substituted coefficients and
-       // expansion point
-       epvector new_seq;
-       new_seq.reserve(seq.size());
-       epvector::const_iterator it = seq.begin(), itend = seq.end();
-       while (it != itend) {
-               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);
+/** Implementation of ex::diff() for a power series.  It treats the series as a
+ *  polynomial.
+ *  @see ex::diff */
+ex pseries::derivative(const symbol & s) const
+{
+    if (s == var) {
+        epvector new_seq;
+        epvector::const_iterator it = seq.begin(), itend = seq.end();
+        
+        // FIXME: coeff might depend on var
+        while (it != itend) {
+            if (is_order_function(it->rest)) {
+                new_seq.push_back(expair(it->rest, it->coeff - 1));
+            } else {
+                ex c = it->rest * it->coeff;
+                if (!c.is_zero())
+                    new_seq.push_back(expair(c, it->coeff - 1));
+            }
+            ++it;
+        }
+        return pseries(relational(var,point), new_seq);
+    } else {
+        return *this;
+    }
 }
 
 
@@ -317,11 +465,18 @@ ex pseries::convert_to_poly(bool no_order) const
                 e += Order(power(var - point, it->coeff));
         } else
             e += it->rest * power(var - point, it->coeff);
-        it++;
+        ++it;
     }
     return e;
 }
 
+/** Returns true if there is no order term, i.e. the series terminates and
+ *  false otherwise. */
+bool pseries::is_terminating(void) const
+{
+    return !is_order_function((seq.end()-1)->rest);
+}
+
 
 /*
  *  Implementation of series expansion
@@ -329,51 +484,57 @@ 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, unsigned options) 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++) {
+    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 = deriv.subs(r);
         if (!coeff.is_zero())
-            seq.push_back(expair(coeff, numeric(n)));
+            seq.push_back(expair(fac.inverse() * coeff, numeric(n)));
     }
     
     // Higher-order terms, if present
-    deriv = deriv.diff(s);
-    if (!deriv.is_zero())
+    deriv = deriv.diff(*s);
+    if (!deriv.expand().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, unsigned options) const
 {
-       epvector seq;
-       if (is_equal(s)) {
-               if (order > 0 && !point.is_zero())
-                       seq.push_back(expair(point, _ex0()));
-               if (order > 1)
-                       seq.push_back(expair(_ex1(), _ex1()));
-               else
-                       seq.push_back(expair(Order(_ex1()), numeric(order)));
-       } else
-               seq.push_back(expair(*this, _ex0()));
-       return pseries(s, point, seq);
+    epvector seq;
+    const ex point = r.rhs();
+    GINAC_ASSERT(is_ex_exactly_of_type(r.lhs(),symbol));
+    const symbol *s = static_cast<symbol *>(r.lhs().bp);
+    
+    if (this->is_equal(*s)) {
+        if (order > 0 && !point.is_zero())
+            seq.push_back(expair(point, _ex0()));
+        if (order > 1)
+            seq.push_back(expair(_ex1(), _ex1()));
+        else
+            seq.push_back(expair(Order(_ex1()), numeric(order)));
+    } else
+        seq.push_back(expair(*this, _ex0()));
+    return pseries(r, seq);
 }
 
 
@@ -389,7 +550,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
@@ -404,7 +565,7 @@ ex pseries::add_series(const pseries &other) const
         if (a == a_end) {
             while (b != b_end) {
                 new_seq.push_back(*b);
-                b++;
+                ++b;
             }
             break;
         } else
@@ -414,7 +575,7 @@ ex pseries::add_series(const pseries &other) const
         if (b == b_end) {
             while (a != a_end) {
                 new_seq.push_back(*a);
-                a++;
+                ++a;
             }
             break;
         } else
@@ -426,13 +587,13 @@ ex pseries::add_series(const pseries &other) const
             new_seq.push_back(*a);
             if (is_order_function((*a).rest))
                 break;
-            a++;
+            ++a;
         } else if (pow_b < pow_a) {
             // b has lesser power, get coefficient from b
             new_seq.push_back(*b);
             if (is_order_function((*b).rest))
                 break;
-            b++;
+            ++b;
         } else {
             // Add coefficient of a and b
             if (is_order_function((*a).rest) || is_order_function((*b).rest)) {
@@ -442,34 +603,34 @@ ex pseries::add_series(const pseries &other) const
                 ex sum = (*a).rest + (*b).rest;
                 if (!(sum.is_zero()))
                     new_seq.push_back(expair(sum, numeric(pow_a)));
-                a++;
-                b++;
+                ++a;
+                ++b;
             }
         }
     }
-    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, unsigned options) const
 {
     ex acc; // Series accumulator
     
     // Get first term from overall_coeff
-    acc = overall_coeff.series(s, point, order);
-
+    acc = overall_coeff.series(r, order, options);
+    
     // Add remaining terms
     epvector::const_iterator it = seq.begin();
     epvector::const_iterator itend = seq.end();
-    for (; it!=itend; it++) {
+    for (; it!=itend; ++it) {
         ex op;
         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, options);
         if (!it->coeff.is_equal(_ex1()))
             op = ex_to_pseries(op).mul_const(ex_to_numeric(it->coeff));
         
@@ -496,9 +657,9 @@ ex pseries::mul_const(const numeric &other) const
             new_seq.push_back(expair(it->rest * other, it->coeff));
         else
             new_seq.push_back(*it);
-        it++;
+        ++it;
     }
-    return pseries(var, point, new_seq);
+    return pseries(relational(var,point), new_seq);
 }
 
 
@@ -514,9 +675,9 @@ 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
     epvector new_seq;
     
@@ -534,42 +695,42 @@ ex pseries::mul_series(const pseries &other) const
         higher_order_a = a_max + b_min;
     if (is_order_function(other.coeff(*s, b_max)))
         higher_order_b = b_max + a_min;
-    int higher_order_c = min(higher_order_a, higher_order_b);
+    int higher_order_c = std::min(higher_order_a, higher_order_b);
     if (cdeg_max >= higher_order_c)
         cdeg_max = higher_order_c - 1;
     
-    for (int cdeg=cdeg_min; cdeg<=cdeg_max; cdeg++) {
+    for (int cdeg=cdeg_min; cdeg<=cdeg_max; ++cdeg) {
         ex co = _ex0();
         // c(i)=a(0)b(i)+...+a(i)b(0)
-        for (int i=a_min; cdeg-i>=b_min; i++) {
+        for (int i=a_min; cdeg-i>=b_min; ++i) {
             ex a_coeff = coeff(*s, i);
             ex b_coeff = other.coeff(*s, cdeg-i);
             if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
-                co += coeff(*s, i) * other.coeff(*s, cdeg-i);
+                co += a_coeff * b_coeff;
         }
         if (!co.is_zero())
             new_seq.push_back(expair(co, numeric(cdeg)));
     }
     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, unsigned options) const
 {
     ex acc; // Series accumulator
     
     // Get first term from overall_coeff
-    acc = overall_coeff.series(s, point, order);
+    acc = overall_coeff.series(r, order, options);
     
     // Multiply with remaining terms
     epvector::const_iterator it = seq.begin();
     epvector::const_iterator itend = seq.end();
-    for (; it!=itend; it++) {
+    for (; it!=itend; ++it) {
         ex op = it->rest;
         if (op.info(info_flags::numeric)) {
             // series * const (special case, faster)
@@ -577,7 +738,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, options);
         if (!it->coeff.is_equal(_ex1()))
             op = ex_to_pseries(op).power_const(ex_to_numeric(it->coeff), order);
 
@@ -604,9 +765,9 @@ ex pseries::power_const(const numeric &p, int deg) const
     ex co0;
     co.push_back(co0 = power(coeff(*s, ldeg), p));
     bool all_sums_zero = true;
-    for (i=1; i<deg; i++) {
+    for (i=1; i<deg; ++i) {
         ex sum = _ex0();
-        for (int j=1; j<=i; j++) {
+        for (int j=1; j<=i; ++j) {
             ex c = coeff(*s, j + ldeg);
             if (is_order_function(c)) {
                 co.push_back(Order(_ex1()));
@@ -622,7 +783,7 @@ ex pseries::power_const(const numeric &p, int deg) const
     // Construct new series (of non-zero coefficients)
     epvector new_seq;
     bool higher_order = false;
-    for (i=0; i<deg; i++) {
+    for (i=0; i<deg; ++i) {
         if (!co[i].is_zero())
             new_seq.push_back(expair(co[i], numeric(i) + p * ldeg));
         if (is_order_function(co[i])) {
@@ -632,27 +793,37 @@ 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);
+}
+
+
+/** Return a new pseries object with the powers shifted by deg. */
+pseries pseries::shift_exponents(int deg) const
+{
+    epvector newseq(seq);
+    for (epvector::iterator i=newseq.begin(); i!=newseq.end(); ++i)
+        i->coeff = i->coeff + deg;
+    return pseries(relational(var, point), newseq);
 }
 
 
 /** 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, unsigned options) 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, options);
         
         // 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, options);
         
         // Singularity encountered, expand basis into series
-        e = basis.series(s, point, order);
+        e = basis.series(r, order, options);
     } else {
         // Basis is a series
         e = basis;
@@ -663,19 +834,63 @@ ex power::series(const symbol & s, const ex & point, int order) const
 }
 
 
+/** Re-expansion of a pseries object. */
+ex pseries::series(const relational & r, int order, unsigned options) const
+{
+    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;
+            epvector::const_iterator it = seq.begin(), itend = seq.end();
+            while (it != itend) {
+                int o = ex_to_numeric(it->coeff).to_int();
+                if (o >= order) {
+                    new_seq.push_back(expair(Order(_ex1()), o));
+                    break;
+                }
+                new_seq.push_back(*it);
+                ++it;
+            }
+            return pseries(r, new_seq);
+        }
+    } else
+        return convert_to_poly().series(r, order, options);
+}
+
+
 /** 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
+ *  @param options  of class series_options
  *  @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, unsigned options) const
 {
     GINAC_ASSERT(bp!=0);
-    return bp->series(s, point, order);
+    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, options);
+    } catch (std::exception &x) {
+        throw (std::logic_error(std::string("unable to compute series (") + x.what() + ")"));
+    }
+    return e;
 }
 
 
@@ -683,6 +898,6 @@ ex ex::series(const symbol &s, const ex &point, int order) const
 const pseries some_pseries;
 const type_info & typeid_pseries = typeid(some_pseries);
 
-#ifndef NO_GINAC_NAMESPACE
+#ifndef NO_NAMESPACE_GINAC
 } // namespace GiNaC
-#endif // ndef NO_GINAC_NAMESPACE
+#endif // ndef NO_NAMESPACE_GINAC