Wyzwanie polega na napisaniu najszybszego możliwego kodu do obliczenia Hafniana matrycy .

Hafnian symetrycznego 2n-by- 2nmatrycę Aokreśla się jako:

Tutaj S _2n reprezentuje zestaw wszystkich permutacji liczb całkowitych od 1do 2n, to znaczy [1, 2n].

Link do wikipedii daje również inną formułę, która może być interesująca (a jeszcze szybsze metody istnieją, jeśli spojrzysz dalej w Internecie). Ta sama strona wiki mówi o macierzach przylegania, ale twój kod powinien również działać dla innych macierzy. Możesz założyć, że wszystkie wartości będą liczbami całkowitymi, ale nie że wszystkie są dodatnie.

Istnieje również szybszy algorytm, ale wydaje się trudny do zrozumienia. a Christian Sievers jako pierwszy wdrożył go (w Haskell).

W tym pytaniu wszystkie macierze są kwadratowe i symetryczne o równych wymiarach.

Implementacja referencyjna (należy pamiętać, że jest to najwolniejsza możliwa metoda).

Oto przykładowy kod python od pana Xcodera.

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([[-1, 1, 1, -1, 0, 0, 1, -1], [1, 0, 1, 0, -1, 0, -1, -1], [1, 1, -1, 1, -1, -1, 0, -1], [-1, 0, 1, -1, -1, 1, -1, 0], [0, -1, -1, -1, -1, 0, 0, -1], [0, 0, -1, 1, 0, 0, 1, 1], [1, -1, 0, -1, 0, 1, 1, 0], [-1, -1, -1, 0, -1, 1, 0, 1]]))
4

M = [[1, 1, 0, 0, 0, 0, 0, 1, 0, 0], [1, 1, -1, 0, -1, 1, 1, 1, 0, -1], [0, -1, -1, -1, 0, -1, -1, 0, -1, 1], [0, 0, -1, 1, -1, 1, -1, 0, 1, -1], [0, -1, 0, -1, -1, -1, -1, 1, -1, 1], [0, 1, -1, 1, -1, 1, -1, -1, 1, -1], [0, 1, -1, -1, -1, -1, 1, 0, 0, 0], [1, 1, 0, 0, 1, -1, 0, 1, 1, -1], [0, 0, -1, 1, -1, 1, 0, 1, 1, 1], [0, -1, 1, -1, 1, -1, 0, -1, 1, 1]]

print(hafnian(M))
-13

M = [[-1, 0, -1, -1, 0, -1, 0, 1, -1, 0, 0, 0], [0, 0, 0, 0, 0, -1, 0, 1, -1, -1, -1, -1], [-1, 0, 0, 1, 0, 0, 0, 1, -1, 1, -1, 0], [-1, 0, 1, -1, 1, -1, -1, -1, 0, -1, -1, -1], [0, 0, 0, 1, 0, 0, 0, 0, 0, 1, -1, 0], [-1, -1, 0, -1, 0, 0, 1, 1, 1, 1, 1, 0], [0, 0, 0, -1, 0, 1, 1, -1, -1, 0, 1, 0], [1, 1, 1, -1, 0, 1, -1, 1, -1, -1, -1, -1], [-1, -1, -1, 0, 0, 1, -1, -1, -1, 1, -1, 0], [0, -1, 1, -1, 1, 1, 0, -1, 1, -1, 1, 1], [0, -1, -1, -1, -1, 1, 1, -1, -1, 1, 0, -1], [0, -1, 0, -1, 0, 0, 0, -1, 0, 1, -1, 1]]

print(hafnian(M))
13

M = [[-1, 1, 0, 1, 0, -1, 0, 0, -1, 1, -1, 1, 0, -1], [1, -1, 1, -1, 1, 1, -1, 0, -1, 1, 1, 0, 0, -1], [0, 1, 1, 1, -1, 1, -1, -1, 0, 0, -1, 0, -1, -1], [1, -1, 1, -1, 1, 0, 1, 1, -1, -1, 0, 0, 1, 1], [0, 1, -1, 1, 0, 1, 0, 1, -1, -1, 1, 1, 0, -1], [-1, 1, 1, 0, 1, 1, -1, 0, 1, -1, -1, -1, 1, -1], [0, -1, -1, 1, 0, -1, -1, -1, 0, 1, -1, 0, 1, -1], [0, 0, -1, 1, 1, 0, -1, 0, 0, -1, 0, 0, 0, 1], [-1, -1, 0, -1, -1, 1, 0, 0, 1, 1, 0, 1, -1, 0], [1, 1, 0, -1, -1, -1, 1, -1, 1, 1, 1, 0, 1, 0], [-1, 1, -1, 0, 1, -1, -1, 0, 0, 1, -1, 0, -1, 0], [1, 0, 0, 0, 1, -1, 0, 0, 1, 0, 0, 1, 1, 1], [0, 0, -1, 1, 0, 1, 1, 0, -1, 1, -1, 1, 1, -1], [-1, -1, -1, 1, -1, -1, -1, 1, 0, 0, 0, 1, -1, -1]]

print(hafnian(M))
83

Zadanie

Powinieneś napisać kod, który podany 2nprzez 2nmacierz, wyświetla swój Hafnian.

Ponieważ będę musiał przetestować kod, pomocne byłoby podanie prostego sposobu na przekazanie macierzy jako danych wejściowych do kodu, na przykład poprzez czytanie ze standardowego wejścia. Przetestuję kod w losowo wybranych matrycach z elementami wybrane z {-1, 0, 1}. Celem takich testów jest zmniejszenie szansy, że Hafnian będzie miał bardzo dużą wartość.

Idealnie twój kod będzie mógł czytać w matrycach dokładnie tak, jak ja je w przykładach w tym pytaniu prosto ze standardowego w. To byłoby wejście, które wyglądałoby [[1,-1],[-1,-1]]na przykład. Jeśli chcesz użyć innego formatu wejściowego, zapytaj, a ja dołożę wszelkich starań, aby to uwzględnić.

Wyniki i krawaty

Przetestuję twój kod na losowych matrycach o rosnącym rozmiarze i zatrzymam się, gdy twój kod po raz pierwszy zajmie więcej niż 1 minutę na moim komputerze. Macierze oceny będą spójne dla wszystkich wniosków, aby zapewnić uczciwość.

Jeśli dwie osoby otrzymają ten sam wynik, zwycięzcą jest ta, która jest najszybsza dla tej wartości n. Jeśli są one w odległości 1 sekundy od siebie, to jest to pierwszy.

Języki i biblioteki

Możesz użyć dowolnego dostępnego języka i bibliotek, które ci się podobają, ale nie ma wcześniejszej funkcji do obliczenia Hafniana. Tam, gdzie to możliwe, dobrze byłoby móc uruchomić kod, więc proszę podać pełne wyjaśnienie, w jaki sposób uruchomić / skompilować kod w systemie Linux, jeśli to w ogóle możliwe.

Moja maszyna Czasy zostaną uruchomione na moim komputerze 64-bitowym. Jest to standardowa instalacja ubuntu z 8 GB pamięci RAM, ośmiordzeniowym procesorem AMD FX-8350 i Radeon HD 4250. Oznacza to również, że muszę mieć możliwość uruchomienia kodu.

Zadzwoń po więcej języków

Byłoby wspaniale uzyskać odpowiedzi w swoim ulubionym, bardzo szybkim języku programowania. Na początek, co powiesz na fortran , nim i rdzę ?

Tabela liderów

• 52 przez mile przy użyciu C ++ . 30 sekund.
• 50 przy użyciu NGN C . 50 sekund.
• 46 Christian Sievers korzystający z Haskell . 40 sekund.
• 40 mil za pomocą Python 2 + pypy . 41 sekund.
• 34 przez ngn przy użyciu Python 3 + pypy . 29 sekund.
• 28 autorstwa Dennisa przy użyciu Pythona 3 . 35 sekund. (Pypy jest wolniejszy)

Haskell

import Control.Parallel.Strategies
import qualified Data.Vector.Unboxed as V
import qualified Data.Vector as VB

type Poly = V.Vector Int

type Matrix = VB.Vector ( VB.Vector Poly )

constpoly :: Int -> Int -> Poly
constpoly n c = V.generate (n+1) (\i -> if i==0 then c else 0)

add :: Poly -> Poly -> Poly
add = V.zipWith (+)

shiftmult :: Poly -> Poly -> Poly
shiftmult a b = V.generate (V.length a) 
                           (\i -> sum [ a!j * b!(i-1-j) | j<-[0..i-1] ])
  where (!) = V.unsafeIndex

x :: Matrix -> Int -> Int -> Int -> Poly -> Int
x  _    0  _ m p = m * V.last p
x mat n c m p =
  let mat' = VB.generate (2*n-2) $ \i ->
             VB.generate i       $ \j ->
                 shiftmult (mat!(2*n-1)!i) (mat!(2*n-2)!j) `add`
                 shiftmult (mat!(2*n-1)!j) (mat!(2*n-2)!i) `add`
                 (mat!i!j)
      p' = p `add` shiftmult (mat!(2*n-1)!(2*n-2)) p
      (!) = VB.unsafeIndex
      r = if c>0 then parTuple2 rseq rseq else r0
      (a,b) = (x mat (n-1) (c-1) m p, x mat' (n-1) (c-1) (-m) p')
              `using` r
  in a+b

haf :: [[Int]] -> Int
haf m = let n=length m `div` 2
        in x (VB.fromList $ map (VB.fromList . map (constpoly n)) m) 
             n  5  ((-1)^n)  (constpoly n 1) 

main = getContents >>= print . haf . read

To implementuje odmianę algorytmu 2 Andreasa Björklunda: liczenie idealnych dopasowań tak szybko jak Ryser .

Kompiluj za ghcpomocą opcji czasu kompilacji -O3 -threadedi używaj opcji czasu wykonywania +RTS -Ndo równoległości. Pobiera dane wejściowe ze standardowego wejścia.

Python 3

from functools import lru_cache

@lru_cache(maxsize = None)
def haf(matrix):
	n = len(matrix)
	if n == 2: return matrix[0][1]
	h = 0
	for j in range(1, n):
		if matrix[0][j] == 0: continue
		copy = list(matrix)
		del copy[:j+1:j]
		copy = list(zip(*copy))
		del copy[:j+1:j]
		h += matrix[0][j] * haf(tuple(copy))
	return h

print(haf(tuple(map(tuple, eval(open(0).read())))))

Wypróbuj online!

C ++ (gcc)

#define T(x) ((x)*((x)-1)/2)
#define S 1
#define J (1<<S)
#define TYPE int

#include <iostream>
#include <vector>
#include <string>
#include <pthread.h>

using namespace std;

struct H {
    int s, w, t;
    TYPE *b, *g;
};

void *solve(void *a);
void hafnian(TYPE *b, int s, TYPE *g, int w, int t);

int n, m, ti = 0;
TYPE r[J] = {0};
pthread_t pool[J];

int main(void) {
    vector<int> a;
    string s;
    getline(cin, s);

    for (int i = 0; i < s.size(); i++)
        if (s[i] == '0' || s[i] == '1')
            a.push_back((s[i-1] == '-' ? -1 : 1)*(s[i] - '0'));

    for (n = 1; 4*n*n < a.size(); n++);
    m = n+1;

    TYPE z[T(2*n)*m] = {0}, g[m] = {0};

    for (int j = 1; j < 2*n; j++)
        for (int k = 0; k < j; k++)
            z[(T(j)+k)*m] = a[j*2*n+k];
    g[0] = 1;

    hafnian(z, 2*n, g, 1, -1);

    TYPE h = 0;
    for (int t = 0; t < ti; t++) {
        pthread_join(pool[t], NULL);
        h += r[t];
    }

    cout << h << endl;

    return 0;
}

void *solve(void *a) {
    H *p = reinterpret_cast<H*>(a);
    hafnian(p->b, p->s, p->g, p->w, p->t);
    delete[] p->b;
    delete[] p->g;
    delete p;
    return NULL;
}

void hafnian(TYPE *b, int s, TYPE *g, int w, int t) {
    if (t == -1 && (n < S || s/2 == n-S)) {
        H *p = new H;
        TYPE *c = new TYPE[T(s)*m], *e = new TYPE[m];
        copy(b, b+T(s)*m, c);
        copy(g, g+m, e);
        p->b = c;
        p->s = s;
        p->g = e;
        p->w = w;
        p->t = ti;
        pthread_create(pool+ti, NULL, solve, p);
        ti++;
    }
    else if (s > 0) {
        TYPE c[T(s-2)*m], e[m];
        copy(b, b+T(s-2)*m, c);
        hafnian(c, s-2, g, -w, t);
        copy(g, g+m, e);

        for (int u = 0; u < n; u++) {
            TYPE *d = e+u+1,
                  p = g[u], *x = b+(T(s)-1)*m;
            for (int v = 0; v < n-u; v++)
                d[v] += p*x[v];
        }

        for (int j = 1; j < s-2; j++)
            for (int k = 0; k < j; k++)
                for (int u = 0; u < n; u++) {
                    TYPE *d = c+(T(j)+k)*m+u+1,
                          p = b[(T(s-2)+j)*m+u], *x = b+(T(s-1)+k)*m,
                          q = b[(T(s-2)+k)*m+u], *y = b+(T(s-1)+j)*m;
                    for (int v = 0; v < n-u; v++)
                        d[v] += p*x[v] + q*y[v];
                }

        hafnian(c, s-2, e, w, t);
    }
    else
        r[t] += w*g[n];
}

Wypróbuj online! (13s dla n = 24)

Na podstawie szybszej implementacji Pythona w moim innym poście. Edytuj drugi wiersz #define S 3na 8-rdzeniowej maszynie i skompiluj z g++ -pthread -march=native -O2 -ftree-vectorize.

Dzieli pracę na pół, więc wartość Spowinna wynosić log2(#threads). Typy można łatwo zmieniać pomiędzy int, long, float, i doublemodyfikując wartość #define TYPE.

Python 3

Oblicza haf (A) jako zapamiętaną sumę (A [i] [j] * haf (A bez wierszy i kolumn i oraz j)).

#!/usr/bin/env python3
import json,sys
a=json.loads(sys.stdin.read())
n=len(a)//2
b={0:1}
def haf(x):
 if x not in b:
  i=0
  while not x&(1<<i):i+=1
  x1=x&~(1<<i)
  b[x]=sum(a[i][j]*haf(x1&~(1<<j))for j in range(2*n)if x1&(1<<j)and a[i][j])
 return b[x]
print(haf((1<<2*n)-1))

do

Kolejny wgląd w pracę Andreasa Björklunda , który jest znacznie łatwiejszy do zrozumienia, jeśli spojrzysz również na kod Haskella Christiana Sieversa . Przez pierwsze kilka poziomów rekurencji rozdziela on wątki typu round-robin na dostępne procesory. Ostatni poziom rekurencji, który stanowi połowę inwokacji, jest optymalizowany ręcznie.

Skompilować z: gcc -O3 -pthread -march=native; dzięki @Dennis za 2x przyspieszenie

n = 24 na 24 sekundy w TIO

#define _GNU_SOURCE
#include<sched.h>
#include<stdio.h>
#include<stdlib.h>
#include<memory.h>
#include<unistd.h>
#include<pthread.h>
#define W while
#define R return
#define S static
#define U (1<<31)
#define T(i)((i)*((i)-1)/2)
typedef int I;typedef long L;typedef char C;typedef void V;
I n,ncpu,icpu;
S V f(I*x,I*y,I*z){I i=n,*z1=z+n;W(i){I s=0,*x2=x,*y2=y+--i;W(y2>=y)s+=*x2++**y2--;*z1--+=s;}}
typedef struct{I m;V*a;V*p;pthread_barrier_t*bar;I r;}A;S V*(h1)(V*);
I h(I m,I a[][n+1],I*p){
 m-=2;I i,j,k=0,u=T(m),v=u+m,b[u][n+1],q[n+1];
 if(!m){I*x=a[v+m],*y=p+n-1,s=0;W(y>=p)s-=*x++**y--;R s;}
 memcpy(b,a,sizeof(b));memcpy(q,p,sizeof(q));f(a[v+m],p,q);
 for(i=1;i<m;i++)for(j=0;j<i;j++){f(a[u+i],a[v+j],b[k]);f(a[u+j],a[v+i],b[k]);k++;}
 if(2*n-m>8)R h(m,a,p)-h(m,b,q);
 pthread_barrier_t bar;pthread_barrier_init(&bar,0,2);pthread_t th;
 cpu_set_t cpus;CPU_ZERO(&cpus);CPU_SET(icpu++%ncpu,&cpus);
 pthread_attr_t attr;pthread_attr_init(&attr);
 pthread_attr_setaffinity_np(&attr,sizeof(cpu_set_t),&cpus);
 A arg={m,a,p,&bar};pthread_create(&th,&attr,h1,&arg);
 I r=h(m,b,q);pthread_barrier_wait(&bar);pthread_join(th,0);pthread_barrier_destroy(&bar);
 R arg.r-r;
}
S V*h1(V*x0){A*x=(A*)x0;x->r=h(x->m,x->a,x->p);pthread_barrier_wait(x->bar);R 0;}
I main(){
 ncpu=sysconf(_SC_NPROCESSORS_ONLN);
 S C s[200000];I i=0,j=0,k,l=0;W((k=read(0,s+l,sizeof(s)-l))>0)l+=k;
 n=1;W(s[i]!=']')n+=s[i++]==',';n/=2;
 I a[T(2*n)][n+1];memset(a,0,sizeof(a));k=0;
 for(i=0;i<2*n;i++)for(j=0;j<2*n;j++){
  W(s[k]!='-'&&(s[k]<'0'||s[k]>'9'))k++;
  I v=0,m=s[k]=='-';k+=m;W(k<l&&('0'<=s[k]&&s[k]<='9'))v=10*v+s[k++]-'0';
  if(i>j)*a[T(i)+j]=v*(1-2*m);
 }
 I p[n+1];memset(p,0,sizeof(p));*p=1;
 printf("%d\n",(1-2*(n&1))*h(2*n,a,p));
 R 0;
}

Algorytm:

Matryca, która jest symetryczna, jest przechowywana w lewym dolnym trójkątnym kształcie. Wskaźniki trójkąta i,jodpowiadają indeksowi liniowemu, dla T(max(i,j))+min(i,j)którego Tjest makrem i*(i-1)/2. Elementy macierzy są wielomianami stopnia n. Wielomian jest reprezentowany jako tablica współczynników uporządkowanych od stałego składnika ( p[0]) do współczynnika x ⁿ ( p[n]). Początkowe wartości macierzy -1,0,1 są najpierw konwertowane na wielomiany const.

Wykonujemy krok rekurencyjny z dwoma argumentami: pół macierzą (tj. Trójkątem) awielomianów i oddzielnym wielomianem p(w artykule określanym jako beta). Zmniejszamy mproblem wielkości (początkowo m=2*n) rekurencyjnie do dwóch problemów wielkości m-2i zwracamy różnicę u ich hafnianów. Jednym z nich jest użycie tego samego abez dwóch ostatnich rzędów i to samo p. Innym jest użycie trójkąta b[i][j] = a[i][j] + shmul(a[m-1][i],a[m-2][j]) + shmul(a[m-1][j],a[m-2][i])(gdzie shmuljest operacja mnożenia z przesunięciem na wielomianach - to jak zwykle iloczyn wielomianowy, dodatkowo pomnożony przez zmienną „x”; moce wyższe niż x ^ n są ignorowane) i oddzielny wielomian q = p + shmul(p,a[m-1][m-2]). Kiedy rekurencja uderza wielkość-0 a, wracamy główną współczynnik p: p[n].

Operacja przełączania i mnożenia jest realizowana w funkcji f(x,y,z). Zmienia zna miejscu. Mówiąc luźno, to prawda z += shmul(x,y). To wydaje się być najbardziej krytyczną częścią wydajności.

Po zakończeniu rekurencji musimy naprawić znak wyniku, mnożąc przez (-1) ⁿ .

Pyton

Jest to w zasadzie prosta referencyjna implementacja algorytmu 2 ze wspomnianego artykułu . Jedynymi zmianami było utrzymanie bieżącej wartości B , opuszczenie wartości p tylko przez aktualizujący g podczas i ∈ X i ściętego wielomian mnożenie tylko przez obliczenie wartości do stopnia 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

Wypróbuj online!

Oto szybsza wersja z kilkoma łatwymi optymalizacjami.

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)

Wypróbuj online!

Dla dodatkowej zabawy, oto implementacja referencyjna w J.

Oktawa

Jest to w zasadzie kopia wpisu Dennisa , ale zoptymalizowana dla Octave. Głównej optymalizacji dokonuje się przy użyciu pełnej macierzy wejściowej (i jej transpozycji) i rekurencji przy użyciu tylko wskaźników macierzowych, zamiast tworzenia zredukowanych macierzy.

Główną zaletą jest ograniczenie kopiowania matryc. Podczas gdy Octave nie ma różnicy między wskaźnikami / referencjami a wartościami, a funkcjonalnie tylko przekazuje wartość, to inna historia za kulisami. Tam stosuje się kopiowanie przy zapisie (leniwa kopia). Oznacza to, że dla kodu a=1;b=a;b=b+1zmienna bjest kopiowana tylko do nowej lokalizacji w ostatniej instrukcji, gdy jest zmieniana. Odmatin i matranspnigdy nie są zmieniane, nigdy nie zostaną skopiowane. Wadą jest to, że funkcja spędza więcej czasu na obliczaniu poprawnych wskaźników. Być może będę musiał wypróbować różne warianty między wskaźnikami liczbowymi i logicznymi, aby to zoptymalizować.

Ważna uwaga: macierz wejściowa powinna być int32! Zapisz funkcję w pliku o nazwiehaf.m

function h=haf(matin,indices,matransp,transp)

    if nargin-4
        indices=int32(1:length(matin));
        matransp=matin';
        transp=false;
    end
    if(transp)
        matrix=matransp;
    else
        matrix=matin;
    end
    ind1=indices(1);
    n=length(indices);
    if n==2
        h=matrix(ind1,indices(2));
        return
    end
    h=0*matrix(1); 
    for j=1:(n-1)
        indj=indices(j+1);
        k=matrix(ind1,indj);
        if logical(k)
            indicestemp=true(n,1);
            indicestemp(1:j:j+1)=false;
            h=h+k.*haf(matin,indices(indicestemp),matransp,~transp);
        end
    end
end

Przykładowy skrypt testowy:

matrix = int32([0 0 1 -1 1 0 -1 -1 -1 0 -1 1 0 1 1 0 0 1 0 0 1 0 1 1;0 0 1 0 0 -1 -1 -1 -1 0 1 1 1 1 0 -1 -1 0 0 1 1 -1 0 0;-1 -1 0 1 0 1 -1 1 -1 1 0 0 1 -1 0 0 0 -1 0 -1 1 0 0 0;1 0 -1 0 1 1 0 1 1 0 0 0 1 0 0 0 1 -1 -1 -1 -1 1 0 -1;-1 0 0 -1 0 0 1 -1 0 1 -1 -1 -1 1 1 0 1 1 1 0 -1 1 -1 -1;0 1 -1 -1 0 0 1 -1 -1 -1 0 -1 1 0 0 0 -1 0 0 1 0 0 0 -1;1 1 1 0 -1 -1 0 -1 -1 0 1 1 -1 0 1 -1 0 0 1 -1 0 0 0 -1;1 1 -1 -1 1 1 1 0 0 1 0 1 0 0 0 0 1 0 1 0 -1 1 0 0;1 1 1 -1 0 1 1 0 0 -1 1 -1 1 1 1 0 -1 -1 -1 -1 0 1 1 -1;0 0 -1 0 -1 1 0 -1 1 0 1 0 0 0 0 0 1 -1 0 0 0 1 -1 -1;1 -1 0 0 1 0 -1 0 -1 -1 0 0 1 0 0 -1 0 -1 -1 -1 -1 -1 1 -1;-1 -1 0 0 1 1 -1 -1 1 0 0 0 -1 0 0 -1 0 -1 -1 0 1 -1 0 0;0 -1 -1 -1 1 -1 1 0 -1 0 -1 1 0 1 -1 -1 1 -1 1 0 1 -1 1 -1;-1 -1 1 0 -1 0 0 0 -1 0 0 0 -1 0 0 -1 1 -1 -1 0 1 0 -1 -1;-1 0 0 0 -1 0 -1 0 -1 0 0 0 1 0 0 1 1 1 1 -1 -1 0 -1 -1;0 1 0 0 0 0 1 0 0 0 1 1 1 1 -1 0 0 1 -1 -1 -1 0 -1 -1;0 1 0 -1 -1 1 0 -1 1 -1 0 0 -1 -1 -1 0 0 -1 1 0 0 -1 -1 1;-1 0 1 1 -1 0 0 0 1 1 1 1 1 1 -1 -1 1 0 1 1 -1 -1 -1 1;0 0 0 1 -1 0 -1 -1 1 0 1 1 -1 1 -1 1 -1 -1 0 1 1 0 0 -1;0 -1 1 1 0 -1 1 0 1 0 1 0 0 0 1 1 0 -1 -1 0 0 0 1 0;-1 -1 -1 1 1 0 0 1 0 0 1 -1 -1 -1 1 1 0 1 -1 0 0 0 0 0;0 1 0 -1 -1 0 0 -1 -1 -1 1 1 1 0 0 0 1 1 0 0 0 0 1 0;-1 0 0 0 1 0 0 0 -1 1 -1 0 -1 1 1 1 1 1 0 -1 0 -1 0 1;-1 0 0 1 1 1 1 0 1 1 1 0 1 1 1 1 -1 -1 1 0 0 0 -1 0])

tic
i=1;
while(toc<60)
    tic
    haf(matrix(1:i,1:i));
    i=i+1;
end

Wypróbowałem to za pomocą TIO i MATLAB (nigdy nie instalowałem Octave). Wyobrażam sobie, że uruchomienie go jest tak proste, jak sudo apt-get install octave. Polecenie octavezaładuje graficzny interfejs użytkownika Octave. Jeśli będzie to bardziej skomplikowane, usunę tę odpowiedź, dopóki nie podam bardziej szczegółowych instrukcji instalacji.

Ostatnio Andreas Bjorklund, Brajesh Gupt i ja opublikowaliśmy nowy algorytm dla Hafnians złożonych matryc: https://arxiv.org/pdf/1805.12498.pdf . Dla macierzy n \ razy n skaluje się jak n ^ 3 2 ^ {n / 2}.

Jeśli dobrze rozumiem oryginalny algorytm Andreasa z https://arxiv.org/pdf/1107.4466.pdf , skaluje się jak n ^ 4 2 ^ {n / 2} lub n ^ 3 log (n) 2 ^ {n / 2} jeśli użyłeś transformacji Fouriera do wykonania mnożenia wielomianowego.
Nasz algorytm jest specjalnie opracowany dla złożonych macierzy, więc nie będzie tak szybki, jak te opracowane tutaj dla macierzy {-1,0,1}. Zastanawiam się jednak, czy można skorzystać z niektórych sztuczek użytych do ulepszenia naszej implementacji? Również, jeśli ludzie są zainteresowani, chciałbym zobaczyć, jak działa ich implementacja, gdy podano liczby zespolone zamiast liczb całkowitych. Wreszcie, wszelkie komentarze, krytyka, ulepszenia, błędy, ulepszenia są mile widziane w naszym repozytorium https://github.com/XanaduAI/hafnian/

Twoje zdrowie!