Trac #27929: division/inversion of lazy power series

URL: https://trac.sagemath.org/27929
Reported by: dkrenn
Ticket author(s): Daniel Krenn
Reviewer(s): Travis Scrimshaw

src/sage/combinat/species/series.py

@@ -1213,6 +1213,97 @@ def __pow__(self, n):
             raise ValueError("n must be a nonnegative integer")
         return prod([self]*n, self.parent().identity_element())
 
+    def __invert__(self):
+        """
+        Return 1 over this power series, i.e. invert this power series.
+
+        EXAMPLES::
+
+            sage: L = LazyPowerSeriesRing(QQ)
+            sage: x = L.gen()
+
+        Geometric series::
+
+            sage: a = ~(1-x); a.compute_coefficients(10); a
+            1 + x + x^2 + x^3 + x^4 + x^5 + x^6 + x^7 + x^8 + x^9 + x^10 + O(x^11)
+
+        (Shifted) Fibonacci numbers::
+
+            sage: b = ~(1-x-x^2); b.compute_coefficients(10); b
+            1 + x + 2*x^2 + 3*x^3 + 5*x^4 + 8*x^5
+            + 13*x^6 + 21*x^7 + 34*x^8 + 55*x^9 + 89*x^10 + O(x^11)
+
+        Series whose constant coefficient is `0` cannot be inverted::
+
+            sage: ~x
+            Traceback (most recent call last):
+            ....
+            ZeroDivisionError: cannot invert x because constant coefficient is 0
+        """
+        if self.get_aorder() > 0:
+            raise ZeroDivisionError(
+                'cannot invert {} because '
+                'constant coefficient is 0'.format(self))
+        return self._new(self._invert_gen, lambda a: 0, self)
+
+    invert = __invert__
+
+    def _invert_gen(self, ao):
+        r"""
+        Return an iterator for the coefficients of 1 over this power series.
+
+        TESTS::
+
+            sage: L = LazyPowerSeriesRing(QQ)
+            sage: f = L([1, -1, 0])
+            sage: g = f._invert_gen(0)
+            sage: [next(g) for i in range(10)]
+            [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
+        """
+        from itertools import count
+
+        assert ao == 0
+
+        ic0 = ~self.coefficient(0)
+        yield ic0
+        if self.order == 0:
+            return
+
+        one = self.parent()(1)
+        base = one - ic0 * self
+        base.coefficient(0)
+        ao_base = base.get_aorder()
+        assert ao_base >= 1
+
+        current = one + base
+        k = 1
+        for n in count(1):
+            while ao_base*k < n:
+                current = one + base * current
+                k += 1
+                current.coefficient(n)  # make sure new current is initialized
+            ao_base = base.get_aorder()  # update this so that while above is faster
+            yield current.coefficient(n) * ic0
+
+    def _div_(self, other):
+        """
+        Divide this power series by ``other``.
+
+        EXAMPLES::
+
+            sage: L = LazyPowerSeriesRing(QQ)
+            sage: x = L.gen()
+
+        Fibonacci numbers::
+
+            sage: b = x / (1-x-x^2); b.compute_coefficients(10); b
+            x + x^2 + 2*x^3 + 3*x^4 + 5*x^5 + 8*x^6
+            + 13*x^7 + 21*x^8 + 34*x^9 + 55*x^10 + O(x^11)
+        """
+        return self * ~other
+
+    div = _div_
+
     def __call__(self, y):
         """
         Return the composition of this power series and the power series y.