Added missing namespace qualifier

[ginac.git] / ginac / pseries.cpp

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index dd4f51ab91ba3becc8007706491b30e4fd6afa01..d622ecc20f6c442c2576bc6836073be7bac27f41 100644 (file)
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -520,7 +520,7 @@ ex basic::series(const relational & r, int order, unsigned options) const
                seq.push_back(expair(coeff, _ex0));
        
        int n;
-       for (n=1; n<order; ++n) {
+       for (n=1; n<=order; ++n) {
                fac = fac.mul(n);
                // We need to test for zero in order to see if the series terminates.
                // The problem is that there is no such thing as a perfect test for
@@ -535,9 +535,34 @@ ex basic::series(const relational & r, int order, unsigned options) const
        }
        
        // Higher-order terms, if present
-       deriv = deriv.diff(s);
-       if (!deriv.expand().is_zero())
-               seq.push_back(expair(Order(_ex1), n));
+       int ldeg;
+       try {
+               ldeg = std::abs(deriv.ldegree(s));
+       }
+       catch (std::runtime_error) {
+               ldeg = 0;
+       }
+       if (ldeg != 0) {
+               // pure polynomial
+               if (!deriv.subs(r, subs_options::no_pattern).is_zero()) {
+                       seq.push_back(expair(Order(_ex1), n));
+               } else {
+               seq.push_back(expair(Order(_ex1), n+ldeg-1));
+               }
+       } else {
+               // something more complicated -> loop until next coefficient is found
+               for (;; ++n) {
+                       deriv = deriv.diff(s).expand();
+                       if (deriv.is_zero()) {
+                               break;
+                       }
+                       if (!deriv.subs(r, subs_options::no_pattern).is_zero()) {
+                               seq.push_back(expair(Order(_ex1), n));
+                               break;
+                       }
+               }
+       }
+       
        return pseries(r, seq);
 }
 
@@ -774,7 +799,7 @@ ex mul::series(const relational & r, int order, unsigned options) const
                        acc = ex_to<pseries>(acc.mul_series(ex_to<pseries>(op)));
        }
 
-       if (negldegree && (seq.size() > 1)) {
+       if (seq.size() > 1) {
 
                // re-calculation of series
 
@@ -857,7 +882,7 @@ ex pseries::power_const(const numeric &p, int deg) const
        co.reserve(deg);
        co.push_back(power(coeff(var, ldeg), p));
        bool all_sums_zero = true;
-       for (int i=1; i<deg; ++i) {
+       for (int i=1; i<=deg; ++i) {
                ex sum = _ex0;
                for (int j=1; j<=i; ++j) {
                        ex c = coeff(var, j + ldeg);
@@ -875,7 +900,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 (int i=0; i<deg; ++i) {
+       for (int i=0; i<=deg; ++i) {
                if (!co[i].is_zero())
                        new_seq.push_back(expair(co[i], p * ldeg + i));
                if (is_order_function(co[i])) {
@@ -884,7 +909,7 @@ ex pseries::power_const(const numeric &p, int deg) const
                }
        }
        if (!higher_order && !all_sums_zero)
-               new_seq.push_back(expair(Order(_ex1), p * ldeg + deg));
+               new_seq.push_back(expair(Order(_ex1), p * ldeg + deg + 1));
 
        return pseries(relational(var,point), new_seq);
 }
@@ -923,7 +948,7 @@ ex power::series(const relational & r, int order, unsigned options) const
        // Is the expression of type something^(-int)?
        if (!must_expand_basis && !exponent.info(info_flags::negint))
                return basic::series(r, order, options);
-               
+
        // Is the expression of type 0^something?
        if (!must_expand_basis && !basis.subs(r, subs_options::no_pattern).is_zero())
                return basic::series(r, order, options);