Ticket #6452: trac_6452-ring-codes.patch

module_list.py

@@ -154 +154 @@
     Extension('sage.coding.binary_code',
               sources = ['sage/coding/binary_code.pyx']),
 
+    Extension('sage.coding.ring_code', 
+              sources = ['sage/coding/ring_code.pyx']), 
+
     ################################
     ## 
     ## sage.combinat

sage/coding/all.py

@@ -55 +55 @@
 
 from sd_codes import self_dual_codes_binary
 
+from ring_code import RingCode 

new file sage/coding/ring_code.pyx

@@ +1 @@
+"""
+Codes over rings of the form ZZ/mZZ.
+
+This module constructs codes over rings of the form ZZ/mZZ, that is, submodules
+of FreeModule(IntegerModRing(m), n).
+
+    AUTHORS:
+        -- Cesar Agustin Garcia-Vazquez
+        -- Carlos A. Lopez-Andrade
+        -- David Joyner
+
+TODO:
+  - Rewrite nearest ngbr decoder completely using pure cython with no GF(q) classes
+    Coerce answer at last step.
+  - Generalize this code to allow for other finite rings.
+
+"""
+
+    #*****************************************************************************
+    #       Copyright (C) 2008 William Stein <wstein@gmail.com>
+    #                     2008 Cesar A. Garcia-Vazquez <cesarnda@gmail.com>
+    #                     2008 Carlos A. Lopez-Andrade <calopez@cs.buap.mx>
+    #
+    #  Distributed under the terms of the GNU General Public License (GPL),
+    #  version 2 or later (at your option).
+    #
+    # This code is distributed in the hope that it will be useful,
+    # but WITHOUT ANY WARRANTY; without even the implied warranty of
+    # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+    # General Public License for more details. 
+    #
+    #  The full text of the GPL is available at:
+    #
+    #                  http://www.gnu.org/licenses/
+    #*****************************************************************************
+
+import sage.modules.free_module as fm
+import sage.modules.module as module
+from sage.rings.integer_mod_ring import IntegerModRing
+
+#cython
+cdef extern from *:
+    void* malloc(int)
+    int memset(void*, int, int)
+    void free(int*)
+    double ceil( double )
+    double floor( double )
+    double pow( double, double)
+
+from sage.modules.free_module_element import vector
+from sage.modules.free_module import FreeModule
+
+cdef class RingCode: #(module.Module):
+    r"""
+    This class generates code over a ring ZZ/mZZ generated by a "generator matrix".
+
+    CALL FORMATS:
+        1. CS = RingCode(matrixGen) 
+
+    INPUT:
+        matrixGen -- matrix to be used as a generator matrix
+
+    OUTPUT:
+        list -- list to be used as a set.
+
+    EXAMPLES:
+        sage: M = Matrix(IntegerModRing(25), [[1, 1, 6],[1, 6, 1],[6, 1, 1]])
+        sage: CS = RingCode(M)
+        sage: CS
+        (3, 625, 1)-code over the Ring of integers modulo 25
+        sage: M = Matrix(IntegerModRing(25), [[6, 6, 6],[6, 6, 6],[6, 6, 6]])
+        sage: CS = RingCode(M)
+        sage: CS
+        (3, 25, 3)-code over the Ring of integers modulo 25
+        sage: M = Matrix(IntegerModRing(9), [[5,4,8,1],[1,5,4,8],[8,1,5,4],[4,8,1,5]])
+        sage: CS = RingCode(M)
+        sage: CS
+        (4, 81, 2)-code over the Ring of integers modulo 9 
+        sage: M = Matrix(IntegerModRing(25), [[1, 0, 4],[0, 4, 1],[4, 1, 0]])
+        sage: CS = RingCode(M)
+        sage: CS
+        (3, 3125, 1)-code over the Ring of integers modulo 25 
+        sage: M = Matrix(IntegerModRing(9), [[5,1,2,1],[1,5,1,2],[2,1,5,1],[1,2,1,5]])
+        sage: CS = RingCode(M)
+        sage: CS
+        (4, 81, 2)-code over the Ring of integers modulo 9 
+        sage: M = Matrix(IntegerModRing(9), [[5,0,0,0,4],[4,5,0,0,0],[0,4,5,0,0],[0,0,4,5,0],[0,0,0,4,5]])
+        sage: CS = RingCode(M)
+        sage: CS
+        (5, 6561, 2)-code over the Ring of integers modulo 9 
+        
+    NOTES:
+        This function uses a d-heap to store the codewords, this way, it is as fast as Magma.
+
+    AUTHORS:
+        -- Cesar A. Garcia-Vazquez
+        -- Carlos A. Lopez-Andrade
+    """
+    cdef int codewords      #Cardinality of the set
+    cdef int size 
+    cdef int sizeCols       #number of columns of the matrix
+    cdef int modulus        
+    cdef int minimum        #minimum Hamming distance
+    cdef int count          #to iterate on self
+    cdef long sizeArray
+    cdef int** heap_matrix
+    cdef _base_ring, ambient_module
+    cdef gen_matrix, gen_mat, minimum_weight_codeword, minimum_distance, mincodeword
+    cdef next, spanning_codewords, characteristic, _length
+    cdef list codeSet
+
+    def __new__(self, matrixGen):
+        cdef int size = <int>matrixGen.nrows()
+        cdef int sizeCols = <int>matrixGen.ncols()
+        #cdef int length = sizeCols
+        cdef int lamb = 2*size
+        self.size = size
+        self.sizeCols = sizeCols
+        self._length = sizeCols
+        self._base_ring = matrixGen.base_ring()
+        self.modulus = <int>matrixGen.base_ring().order()
+        if sizeCols <= size + lamb:
+            self._use_heap_matrix(matrixGen)
+        self.gen_matrix = matrixGen
+
+    cdef void _use_heap_matrix(self, matrixGen):
+        cdef int modulus = self.modulus
+        cdef int size = self.size
+        cdef int sizeCols = self.sizeCols
+        cdef int* vectorGen =  <int*>malloc(sizeof(int) * size)
+        cdef int* newVector =  <int*>malloc(sizeof(int) * sizeCols)
+        
+        #Number of columns of the matrix
+        cdef long sizeArray = (<int>pow(modulus,sizeCols) - 1)/(modulus - 1) - 1    
+        self.sizeArray = sizeArray
+        self.count = 0
+        #Limit used to determinate how many codewords are going to be generated.
+        cdef int limit = <int>pow(modulus,size)
+        
+        cdef int localLimit     #Local limit to modify the values of the vector that multiplies the matrixGen.
+        cdef int localIndex     #Local index to determinate which index is going to change.
+        cdef int value      #Takes the first element of the vector to check if the second has to change.
+        cdef int element        #Element that must be in the matrixGen.
+        cdef int head       #Second element of the vector generated by multiplying the matrixGen and the vector.
+        cdef int index1     #Index used in some fors.
+        cdef int index2     #Another index used in some fors
+        cdef int index3
+        cdef int mul
+             
+        cdef long index
+        cdef int** heap_matrix = <int**>malloc(sizeof(int*)*modulus)
+        
+        for index1 from 0 <= index1 < modulus:
+            heap_matrix[index1]= <int*>malloc(sizeof(int)*sizeArray)
+            memset(heap_matrix[index1],0,sizeof(int)*sizeArray)
+        
+        cdef int** matrixGenerator = <int**>malloc(sizeof(int*)*size)
+        
+        
+        for index1 from 0 <= index1 < size:
+            matrixGenerator[index1]= <int*>malloc(sizeof(int)*sizeCols)
+            memset(matrixGenerator[index1],0,sizeof(int)*sizeCols)
+        
+        for index1 from 0 <= index1 < size:
+            for index2 from 0 <= index2 < sizeCols:
+                matrixGenerator[index1][index2] = matrixGen[index1][index2]
+        
+        for index1 from 0 <= index1 < modulus:
+            for index2 from 0 <= index2 < sizeArray:
+                heap_matrix[index1][index2] = -1
+        
+        for index1 from 0 <= index1 < size:
+            vectorGen[index1] = 0
+        
+        vectorGen[0] = -1
+        
+        #main for for obtaining the whole set
+        for index from 0 <= index < limit:          
+            localLimit = <int>ceil( <double>index / <double>modulus )
+            value = vectorGen[0] + 1
+            vectorGen[0] = value % modulus
+            localIndex = 1
+        
+            for index2 from 0 <= index2 < localLimit:
+                if value == modulus:
+                    value = vectorGen[localIndex] + 1
+                    vectorGen[localIndex] = value % modulus
+                    localIndex = localIndex + 1
+                else:
+                    break
+        
+            for index1 from 0 <= index1 < sizeCols:
+                newVector[index1] = 0
+            # multiplying vector and matrix
+            for index1 from  0 <= index1 < sizeCols:
+                for index2  from 0 <= index2 < size:
+                    mul = matrixGenerator[index2][index1] * vectorGen[index2]
+                    newVector[index1] = newVector[index1] + mul
+                newVector[index1] = newVector[index1] % modulus
+            head = newVector[0]
+            element = -1    #To start filling
+            index1 = 1
+        
+            #Filling While
+            while index1 < sizeCols:
+                element = (element + 1)* modulus + newVector[index1]
+                heap_matrix[head][element] = newVector[index1]
+                index1 = index1 + 1
+        #End of main for
+        
+        self.heap_matrix = heap_matrix
+        free(vectorGen)
+        self._remove_zero()
+        self._generate_codewords()
+        
+    cdef _remove_zero(self):
+        cdef int sizeCols = self.sizeCols
+        cdef int modulus = self.modulus
+        cdef int* vectorGen =  <int*>malloc(sizeof(int) * sizeCols)
+        cdef int index
+        cdef int element = -1
+        cdef int head = 0
+        index = 1
+        while index < sizeCols:
+            element = (element + 1)* modulus
+            self.heap_matrix[head][element] = -1
+            index = index + 1
+        
+    cdef _generate_codewords(self):
+        cdef int index1
+        cdef int index2
+        cdef int index3
+        cdef int index
+        cdef int head
+        cdef int codewords = 0
+        cdef int sizeCols = self.sizeCols
+        cdef int modulus = self.modulus
+        cdef int limit = sizeCols - 1
+        cdef long sizeArray = self.sizeArray
+        cdef long beginArray = sizeArray - <int>pow(modulus,(sizeCols - 1))
+        cdef long endArray = <long>pow(modulus,(sizeCols - 2))
+        cdef long index2Start = beginArray
+        cdef int minimum = sizeCols #minimum distance
+        cdef int localMin    #temp minium distance
+        cdef int* vectorGen =  <int*>malloc(sizeof(int) * sizeCols) 
+        cdef list codeSet = []
+        #cdef int** heap_matrix = self.heap_matrix
+        for index1 from 0 <= index1 < modulus:
+            for index2 from 0 <= index2 < endArray:
+                for index from 0 <= index < modulus:
+                    if self.heap_matrix[index1][beginArray + index] != -1:
+                        addVector = ([])
+                        #mincodeword = ([])
+                        codewords = codewords + 1
+                        vectorGen[limit] = index
+                        localMin = 0
+                        if index != 0:
+                            localMin = localMin + 1
+                        head = beginArray
+                        #print "Before searching for fathers"
+                        for index3 from limit > index3 >= 1:
+                            head = <int> floor( (head - 1) / modulus)
+                            vectorGen[index3] = self.heap_matrix[index1][ head ]
+                            if vectorGen[index3] != 0:
+                                localMin = localMin + 1
+                            head = head - vectorGen[index3]
+                        vectorGen[0] = index1
+                        if index1 != 0:
+                            localMin = localMin + 1
+                        #print "Every father found"
+                        for index3 from 0 <= index3 < sizeCols:
+                            addVector.append(vectorGen[index3])
+                        addVector = vector(IntegerModRing(modulus),addVector)
+                        codeSet.append(addVector)
+                        #print "Every vector in the set"
+                        if localMin < minimum:              # Compares the min dist with a probable new min dist
+                            minimum = localMin              # A new minimum distance has been found
+                            #print "A new minimum distance has been found", minimum, addVector
+                            mincodeword = addVector
+                beginArray = beginArray + modulus
+            beginArray = index2Start
+        #free(vectorGen)
+        addVector = ([])
+        for index3 from 0 <= index3 < sizeCols:
+            addVector.append(0) 
+        codeSet.append(vector(IntegerModRing(modulus), addVector))
+        self.minimum = minimum
+        self.mincodeword = mincodeword
+        self.codewords = codewords + 1
+        self.codeSet = codeSet
+
+    def base_ring(self):
+        """
+        Returns base ring.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[1, 1, 6],[1, 6, 1],[6, 1, 1]])
+            sage: C = RingCode(M)
+            sage: C.base_ring()
+            Ring of integers modulo 12
+
+        """
+        return self._base_ring
+
+    def ambient_module(self):
+        """
+        Returns ambient module.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[1, 1, 6],[1, 6, 1],[6, 1, 1]])
+            sage: C = RingCode(M)
+            sage: C.ambient_module()
+            Ambient free module of rank 3 over Ring of integers modulo 12
+        """
+        return FreeModule(self.base_ring(),self._length)
+
+    def characteristic(self):
+        """
+        Returns characteristic of base ring.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[1, 1, 6],[1, 6, 1],[6, 1, 1]])
+            sage: C = RingCode(M)
+            sage: C.characteristic()
+            12
+        """
+        return (self.base_ring()).characteristic()
+
+    def __contains__(self,v):
+        C = self.codeSet
+        return C.__contains__(v)
+
+    def gen_mat(self):
+        """
+        Returns a generator matrix.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[1, 1, 6, -1],[1, 6, 1, 2],[6, 1, 1, 0]])
+            sage: C = RingCode(M)
+            sage: C.gen_mat()
+            [ 1  1  6 11]
+            [ 1  6  1  2]
+            [ 6  1  1  0]
+
+        """
+        return self.gen_matrix
+
+    def __getitem__(self, item):
+        """
+        Iterator method.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(6), [[2, 2, 2],[2, 2, 2],[2, 2, 2]])
+            sage: CS = RingCode(M); CS
+            (3, 3, 3)-code over the Ring of integers modulo 6
+            sage: [c for c in CS]
+            [(2, 2, 2), (4, 4, 4), (0, 0, 0)]
+        """
+        return self.codeSet[item]
+    
+    def __latex__(self):
+        s = 'Linear Code over the %s \n(%d, %d, %d)'%(self.base_ring, self.sizeCols, self.codewords, self.minimum)
+        return s
+
+    def length(self):
+        """
+        Returns the usual length of code.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[1, 1, 6],[1, 6, 1],[6, 1, 1]])
+            sage: C = RingCode(M)
+            sage: C.length()
+            3
+        """
+        return self._length
+
+    def list(self):
+        """
+        Returns list of all codewords.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(2), [[1, 1, 0],[1, 0, 1],[0, 1, 1]])
+            sage: C = RingCode(M)
+            sage: C
+            (3, 4, 2)-code over the Ring of integers modulo 2
+            sage: C.list()
+            [(0, 1, 1), (1, 0, 1), (1, 1, 0), (0, 0, 0)]
+        """
+        return self.codeSet
+
+    def minimum_distance(self):
+        """
+        Returns the usual minimum distance.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[0, 1, 6, -1],[1, 6, 1, 2],[6, 1, 1, 0]])
+            sage: C = RingCode(M)
+            sage: C
+            (4, 1728, 2)-code over the Ring of integers modulo 12
+            sage: C.minimum_distance()
+            2
+
+        """
+        return self.minimum
+
+    def minimum_weight_codeword(self):
+        """
+        Returns a non-zero codeword of minimum weight.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[0, 1, 6, -1],[1, 6, 1, 2],[6, 1, 1, 0]])
+            sage: C = RingCode(M)
+            sage: C
+            (4, 1728, 2)-code over the Ring of integers modulo 12
+            sage: c = C.minimum_weight_codeword(); c
+            (0, 1, 0, 5)
+            sage: c in C
+            True
+
+        """
+        return self.mincodeword
+
+    def next(self):
+        """
+        Returns an iterator fr the codewords.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[0, 1, 6, -1],[1, 6, 1, 2],[6, 1, 1, 0]])
+            sage: C = RingCode(M)
+            sage: C[0]
+            (0, 11, 1, 1)
+            sage: C.next()
+            (0, 11, 1, 1)
+            sage: C[1]
+            (0, 11, 2, 2)
+            sage: C.next()
+            (0, 11, 2, 2)
+
+        """
+        cdef int c
+        if self.count >= len(self.codeSet):
+            self.count = 0
+            raise StopIteration
+        else:
+            c = self.count
+            self.count = c + 1
+            return self.codeSet[c]
+
+    def __repr__(self):
+        s = '(%d, %d, %d)-code over the %s'%(self.length(), self.codewords, self.minimum, self.base_ring())
+        return s
+        
+    def spanning_codewords(self):
+        """
+        Returns a list of codewords which span the code.
+
+        EXAMPLES:
+            sage: M = Matrix(IntegerModRing(12), [[0, 1, 6, -1],[1, 6, 1, 2],[6, 1, 1, 0]])
+            sage: C = RingCode(M)
+            sage: C.spanning_codewords()
+            [(0, 1, 6, 11), (1, 6, 1, 2), (6, 1, 1, 0)]
+
+        """
+        return self.gen_matrix.rows()
+ No newline at end of file