# Codegolf the Hafnian

The challenge is to write codegolf for the Hafnian of a matrix. The Hafnian of an 2n-by-2n symmetric matrix A is defined as: Here S2n represents the set of all permutations of the integers from 1 to 2n, that is [1, 2n].

For those interested in applications of the Hafnian, the mathoverflow 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.

• I think this would be easier to digest with an informal explanation of the hafnian. Something like, take all subsets of n matrix entries where their n row indices and n column indices form a partition of 1..2n, take the product of each, add them, and scale the sum. – xnor Mar 1 '18 at 17:11

# R, 150142127 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 <- 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
}


• Apply is cheaper by 2 bytes, which allows for additional 4 bytes of saving by cramming the other lines together. tio.run/##PY6xDoIwEIZ3nsLxzpxiS4ymkYEXYHIjDFDEEKBtSokS47PX4sDw5/… It is also quite interesting how base R lacks a permutation function for a statistical programming language – Vlo Mar 1 '18 at 21:36
• @Vlo very nice! we can move N and k into function arguments to get it to one statement, removing the {} and saving another two bytes. – Giuseppe Mar 1 '18 at 21:41
• @Giuseppe Darn keeps on forgetting that you can define them in the function args. Spent a few minutes trying to ember those variables... – Vlo Mar 1 '18 at 21:42

# Pyth, 24 bytes

sm*Fmc@@Qhkek2d{mScd2.pU


Try it here!

### Old version, 35 bytes

*c1**FK/lQ2^2Ksm*Fm@@Q@dtyk@dykK.pU


Try it here!

• Currently the lead but you have to fear the Jelly answers to come.... :) – user9206 Mar 1 '18 at 17:25
• Eh Jelly will definitely beat my by about 10 bytes. Pyth ain’t the best tool for the job – Mr. Xcoder Mar 1 '18 at 17:26
• 05AB1E is looking like it may even tie Pyth (believe it or not, finally a matrix challenge where a[b] is enough to compete). – Magic Octopus Urn Mar 1 '18 at 17:31
• @MagicOctopusUrn I already have an 05AB1E solution that beats Pyth :-) Not going to post it (for now, at least) – Mr. Xcoder Mar 1 '18 at 17:32
• Is it something along the lines of xÍysè<¹sès·<ysè<è lmao? P.S. Mine is 40 bytes and isn't working so well hah, so feel free to post it, unsure I'm able to finish before I have to go home. – Magic Octopus Urn Mar 1 '18 at 17:38

# Stax, 232219 17 bytes

ü;Y╙◘▌Φq↓ê²╧▐å↑┌C


Run and debug it online

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

• Very impressive! – user9206 Mar 1 '18 at 20:11

# 05AB1E, 21 bytes

ā<œε2ô{}ÙεεIsèsè;]PO


Try it online!

### Old version, 32 bytes

āœvIg;©Lε·UIyX<èèyXèè}P}Oθ®!/®o/


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

• No kidding èsè, hah... haha... I'm punny. – Magic Octopus Urn Mar 2 '18 at 17:13
• @MagicOctopusUrn Fixed... I forgot 05AB1E is 0-indexed >_< – Mr. Xcoder Mar 2 '18 at 17:21

# Jelly, 19 bytes

LŒ!s€2Ṣ€QḅL_LịFHP€S


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!

### 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.


# C (gcc), 288285282293292272 271 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; 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);.
float S,p,F;j,i;s(A,n,P,l,o,k)float*A;int*P;{if(k=j=0,o-l)for(;k<l;s(A,n,P,l,o+1))P[o]=k++;else{for(p=-l;j<l;j++)for(i=0;i<l;)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!

# Explanation

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
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!

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!

• Great to have a C answer but, as you suggest, it is non compliant currently. – user9206 Mar 2 '18 at 7:27
• @Lembik Fixed. Now supports float matrices. – Jonathan Frech Mar 2 '18 at 15:51
• 2*j+++1 is equivalent to j+j+++1, which is the same as j-(-j++-1), so we can use the bitwise complement efficiently to save a byte: j-~j++ (Try it online) – Mr. Xcoder Mar 2 '18 at 18:51

# R, 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!

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, 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
}

• Very nice answer. -1 for calling If with parenthesis, -4 for using F as initialized variable, -1 for assigning n within the if. tio.run/##XU/LCsIwELz7FcFTVtOQl1pf1/… – Vlo Mar 27 '18 at 21:35
• neat! I'd say post it in the speed challenge, but there are probably some more optimizations (like threading) that can be made, and while R isn't exactly known for its speed, it'd be good to have it there for reference. – Giuseppe Mar 28 '18 at 12:44
• Do it for benchmark purposes! – Vlo Mar 28 '18 at 15:32
• I actually tried testing this for speed, but got quickly discouraged by the results. The slowest Python submission in the speed challenge using the same exact algorithm crunches 24x24 matrix in a few seconds on TIO, but R times out. On my local machine it also didn't respond in a reasonable time, even when aided with memoization from package 'memo'... – Kirill L. Mar 29 '18 at 8:48

# Jelly, 29 bytes

LHµ2*×!
LŒ!s€2;@€€Wị@/€€P€S÷Ç


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.

• Identical to my solution (except the J) before golfed down. Jelly minds think alike? source – user202729 Mar 2 '18 at 7:17
• I was able to reduce it a bit more by refactoring the part that you mentioned as well as the division by 2 and the factorial. LḶŒ!s€2ḅL‘ịFZµPS÷JḤ$P$ TIO – miles Mar 2 '18 at 8:59
• @user202729 haha nice – dylnan Mar 2 '18 at 15:39
• @miles Wow that's a lot of savings. I'll edit it into my answer but it's pretty different so feel free to submit your own answer if you want – dylnan Mar 2 '18 at 15:40

-14 bytes thanks to ovs.

import Data.List
f m|n<-length mdiv2=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!

Ugh...

# MATL, 2924 22 bytes

Zy:Y@!"G@2eZ{)tn:E/pvs


Try it online! Or verify all test cases: 1, 2, 3.

### 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)


# Wolfram Language (Mathematica), 85 bytes

#[[s[[j]],s[[n+j]]]]~Product~{j,n}~Sum~{s,n=Tr[1^#]/2;Permutations@Range[2n]}/n!/2^n&


Try it online!

• Very slightly disappointed Mathematica doesn't have a function for this :) – user9206 Mar 2 '18 at 8:27

# Coconut, 165145128 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!

## Perl 6, 86 bytes

{my \n=$^m/2;^$m .permutations.map({[*] .map(->\a,\b{\$m[a][b]})}).sum/(2**n*[*] 1..n)}
`