cache pseries coeff accesses in pseries::mul_series.

[ginac.git] / ginac / pseries.cpp

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index 31a92002958399818659b65e83094b498a66dbc4..8f267e8de4d49c0062d5d1ed910aea1c61f53cff 100644 (file)
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -177,7 +177,7 @@ void pseries::print_series(const print_context & c, const char *openbrace, const
                                }
                        }
                } else
-                       Order(power(var-point,i->coeff)).print(c);
+                       Order(pow(var - point, i->coeff)).print(c);
                ++i;
        }
 
@@ -281,8 +281,8 @@ ex pseries::op(size_t i) const
                throw (std::out_of_range("op() out of range"));
 
        if (is_order_function(seq[i].rest))
-               return Order(power(var-point, seq[i].coeff));
-       return seq[i].rest * power(var - point, seq[i].coeff);
+               return Order(pow(var-point, seq[i].coeff));
+       return seq[i].rest * pow(var - point, seq[i].coeff);
 }
 
 /** Return degree of highest power of the series.  This is usually the exponent
@@ -500,8 +500,7 @@ ex pseries::evalm() const
                        ex newcoeff = i->rest.evalm();
                        if (!newcoeff.is_zero())
                                newseq.push_back(expair(newcoeff, i->coeff));
-               }
-               else {
+               } else {
                        ex newcoeff = i->rest.evalm();
                        if (!are_ex_trivially_equal(newcoeff, i->rest)) {
                                something_changed = true;
@@ -589,9 +588,9 @@ ex pseries::convert_to_poly(bool no_order) const
        for (auto & it : seq) {
                if (is_order_function(it.rest)) {
                        if (!no_order)
-                               e += Order(power(var - point, it.coeff));
+                               e += Order(pow(var - point, it.coeff));
                } else
-                       e += it.rest * power(var - point, it.coeff);
+                       e += it.rest * pow(var - point, it.coeff);
        }
        return e;
 }
@@ -644,7 +643,7 @@ ex basic::series(const relational & r, int order, unsigned options) const
 
        int n;
        for (n=1; n<order; ++n) {
-               fac = fac.mul(n);
+               fac = fac.div(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
                // zero.  Expanding the term occasionally helps a little...
@@ -654,7 +653,7 @@ ex basic::series(const relational & r, int order, unsigned options) const
 
                coeff = deriv.subs(r, subs_options::no_pattern);
                if (!coeff.is_zero())
-                       seq.push_back(expair(fac.inverse() * coeff, n));
+                       seq.push_back(expair(fac * coeff, n));
        }
        
        // Higher-order terms, if present
@@ -824,11 +823,11 @@ ex pseries::mul_series(const pseries &other) const
        
        // Series multiplication
        epvector new_seq;
-       int a_max = degree(var);
-       int b_max = other.degree(var);
-       int a_min = ldegree(var);
-       int b_min = other.ldegree(var);
-       int cdeg_min = a_min + b_min;
+       const int a_max = degree(var);
+       const int b_max = other.degree(var);
+       const int a_min = ldegree(var);
+       const int b_min = other.ldegree(var);
+       const int cdeg_min = a_min + b_min;
        int cdeg_max = a_max + b_max;
        
        int higher_order_a = std::numeric_limits<int>::max();
@@ -837,18 +836,30 @@ ex pseries::mul_series(const pseries &other) const
                higher_order_a = a_max + b_min;
        if (is_order_function(other.coeff(var, b_max)))
                higher_order_b = b_max + a_min;
-       int higher_order_c = std::min(higher_order_a, higher_order_b);
+       const int higher_order_c = std::min(higher_order_a, higher_order_b);
        if (cdeg_max >= higher_order_c)
                cdeg_max = higher_order_c - 1;
-       
+
+       std::map<int, ex> rest_map_a, rest_map_b;
+       for (const auto& it : seq)
+               rest_map_a[ex_to<numeric>(it.coeff).to_int()] = it.rest;
+
+       if (other.var.is_equal(var))
+               for (const auto& it : other.seq)
+                       rest_map_b[ex_to<numeric>(it.coeff).to_int()] = it.rest;
+
        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) {
-                       ex a_coeff = coeff(var, i);
-                       ex b_coeff = other.coeff(var, cdeg-i);
-                       if (!is_order_function(a_coeff) && !is_order_function(b_coeff))
-                               co += a_coeff * b_coeff;
+                       const auto& ita = rest_map_a.find(i);
+                       if (ita == rest_map_a.end())
+                               continue;
+                       const auto& itb = rest_map_b.find(cdeg-i);
+                       if (itb == rest_map_b.end())
+                               continue;
+                       if (!is_order_function(ita->second) && !is_order_function(itb->second))
+                               co += ita->second * itb->second;
                }
                if (!co.is_zero())
                        new_seq.push_back(expair(co, numeric(cdeg)));
@@ -1027,7 +1038,7 @@ ex pseries::power_const(const numeric &p, int deg) const
        // Compute coefficients of the powered series
        exvector co;
        co.reserve(numcoeff);
-       co.push_back(power(coeff(var, ldeg), p));
+       co.push_back(pow(coeff(var, ldeg), p));
        for (int i=1; i<numcoeff; ++i) {
                ex sum = _ex0;
                for (int j=1; j<=i; ++j) {