- pseries::coeff() now uses a binary instead of a linear search algorithm

[ginac.git] / ginac / pseries.cpp

diff --git a/ginac/pseries.cpp b/ginac/pseries.cpp
index 459b46920b59793071185925a9efd61aa502c107..87d0e4377a6ccbb3023723759b5c6be3773b45d8 100644 (file)
--- a/ginac/pseries.cpp
+++ b/ginac/pseries.cpp
@@ -242,20 +242,42 @@ 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++;
-        }
-        return _ex0();
+               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
+{
+       if (var.is_equal(s))
+               return convert_to_poly();
+       else
+               return inherited::collect(s);
+}
+
 ex pseries::eval(int level) const
 {
     if (level == 1)
@@ -298,6 +320,32 @@ ex pseries::subs(const lst & ls, const lst & lr) const
     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(var, point, new_seq);
+    } else {
+        return *this;
+    }
+}
+
 
 /*
  *  Construct ordinary polynomial out of series
@@ -663,6 +711,31 @@ 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
+{
+       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(var, point, new_seq);
+               }
+       } else
+               return convert_to_poly().series(s, p, 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