# A000794: 1, 2, 24, 3852, 18534400, ....
#
# For any integer q, this Sage code computes the incidence matrix I
# of a projective plane over a finite field of order q, using a method
# described in the paper "Incidence Matrices of Projective Planes and
# of Some Regular Bipartite Graphs of Girth 6 with Few Vertices" by
# Balbuena, see http://dx.doi.org/10.1137/070688225 .
# 
# The Sage function for computing the permanent quickly runs out of
# memory, so we provide a custom function for computing the permanent
# of this incidence matrix, which is slower but uses less memory.

%cython
from sage.matrix.matrix2 import _binomial

def my_combinations_iterator(mset, n):
    items = mset
    if n is None:
        n = len(items)
    for i in xrange(len(items)):
        v = items[i:i+1]
        if n == 1:
            yield v
        else:
            rest = items[i+1:]
            for c in my_combinations_iterator(rest, n-1):
                yield v + c

def my_permanent(M):
    cdef Py_ssize_t m, n, r
    cdef int sn

    perm = 0
    m = M.nrows()
    n = M.ncols()

    for r from 1 <= r < m+1:
        print r,

        s = 0
        for cols in my_combinations_iterator(range(n), r):
            s += M.prod_of_row_sums(cols)

        if (m - r) % 2 == 0:
            sn = 1
        else:
            sn = -1
            
        perm = perm + sn * _binomial(n-r, m-r) * s

    return perm
    
%sage
def pos_mat(M, Lk):
    return [ flatten([[M[i][j].count(k) for j in range(len(M[0]))] for k in Lk]) for i in range(len(M))]

def to_int(expr):
    if K.is_prime_field():
        return ZZ(expr)
    else:
        return int(expr.int_repr())

def is_proj_plane(I):
    flag = True
    for (i,j) in combinations_iterator(range(I.nrows()), 2):
        flag = flag * [ int(I[i][k] == I[j][k] == 1) for k in range(I.ncols()) ].count(1) == 1
        flag = flag * [ int(I[k][i] == I[k][j] == 1) for k in range(I.nrows()) ].count(1) == 1

    return flag

for q in [2,3,4]:
    t0 = cputime()
    K = GF(q,name='a')
    LK = [x for x in K]

    LSigma = []
    for u0 in LK:
        if u0 != K(0):
            LSigma.append([[ [ to_int(u0*j + i - K(1)) + 1 ] for j in LK] for i in LK])

    qIq = [[ (i == j)*range(1, q+1) for j in range(q)] for i in range(q)]
    Oq = [[ [i+1]  for j in range(q)] for i in range(q)]

    A = []
    for M0 in LSigma + [qIq, Oq]:
        A += list(pos_mat(M0, [ to_int(x - K(1)) + 1 for x in LK ]))

    Oqq1  = [[ [i+1]  for j in range(q)] for i in range(q+1)]
    Oqq1T = list(matrix(pos_mat(Oqq1, range(1,q+2))).transpose())
    I = matrix([A[i] + list(Oqq1T[i]) for i in range(len(A))] + [(q^2)*[0] + (q+1)*[1]])
    print I.str()
    print "Does it correspond to a projective plane? ", is_proj_plane(I)
    
    p = I.permanent()
    #p = my_permanent(I)
    print "\nThis matrix has permanent", p, ", which was computed in", cputime(t0), "seconds."