Codegolf the Hafnian
R, 150 142 127 119 bytes



function(A,N=nrow(A),k=1:(N/2)*2)sum(apply(gtools::permutations(N,N),1,function(r)prod(A[cbind(r[k-1],r[k])])))/prod(k)


Try it online!

Uses the same trick I discovered golfing down this answer to index the matrix P, and @Vlo suggested an approach to entirely remove the for loop for -6 bytes!

To create a new test case, you can do matrix(c(values,separated,by,commas,going,across,rows),nrow=2n,ncol=2n,byrow=T).

Explanation: (the code is the same; it uses an apply rather than a for loop but the logic is otherwise identical).

function(A){
N &lt;- nrow(A)                   #N = 2*n
k &lt;- 1:(N/2) * 2               #k = c(2,4,...,N) -- aka 2*j in the formula
P &lt;- gtools::permutations(N,N) #all N-length permutations of 1:N
for(i in 1:nrow(P))
 F &lt;- F + prod(A[cbind(P[i,k-1],P[i,k])]) # takes the product of all A_sigma(2j-1)sigma(2j) for fixed i and adds it to F (initialized to 0)
F / prod(k)                    #return value; prod(k) == n! * 2^n
}