gh-35894: Regular sequences: implement convolution / ring structure
    
### 📚 Description

Add/implement convolution of two regular sequences, therefore making the
parent a ring.

Details of the algorithm: [convolution-sagemath-pr-
35894.pdf](https://github.com/sagemath/sage/files/12216984/convolution-
sagemath-pr-35894.pdf)

See also meta issue #21202.

### 📝 Checklist

<!-- Put an `x` in all the boxes that apply. It should be `[x]` not `[x
]`. -->

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies

No formal dependencies, but there are trivial merge conflicts---in both
branches, methods are inserted at the same position in the code---
between:
- #21343: deal with mu[0]*right != right in k-regular sequences
- #35894: Regular sequences: implement convolution / ring structure
    
URL: #35894
Reported by: Daniel Krenn
Reviewer(s): cheuberg, Daniel Krenn

Author: Release Manager

src/sage/combinat/k_regular_sequence.py

@@ -1037,6 +1037,115 @@ def forward_differences(self, **kwds):
         """
         return self.subsequence(1, {1: 1, 0: -1}, **kwds)
 
+    @minimize_result
+    def _mul_(self, other):
+        r"""
+        Return the product of this `k`-regular sequence with ``other``,
+        where the multiplication is convolution of power series.
+
+        The operator `*` is mapped to :meth:`convolution`.
+
+        INPUT:
+
+        - ``other`` -- a :class:`kRegularSequence`
+
+        - ``minimize`` -- (default: ``None``) a boolean or ``None``.
+          If ``True``, then :meth:`~RecognizableSeries.minimized` is called after the operation,
+          if ``False``, then not. If this argument is ``None``, then
+          the default specified by the parent's ``minimize_results`` is used.
+
+        OUTPUT:
+
+        A :class:`kRegularSequence`
+
+        ALGORITHM:
+
+        See pdf attached to
+        `github pull request #35894 <https://github.com/sagemath/sage/pull/35894>`_
+        which contains a draft describing the details of the used algorithm.
+
+        EXAMPLES::
+
+            sage: Seq2 = kRegularSequenceSpace(2, ZZ)
+            sage: E = Seq2((Matrix([[0, 1], [0, 1]]), Matrix([[0, 0], [0, 1]])),
+            ....:          vector([1, 0]), vector([1, 1]))
+            sage: E
+            2-regular sequence 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, ...
+
+        We can build the convolution (in the sense of power-series) of `E` by
+        itself via::
+
+            sage: E.convolution(E)
+            2-regular sequence 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, ...
+
+        This is the same as using multiplication operator::
+
+            sage: E * E
+            2-regular sequence 1, 0, 2, 0, 3, 0, 4, 0, 5, 0, ...
+
+        Building :meth:`partial_sums` can also be seen as a convolution::
+
+            sage: o = Seq2.one_hadamard()
+            sage: E * o
+            2-regular sequence 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, ...
+            sage: E * o == E.partial_sums(include_n=True)
+            True
+
+        TESTS::
+
+            sage: E * o == o * E
+            True
+        """
+        from sage.arith.srange import srange
+        from sage.matrix.constructor import Matrix
+        from sage.matrix.special import zero_matrix
+        from sage.modules.free_module_element import vector
+
+        P = self.parent()
+        k = P.k
+
+        def tensor_product(left, right):
+            T = left.tensor_product(right)
+            T.subdivide()
+            return T
+
+        matrices_0 = {r: sum(tensor_product(self.mu[s], other.mu[r-s])
+                             for s in srange(0, r+1))
+                      for r in P.alphabet()}
+        matrices_1 = {r: sum(tensor_product(self.mu[s], other.mu[k+r-s])
+                             for s in srange(r+1, k))
+                      for r in P.alphabet()}
+        left = vector(tensor_product(Matrix(self.left), Matrix(other.left)))
+        right = vector(tensor_product(Matrix(self.right), Matrix(other.right)))
+
+        def linear_representation_morphism_recurrence_order_1(C, D):
+            r"""
+                Return the morphism of a linear representation
+                for the sequence `z_n` satisfying
+                `z_{kn+r} = C_r z_n + D_r z_{n-1}`.
+            """
+            Z = zero_matrix(C[0].dimensions()[0])
+
+            def blocks(r):
+                upper = list([C[s], D[s], Z]
+                             for s in reversed(srange(max(0, r-2), r+1)))
+                lower = list([Z, C[s], D[s]]
+                             for s in reversed(srange(k-3+len(upper), k)))
+                return upper + lower
+
+            return {r: Matrix.block(blocks(r)) for r in P.alphabet()}
+
+        result = P.element_class(
+            P,
+            linear_representation_morphism_recurrence_order_1(matrices_0,
+                                                              matrices_1),
+            vector(list(left) + (2*len(list(left)))*[0]),
+            vector(list(right) + (2*len(list(right)))*[0]))
+
+        return result
+
+    convolution = _mul_
+
     @minimize_result
     def partial_sums(self, include_n=False):
         r"""
@@ -1223,7 +1332,10 @@ class kRegularSequenceSpace(RecognizableSeriesSpace):
     Element = kRegularSequence
 
     @classmethod
-    def __normalize__(cls, k, coefficient_ring, **kwds):
+    def __normalize__(cls, k,
+                      coefficient_ring,
+                      category=None,
+                      **kwds):
         r"""
         Normalizes the input in order to ensure a unique
         representation.
@@ -1234,13 +1346,17 @@ def __normalize__(cls, k, coefficient_ring, **kwds):
 
             sage: Seq2 = kRegularSequenceSpace(2, ZZ)
             sage: Seq2.category()
-            Category of modules over Integer Ring
+            Category of algebras over Integer Ring
             sage: Seq2.alphabet()
             {0, 1}
         """
         from sage.arith.srange import srange
+        from sage.categories.algebras import Algebras
+        category = category or Algebras(coefficient_ring)
         nargs = super().__normalize__(coefficient_ring,
-                                      alphabet=srange(k), **kwds)
+                                      alphabet=srange(k),
+                                      category=category,
+                                      **kwds)
         return (k,) + nargs
 
     def __init__(self, k, *args, **kwds):
@@ -1331,6 +1447,39 @@ def _n_to_index_(self, n):
         except OverflowError:
             raise ValueError('value {} of index is negative'.format(n)) from None
 
+    @cached_method
+    def one(self):
+        r"""
+        Return the one element of this :class:`kRegularSequenceSpace`,
+        i.e. the unique neutral element for `*` and also
+        the embedding of the one of the coefficient ring into
+        this :class:`kRegularSequenceSpace`.
+
+        EXAMPLES::
+
+            sage: Seq2 = kRegularSequenceSpace(2, ZZ)
+            sage: O = Seq2.one(); O
+            2-regular sequence 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
+            sage: O.linear_representation()
+            ((1), Finite family {0: [1], 1: [0]}, (1))
+
+        TESTS::
+
+            sage: Seq2.one() is Seq2.one()
+            True
+        """
+        from sage.matrix.constructor import Matrix
+        from sage.modules.free_module_element import vector
+
+        R = self.coefficient_ring()
+        one = R.one()
+        zero = R.zero()
+        return self.element_class(self,
+                                  [Matrix([[one]])]
+                                  + (self.k-1)*[Matrix([[zero]])],
+                                  vector([one]),
+                                  vector([one]))
+
     def some_elements(self):
         r"""
         Return some elements of this `k`-regular sequence.

src/sage/combinat/recognizable_series.py

@@ -1140,6 +1140,12 @@ def minimized(self):
             sage: all(c == d and v == w
             ....:     for (c, v), (d, w) in islice(zip(iter(S), iter(M)), 20))
             True
+
+        TESTS::
+
+            sage: Rec((Matrix([[0]]), Matrix([[0]])),
+            ....:     vector([1]), vector([0])).minimized().linear_representation()
+            ((), Finite family {0: [], 1: []}, ())
         """
         return self._minimized_right_()._minimized_left_()
 
@@ -1314,7 +1320,7 @@ def _add_(self, other):
 
         result = P.element_class(
             P,
-            dict((a, self.mu[a].block_sum(other.mu[a])) for a in P.alphabet()),
+            {a: self.mu[a].block_sum(other.mu[a]) for a in P.alphabet()},
             vector(tuple(self.left) + tuple(other.left)),
             vector(tuple(self.right) + tuple(other.right)))
 
@@ -1376,6 +1382,11 @@ def _rmul_(self, other):
             sage: 1 * E is E
             True
 
+        ::
+
+            sage: 0 * E is Seq2.zero()
+            True
+
         We test that ``_rmul_`` and ``_lmul_`` are actually called::
 
             sage: def print_name(f):
@@ -1394,9 +1405,11 @@ def _rmul_(self, other):
             _lmul_
             2-regular sequence 2, 0, 2, 0, 2, 0, 2, 0, 2, 0, ...
         """
+        P = self.parent()
+        if other.is_zero():
+            return P._zero_()
         if other.is_one():
             return self
-        P = self.parent()
         return P.element_class(P, self.mu, other * self.left, self.right)
 
     def _lmul_(self, other):
@@ -1433,6 +1446,11 @@ def _lmul_(self, other):
             sage: E * 1 is E
             True
 
+        ::
+
+            sage: E * 0 is Seq2.zero()
+            True
+
         The following is not tested, as `MS^i` for integers `i` does
         not work, thus ``vector([m])`` fails. (See :trac:`21317` for
         details.)
@@ -1449,9 +1467,11 @@ def _lmul_(self, other):
             sage: M  # not tested
             sage: M.linear_representation()  # not tested
         """
+        P = self.parent()
+        if other.is_zero():
+            return P._zero_()
         if other.is_one():
             return self
-        P = self.parent()
         return P.element_class(P, self.mu, self.left, self.right * other)
 
     @minimize_result
@@ -1557,8 +1577,7 @@ def tensor_product(left, right):
             return T
         result = P.element_class(
             P,
-            dict((a, tensor_product(self.mu[a], other.mu[a]))
-                 for a in P.alphabet()),
+            {a: tensor_product(self.mu[a], other.mu[a]) for a in P.alphabet()},
             vector(tensor_product(Matrix(self.left), Matrix(other.left))),
             vector(tensor_product(Matrix(self.right), Matrix(other.right))))
 
@@ -1952,6 +1971,59 @@ def some_elements(self, **kwds):
                     break
                 yield self(mu, *LR, **kwds)
 
+    @cached_method
+    def _zero_(self):
+        r"""
+        Return the zero element of this :class:`RecognizableSeriesSpace`,
+        i.e. the unique neutral element for `+`.
+
+        TESTS::
+
+            sage: Rec = RecognizableSeriesSpace(ZZ, [0, 1])
+            sage: Z = Rec._zero_(); Z
+            0
+            sage: Z.linear_representation()
+            ((), Finite family {0: [], 1: []}, ())
+        """
+        from sage.matrix.constructor import Matrix
+        from sage.modules.free_module_element import vector
+        from sage.sets.family import Family
+
+        return self.element_class(
+            self, Family(self.alphabet(), lambda a: Matrix()),
+            vector([]), vector([]))
+
+    @cached_method
+    def one(self):
+        r"""
+        Return the one element of this :class:`RecognizableSeriesSpace`,
+        i.e. the embedding of the one of the coefficient ring into
+        this :class:`RecognizableSeriesSpace`.
+
+        EXAMPLES::
+
+            sage: Rec = RecognizableSeriesSpace(ZZ, [0, 1])
+            sage: O = Rec.one(); O
+            [] + ...
+            sage: O.linear_representation()
+            ((1), Finite family {0: [0], 1: [0]}, (1))
+
+        TESTS::
+
+            sage: Rec.one() is Rec.one()
+            True
+        """
+        from sage.matrix.constructor import Matrix
+        from sage.modules.free_module_element import vector
+
+        R = self.coefficient_ring()
+        one = R.one()
+        zero = R.zero()
+        return self.element_class(self,
+                                  len(self.alphabet())*[Matrix([[zero]])],
+                                  vector([one]),
+                                  vector([one]))
+
     @cached_method
     def one_hadamard(self):
         r"""
@@ -1979,7 +2051,7 @@ def one_hadamard(self):
         from sage.modules.free_module_element import vector
 
         one = self.coefficient_ring()(1)
-        return self(dict((a, Matrix([[one]])) for a in self.alphabet()),
+        return self({a: Matrix([[one]]) for a in self.alphabet()},
                     vector([one]), vector([one]))
 
     def _element_constructor_(self, data,
@@ -2006,6 +2078,19 @@ def _element_constructor_(self, data,
             sage: Rec(S) is S
             True
 
+        ::
+
+            sage: A = Rec(42); A
+            42*[] + ...
+            sage: A.linear_representation()
+            ((42), Finite family {0: [0], 1: [0]}, (1))
+            sage: Z = Rec(0); Z
+            0
+            sage: Z.linear_representation()
+            ((), Finite family {0: [], 1: []}, ())
+
+        ::
+
             sage: Rec((M0, M1))
             Traceback (most recent call last):
             ...
@@ -2024,20 +2109,18 @@ def _element_constructor_(self, data,
             ValueError: left or right vector is None
         """
         if isinstance(data, int) and data == 0:
-            from sage.matrix.constructor import Matrix
-            from sage.modules.free_module_element import vector
-            from sage.sets.family import Family
-
-            return self.element_class(
-                self, Family(self.alphabet(), lambda a: Matrix()),
-                vector([]), vector([]))
+            return self._zero_()
 
         if type(data) == self.element_class and data.parent() == self:
             element = data
 
         elif isinstance(data, RecognizableSeries):
             element = self.element_class(self, data.mu, data.left, data.right)
 
+        elif data in self.coefficient_ring():
+            c = self.coefficient_ring()(data)
+            return c * self.one()
+
         else:
             mu = data
             if left is None or right is None: