Die Herausforderung besteht darin, den schnellstmöglichen Code für die Berechnung des Hafnian einer Matrix zu schreiben .

Die Hafnian einer symmetrischen 2n-by- 2nMatrix Aist definiert als:

Hier repräsentiert S _2n die Menge aller Permutationen der ganzen Zahlen von 1bis 2n, das heißt [1, 2n].

Der Wikipedia-Link gibt auch eine anders aussehende Formel an, die von Interesse sein kann (und es gibt noch schnellere Methoden, wenn Sie weiter im Web suchen). Die gleiche Wiki-Seite behandelt Adjazenzmatrizen, aber Ihr Code sollte auch für andere Matrizen funktionieren. Sie können davon ausgehen, dass alle Werte Ganzzahlen sind, aber nicht, dass sie alle positiv sind.

Es gibt auch einen schnelleren Algorithmus, der jedoch schwer zu verstehen ist. und Christian Sievers war der erste, der es umsetzte (in Haskell).

In dieser Frage sind die Matrizen alle quadratisch und symmetrisch mit einer gleichmäßigen Dimension.

Referenzimplementierung (beachten Sie, dass dies die langsamstmögliche Methode verwendet).

Hier ist ein Beispiel für Python-Code von 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([[-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

Die Aufgabe

Sie sollten Code schreiben, der, gegeben 2ndurch 2nMatrix, sein Hafnian ausgibt.

Da ich Ihren Code testen muss, wäre es hilfreich, wenn Sie mir eine einfache Möglichkeit geben könnten, eine Matrix als Eingabe für Ihren Code zu verwenden, z. B. durch Einlesen des Standards. Ich werde Ihren Code in zufällig ausgewählten Matrizen mit Elementen testen ausgewählt aus {-1, 0, 1}. Der Zweck solcher Tests ist es, die Wahrscheinlichkeit zu verringern, dass der Hafnianer einen sehr großen Wert hat.

Idealerweise kann Ihr Code Matrizen genau so einlesen, wie ich sie in den Beispielen in dieser Frage angegeben habe. So würde die Eingabe [[1,-1],[-1,-1]]beispielsweise aussehen . Wenn Sie ein anderes Eingabeformat verwenden möchten, fragen Sie bitte und ich werde mein Bestes tun, um dies zu berücksichtigen.

Partituren und Krawatten

Ich werde Ihren Code auf zufälligen Matrizen von zunehmender Größe testen und stoppen, wenn Ihr Code zum ersten Mal länger als 1 Minute auf meinem Computer dauert. Die Bewertungsmatrizen werden für alle Einreichungen konsistent sein, um die Fairness zu gewährleisten.

Wenn zwei Personen die gleiche Punktzahl erzielen, ist der Gewinner derjenige, der für diesen Wert von am schnellsten ist n. Wenn diese innerhalb einer Sekunde voneinander entfernt sind, ist dies die zuerst veröffentlichte.

Sprachen und Bibliotheken

Sie können jede verfügbare Sprache und Bibliothek verwenden, die Sie möchten, aber keine bereits vorhandene Funktion, um das Hafnian zu berechnen. Wo immer möglich, wäre es gut, wenn Sie Ihren Code ausführen könnten. Fügen Sie daher bitte eine vollständige Erklärung dazu bei, wie Sie Ihren Code unter Linux ausführen / kompilieren, wenn dies überhaupt möglich ist. "

Mein Computer Die Timings werden auf meinem 64-Bit-Computer ausgeführt. Dies ist eine Ubuntu-Standardinstallation mit 8 GB RAM, AMD FX-8350 Eight-Core-Prozessor und Radeon HD 4250. Dies bedeutet auch, dass ich in der Lage sein muss, Ihren Code auszuführen.

Fordern Sie Antworten in mehr Sprachen

Es wäre toll, Antworten in Ihrer bevorzugten superschnellen Programmiersprache zu erhalten. So starten Sie die Dinge aus, wie etwa Fortran , nim und Rost ?

Bestenliste

• 52 Meilen mit C ++ . 30 Sekunden.
• 50 durch NGN Verwendung C . 50 Sekunden.
• 46 von Christian Sievers mit Haskell . 40 Sekunden.
• 40 Meilen mit Python 2 + Pypy . 41 Sekunden.
• 34 von ngn mit Python 3 + pypy . 29 Sekunden.
• 28 von Dennis mit Python 3 . 35 Sekunden. (Pypy ist langsamer)

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

Dies implementiert eine Variation von Algorithmus 2 von Andreas Björklund: Das Zählen perfekter Übereinstimmungen so schnell wie Ryser .

Kompilieren Sie ghcmit den -O3 -threadedOptionen +RTS -Nfür die Kompilierungszeit und verwenden Sie die Laufzeitoptionen für die Parallelisierung. Übernimmt die Eingabe von stdin.

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

Probieren Sie es 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];
}

Probieren Sie es online! (13s für n = 24)

Basiert auf der schnelleren Python- Implementierung in meinem anderen Beitrag. Bearbeiten Sie die zweite Zeile zu#define S 3 auf Ihrem 8-Core-Rechner und kompilieren Sie mit g++ -pthread -march=native -O2 -ftree-vectorize.

Das halbiert die Arbeit, also sollte der Wert von Ssein log2(#threads). Die Typen können leicht zwischen geändert werden int, long, floatund doubledurch den Wert des Modifizieren #define TYPE.

Python 3

Dies berechnet haf (A) als memoisierte Summe (A [i] [j] * haf (A ohne Zeilen und Spalten i und 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))

C

Ein weiteres Werkzeug aus Andreas Björklunds Arbeit , das viel einfacher zu verstehen ist, wenn man sich auch Christian Sievers 'Haskell-Code ansieht . In den ersten Stufen der Rekursion werden Round-Robin-Threads auf die verfügbaren CPUs verteilt. Die letzte Stufe der Rekursion, die die Hälfte der Aufrufe ausmacht, wird von Hand optimiert.

Kompilieren mit: gcc -O3 -pthread -march=native; danke @Dennis für ein 2x Speed-up

n = 24 in 24s bei 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;
}

Algorithmus:

Die symmetrische Matrix wird in Form eines Dreiecks unten links gespeichert. Dreieck-Indizes i,jentsprechen dem linearen Index, T(max(i,j))+min(i,j)für den Tein Makro steht i*(i-1)/2. Matrixelemente sind Polynome vom Grad n. Ein Polynom wird dargestellt als eine Anordnung von Koeffizienten aus dem konstanten Term (geordnete p[0]) bis x ⁿ ‚s Koeffizienten ( p[n]). Die anfänglichen -1,0,1-Matrixwerte werden zuerst in konstante Polynome umgewandelt.

Wir führen einen rekursiven Schritt mit zwei Argumenten durch: der Halbmatrix (dh dem Dreieck) avon Polynomen und einem separaten Polynom p(im Artikel als Beta bezeichnet). Wir reduzieren das Größenproblem m(anfangs m=2*n) rekursiv auf zwei Größenprobleme m-2und geben die Differenz ihrer Hafnier zurück. Eine davon ist, dass Sie dasselbe aohne die letzten beiden Zeilen und dasselbe verwenden p. Eine andere Möglichkeit ist die Verwendung des Dreiecks b[i][j] = a[i][j] + shmul(a[m-1][i],a[m-2][j]) + shmul(a[m-1][j],a[m-2][i])(bei der shmuldie Verschiebungs-Multiplikations-Operation für Polynome ausgeführt wird - es ist wie üblich ein Polynomprodukt, zusätzlich multipliziert mit der Variablen "x"; Potenzen über x ^ n werden ignoriert) und des separaten Polynoms q = p + shmul(p,a[m-1][m-2]). Wenn die Rekursion eine Größe von 0 erreicht a, geben wir den Hauptkoeffizienten von p: zurück p[n].

Die Verschiebe- und Multiplikationsoperation ist in Funktion implementiert f(x,y,z). Es ändert sich zan Ort und Stelle. Im Grunde genommen schon z += shmul(x,y). Dies scheint der leistungskritischste Teil zu sein.

Nachdem die Rekursion beendet ist, müssen wir das Vorzeichen des Ergebnisses durch Multiplikation mit (-1) ^{n korrigieren} .

Python

Dies ist so ziemlich eine direkte Referenzimplementierung von Algorithmus 2 aus dem erwähnten Artikel . Die einzigen Änderungen bestanden darin, nur den aktuellen Wert von B beizubehalten, die Werte von β zu löschen, indem nur g aktualisiert wurde, wenn i ∈ X ist , und die Multiplikation mit abgeschnittenen Polynomen durchzuführen, indem nur die Werte bis zum Grad n berechnet wurden .

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

Probieren Sie es online!

Hier ist eine schnellere Version mit einigen der einfachen Optimierungen.

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)

Probieren Sie es online!

Für noch mehr Spaß, hier ist eine Referenz Implementierung in J.

Oktave

Dies ist im Grunde eine Kopie von Dennis 'Eintrag , aber optimiert für Octave. Die Hauptoptimierung erfolgt unter Verwendung der vollständigen Eingabematrix (und ihrer Transponierung) und Rekursion, wobei nur Matrixindizes verwendet werden, anstatt reduzierte Matrizen zu erstellen.

Der Hauptvorteil ist das reduzierte Kopieren von Matrizen. Während Octave keinen Unterschied zwischen Zeigern / Referenzen und Werten hat und funktional nur einen vorübergehenden Wert hat, ist es eine andere Geschichte hinter den Kulissen. Dort wird Copy-on-Write (Lazy Copy) verwendet. Das heißt, für den Code a=1;b=a;b=b+1wird die Variable berst bei der letzten Anweisung an eine neue Position kopiert, wenn sie geändert wird. Da matinund matranspwerden nie geändert, werden sie nie kopiert. Der Nachteil ist, dass die Funktion mehr Zeit für die Berechnung der korrekten Indizes benötigt. Möglicherweise muss ich verschiedene Variationen zwischen numerischen und logischen Indizes ausprobieren, um dies zu optimieren.

Wichtiger Hinweis: Eingabematrix sollte sein int32! Speichern Sie die Funktion in einer Datei mit dem Namenhaf.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

Beispiel-Testskript:

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

Ich habe das mit TIO und MATLAB ausprobiert (Octave habe ich eigentlich nie installiert). Ich kann mir vorstellen, dass es so einfach ist, es zum Laufen zu bringen sudo apt-get install octave. Der Befehl octavelädt die Octave-GUI. Wenn es komplizierter ist, werde ich diese Antwort löschen, bis ich eine ausführlichere Installationsanleitung zur Verfügung gestellt habe.

Kürzlich haben Andreas Bjorklund, Brajesh Gupt und ich einen neuen Algorithmus für Hafnianer komplexer Matrizen veröffentlicht: https://arxiv.org/pdf/1805.12498.pdf . Für eine n \ times n-Matrix skaliert sie wie n ^ 3 2 ^ {n / 2}.

Wenn ich den ursprünglichen Algorithmus von Andreas von https://arxiv.org/pdf/1107.4466.pdf richtig verstehe, wird er wie folgt skaliert: n ^ 4 2 ^ {n / 2} oder n ^ 3 log (n) 2 ^ {n / 2} Wenn Sie Fouriertransformationen für Polynommultiplikationen verwendet haben.
Unser Algorithmus ist speziell auf komplexe Matrizen zugeschnitten, sodass er nicht so schnell ist wie die hier für {-1,0,1} Matrizen entwickelten. Ich frage mich jedoch, ob man einige der Tricks anwenden kann, mit denen Sie unsere Implementierung verbessert haben. Auch wenn die Leute interessiert sind, würde ich gerne sehen, wie sich ihre Implementierung verhält, wenn komplexe Zahlen anstelle von ganzen Zahlen verwendet werden. Schließlich sind Kommentare, Kritik, Verbesserungen, Bugs und Verbesserungen in unserem Repository https://github.com/XanaduAI/hafnian/ willkommen.

Prost!