[GSoC] Implemented convolution operation of two formal power series

sympy/series/formal.py

@@ -16,14 +16,15 @@
 from sympy.core.singleton import S
 from sympy.core.symbol import Wild, Dummy, symbols, Symbol
 from sympy.core.sympify import sympify
+from sympy.discrete.convolutions import convolution
 from sympy.functions.combinatorial.factorials import binomial, factorial, rf
 from sympy.functions.elementary.integers import floor, frac, ceiling
 from sympy.functions.elementary.miscellaneous import Min, Max
 from sympy.functions.elementary.piecewise import Piecewise
 from sympy.series.limits import Limit
 from sympy.series.order import Order
 from sympy.simplify.powsimp import powsimp
-from sympy.series.sequences import sequence
+from sympy.series.sequences import sequence, SeqMul
 from sympy.series.series_class import SeriesBase
 
 
@@ -1143,6 +1144,98 @@ def integrate(self, x=None, **kwargs):
 
         return self.func(f, self.x, self.x0, self.dir, (ak, self.xk, ind))
 
+    def convolve(self, other, x=None, n=6, cycle=0, method=None):
+        """Convolute two Formal Power Series and return the truncated terms
+        upto specified order.
+
+        Parameters
+        ==========
+
+        n : Number, optional
+            Specifies the order of the term up to which the polynomial should
+            be truncated.
+        cycle : Integer
+            Specifies the length for doing cyclic convolution.
+        method : 'dyadic', 'subset' or None
+            Set dyadic to identify the convolution type as dyadic (*bitwise-XOR*)
+            convolution, which is performed using **FWHT**.
+            Set subset to identify the convolution type as subset convolution.
+            Else, set None to interpret it as simple linear convolution. If
+            cycle is non-zero, then it will be interpreted as cyclic convolution.
+
+        Examples
+        ========
+
+        >>> from sympy import fps, sin, exp, convolution
+        >>> from sympy.abc import x
+        >>> f1 = fps(sin(x))
+        >>> f2 = fps(exp(x))
+
+        >>> f1.convolve(f2, x, 4)
+        x + x**2 + x**3/3 + O(x**4)
+
+        >>> f1.convolve(f2, x, 4, cycle=4)
+        11*x/12 + 35*x**2/36 + x**3/3 + O(x**4)
+
+        >>> f1.convolve(f2, x, 6, method='dyadic')
+        4667/4800 + 2641*x/2880 + x**3/3 + x**4/60 + x**5/20 + O(x**6)
+
+        >>> f1.convolve(f2, x, 6, method='subset')
+        x + x**3/3 + x**5/20 + O(x**6)
+
+        See Also
+        ========
+
+        sympy.discrete.convolutions
+        """
+        if x is None:
+            x = self.x
+        if n is None:
+            return iter(self)
+
+        other = sympify(other)
+
+        if  not isinstance(other, FormalPowerSeries):
+            raise ValueError("Both series should be an instance of FormalPowerSeries"
+                             " class.")
+
+        if self.dir != other.dir:
+            raise ValueError("Both series should be calculated from the"
+                             " same direction.")
+        elif self.x0 != other.x0:
+            raise ValueError("Both series should be calculated about the"
+                             " same point.")
+
+        elif self.x != other.x:
+            raise ValueError("Both series should have the same symbol.")
+
+        if cycle > n:
+            raise ValueError("Value of cycle should be less than or equal to n.")
+
+        if method not in ['dyadic', 'subset', None]:
+            raise ValueError("Method of convolution should either be dyadic, subset or None.")
+
+        k = self.ak.variables[0]
+        coeff1 = sequence(self.ak.formula, (k, 0, n-1))
+
+        k = other.ak.variables[0]
+        coeff2 = sequence(other.ak.formula, (k, 0, n-1))
+
+        conv_coeff =[]
+        if method == 'dyadic':
+            conv_coeff = convolution(coeff1, coeff2, cycle, None, None, True, False)
+        elif method == 'subset':
+            conv_coeff = convolution(coeff1, coeff2, cycle, None, None, False, True)
+        elif method is None:
+            conv_coeff = convolution(coeff1, coeff2, cycle, None, None, False, False)
+
+        conv_seq = sequence(tuple(conv_coeff), (k, 0, n-1))
+        k = self.xk.variables[0]
+        xk_seq = sequence(self.xk.formula, (k, 0, n-1))
+        terms = xk_seq * conv_seq
+
+        return Add(*terms) + Order(self.xk.coeff(n), (self.x, self.x0))
+
     def __add__(self, other):
         other = sympify(other)
 

sympy/series/tests/test_formal.py

@@ -507,3 +507,29 @@ def test_fps__operations():
     fi = f2.integrate(x)
     assert fi.function == sin(x)
     assert fi.truncate() == x - x**3/6 + x**5/120 + O(x**6)
+
+def test_fps__convolution():
+    f1, f2, f3 = fps(sin(x)), fps(exp(x)), fps(cos(x))
+
+    raises(ValueError, lambda: f1.convolve(exp(x), x))
+    raises(ValueError, lambda: f1.convolve(fps(exp(x), dir=-1), x, 4))
+    raises(ValueError, lambda: f1.convolve(fps(exp(x), x0=1), x, 4))
+    raises(ValueError, lambda: f1.convolve(fps(exp(y)), x, 4))
+
+    assert f1.convolve(f2, x, 3) == x + x**2 + O(x**3)
+    assert f1.convolve(f2, x, 4) == x + x**2 + x**3/3 + O(x**4)
+    assert f1.convolve(f3, x, 4) == x - 2*x**3/3 + O(x**4)
+
+    assert f1.convolve(f2, x, 4, cycle=4) == 11*x/12 + 35*x**2/36 + x**3/3 + O(x**4)
+    assert f1.convolve(f3, x, 4, cycle=3) == -S(2)/3 + x + x**2/12 - 2*x**3/3 + O(x**4)
+    raises(ValueError, lambda: f1.convolve(f2, x, 4, cycle=6))
+
+    assert f1.convolve(f2, x, 6, method='dyadic') == \
+        S(4667)/4800 + 2641*x/2880 + x**3/3 + x**4/60 + x**5/20 + O(x**6)
+    assert f1.convolve(f3, x, 5, cycle=4, method='dyadic') == 9*x/8 - 97*x**3/144 + O(x**5)
+
+    raises(ValueError, lambda: f1.convolve(f2, x, 6, method='linear'))
+
+    assert f1.convolve(f2, x, 6, method='subset') == x + x**3/3 + x**5/20 + O(x**6)
+    assert f1.convolve(f3, x, 5, cycle=5, method='subset') == \
+        S(1)/24 + x - x**2/144 - 2*x**3/3 + O(x**5)