code for symmetric matrices

src/sage/matrix/matrix_space.py

@@ -1463,6 +1463,119 @@ def __iter__(self):
                 for iv in sage.combinat.integer_vector.IntegerVectors(weight, number_of_entries, max_part=(order-1)):
                    yield self(entries=[base_elements[i] for i in iv])
 
+
+    def symmetric_matrices(self, f, g=None ):
+        r"""
+        Return a generator of the matrices in this matrix space that satisfy:
+        `A[j,i] = f(A[i,j])` for `i != j`
+        `A[i,i] = g(A[i,i])`
+        
+        If the matrix space doesn't contains square matrices, then a
+        `ValueError` is raised.
+
+        INPUT:
+
+        - `f` -- function
+
+        - `g` -- (optional) funcition; if it is None, then we assume `g=f`,
+          default value: `None`.
+
+        EXAMPLES:
+
+        
+        TESTS::
+        
+        """
+        if self.__nrows != self.__ncols:
+            raise ValueError("can't have symmetric matrices if they are not square")
+        if g == None:
+            g = f
+
+        def make_symmetric( M ):
+            for i in range(M.nrows()):
+                for j in range(i+1, M.nrows()):
+                    M[j,i] = f(M[i,j])
+
+        #Make sure that we can iterate over the base ring
+        base_ring = self.base_ring()
+        base_iter = iter(base_ring)
+
+        nrows = self.__nrows
+        number_of_entries = nrows**2
+        entries_in_upper_half = (nrows*(nrows-1))//2
+
+        #If the number of entries is zero, then just
+        #yield the empty matrix in that case and return
+        if number_of_entries == 0:
+            yield self(0)
+            return
+
+        import sage.combinat.integer_vector
+
+        if not base_ring.is_finite():
+            #When the base ring is not finite, then we should go
+            #through and yield the matrices by "weight", which is
+            #the total number of iterations that need to be done
+            #on the base ring to reach the matrix.
+            base_elements = [ next(base_iter) ]
+            weight = 0
+            while True:
+                for iv in sage.combinat.integer_vector.IntegerVectors(weight, entries_in_upper_half+nrows):
+                    #Now we need to contruct the entries of the matrix so that is symmetric
+                    #with respect to f and g
+                    matrix_entries = [] #the upper half (with diagonal) entries of the matrix
+                    length_of_row = nrows # == self.__ncols
+                    valid_diagonal = True#if false, then we need a new integer vector
+                    for _ in range(nrows):
+                        #check diagonal
+                        if base_elements[iv[0]] != g( base_elements[iv[0]] ):
+                            valid_diagonal = False
+                            break
+                        #now construct the row
+                        zeros = [0]*(nrows - length_of_row)
+                        row = zeros + [base_elements[i] for i in iv[:length_of_row]]
+                        matrix_entries.extend(row) #append row
+
+                        iv = iv[length_of_row:]
+                        length_of_row -= 1
+
+                    if valid_diagonal:
+                        M = self(entries=matrix_entries)
+                        #now make M symmetric with respect to f
+                        make_symmetric(M)
+                        yield M
+                    #else iterate
+                weight += 1
+                base_elements.append( next(base_iter) )
+        else:
+            #In the finite case, we do a similar thing except that
+            #instead of checking if the diagonal is correct after creating the vector
+            #we can select all possible diagonal elements a priori
+            order = base_ring.order()
+            base_elements = list(base_ring)
+            diagonal_elements = [ x for x in base_elements if g(x) == x ]
+            number_diagonal_elements = len(diagonal_elements)
+            for weight1 in range((order-1)*entries_in_upper_half+1):
+                for iv2 in sage.combinat.integer_vector.IntegerVectors(weight1, entries_in_upper_half, max_part=(order-1)):
+                    for weight2 in range((number_diagonal_elements-1)*nrows+1):
+                        for dia in sage.combinat.integer_vector.IntegerVectors(weight2, nrows, max_part=(number_diagonal_elements-1)):
+                            iv = iv2.clone() # iv is going to be changed within the next loop, so we keep a copy iv2
+                            #construct upper half matrix
+                            matrix_entries = [] #entries of matrix with lower half 0
+                            length_of_row = nrows-1
+                            for r in range(nrows):
+                                zeros = [0]*(nrows - length_of_row - 1)
+                                row = zeros + [diagonal_elements[dia[r]]] + [base_elements[i] for i in iv[:length_of_row]]
+                                matrix_entries.extend(row)
+
+                                iv = iv[length_of_row:]
+                                length_of_row -= 1
+
+                            M = self(entries=matrix_entries)
+                            #make M symmetric
+                            make_symmetric(M)
+                            yield M
+
     def __getitem__(self, x):
         """
         Return a polynomial ring over this ring or the `n`-th element of this ring.