(PARI) mj(j, n) = (-1)^j*sum(k=0, n\3, binomial(n, 3*k+j));
a(n) = {m = matrix(3, 3); for (j=1, 3, m[1, j] = mj(j-1, n)); for (j=2, 3, m[2, j] = m[1, j-1]); m[2, 1] = m[1, 3]; for (j=2, 3, m[3, j] = m[2, j-1]); m[3, 1] = m[2, 3]; matdet(m); } \\
Michel Marcus, Jul 26 2017
(Python)
from sympy.matrices import Matrix
from sympy import binomial
def mj(j, n):
return (-1)**j*sum(binomial(n, 3*k + j) for k in range(n//3 + 1))
def a(n):
m=Matrix(3, 3, [0]*9)
for j in range(3):m[0, j]=mj(j, n)
for j in range(1, 3):m[1, j]=m[0, j - 1]
m[1, 0]=m[0, 2]
for j in range(1, 3):m[2, j] = m[1, j - 1]
m[2, 0]=m[1, 2]
return m.det()