OEIS

On this page I collect some Sage code that I used to compute entries in the Online Encyclopedia of Integer Sequences. In particular this includes some code that is too long to be included in the OEIS.

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 "\nThe matrix has permanent", p, ", computed in", cputime(t0), "seconds."

# Output

[1 0 0 1 1 0 0]

[0 1 1 0 1 0 0]

[1 0 1 0 0 1 0]

[0 1 0 1 0 1 0]

[0 0 1 1 0 0 1]

[1 1 0 0 0 0 1]

[0 0 0 0 1 1 1]

Does it correspond to a projective plane? True

This matrix has permanent 24, which was computed in 0.080005 seconds.

[1 0 0 0 1 0 0 0 1 1 0 0 0]

[0 0 1 1 0 0 0 1 0 1 0 0 0]

[0 1 0 0 0 1 1 0 0 1 0 0 0]

[1 0 0 0 0 1 0 1 0 0 1 0 0]

[0 1 0 1 0 0 0 0 1 0 1 0 0]

[0 0 1 0 1 0 1 0 0 0 1 0 0]

[1 0 0 1 0 0 1 0 0 0 0 1 0]

[0 1 0 0 1 0 0 1 0 0 0 1 0]

[0 0 1 0 0 1 0 0 1 0 0 1 0]

[0 0 0 1 1 1 0 0 0 0 0 0 1]

[0 0 0 0 0 0 1 1 1 0 0 0 1]

[1 1 1 0 0 0 0 0 0 0 0 0 1]

[0 0 0 0 0 0 0 0 0 1 1 1 1]

Does it correspond to a projective plane? True

This matrix has permanent 3852, which was computed in 0.072004 seconds.

[1 0 0 0 0 0 0 1 0 1 0 0 0 0 1 0 1 0 0 0 0]

[0 0 0 1 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0]

[0 1 0 0 0 0 1 0 1 0 0 0 0 0 0 1 1 0 0 0 0]

[0 0 1 0 0 1 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0]

[1 0 0 0 0 0 1 0 0 0 0 1 0 1 0 0 0 1 0 0 0]

[0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0]

[0 0 0 1 0 1 0 0 1 0 0 0 0 0 1 0 0 1 0 0 0]

[0 1 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 1 0 0 0]

[1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 1 0 0]

[0 1 0 0 1 0 0 0 0 0 0 1 0 0 1 0 0 0 1 0 0]

[0 0 1 0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 1 0 0]

[0 0 0 1 0 0 1 0 0 1 0 0 1 0 0 0 0 0 1 0 0]

[1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 0]

[0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 1 0]

[0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 1 0]

[0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 1 0]

[0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1]

[1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]

[0 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 1]

[0 0 0 0 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 1]

[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1]

Does it correspond to a projective plane? True

This matrix has permanent 18534400, which was computed in 13.964873 seconds.

A007878: 1, 2, 5, 16, 59, ...

Number of terms in the disriminant of a generic polynomial of degree n.

for N in range(1, 9):

Lvars = ['x'] + ['a'+str(i) for i in range(N)]

R = PolynomialRing(QQ, Lvars)

R.inject_variables(verbose=False)

P = sum([globals()["a"+str(i)]*x^i for i in range(N)]) + x^N

M = P.sylvester_matrix(diff(P, x), x)

print len(M.determinant().coefficients())

A055547, A055548, A055549

The code below computes the higher entries in these sequences quickly, by recursively adding rows and columns. Rows r_i, r_j and columns c_i, c_j for which the inner products <r_i, r_j> and <c_i, c_j> are not equal, can then be discarded at an early stage. To speed things up, the top three functions can be precompiled using Cython. This code works for each of the three sequences, depending on which of the lines starting "entries" is uncommented.

%cython

def ip(v1, v2):

return sum([v1[i]*v2[i] for i in range(len(v1))])

def ips(Lr, Lc):

flag = True

i = 0

while flag and i <= len(Lr) - 1:

if ip(Lr[i], Lr[-1]) != ip(Lc[i], Lc[-1]):

return False

i += 1

return True

def rc(Ltup, n):

m = len(Ltup)

r = [Ltup[i][n-m] for i in range(m-1)] + Ltup[-1][n-m:]

c = Ltup[-1][:n-m+1] + [ Ltup[m-2-i][n+2-m+2*i] for i in range(m-1) ]

c.reverse()

return r, c

%sage

entries = [ 0, 1] # A055547

#entries = [-1, 1] # A055548

#entries = [-1, 0, 1] # A055549

def NrMat(Ltup, Lr, Lc, n):

if ips(Lr, Lc):

m = len(Lc)

if m == n:

return 1

else:

N = 0

L0 = [entries for i in range(2*n - 1 - 2*m)]

for tup in cartesian_product_iterator(L0):

Ltup1 = Ltup + [list(tup)]

r, c = rc(Ltup1, n)

N += NrMat(Ltup1, Lr + [r], Lc + [c], n)

return N

return 0

print [NrMat([], [], [], n) for n in range(1,4)]