Return to Answer

added 1863 characters in body

edited Mar 12, 2018 at 23:57

17.3k
2
31
95

Here is a faster version with some of the easy optimizations.

def hafnian(m): n = len(m)/2 z = [[0]*(n+1) for _ in range(n*(2*n-1))] for j in range(1, 2*n): for k in range(j): z[j*(j-1)/2+k][0] = m[j][k] return solve(z, 2*n, 1, [1] + [0]*n, n) def solve(b, s, w, g, n): if s == 0: return w*g[n] c = [b[(j+1)*(j+2)/2+k+2][:] for j in range(1, s-2) for k in range(j)] h = solve(c, s-2, -w, g, n) e = g[:] for u in range(n): for v in range(n-u): e[u+v+1] += g[u]*b[0][v] for j in range(1, s-2): for k in range(j): for u in range(n): for v in range(n-u): c[j*(j-1)/2+k][u+v+1] += b[(j+1)*(j+2)/2][u]*b[(k+1)*(k+2)/2+1][v] + b[(k+1)*(k+2)/2][u]*b[(j+1)*(j+2)/2+1][v] return h + solve(c, s-2, w, e, n)

Try it online!

For added fun, here is a reference implementation in J.

Here is a faster version with some of the easy optimizations.

def hafnian(m): n = len(m)/2 z = [[0]*(n+1) for _ in range(n*(2*n-1))] for j in range(1, 2*n): for k in range(j): z[j*(j-1)/2+k][0] = m[j][k] return solve(z, 2*n, 1, [1] + [0]*n, n) def solve(b, s, w, g, n): if s == 0: return w*g[n] c = [b[(j+1)*(j+2)/2+k+2][:] for j in range(1, s-2) for k in range(j)] h = solve(c, s-2, -w, g, n) e = g[:] for u in range(n): for v in range(n-u): e[u+v+1] += g[u]*b[0][v] for j in range(1, s-2): for k in range(j): for u in range(n): for v in range(n-u): c[j*(j-1)/2+k][u+v+1] += b[(j+1)*(j+2)/2][u]*b[(k+1)*(k+2)/2+1][v] + b[(k+1)*(k+2)/2][u]*b[(j+1)*(j+2)/2+1][v] return h + solve(c, s-2, w, e, n)

Try it online!

For added fun, here is a reference implementation in J.

Source Link

answered Mar 7, 2018 at 6:10

miles

17.3k
2
31
95

Python

This is pretty much a straight-forward, reference implementation of Algorithm 2 from the mentioned paper. The only changes were to only keep the current value of B, dropping the values of β by only updating g when i ∈ X, and truncated polynomial multiplication by only calculating the values up to degree n.

from itertools import chain,combinations def powerset(s): return chain.from_iterable(combinations(s, k) for k in range(len(s)+1)) def padd(a, b): return [a[i]+b[i] for i in range(len(a))] def pmul(a, b): n = len(a) c = [0]*n for i in range(n): for j in range(n): if i+j < n: c[i+j] += a[i]*b[j] return c def hafnian(m): n = len(m) / 2 z = [[[c]+[0]*n for c in r] for r in m] h = 0 for x in powerset(range(1, n+1)): b = z g = [1] + [0]*n for i in range(1, n+1): if i in x: g = pmul(g, [1] + b[0][1][:n]) b = [[padd(b[j+2][k+2], [0] + padd(pmul(b[0][j+2], b[1][k+2]), pmul(b[0][k+2], b[1][j+2]))[:n]) if j != k else 0 for k in range(2*n-2*i)] for j in range(2*n-2*i)] else: b = [r[2:] for r in b[2:]] h += (-1)**(n - len(x)) * g[n] return h

Try it online!

For added fun, here is a reference implementation in J.