"""
Build as follows:

# pyrex mult2.pyx
# ls -lh mult2.c
-rw-r--r--  1 was was 28K Nov 11 14:09 mult2.c
# gcc -shared -O3 -fPIC -I/home/was/local/include/python2.4 mult2.c -o mult2.so
-rwxr-xr-x  1 was was 21K Nov 11 14:13 mult2.so

(Your Python include directory is probably different than mine.)
"""

cdef extern from "Python.h":
    void* PyMem_Malloc(int)
    void PyMem_Free(void *p)

cdef class Matrix:
    cdef int *entries
    cdef int p, n

    def __new__(self, int p, int n, entries=None):
        self.n = n
        self.p = p
        self.entries = <int*> PyMem_Malloc(sizeof(int)*n*n)

    def __dealloc__(self):
        PyMem_Free(self.entries)

    def __init__(self, int p, int n, entries=None):
        """
        Create a matrix over a finite field.
        INPUT:
            p -- prime
            n -- positive integer
            entries -- entries of the matrix (defaults to None, which means 0 matrix).
        """
        cdef int i, j, k, x
        if entries != None:
            for i from 0 <= i < self.n:
                for j from 0 <= j < self.n:
                    k = i*self.n + j
                    x = entries[k] % p
                    if x<0: x = x + p
                    self.entries[k] = x 
        else:
            for i from 0 <= i < self.n:
                for j from 0 <= j < self.n:
                    self.entries[i*self.n + j] = 0
            
    def __repr__(self):
        cdef int i, j, n
        s = ""
        n = self.n
        for i from 0 <= i < n:
            for j from 0 <= j < n:
                s = s + "%s, "%self.entries[n*i + j]
            s = s + "\n"
        return s
    
    cdef Matrix __mul__(Matrix self, Matrix B):
        cdef int s, i, j, k, n
        cdef Matrix ans
        ans = Matrix(self.p, self.n)
        n = self.n
        for i from 0 <= i < n:
            for j from 0 <= j < n:
                # The i,j entry of the product
                s = 0
                for k from 0 <= k < n:
                    s = (s + (self.entries[i*n+k] * B.entries[k*n+j])) % self.p
                    if s < 0: s = s + self.p
                    ans.entries[i*n+j] = s
        return ans