Calculez le Hafnian le plus rapidement possible

12

Le défi est d'écrire le code le plus rapide possible pour calculer le Hafnien d'une matrice .

Le Hafnian d'une matrice symétrique 2n-par- est défini comme:2nA

Ici S 2n représente l'ensemble de toutes les permutations des entiers de 1à 2n, c'est-à-dire [1, 2n].

Le lien wikipedia donne également une formule différente qui peut être intéressante (et des méthodes encore plus rapides existent si vous regardez plus loin sur le web). La même page wiki parle des matrices de contiguïté mais votre code devrait également fonctionner pour d'autres matrices. Vous pouvez supposer que les valeurs seront toutes des entiers mais pas qu'elles soient toutes positives.

Il existe également un algorithme plus rapide mais il semble difficile à comprendre. et Christian Sievers a été le premier à le mettre en œuvre (à Haskell).

Dans cette question, les matrices sont toutes carrées et symétriques avec une dimension paire.

Implémentation de référence (notez que cela utilise la méthode la plus lente possible).

Voici un exemple de code python de M. 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

La tâche

Vous devez écrire du code qui, donné 2npar une 2nmatrice, génère son Hafnian.

Comme je devrai tester votre code, il serait utile que vous me donniez un moyen simple de donner une matrice en entrée à votre code, par exemple en lisant à partir de standard in. Je testerai votre code dans des matrices choisies au hasard avec des éléments sélectionné parmi {-1, 0, 1}. Le but de tests comme celui-ci est de réduire les chances que le Hafnian soit d'une très grande valeur.

Idéalement, votre code sera capable de lire dans les matrices exactement comme je les ai dans les exemples de cette question directement à partir de la norme en. C'est-à-dire que l'entrée ressemblerait [[1,-1],[-1,-1]]par exemple. Si vous souhaitez utiliser un autre format d'entrée, veuillez demander et je ferai de mon mieux pour s'adapter.

Scores et égalités

Je vais tester votre code sur des matrices aléatoires de taille croissante et arrêter la première fois que votre code prend plus d'une minute sur mon ordinateur. Les matrices de notation seront cohérentes pour toutes les soumissions afin d'assurer l'équité.

Si deux personnes obtiennent le même score, le gagnant est celui qui est le plus rapide pour cette valeur de n. Si ceux-ci sont à moins d'une seconde les uns des autres, c'est celui affiché en premier.

Langues et bibliothèques

Vous pouvez utiliser n'importe quelle langue et bibliothèque disponibles, mais aucune fonction préexistante pour calculer le hafnien. Dans la mesure du possible, il serait bon de pouvoir exécuter votre code, veuillez donc inclure une explication complète sur la façon d'exécuter / compiler votre code sous Linux si possible. »

Ma machine Les temporisations seront exécutées sur ma machine 64 bits. Il s'agit d'une installation Ubuntu standard avec 8 Go de RAM, processeur AMD FX-8350 à huit cœurs et Radeon HD 4250. Cela signifie également que je dois pouvoir exécuter votre code.

Appel à réponses dans plus de langues

Ce serait formidable d'obtenir des réponses dans votre langage de programmation super rapide préféré. Pour commencer, que diriez-vous du fortran , du nim et de la rouille ?

Classement

  • 52 par miles en utilisant C ++ . 30 secondes.
  • 50 en utilisant ngn C . 50 secondes.
  • 46 par Christian Sievers utilisant Haskell . 40 secondes.
  • 40 miles en utilisant Python 2 + pypy . 41 secondes.
  • 34 par ngn en utilisant Python 3 + pypy . 29 secondes.
  • 28 par Dennis en utilisant Python 3 . 35 secondes. (Pypy est plus lent)

la source
Y a-t-il une limite pour les valeurs absolues des entrées de la matrice? Pouvons-nous renvoyer une approximation en virgule flottante? Faut-il utiliser des entiers de précision arbitraires?
Dennis
@Dennis En pratique, je n'utiliserai que -1,0,1 pour tester (choisi au hasard). Je ne veux pas que ce soit un gros défi. En toute honnêteté, je ne sais pas si nous atteindrons les limites de 64 bits avant que le code ne soit trop lent à exécuter, mais je suppose que nous ne le ferons pas. Actuellement, nous sommes loin de cela.
Si les entrées sont limitées à -1,0,1 , cela doit être mentionné sur la question. Notre code doit-il fonctionner du tout pour d'autres matrices?
Dennis
@Dennis Une ancienne version disait cela, mais je dois avoir écrit dessus. Je préférerais que le code ne soit pas spécialisé pour les entrées -1,0,1 mais je suppose que je ne peux pas arrêter ça.
Avez-vous plus de cas de test? peut-être pour un n plus grand ?
miles

Réponses:

14

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

Cela implémente une variante de l'algorithme 2 d' Andreas Björklund: compter les correspondances parfaites aussi vite que Ryser .

Compilez à l'aide des ghcoptions de compilation -O3 -threadedet utilisez les options +RTS -Nd' exécution pour la parallélisation. Prend l'entrée de stdin.

Christian Sievers
la source
2
Peut-être noter cela parallelet vectordoit être installé?
H.PWiz
@ H.PWiz Personne ne s'est plaint ici , mais bien sûr, notant que cela ne fera pas de mal. Eh bien, maintenant vous l'avez fait.
Christian Sievers
@ChristianSievers Je ne pense pas qu'ils se plaignent. L'OP peut ne pas être familier avec Haskell, il est donc judicieux d'indiquer explicitement ce qui doit être installé pour pouvoir chronométrer le code.
Dennis
@Dennis Je ne voulais pas dire "vous vous êtes plaint" mais "vous l'avez noté". Et je n'ai pas pensé à me plaindre comme une chose négative. Le PO est le même que dans le défi auquel j'ai lié, donc il ne devrait pas y avoir de problème.
Christian Sievers
N = 40 en 7,5 secondes sur TIO ... Mec, c'est rapide!
Dennis
6

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

Essayez-le en ligne!

Dennis
la source
6

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];
}

Essayez-le en ligne! (13s pour n = 24)

Basé sur la mise en œuvre plus rapide de Python dans mon autre article. Modifiez la deuxième ligne #define S 3sur votre machine à 8 cœurs et compilez avec g++ -pthread -march=native -O2 -ftree-vectorize.

Le divise le travail en deux, donc la valeur de Sdevrait être log2(#threads). Les types peuvent facilement être changés entre int, long, floatet doubleen modifiant la valeur #define TYPE.

miles
la source
C'est la principale réponse jusqu'à présent. Votre code ne lit pas réellement dans l'entrée comme spécifié car il ne peut pas gérer les espaces. Je devais faire par exempletr -d \ < matrix52.txt > matrix52s.txt
@Lembik Désolé, seulement utilisé contre la matrice sans espace de taille 24. Correction maintenant pour travailler avec les espaces.
miles
4

Python 3

Ceci calcule haf (A) comme une somme mémorisée (A [i] [j] * haf (A sans lignes et cols i et 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))
ngn
la source
3

C

Un autre élément du document d' Andreas Björklund , qui est beaucoup plus facile à comprendre si vous regardez également le code Haskell de Christian Sievers . Pour les premiers niveaux de la récursivité, il distribue les threads à tour de rôle sur les processeurs disponibles. Le dernier niveau de la récursivité, qui représente la moitié des invocations, est optimisé à la main.

Compiler avec: gcc -O3 -pthread -march=native; merci @Dennis pour une accélération 2x

n = 24 en 24s sur 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;
}

Algorithme:

La matrice, qui est symétrique, est stockée sous forme triangulaire en bas à gauche. Les indices triangulaires i,jcorrespondent à l'indice linéaire T(max(i,j))+min(i,j)Test une macro pour i*(i-1)/2. Les éléments matriciels sont des polynômes de degré n. Un polynôme est représenté comme un tableau de coefficients ordonnés du terme constant ( p[0]) au coefficient de x n (p[n] ). Les valeurs initiales de la matrice -1,0,1 sont d'abord converties en polynômes const.

Nous effectuons une étape récursive avec deux arguments: la demi-matrice (c'est-à-dire le triangle) ades polynômes et un polynôme séparé p(appelé bêta dans l'article). Nous réduisons le mproblème de taille (initialement m=2*n) récursivement à deux problèmes de taille m-2et retournons la différence de leurs hafniens. L'un d'eux est d'utiliser le même asans ses deux dernières lignes, et le même p. Une autre consiste à utiliser le triangle b[i][j] = a[i][j] + shmul(a[m-1][i],a[m-2][j]) + shmul(a[m-1][j],a[m-2][i])(où se shmultrouve l'opération de décalage-multiplication sur les polynômes - c'est comme le produit polynomial comme d'habitude, en plus multiplié par la variable "x"; les puissances supérieures à x ^ n sont ignorées), et le polynôme séparé q = p + shmul(p,a[m-1][m-2]). Lorsque la récursivité frappe une taille 0 a, nous revenons le grand coefficient de p: p[n].

L'opération shift-and-multiply est implémentée en fonction f(x,y,z). Il modifie zsur place. En gros, c'est le cas z += shmul(x,y). Cela semble être la partie la plus critique en termes de performances.

Une fois la récursivité terminée, nous devons fixer le signe du résultat en multipliant par (-1) n .

ngn
la source
Pourriez-vous montrer un exemple explicite de l'entrée que votre code accepte? Disons pour une matrice 2 par 2. De plus, vous semblez avoir joué à votre code! (C'est un défi de code le plus rapide, pas un défi de golf.)
@Lembik Pour mémoire, comme je l'ai dit dans le chat, l'entrée est au même format que les exemples - json (en fait, elle ne lit que les chiffres et utilise n = sqrt (len (entrée)) / 2). J'écris généralement du code court, même lorsque le golf n'est pas obligatoire.
2018
Quelle est la plus grande matrice de taille que ce nouveau code devrait prendre en charge?
1
-march=nativefera une grande différence ici. Au moins sur TIO, il réduit presque de moitié le temps de mur.
Dennis
1
De plus, au moins sur TIO, l'exécutable produit par gcc sera encore plus rapide.
Dennis
3

Python

Il s'agit d'une implémentation de référence simple de l'algorithme 2 du document mentionné . Les seules modifications consistaient à ne conserver que la valeur actuelle de B , à supprimer les valeurs de β en ne mettant à jour g que lorsque iX , et à multiplier les polynômes tronqués en calculant uniquement les valeurs jusqu'au degré 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

Essayez-le en ligne!

Voici une version plus rapide avec certaines des optimisations faciles.

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)

Essayez-le en ligne!

Pour plus de plaisir, voici une implémentation de référence en J.

miles
la source
C'est assez lent dans toutes les compréhensions de liste et dans le calcul de valeurs équivalentes sur la diagonale, il n'est donc pas nécessaire de comparer cela.
miles
Plutôt génial!
Très agréable! J'ai essayé une chose similaire avec sympy qui était incroyablement lente et, tout en étant correcte pour les petits exemples, a renvoyé - après un long moment - un mauvais résultat pour la matrice 24 * 24. Je n'ai aucune idée de ce qui se passe là-bas. - Par la description ci-dessus de l'algorithme 2, la multiplication polynomiale y est déjà destinée à être tronquée.
Christian Sievers
2
Dans pmul, utilisez for j in range(n-i):et évitez leif
Christian Sievers
1
@Lembik Il calcule la matrice entière; pour un facteur d'environ deux remplacer j != kpar j < k. Il copie une sous-matrice dans le cas contraire, ce qui peut être évité lorsque nous traitons et supprimons les deux derniers au lieu des deux premières lignes et colonnes. Et quand il calcule avec x={1,2,4}et plus tard avec x={1,2,4,6}alors il répète ses calculs jusqu'à i=5. J'ai remplacé la structure des deux boucles externes avec un premier bouclage ipuis en supposant récursivement i in Xet i not in X. - BTW, Il pourrait être intéressant de regarder la croissance du temps nécessaire par rapport aux autres programmes plus lents.
Christian Sievers
1

Octave

Il s'agit essentiellement d'une copie de l'entrée de Dennis , mais optimisée pour Octave. L'optimisation principale se fait en utilisant la matrice d'entrée complète (et sa transposition) et la récursivité en utilisant uniquement des indices de matrice, plutôt que de créer des matrices réduites.

Le principal avantage est une copie réduite des matrices. Bien qu'Octave n'ait pas de différence entre les pointeurs / références et les valeurs et ne fonctionne que par valeur, c'est une histoire différente dans les coulisses. Là, la copie sur écriture (copie paresseuse) est utilisée. Cela signifie que, pour le code a=1;b=a;b=b+1, la variable bn'est copiée vers un nouvel emplacement qu'à la dernière instruction, lorsqu'elle est modifiée. Depuis matinet matranspne sont jamais modifiés, ils ne seront jamais copiés. L'inconvénient est que la fonction passe plus de temps à calculer les bons indices. Je devrai peut-être essayer différentes variations entre les indices numériques et logiques pour optimiser cela.

Remarque importante: la matrice d'entrée devrait l'être int32! Enregistrez la fonction dans un fichier appeléhaf.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

Exemple de script de test:

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

J'ai essayé cela en utilisant TIO et MATLAB (je n'ai en fait jamais installé Octave). J'imagine que le faire fonctionner est aussi simple que sudo apt-get install octave. La commande octavechargera l'interface graphique d'Octave. Si c'est plus compliqué que cela, je supprimerai cette réponse jusqu'à ce que j'aie fourni des instructions d'installation plus détaillées.

Sanchises
la source
0

Récemment, Andreas Bjorklund, Brajesh Gupt et moi-même avons publié un nouvel algorithme pour les Hafniens de matrices complexes: https://arxiv.org/pdf/1805.12498.pdf . Pour une matrice n \ fois n, elle évolue comme n ^ 3 2 ^ {n / 2}.

Si je comprends bien l'algorithme original d'Andreas de https://arxiv.org/pdf/1107.4466.pdf, il évolue comme n ^ 4 2 ^ {n / 2} ou n ^ 3 log (n) 2 ^ {n / 2} si vous avez utilisé des transformées de Fourier pour effectuer des multiplications polynomiales.
Notre algorithme est spécifiquement adapté aux matrices complexes, il ne sera donc pas aussi rapide que ceux développés ici pour les matrices {-1,0,1}. Je me demande cependant si on peut utiliser certaines des astuces que vous avez utilisées pour améliorer notre implémentation? De plus, si les gens sont intéressés, j'aimerais voir comment leur implémentation fonctionne quand on leur donne des nombres complexes au lieu d'entiers. Enfin, tous les commentaires, critiques, améliorations, bugs, améliorations sont les bienvenus dans notre référentiel https://github.com/XanaduAI/hafnian/

À votre santé!

Nicolás Quesada
la source
Bienvenue sur le site! Cependant, les réponses à cette question doivent contenir du code. Ce serait mieux de laisser un commentaire, (que malheureusement vous n'avez pas de représentant à faire).
Ad Hoc Garf Hunter
Bienvenue chez PPCG. Bien que votre réponse puisse faire un bon commentaire, le site n'est pas destiné à l'assurance qualité. Il s'agit d'un site de défi et la réponse à un défi doit avoir du code et non une explication d'autre chose. Veuillez mettre à jour ou supprimer (si vous ne le faites pas, les mods le feront)
Muhammad Salman
Eh bien, le code est sur github, mais je suppose qu'il est absurde de simplement le copier-coller ici.
Nicolás Quesada
2
Si cela correspond à une réponse, surtout si vous êtes l'un des auteurs, je ne pense pas qu'il y ait quelque chose de mal à publier une solution compétitive correctement attribuée qui a été publiée ailleurs.
Dennis
@ NicolásQuesada Les réponses sur ce site doivent être autonomes si possible, ce qui signifie que nous ne devrions pas avoir à aller sur un autre site pour voir votre réponse / code.
mbomb007