Fix series of positive powers of polynomials.

Before calling pseries::power_const(), make sure the basis series has
sufficient terms. In particular, do not shrink the order of expansion
- only grow it.

Fix various hairy problems in pseries::power_const() when a polynomial
is raised to a positive integer power.

(The special cases here can actually also be computed by Taylor
expansion but the fixes should be more general.)

Also add test cases.

Thanks to Vitaly Magerya <vmagerya@gmail.com> for reporting this.

diff --git a/check/exam_pseries.cpp b/check/exam_pseries.cpp
index f0fead05232ccd13fd4c6f5559fc28b7c7f41afd..671d57512205d0476eb2ba869e7e99e3e429315b 100644 (file)
--- a/check/exam_pseries.cpp
+++ b/check/exam_pseries.cpp
@@ -390,6 +390,24 @@ static unsigned exam_series14()
        return result;
 }
 
+// Test expansion of powers of polynomials.
+static unsigned exam_series15()
+{
+       unsigned result = 0;
+
+       ex e = pow(x + pow(x,2), 2);
+
+       result += check_series(e, 0, Order(1), 0);
+       result += check_series(e, 0, Order(x), 1);
+       result += check_series(e, 0, Order(pow(x,2)), 2);
+       result += check_series(e, 0, pow(x,2) + Order(pow(x,3)), 3);
+       result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + Order(pow(x,4)), 4);
+       result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + pow(x,4), 5);
+       result += check_series(e, 0, pow(x,2) + 2*pow(x,3) + pow(x,4), 6);
+
+       return result;
+}
+
 unsigned exam_pseries()
 {
        unsigned result = 0;
@@ -410,6 +428,7 @@ unsigned exam_pseries()
        result += exam_series12();  cout << '.' << flush;
        result += exam_series13();  cout << '.' << flush;
        result += exam_series14();  cout << '.' << flush;
+       result += exam_series15();  cout << '.' << flush;
        
        return result;
 }

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index 00188f0147d599d2f1d15ebf52f8362cf6242e74..af26a6dd5250326845b5832891d1380c0d3191f0 100644 (file)
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -992,8 +992,8 @@ ex pseries::power_const(const numeric &p, int deg) const
        // which can easily be solved given the starting value c_0 = (a_0)^p.
        // For the more general case where the leading coefficient of A(x) is not
        // a constant, just consider A2(x) = A(x)*x^m, with some integer m and
-       // repeat the above derivation.  The leading power of C2(x) = A2(x)^2 is
-       // then of course x^(p*m) but the recurrence formula still holds.
+       // repeat the above derivation.  The leading power of C2(x) = A2(x)^p is
+       // then of course a_0^p*x^(p*m) but the recurrence formula still holds.
        
        if (seq.empty()) {
                // as a special case, handle the empty (zero) series honoring the
@@ -1005,16 +1005,27 @@ ex pseries::power_const(const numeric &p, int deg) const
                else
                        return *this;
        }
-       
-       const int ldeg = ldegree(var);
-       if (!(p*ldeg).is_integer())
+
+       const int base_ldeg = ldegree(var);
+       if (!(p*base_ldeg).is_integer())
                throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
+       int new_ldeg = (p*base_ldeg).to_int();
+
+       const int base_deg = degree(var);
+       int new_deg = deg;
+       if (p.is_pos_integer()) {
+               // No need to compute beyond p*base_deg.
+               new_deg = std::min((p*base_deg).to_int(), deg);
+       }
 
        // adjust number of coefficients
-       int numcoeff = deg - (p*ldeg).to_int();
+       int numcoeff = new_deg - new_ldeg;
+       if (new_deg < deg)
+               numcoeff += 1;
+
        if (numcoeff <= 0) {
-               epvector epv { expair(Order(_ex1), deg) };
-               return dynallocate<pseries>(relational(var,point), std::move(epv));
+               return dynallocate<pseries>(relational(var, point),
+                                           epvector{{Order(_ex1), deg}});
        }
        
        // O(x^n)^(-m) is undefined
@@ -1024,33 +1035,35 @@ ex pseries::power_const(const numeric &p, int deg) const
        // Compute coefficients of the powered series
        exvector co;
        co.reserve(numcoeff);
-       co.push_back(pow(coeff(var, ldeg), p));
+       co.push_back(pow(coeff(var, base_ldeg), p));
        for (int i=1; i<numcoeff; ++i) {
                ex sum = _ex0;
                for (int j=1; j<=i; ++j) {
-                       ex c = coeff(var, j + ldeg);
+                       ex c = coeff(var, j + base_ldeg);
                        if (is_order_function(c)) {
                                co.push_back(Order(_ex1));
                                break;
                        } else
                                sum += (p * j - (i - j)) * co[i - j] * c;
                }
-               co.push_back(sum / coeff(var, ldeg) / i);
+               co.push_back(sum / coeff(var, base_ldeg) / i);
        }
        
        // Construct new series (of non-zero coefficients)
        epvector new_seq;
        bool higher_order = false;
        for (int i=0; i<numcoeff; ++i) {
-               if (!co[i].is_zero())
-                       new_seq.emplace_back(expair(co[i], p * ldeg + i));
+               if (!co[i].is_zero()) {
+                       new_seq.emplace_back(expair(co[i], new_ldeg + i));
+               }
                if (is_order_function(co[i])) {
                        higher_order = true;
                        break;
                }
        }
-       if (!higher_order)
-               new_seq.emplace_back(expair(Order(_ex1), p * ldeg + numcoeff));
+       if (!higher_order && new_deg == deg) {
+               new_seq.emplace_back(expair{Order(_ex1), new_deg});
+       }
 
        return pseries(relational(var,point), std::move(new_seq));
 }
@@ -1150,7 +1163,8 @@ ex power::series(const relational & r, int order, unsigned options) const
 
        if (!(real_ldegree*numexp).is_integer())
                throw std::runtime_error("pseries::power_const(): trying to assemble a Puiseux series");
-       ex e = basis.series(r, (order + real_ldegree*(1-numexp)).to_int(), options);
+       int extra_terms = (real_ldegree*(1-numexp)).to_int();
+       ex e = basis.series(r, order + std::max(0, extra_terms), options);
        
        ex result;
        try {