[Question]
  [
The challenge is to write codegolf for the [Hafnian of a matrix](https://en.wikipedia.org/wiki/Hafnian).
The Hafnian of an \$2n \times 2n\$ symmetric matrix \$A\$ is defined as:
$$\operatorname{haf}(A) = \frac 1 {n!2^n} \sum\_{\sigma \in S\_{2n}} \prod^n\_{j=1} A\_{\sigma(2j-1),\sigma(2j)}$$
Here \$S\_{2n}\$ represents the set of all permutations of the integers from \$1\$ to \$2n\$, that is \$[1, 2n]\$.
The Wikipedia link talks about adjacency matrices but your code should work for any real valued symmetric input matrices.
For those interested in applications of the Hafnian, the [mathoverflow](https://mathoverflow.net/questions/273324/applications-of-hafnians/273325) link discusses some more.
Your code can take input however it wishes and give output in any sensible format but please include in your answer a full worked example including clear instructions for how to supply input to your code.
The input matrix is always square and will be at most 16 by 16. There is no need to be able to handle the empty matrix or matrices of odd dimension.
## Reference implementation
Here is some example python code from Mr. Xcoder.
```
from itertools import permutations
from math import factorial
def hafnian(matrix):
    my_sum = 0
    n = len(matrix) // 2
    for sigma in permutations(range(n*2)):
        prod = 1
        for j in range(n):
            prod *= matrix[sigma[2*j]][sigma[2*j+1]]
        my_sum += prod
    return my_sum / (factorial(n) * 2 ** n)
print(hafnian([[0, 4.5], [4.5, 0]]))
4.5
print(hafnian([[0, 4.7, 4.6, 4.5], [4.7, 0, 2.1, 0.4], [4.6, 2.1, 0, 1.2], [4.5, 0.4, 1.2, 0]])
16.93
print(hafnian([[1.3, 4.1, 1.2, 0.0, 0.9, 4.4], [4.1, 4.2, 2.7, 1.2, 0.4, 1.7], [1.2, 2.7, 4.9, 4.7, 4.0, 3.7], [0.0, 1.2, 4.7, 2.2, 3.3, 1.8], [0.9, 0.4, 4.0, 3.3, 0.5, 4.4], [4.4, 1.7, 3.7, 1.8, 4.4, 3.2]])
262.458
```
---
The wiki page has now (March 2 2018) been updated by ShreevatsaR to include a different way of calculating the Hafnian. It would be very interesting to see this golfed.
  
[Answer]
# [R](https://www.r-project.org/), ~~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!](https://tio.run/##PY6xDoIwGIR3n8Lx/82P2BKjaWTgBTq5EQYoYgjQNqVEefpaHBguXy6Xu5wL3fGRhG7RyvdGQ0Ey1858oEAaciZAphxPHOdlgtracYW3N2achbAvNy2@3lozSJJIjPYZh9aZFopSNb1uwZVDwiqKqLBCxPSfDhg6mGrv@i8oyOhyvkbdKEPixBEPe8rEnbZXeYQyY0SzbvaJGH4 "R – Try It Online")
Uses the same trick I discovered golfing down [this answer](https://codegolf.stackexchange.com/a/156888/67312) to index the matrix `P`, and [@Vlo](https://codegolf.stackexchange.com/users/30693/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 <- nrow(A)                   #N = 2*n
k <- 1:(N/2) * 2               #k = c(2,4,...,N) -- aka 2*j in the formula
P <- gtools::permutations(N,N) #all N-length permutations of 1:N
for(i in 1:nrow(P))
 F <- 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
}
```
[Answer]
# [Pyth](https://pyth.readthedocs.io), 24 bytes
```
sm*Fmc@@Qhkek2d{mScd2.pU
```
**[Try it here!](http://pyth.herokuapp.com/?code=sm%2aFmc%40%40Qhkek2d%7BmScd2.pU&input=%5B%5B0%2C+4.7%2C+4.6%2C+4.5%5D%2C+%5B4.7%2C+0%2C+2.1%2C+0.4%5D%2C+%5B4.6%2C+2.1%2C+0%2C+1.2%5D%2C+%5B4.5%2C+0.4%2C+1.2%2C+0%5D%5D&debug=0)**
---
### Old version, 35 bytes
```
*c1**FK/lQ2^2Ksm*Fm@@Q@[[email protected]](/cdn-cgi/l/email-protection)
```
[Try it here!](http://pyth.herokuapp.com/?code=%2ac1%2a%2aFK%2FlQ2%5E2Ksm%2aFm%40%40Q%40dtyk%40dykK.pU&input=%5B%5B3%2C0.7%5D%2C%5B0.5%2C3%5D%5D&debug=0)
[Answer]
# [Stax](https://github.com/tomtheisen/stax), 23 22 19 17 [bytes](https://github.com/tomtheisen/stax/blob/master/docs/packed.md#packed-stax)
```
ü;Y╙◘▌Φq↓ê²╧▐å↑┌C
```
[Run and debug it online](https://staxlang.xyz/#c=%C3%BC%3BY%E2%95%99%E2%97%98%E2%96%8C%CE%A6q%E2%86%93%C3%AA%C2%B2%E2%95%A7%E2%96%90%C3%A5%E2%86%91%E2%94%8CC&i=%5B%5B3.0%2C0.7%5D%2C%5B0.5%2C3.0%5D%5D%0A%5B%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%2C%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%2C%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%2C%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%2C+%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%2C%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%2C%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%2C%5B1.0%2C2.0%2C3.0%2C4.0%2C5.0%2C6.0%2C7.0%2C8.0%5D%5D&a=1&m=2)
The corresponding ascii representation of the same program is this.
```
%r|TF2/{xsE@i^H/m:*+
```
The program suffers from some floating point rounding error. In particular, it reports `33673.5000000011` instead of `33673.5`. But I think the accuracy is acceptable given that this program operates on floating point values. It's also very slow, taking almost a minute for the example inputs on this machine.
```
%                             get size of matrix
 r|T                          get all permutations of [0 ... size-1]
    F                         for each, execute the rest of the program
     2/                       get consecutive pairs
       {        m             map each pair... 
        xsE@                      the matrix element at that location
            i^H/                  divided by 2*(i+1) where i=iteration index
                 :*           product of array
                   +          add to running total
```
[Answer]
# [05AB1E](https://github.com/Adriandmen/05AB1E), 21 [bytes](https://github.com/Adriandmen/05AB1E/wiki/Codepage)
```
ā<œε2ô{}Ùεε`Isèsè;]PO
```
[Try it online!](https://tio.run/##MzBNTDJM/f//SKPN0cnnthod3lJde3jmua3ntiZ4Fh9eAUTWsQH@//9HRxvoKJjomYMIMxBhGqujEA0WAEoY6RkCaT0TiJgZTEBHwVDPCCJmCpYHCwBZsbEA "05AB1E – Try It Online")
---
### Old version, 32 [bytes](https://github.com/Adriandmen/05AB1E/wiki/Codepage)
```
āœvIg;©Lε·UIyX<èèyXèè}P}Oθ®!/®o/
```
[Try it online!](https://tio.run/##MzBNTDJM/f//SOPRyWWe6daHVvqc23poe6hnZYTN4RWHV1RGgMjagFr/czsOrVPUP7QuX////@hoYx0DPfNYnWgDPVMd49hYAA "05AB1E – Try It Online")
**How it works?**
```
āœvIg;©Lε·UIyX<èèyXèè}P}Oθ®!/®o/ – Full program. Argument: A matrix M.
ā                                – The range [1 ... len(M)].
 œ                               – Permutations.
  v                    }         – Iterate over the above with a variable y.
   Ig;©                          – Push len(M) / 2 and also store it in register c.
       Lε            }           – For each integer in the range [1 ... ^]:
         ·U                      – Double it and store it in a variable X.
            yX<                  – Push the element of y at index X-1.
           I   è                 – And index with the result into M.
                yXè              – Push the element of y at index X.
                   è             – And index with the result into ^^.
                      P          – Take the product of the resulting list.
                        O        – Sum the result of the mapping.
                         θ       – And take the last element*.
                          ®!     – Take the factorial of the last item in register c.
                             ®o  – Raise 2 to the power of the last item in register c.
                            /  / – And divide the sum of the mapping accordingly.
* – Yeah, this is needed because I mess up the stack when pushing so many values in the loop and not popping correctly ;P
```
[Answer]
# [R](https://www.r-project.org/), ~~84~~ 78 bytes
```
h=function(m)"if"(n<-nrow(m),{for(j in 2:n)F=F+m[1,j]*h(m[v<--c(1,j),v]);F},1)
```
[Try it online!](https://tio.run/##XU/LCsIwELz7FcFTVtOQl1pf1/6EeJBCsUJTKFoF8dvrbrtVFMKG2RlmZpuuO@@LW8yvZR1lBdOymMq4S2JT3xGqZ1E38iLKKNwmQrbP5tXBqstxdpbVod0lSS4RgmqPsM1eykKHxOnalA@ZS6NE0AseBpRw@AAm/5IVjeVHjBDXTlv8dWBugEpY7UZD4npI1qT7tbbak9KOIm1orGkXmAhEOApkSe84QjcUW38qooHv2@lvkRXpHBGedilncHHDhBkuCzz6DD8Epbzz2uEZeOkSoHsD "R – Try It Online")
Edit: Thanks to Vlo for -6 bytes.
It seems that everybody here is implementing the standard reference algorithm with permutations, but I tried to take advantage of the community knowledge gained in the [related challenge](https://codegolf.stackexchange.com/questions/157049/calculate-the-hafnian-as-quickly-as-possible), which is basically the same task targeted for fastest code instead of golf.
It turns out that for a language that is good at slicing matrices (like R), the recursive algorithm: `hafnian(m) = sum(m[i,j] * hafnian(m[-rows and columns at i,j])` is not only faster, but also quite golfy. Here is the ungolfed code:
```
hafnian<-function(m)
{
    n=nrow(m)
    #Exits one step earlier than golfed version
    if(n == 2) return(m[1,2])
    h = 0
    for(j in 2:n) {
        if(m[1,j] == 0) next
        h = h + m[1,j] * hafnian(m[c(-1,-j),c(-1,-j)])
    }
    h
}
```
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 19 bytes
```
LŒ!s€2Ṣ€QḅL_LịFHP€S
```
[Try it online!](https://tio.run/##y0rNyan8/9/n6CTF4kdNa4we7lwEpAIf7mj1ifd5uLvbzSMAyA/@f7gdSEX@/x/NpRAdbaCjYKJnGqujEA2kdBQMYmN14MLmIMIMSQFQAChhpGcIpPVMIGJmMAEdBUM9I4RBeiZgAYSRhnrGIKMMYcJ6BiDCEiQGNcoQxDYCmWcOVwQ2xRwkbwiXMoHoArOAhhhD5MEGghWBpYxALGOQnYZ6FhB5S6iBUF3GIK4pkv0Qu8AGgnWBpUBcI5AXYgE "Jelly – Try It Online")
### Alternate version, 15 bytes, postdates challenge
```
LŒ!s€2Ṣ€QœịHP€S
```
Jelly *finally* got n-dimensional array indexing.
[Try it online!](https://tio.run/##TY89DsIwDIV3ThH2yKqTtKU3YGAAMaEoIwvqxsTKwsBWbsDMBRAj6kHKRYLt9G@Jn9@zPyunY11fYty0zfL8u75M935S2bWP7nNfb0nu4/dG5RCjXyjvM60c5EErT0WrLAQ92iU/xWyADAoMIFVwySsGQysEM4HAiTEhESyjcLAh46dir0cha8O8chwSSsk5jpFLW6IIYlMuQBmSyLCyfBNhlfKqB/Zbltt8dj/dEqBsScSt4S@EPw "Jelly – Try It Online")
### How it works
```
LŒ!s€2Ṣ€QœiHP€S  Main link. Argument: M (matrix / 2D array)
L                Take the length, yielding 2n.
 Œ!              Generate all permutations of [1, ..., 2n].
   s€2           Split each permutation into pairs.
      Ṣ€         Sort the pair arrays.
        Q        Unique; deduplicate the array of pair arrays.
                 This avoids dividing by n! at the end.
           H     Halve; yield M, with all of its elements divided by 2.
                 This avoids dividing by 2**n at the end.
         œị      At-index (n-dimensional); take each pair of indices [i, j] and
                 yield M[i][j].
            P€   Take the product the results corresponding the same permutation.
              S  Take the sum of the products.
```
The 19-byte version works in a similar fashion; it just has to implement `œị` itself.
```
...ḅL_LịFH...    Return value: Array of arrays of index pairs. Argument: M
    L            Length; yield 2n.
   ḅ             Convert each pair of indices [i, j] from base 2n to integer,
                 yielding ((2n)i + j).
     _L          Subtract 2n, yielding ((2n)(i - 1) + j).
                 This is necessary because indexing is 1-based in Jelly, so the
                 index pair [1, 1] must map to index 1.
        F        Yield M, flattened.
       ị         Take the indices to the left and get the element at these indices
                 from the array to the right.
         H       Halve; divide all retrieved elements by 2.
```
[Answer]
# [C (gcc)](https://gcc.gnu.org/), ~~288~~ ~~285~~ ~~282~~ ~~293~~ ~~292~~ ~~272~~ ~~271~~ 263 bytes
* Saved three bytes by fiddling with two post-increments and for loop placement.
* Saved three bytes by fiddling with anothter post-increment, moving both variable initializations before the branch -- golfed `if(...)...k=0...else...,j=0...` to `if(k=j=0,...)...else...` -- and performed an index shift.
* Required eleven bytes by supporting `float` matrices.
* Saved a byte thanks to [Mr. Xcoder](https://codegolf.stackexchange.com/users/59487/mr-xcoder); golfing `2*j+++1` to `j-~j++`.
* Saved twenty bytes by removing a superfluous `int` variable type declaration and not using a factorial function but instead calculating the factorial value by using an already existing for loop.
* Saved a byte by golfing `S=S/F/(1<<n);` to `S/=F*(1<<n);`.
* Saved eight bytes thanks to [ceilingcat](https://codegolf.stackexchange.com/users/52904/ceilingcat).
```
float S,p,F;j,i;s(A,n,P,l,o,k)float*A;int*P;{if(o-l)for(k=0;k<l;s(A,n,P,l,o+1))P[o]=k++;else{for(p=j=l;j--;)for(i=l;i--;)p-=P[j]==P[i];if(!p){for(F=p=1,j=0;j<n;F*=j)p*=A[P[2*j]*2*n+P[j+++j]];S+=p;}}}float h(A,n)float*A;{int P[j=2*n];S=0;s(A,n,P,j,0);S/=F*(1<<n);}
```
[Try it online!](https://tio.run/##vY89a8MwFEXn9le4gYI@rlvLsdMWWUMWz4aMRkMIcWPZtUWSzei3u7KblNKh2QriwUPnXs7bhe@73ThWbb89BxtY5NKglieyRocCLXo0dP5la1l3Z1bIoa5IH7a06o@kUZFssvYnzwWlRdlr1XAu9@1pP0ygVUa10oShnHO1X@ppsaEqSqOVn7WWvvnB0jmQK6sEjO83WSdzpgy1TK3LooyZ0SxmHfdBzrnRWm64stI593XGYZL5lh68deBR5SOe9IVXWYOIys2zyhkRWdZR6caPbd0ROtzf2aOPVWTxWCFY4EDIXFdqOgjEWCJBihVe8Io3iAhCQMQQS4gEIoVYOcSUyr960glPEAoHcQMNBa7vNvzb7593h@SGYIQgeUovI7oc5MZP "C (gcc) – Try It Online")
# Explanation (271 bytes version)
```
float S,p,F;                    // global float variables: total sum, temporary, factorial
j,i;                            // global integer variables: indices
s(A,n,P,l,o,k)float*A;int*P;{   // recursively look at every permutation in S_n
 if(k=j=0,o-l)                  // initialize k and j, check if o != l (possible  permutation not yet fully generated)
  for(;k<l;s(A,n,P,l,o+1))      // loop through possible values for current possible  permuation position
   P[o]=k++;                    // set possible  permutation, recursively call (golfed into the for loop)
 else{for(p=-l;j<l;j++)         // there exists a possible permutation fully generated
  for(i=0;i<l;)                 // test if the possible permutation is a bijection
   p+=P[j]==P[i++];             // check for unique elements
  if(!p){                       // indeed, it is a permutation
   for(F=p=1,j=0;j<n;F*=j)      // Hafnian product loop and calculate the factorial (over and over to save bytes)
    p*=A[P[2*j]*2*n+P[j-~j++]]; // Hafnian product
   S+=p;}}}                     // add to sum
float h(A,n)float*A;{           // Hafnian function
 int P[j=2*n];S=0;              // allocate permutation memory, initialize sum
 s(A,n,P,j,0);                  // calculate Hafnian sum
 S/=F*(1<<n);}                  // calculate Hafnian
```
[Try it online!](https://tio.run/##vVRLb9pAED6XXzGpVMmPocEE0ofjQy6oRySOCFWLvYY1613XXqNS5P51OoshIUCanoqQZVsz32tmHXcXcbzbpVIzAxMscBTCld/tLSyknjMJbeWalYLNJa@@gtGGXld1jmB4XuiSlRuElMVGU43sZCiuQl5CC2X4gpen4EIlIuZVp3IeUeEYJWpcuXsN3mNIDd443LYgJY/rshJrLjcgtV4ByeRrXm6g4GVeG2aEVgQIk@@qAyJ1VlEW9VB3pXtVlFDCkH7xiwNhqQQyhHjJ4xX1goabCCQ4ha4qQVLhBYnSBjbcQFpLErPgipfM8MTtAKS6dMLVgwxPHPmB6z7RkvQCzLLU9WIJT/BrJmte2W4glyVXBs6oW2Z6KewNMcF4qmfRyvdfG2jFL0Ba/fgiy5hJMrrQMuWJHZEmdXyvxEolT1xWfGt9FVFXhhl5y3zfPWWihpID/ykqUwF7Jj3N7CyrQ1Qi6oWCEN1rBgyvjB2G1XMVU1i2uch4fMyk8KPxNJtFdBW@PwvPIdsBW3O1Ej9qEi15TmlX1Ewrc1O429fXmJaV8wRBmJb5RIrltn5GUREFSHtHOalw5EXZ8@S/sVQJRjMsdVLHpt0Eu3g0gbiWlEqb/PFkgaNpvfcV@xuaTMXWHOYbCsbuGtn1osfpeNr3spnX95RP3ru/aTozcn7JaFsmflSETdO85pElyZ6ozjvtl2Bp9/jpRG5fVh8J0lodRkALRIuZRaRmFk4oh0sGKXVszZ5OMue5tp@Vk0NpJcDxFGXYc8Orgp/DO4rZN05uo5HnBA8Pyg2bf2vc5Uwox9123hUl2Uid9x9ShPe4dJy9/enM3QbYxzsc4BDv8RN@xi8Y9DAIMOhjcIfBAIMhBvcN9l03/BvO0JYPsBs0GLxR2g3w@H@7@Fzff35ucPCGwB7C4OPwcOkdDDW7Pw "C (gcc) – Try It Online")
At the program's core is the following permutation generator which loops through `S_n`. All Hafnian computation is simply build upon it -- and further golfed.
```
j,i,p;Sn(A,l,o,k)int*A;{          // compute every element in S_n
 if(o-l)                          // o!=l, the permutation has not fully been generated
  for(k=0;k<l;k++)                // loop through the integers [0, n)
   A[o]=k,Sn(A,l,o+1);            // modify permutation, call recursively
 else{                            // possible permutation has been generated
  for(p=-l,j=0;j<l;j++)           // look at the entire possible permutation
   for(i=0;i<l;i++)p+=A[j]==A[i]; // check that all elements appear uniquely
  if(!p)                          // no duplicat elements, it is indeed a permutation
   for(printf("["),j=0;j<l        // print
   ||printf("]\n")*0;)            //  the
    printf("%d, ",A[j++]);}}      //   permutation
main(){int l=4,A[l];Sn(A,l,0);}   // all permutations in S_4
```
[Try it online!](https://tio.run/##fZLBTuMwEIbvfYqh0koxcUWRuGVz6DNw7EbIJJN2Esf22g5SVXj2Mm4JFCj4kEv@@TTfzNSLTV0fDp0k6Yp7k62kllb2gky8XhV7eH83N1DbwY0RAZ/Q7wA1DmgikIH7BzMDajO70AJ@fEywV6WWELcIDv0wRhXJGtiqAMZGaEetd/CIaGCDBr2K2MwAWuuzvlwW/V9d9HkuLnC1tY6x3o6b7RHP7eMGfYD1UoIRTIHV2lZlLyfH/FYUXyiDbajdnbcmoVZag8d69IGeUO9m7B1wD/Crp7Mh0KP@rnlRzpULLTs27Niw@2x4kutBxaMXD5w8XuQnx0QjBhGDiEEuL1frrir5S1Vx3OEW655RzEtmb0sMoJxD5WE09H88aqZ9Xrnf12ksNKPTVDNtIkkgvonAG2gQG1CXWnSe99Nm8/VcTN7nw0t/U/T5eQpW/8xcXC8L8aWBNJOUhCn4p5Ewl@yc55UoXl4@kp/6GBSZTOy5BnR5x3ldTce/5LJTSZrPWVE4Xfrd4fAK "C (gcc) – Try It Online")
[Answer]
# [Jelly](https://github.com/DennisMitchell/jelly), 29 bytes
```
LHµ2*×!
LŒ!s€2;@€€Wị@/€€P€S÷Ç
```
[Try it online!](https://tio.run/##y0rNyan8/9/H49BWI63D0xW5fI5OUix@1LTGyNoBSAJR@MPd3Q76EHYAEAcf3n64/f///9HRhjpGOsY6JjqmOmY65joWsTrkiShQyaBYAA "Jelly – Try It Online")
I think the `;@€€Wị@/€€P€` part can likely be golfed down. Need to come back later to check and add an explanation.
[Answer]
# [Haskell](https://www.haskell.org/), 136 bytes
*-14 bytes thanks to ovs.*
```
import Data.List
f m|n<-length m`div`2=sum[product[m!!(s!!(2*j-2)-1)!!(s!!(2*j-1)-1)/realToFrac(2*j)|j<-[1..n]]|s<-permutations[1..n*2]]
```
[Try it online!](https://tio.run/##3ZBBboMwEEX3nGKyCxG42JDQSGRXddXuukOWYiWkIcWAwHSVu9MZh0QQblAsLM@fr/fHPqv2JyuKvs91XTUG3pRR7CNvjXMCfS0Tv8jKb3MGvT/mv3uxazud1k117A4m1YvFssVfrC6@cH3ujmpO9UuTqeKrem/UgUT3ekn8lDNWSnltE7/OGt0ZZfKqbK28ElL2WuUl7ECr@hOWdZOXBhicXAfwS2kFHkRsLa0A4EGKlQeBBOk4gwJpiAqLR6aATOHUZEkxbZs5E3XsC8aJFE1am7vuAWdiNgmLrD6bibOQYvi9ywLatqRN@ZwkQSHxw2uZ4wvxhyO6MewJkeHTvYOBYR2CTiENwtnrxLYdUgaGfcD1fLbbHDbFMqyDSvF0WcxCGdseIAafDP3jxH9tAFz9Hw "Haskell – Try It Online")
Ugh...
[Answer]
# [MATL](https://github.com/lmendo/MATL), ~~29~~ ~~24~~ 22 bytes
```
Zy:Y@!"G@2eZ{)tn:E/pvs
```
[Try it online!](https://tio.run/##y00syfn/P6rSKtJBUcndwSg1qlqzJM/KVb@grPj//2hjBQM9c2sgYapgHAsA) Or verify all test cases: [1](https://tio.run/##y00syfn/P6rSKtJBUcndwSg1qlqzJM/KVb@grPj//2hjBQM9c2sgYapgHAsA), [2](https://tio.run/##y00syfn/P6rSKtJBUcndwSg1qlqzJM/KVb@grPj//2gDHQUTPXMQYQYiTK0hXKCwkZ4hkNYzsYbIQbg6CoZ6RiARU7AcmAtkxQIA), [3](https://tio.run/##y00syfn/P6rSKtJBUcndwSg1qlqzJM/KVb@grPj//2hDPWMdBRM9Qx0FQz0jHQUDPQMQYQkSM7GGSJiAJIz0zOFKTEAsc2sIFyxhAtEBZgENMAbJgo0CKwFLGIFYxiDbDPUsrCF2gI2C6jAGcU3h9kLsABsF1qEDETPWM4oFAA).
### How it works
```
Zy       % Size of (implicit) input: pushes [2*n 2*n], where the
         % input is a 2*n × 2*n matrix. 
:        % Range: gives row vector [1 2 ... 2*n]
Y@       % All permutation of that vector as rows of a matrix
!"       % For each permutation 
  G      %   Push input matrix
  @      %   Push current permutation
  2e     %   Reshape as a 2-row array
  Z{     %   Split rows into a cell array of size 2
  )      %   Reference indexing. With a cell array as index this
         %   applies element-wise indexing (similar to sub2ind).
         %   Gives a row vector with the n matrix entries selected
         %   by the current permutation
  t      %   Duplicate
  n:     %   Number of elements, range: this gives [1 2 ... n]
  E      %   Double, element-wise: gives [2 4 ... 2*n]
  /      %   Divide, element-wise
  p      %   Product
  vs     %   Vertically concatenate and sum
         % End (implicit). Display (implicit)
```
[Answer]
# [JavaScript (Node.js)](https://nodejs.org), ~~121~~ ~~119~~ ~~117~~ ~~115~~ 101 bytes
```
A=>(F=(a,t,r)=>{a--?A[t].map((x,i)=>F[i]=F[i]||F(F[i]=a,i,a%2?r:r*x/(a+2))):s+=r})(A.length,s=0,1)||s
```
[Try it online!](https://tio.run/##bY/LboMwEEX3@Qo2lezGDNi8QiQnYsNPIBZWShNX1CCDqkil305tQx5SuxnP3DtzD3yILzGctOxHX3VvzdxrqUZ@6tTQtQ203XlzEe9KCsXngh9QyZEgI9GYH76F7x@LaqzhU/QIXYk0YlnJmtsyTSVygyCSiBd21Hv9eg2Q2DKM8X7Ycv2DUQFto87jhQw8JBRP07B8AFqhqKpC4sWQ1MSrzEO8sK4x3gSBGTb/rma2pE9HRjAGA2peiBctvQnEo8Ae4RA74YGhKeTRHxCFyALobRlCW3KrrQBqe2Yp2X3JZWfWp3crXq5cZ0KixXeBbslZzHaRZVLYLX6@Bq5XkR2TJ/7CcoHuyll2ZOuPsZRBnOzmXw "JavaScript (Node.js) – Try It Online")
## Explanation & Ungolfed
```
function hafnian(A) {            // Main function taking an array
 return (
  F = function(                  //  Helper function, and also a temporary storage
   a,                            //   how many numbers left
   t,                            //   previous number
   r                             //   the product
  ) {                            //   The braces are to avoid returning value
   a-- ?                         //   If a > 0:
    A[t].map((x, i) =>           //    Loop from i = 2n - 1 downto 0:
     F[i] = F[i] || F(           //     Skip if F[i] is truthy, i.e. i used in the permutation    
      F[i] = a,                  //     Otherwise: Set F[i] to truthy* & start a new recursion
      i,                         //      with i as the previous number
      a % 2                      //      If a - 1 is odd:
       ? r                       //       r does not change
       : r * x / (a + 2)         //      Otherwise: multiply A[t][i] to r, divide by a + 1
     )                           //      Set back F[i] to undefined after the recursion
    )
   : s += r;                     //   Otherwise: add r to s. The factor 1/(n!*2^n) = 1/(2n)!!
  }                              //   has been included in the recursion
 )(
  A.length,                      //  2n = size of the array
  s = 0,                         //  Initial sum s = 0
  1                              //  Initial product r = 1
 ) || s;                         //  Return s
}
/* Note that when a = 1 F[i] will be falsy, but the next recursion will just 
   skip this part anyway. */ 
```
[Answer]
# [Coconut](http://coconut-lang.org/), ~~165~~ ~~145~~ ~~128~~ 127 bytes
-1 byte thanks to Mr. Xcoder
```
m->sum(reduce((o,j)->o*m[p[2*j]][p[j-~j]]/-~j/2,range(len(m)//2),1)for p in permutations(range(len(m))))
from itertools import*
```
[Try it online!](https://tio.run/##rU7NDsIgGLvvKTjC8m3sR50X9yLIYZnMsAw@wuBmfHXEm/Fo7KFtmqbpjDPaGNJyuSZTjXs01KtbnBWlCCurRiyNcKIrVymzrtUzG56Zd@Ane1d0U5YaxnnHoGULeuKItsQpb2KYgka7089iRrF4NEQH5QPithNtHPpQJiF6aOpBgmjqI/RSksdIljc5r20oCiFa6KCHAxzhBAOcc/WnhPxp6OthegE "Coconut – Try It Online")
[Answer]
## Perl 6, 86 bytes
```
{my \n=$^m/2;^$m .permutations.map({[*] .map(->\a,\b{$m[a][b]})}).sum/(2**n*[*] 1..n)}
```
[Answer]
# [JavaScript (Node.js)](https://nodejs.org), ~~87~~ 86 bytes
```
f=A=>A[0]?A[0].reduce((y,x,i)=>i&&y+x*f((g=B=>B.filter((_,j)=>j&&i-j))(A).map(g)),0):1
```
[Try it online!](https://tio.run/##bY/dboMwDIXveQquULKlLgl/ZVKY6GsgNCEKLIgBomxqn57Fgf5I241jn2OfT2mLn@JcTmqcd/1wqpZxUv0sy6E/D10F3dBYn0Xdq6KXSy1TmaSZm79jgak6fZcVIVd2YYrKRDnO9fXyUhPSyKNMjlCrbq4mQj5Yq@3WcdSupZSkFL6KkTSUMpe@8ZVINgrJMpfZPgQ5szP9MNvNc0qt/V4P1r@rEZbw6UgL2hDA9Qv@qoU3gdkcxCMcfCM8MDyE2PsD4uAhgN@WwcUSo7YBOPYCKdF9yWRH6PO75a9XptMh3uqbQLNkLIGdh0wOh9WPt8DtysMxeOKvLBNoroyFo9g@JkIBfnBYfgE "JavaScript (Node.js) – Try It Online")
This is a port of Kirill's algorithm, and this turns out to be shorter than the naive approach! However I'm retaining my original answer and posting this as a separate one.
### Explanation
```
f = A =>                               // A function receiving an array as input
 A[0]                                  //  If A is not empty:
  ? A[0].reduce((y, x, i) =>           //   Calculate the sum
   i                                   //   for i = 1 to 2n - 1:
   && y + x * f(                       //    A[0][i] * the hafnian of
    (g = B =>
     B.filter((_, j) => j && i - j)    //    array A removing rows 0 and i
    )(A).map(g)                        //    as well as columns 0 and i
   ),
   0
  )
  : 1                                  //  If A is empty: return 1
```
[Answer]
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 84 bytes
```
aTr[1##&@@(a[[##]]&@@@#~Partition~2)&/@Permutations@Range[n=Tr[1^a]]]/(n/=2)!/2^n
```
[Try it online!](https://tio.run/##TY1NasMwEIX3OcUUgUlAlSzZiZOFi45gSnfCAVGc1gsr4KoroxyiN@gJewRXIzk/m/l5b943g3Gf3WBc/27mE9Sz@fv5fRu1ICRTam20JqRtw6jIpTGj611/the5ybhqunH4dgaFL/Vq7EenbY3Ro2nblq8tr@XmicujnZuxt04TeH6Bkw48yICr1TRNOYWSbT2FKTQKufd0BYtcYdk9HAQhGJKJ0FmZtN1VoCCYvINYGYU7UrACUeIqsxzLAbUFJXCWyKtuR5FSoS9uVplScQqQIvkRGI@iJXEq8Kdg@@QfFuCSKnDdPvxPvyIwpqKFq/Tez/8 "Wolfram Language (Mathematica) – Try It Online")
---
# [Wolfram Language (Mathematica)](https://www.wolfram.com/wolframscript/), 68 bytes
```
If[#=={},1,Sum[d=Rest@*Delete[i];#[[1,i]]#0@d[d/@#],{i,2,Tr[1^#]}]]&
```
[Try it online!](https://tio.run/##TY09T8MwEIb3/oqTLHVAh5tz0oYKBWVgYUPAZrlSRFxhiTAUM1n@7cFnpx/Lfbzvvc9Ng/@y0@Dd5zAfoZtfjlp0XYhI@P436bF7s7@@v3u239Zb7cyj0JrQGSOqftTjphcGg0OFHydNB2GiMev59eR@vBZw/wQJZwysYdOvQggVQiO3ESGkhlDFiCtY5JbL7uYgCclQklKXTdF2ZwGBpLqCZJOFK5JkzSg6y7LismdtQRHPinnt5ShTWvbpYjUllacEqYufgfkoW4qnmn@SfCj@fgEuqZrX7c3/8isDcypbvKoY4/wP "Wolfram Language (Mathematica) – Try It Online")
[Answer]
# [Pari/GP](http://pari.math.u-bordeaux.fr/), 51 bytes
```
h(m)=if(#m,sum(n=2,#m,m[1,n]*h(m[^n,^n][^1,^1])),1)
```
[Try it online!](https://tio.run/##PY7dDsIgDIVfhcybYbpmhek0y3wRwhJv1CVCiD8XPv2kRXdTes6hX5vOj7m5pmW51UGP86XeBHi@Qx1HA7kNjiD6bQ7dFGGK3k0EE3mtgfRyTun@qYNqTio95vjKbcWiUkzToJxrQXW4G7iAan22xOm57Ncsy2wbpPxiN5SsSFCE5j@PncgfidAygf4etlyO7AmDuDMM6tcvAuiHIk2547helAGWU0HJFwkMd5a3ER6GskNQvwnLcrfuLTsEJRNQPIvGe718AQ "Pari/GP – Try It Online")
[Answer]
# [Python 3.8 (pre-release)](https://docs.python.org/3.8/), ~~185~~ 181 bytes
```
def f(m):
 a=0;n=len(m)//2
 for x in itertools.permutations(range(n*2)):p=1;[p:=p*m[x[2*j]][x[2*j+1]]for j in range(n)];a+=p
 return a/(math.factorial(n)*2**n)
import math,itertools
```
[Try it online!](https://tio.run/##ZY7NbsMgEITvfoo9GgcRg8m//CSIA2rtxlEMiBApfXp3107Tqj3AwMzsB/Ezn4Nv9jFN03vXQ1@O7FiAa@uTb6@dx@t6rQroQ4IHDB6G3KUcwvUmYpfGe3Z5CP5WJuc/utJXirFjbOXJxGMbq9E8jKou1i66ktYS6EKg5wSzJ7dqYwGpy/fkwa3L0eWz6N1bDmlwV6xUqqo8K4YxhpSBYv76xhTT4HPZl8bUHLTYWA4GhUNtLWPFn3hH2/ZXEQ0MlJCoQi/e9tvgIIX6AQo9G//RUjRc4wSFtahxHfCONEOuRlfhQ0uqUXeYyKer5y5pzZs5IQKl5CrUBulS7OfkMBOWboPnzeudmUsE6pKLZ0X/nL4A "Python 3.8 (pre-release) – Try It Online")
Golfed reference implementation provided in the question.
]