Hitung Hafnian secepat mungkin

12

Tantangannya adalah untuk menulis kode tercepat yang mungkin untuk menghitung Hafnian dari sebuah matriks .

The Hafnian dari simetris 2n-by- 2nmatriks Adidefinisikan sebagai:

Di sini S 2n mewakili himpunan semua permutasi bilangan bulat dari 1ke 2n, yaitu [1, 2n].

Tautan wikipedia juga memberikan rumus tampilan berbeda yang mungkin menarik (dan bahkan metode yang lebih cepat ada jika Anda melihat lebih jauh di web). Halaman wiki yang sama berbicara tentang matriks kedekatan tetapi kode Anda juga bisa digunakan untuk matriks lain. Anda dapat mengasumsikan nilai-nilai semua akan menjadi bilangan bulat tetapi bukan berarti semuanya positif.

Ada juga algoritma yang lebih cepat tetapi tampaknya sulit untuk dipahami. dan Christian Sievers adalah yang pertama mengimplementasikannya (dalam Haskell).

Dalam pertanyaan ini semua matriks berbentuk bujur sangkar dan simetris dengan dimensi genap.

Implementasi referensi (perhatikan ini menggunakan metode paling lambat yang mungkin).

Berikut ini beberapa contoh kode python dari 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

Tugas

Anda harus menulis kode itu, diberikan 2noleh 2nmatriks, menampilkan Hafniannya.

Karena saya perlu menguji kode Anda, akan sangat membantu jika Anda dapat memberikan cara sederhana bagi saya untuk memberikan matriks sebagai input ke kode Anda, misalnya dengan membaca dari standar. Saya akan menguji kode Anda dalam matriks yang dipilih secara acak dengan elemen. dipilih dari {-1, 0, 1}. Tujuan pengujian seperti ini adalah untuk mengurangi peluang Hafnian akan menjadi nilai yang sangat besar.

Idealnya kode Anda akan dapat dibaca dalam matriks persis seperti yang saya miliki dalam contoh-contoh dalam pertanyaan ini langsung dari standar masuk. Itu adalah input akan terlihat seperti [[1,-1],[-1,-1]]misalnya. Jika Anda ingin menggunakan format input lain, tanyakan dan saya akan melakukan yang terbaik untuk mengakomodasi.

Skor dan dasi

Saya akan menguji kode Anda pada matriks acak dengan ukuran yang meningkat dan menghentikan pertama kali kode Anda membutuhkan lebih dari 1 menit di komputer saya. Matriks penilaian akan konsisten untuk semua pengiriman untuk memastikan keadilan.

Jika dua orang mendapatkan skor yang sama maka pemenangnya adalah yang tercepat untuk nilai tersebut n. Jika itu adalah dalam 1 detik satu sama lain maka itu adalah yang diposting pertama.

Bahasa dan perpustakaan

Anda dapat menggunakan bahasa dan pustaka yang tersedia yang Anda suka tetapi tidak ada fungsi yang sudah ada sebelumnya untuk menghitung Hafnian. Jika memungkinkan, alangkah baiknya untuk dapat menjalankan kode Anda jadi harap sertakan penjelasan lengkap tentang cara menjalankan / kompilasi kode Anda di Linux jika memungkinkan. `

Mesin Saya Pengaturan waktu akan dijalankan pada mesin 64-bit saya. Ini adalah instalasi ubuntu standar dengan RAM 8GB, Prosesor Delapan-Core AMD FX-8350 dan Radeon HD 4250. Ini juga berarti saya harus dapat menjalankan kode Anda.

Panggil jawaban dalam lebih banyak bahasa

Akan sangat bagus untuk mendapatkan jawaban dalam bahasa pemrograman super cepat favorit Anda. Untuk memulai, bagaimana dengan fortran , nim dan rust ?

Papan peringkat

  • 52 oleh miles menggunakan C ++ . 30 detik.
  • 50 oleh ngn menggunakan C . 50 detik.
  • 46 oleh Christian Sievers menggunakan Haskell . 40 detik.
  • 40 oleh miles menggunakan Python 2 + pypy . 41 detik.
  • 34 oleh ngn menggunakan Python 3 + pypy . 29 detik.
  • 28 oleh Dennis menggunakan Python 3 . 35 detik. (Pypy lebih lambat)

sumber
Apakah ada batasan untuk nilai absolut dari entri matriks? Bisakah kita mengembalikan perkiraan floating point? Apakah kita harus menggunakan bilangan bulat presisi acak?
Dennis
@Dennis Dalam prakteknya saya hanya akan menggunakan -1,0,1 untuk menguji (dipilih secara acak). Saya tidak ingin ini menjadi tantangan besar. Dalam semua kejujuran saya tidak tahu apakah kita akan mencapai batas int 64 bit sebelum kode terlalu lambat untuk berjalan tetapi saya duga adalah bahwa kita tidak akan melakukannya. Saat ini kami tidak berada di dekat itu.
Jika entri dibatasi hingga -1,0,1 , yang harus disebutkan pada pertanyaan. Apakah kode kita harus bekerja sama sekali untuk matriks lain?
Dennis
@ Dennis Versi lama yang digunakan untuk mengatakan itu tetapi saya pasti telah menulisnya. Saya lebih suka jika kode tidak khusus untuk entri -1,0,1 tapi saya kira saya tidak bisa menghentikannya.
Apakah Anda memiliki lebih banyak test case? mungkin untuk n lebih besar ?
mil

Jawaban:

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

Ini mengimplementasikan variasi Algoritma 2 dari Andreas Björklund: Menghitung Pencocokan Sempurna secepat Ryser .

Kompilasi menggunakan ghcdengan opsi waktu kompilasi -O3 -threadeddan gunakan opsi run time +RTS -Nuntuk paralelisasi. Mengambil input dari stdin.

Sievers Kristen
sumber
2
Mungkin perhatikan itu paralleldan vectorharus diinstal?
H.PWiz
@ H.PWiz Tidak ada yang mengeluh di sini , tapi tentu saja, mencatat itu tidak akan menyakitkan. Nah, sekarang sudah.
Christian Sievers
@ChristianSievers Saya tidak berpikir mereka mengeluh. OP mungkin tidak terbiasa dengan Haskell, jadi secara eksplisit menyatakan apa yang harus diinstal untuk dapat menentukan waktu kode adalah ide yang baik.
Dennis
@ Dennis Saya tidak bermaksud "Anda memang mengeluh" tetapi "Anda memang mencatatnya". Dan saya tidak berpikir untuk mengeluh sebagai hal negatif. OP sama dengan tantangan yang saya tautkan, jadi seharusnya tidak ada masalah.
Christian Sievers
N = 40 dalam 7,5 detik pada TIO ... Ya ampun, ini cepat!
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())))))

Cobalah online!

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

Cobalah online! (13d untuk n = 24)

Berdasarkan implementasi Python yang lebih cepat di posting saya yang lain. Edit baris kedua ke #define S 3pada mesin 8-core Anda dan kompilasi dengan g++ -pthread -march=native -O2 -ftree-vectorize.

Membagi pekerjaan menjadi dua, jadi nilai Sseharusnya log2(#threads). Jenis dapat dengan mudah diubah antara int, long, float, dan doubledengan memodifikasi nilai #define TYPE.

mil
sumber
Ini adalah jawaban utama sejauh ini. Kode Anda tidak benar-benar membaca input yang ditentukan karena tidak dapat mengatasi spasi. Saya harus melakukan misalnyatr -d \ < matrix52.txt > matrix52s.txt
@ Lembik Maaf, hanya menggunakannya terhadap matriks tanpa spasi ukuran 24. Tetap sekarang untuk bekerja dengan spasi.
mil
4

Python 3

Ini menghitung haf (A) sebagai jumlah memoised (A [i] [j] * haf (A tanpa baris dan cols i dan 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
sumber
3

C

Impl lain dari makalah Andreas Björklund , yang jauh lebih mudah dipahami jika Anda juga melihat kode Haskell Christian Sievers . Untuk beberapa level pertama rekursi, ia mendistribusikan putaran-robin di atas CPU yang tersedia. Level terakhir dari rekursi, yang merupakan setengah dari doa, dioptimalkan dengan tangan.

Mengkompilasi dengan: gcc -O3 -pthread -march=native; terima kasih @Dennis untuk 2x percepatan

n = 24 dalam 24 detik pada 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;
}

Algoritma:

Matriks, yang simetris, disimpan dalam bentuk segitiga kiri bawah. Indeks segitiga i,jsesuai dengan indeks linier di T(max(i,j))+min(i,j)mana Tmerupakan makro untuk i*(i-1)/2. Elemen matriks adalah polinomial derajat n. Polinomial direpresentasikan sebagai larik koefisien yang diurutkan dari suku konstan ( p[0]) ke koefisien x n ( p[n]). Nilai matriks -1,0,1 awal dikonversi terlebih dahulu ke const polinomials.

Kami melakukan langkah rekursif dengan dua argumen: setengah-matriks (yaitu segitiga) adari polinomial dan polinomial terpisah p(disebut beta dalam makalah). Kami mengurangi mmasalah ukuran (awalnya m=2*n) secara rekursif menjadi dua masalah ukuran m-2dan mengembalikan perbedaan hafnianya. Salah satunya adalah menggunakan yang sama atanpa dua baris terakhir, dan yang sama p. Lain adalah menggunakan segitiga b[i][j] = a[i][j] + shmul(a[m-1][i],a[m-2][j]) + shmul(a[m-1][j],a[m-2][i])(di mana shmuladalah operasi shift-multiply pada polinomial - itu seperti produk polinom seperti biasa, juga dikalikan dengan variabel "x"; kekuatan yang lebih tinggi dari x ^ n diabaikan), dan polinomial yang terpisah q = p + shmul(p,a[m-1][m-2]). Ketika rekursi hits ukuran-0 a, kita kembali koefisien utama p: p[n].

Operasi shift-and-multiply diimplementasikan dalam fungsi f(x,y,z). Itu memodifikasi zdi tempat. Secara longgar, itu benar z += shmul(x,y). Ini tampaknya menjadi bagian yang paling kritis terhadap kinerja.

Setelah rekursi selesai, kita perlu memperbaiki tanda hasil dengan mengalikannya dengan (-1) n .

ngn
sumber
Bisakah Anda menunjukkan contoh eksplisit dari input yang diterima oleh kode Anda? Katakanlah untuk matriks 2 dengan 2. Juga, Anda tampaknya telah memodernisasi kode Anda! (Ini tantangan kode tercepat, bukan tantangan golf.)
@Lembik Sebagai catatan, seperti yang saya katakan dalam obrolan, inputnya dalam format yang sama dengan contoh - json (sebenarnya, ia hanya membaca angka dan menggunakan n = sqrt (len (input)) / 2). Saya biasanya menulis kode pendek, bahkan ketika bermain golf bukan keharusan.
ngn
Apa ukuran matriks terbesar yang harus didukung kode baru ini?
1
-march=nativeakan membuat perbedaan besar di sini. Setidaknya pada TIO, hampir memotong waktu dinding menjadi dua.
Dennis
1
Juga, setidaknya pada TIO, executable yang diproduksi oleh gcc akan lebih cepat.
Dennis
3

Python

Ini cukup mudah, implementasi referensi Algoritma 2 dari makalah yang disebutkan . Satu-satunya perubahan adalah hanya menjaga nilai B saat ini , menjatuhkan nilai β dengan hanya memperbarui g ketika iX , dan memotong perkalian polinomial dengan hanya menghitung nilai hingga derajat 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

Cobalah online!

Ini adalah versi yang lebih cepat dengan beberapa optimasi yang mudah.

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)

Cobalah online!

Untuk tambahan kesenangan, berikut ini adalah implementasi referensi dalam J.

mil
sumber
Ini cukup lambat dari semua pemahaman daftar dan dari penghitungan nilai yang setara di seluruh diagonal, jadi tidak perlu membandingkan ini.
mil
Cukup mengagumkan!
Sangat bagus! Saya mencoba hal serupa dengan sympy yang sangat lambat dan, sementara benar untuk contoh-contoh kecil, kembali - setelah waktu yang lama - hasil yang salah untuk matriks 24 * 24. Saya tidak tahu apa yang terjadi di sana. - Dengan uraian di atas Algoritma 2, perkalian polinomial sudah ada yang dimaksudkan untuk dipotong.
Christian Sievers
2
Masuk pmul, gunakan for j in range(n-i):dan hindariif
Christian Sievers
1
@ Lembik Ini menghitung seluruh matriks; untuk faktor sekitar dua ganti j != kdengan j < k. Ini menyalin submatrix dalam kasus lain yang dapat dihindari ketika kita menangani dan menghapus dua terakhir bukan dua baris dan kolom pertama. Dan ketika ia menghitung dengan x={1,2,4}dan kemudian dengan x={1,2,4,6}itu mengulangi perhitungannya hingga i=5. Saya mengganti struktur dua loop luar dengan looping pertama idan kemudian secara rekursif mengasumsikan i in Xdan i not in X. - BTW, Mungkin menarik untuk melihat pertumbuhan waktu yang dibutuhkan dibandingkan dengan program lain yang lebih lambat.
Christian Sievers
1

Oktaf

Ini pada dasarnya adalah salinan entri Dennis , tetapi dioptimalkan untuk Oktaf. Optimalisasi utama dilakukan dengan menggunakan matriks input penuh (dan transposinya) dan rekursi hanya menggunakan indeks matriks, daripada membuat matriks yang dikurangi.

Keuntungan utama adalah mengurangi penyalinan matriks. Sementara Oktaf tidak memiliki perbedaan antara pointer / referensi dan nilai-nilai dan secara fungsional hanya melewati nilai, itu cerita yang berbeda di belakang layar. Di sana, copy-on-write (lazy copy) digunakan. Itu berarti, untuk kode a=1;b=a;b=b+1, variabel bhanya disalin ke lokasi baru di pernyataan terakhir, ketika diubah. Karena matindan matransptidak pernah berubah, mereka tidak akan pernah dengan cara disalin. Kerugiannya adalah bahwa fungsi menghabiskan lebih banyak waktu menghitung indeks yang benar. Saya mungkin harus mencoba variasi yang berbeda antara indeks numerik dan logis untuk mengoptimalkan ini.

Catatan penting: matriks input seharusnya int32! Simpan fungsi dalam file yang disebuthaf.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

Contoh skrip pengujian:

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

Saya sudah mencoba menggunakan TIO dan MATLAB (sebenarnya saya tidak pernah menginstal Octave). Saya membayangkan membuatnya bekerja sesederhana itu sudo apt-get install octave. Perintah octaveakan memuat GUI Oktaf. Jika lebih rumit dari ini, saya akan menghapus jawaban ini sampai saya memberikan instruksi instalasi yang lebih rinci.

Sanchises
sumber
0

Baru-baru ini Andreas Bjorklund, Brajesh Gupt dan saya sendiri menerbitkan algoritma baru untuk Hafnians dari matriks kompleks: https://arxiv.org/pdf/1805.12498.pdf . Untuk matriks n \ kali n, matriks akan berskala seperti n ^ 3 2 ^ {n / 2}.

Jika saya mengerti dengan benar, algoritma asli Andreas dari https://arxiv.org/pdf/1107.4466.pdf berskala seperti n ^ 4 2 ^ {n / 2} atau n ^ 3 log (n) 2 ^ {n / 2} jika Anda menggunakan transformasi Fourier untuk melakukan perkalian polinomial.
Algoritme kami secara khusus disoroti untuk matriks kompleks sehingga tidak akan secepat yang dikembangkan di sini untuk matriks {-1,0,1}. Namun saya ingin tahu apakah seseorang dapat menggunakan beberapa trik yang Anda gunakan untuk meningkatkan implementasi kami? Juga jika orang tertarik saya ingin melihat bagaimana implementasi mereka lakukan ketika diberi bilangan kompleks bukan bilangan bulat. Akhirnya, setiap komentar, kritik, perbaikan, bug, perbaikan disambut di repositori kami https://github.com/XanaduAI/hafnian/

Bersulang!

Nicolás Quesada
sumber
Selamat datang di situs ini! Namun jawaban untuk pertanyaan ini harus mengandung kode. Ini akan lebih baik dibiarkan sebagai komentar, (Yang sayangnya Anda tidak memiliki perwakilan untuk membuat).
Ad Hoc Garf Hunter
Selamat datang di PPCG. Meskipun jawaban Anda dapat memberikan komentar yang bagus, situs ini bukan untuk QA. Ini adalah situs tantangan dan jawaban untuk tantangan harus memiliki kode dan bukan penjelasan tentang hal lain. Silakan perbarui atau hapus (Jika tidak, mod akan)
Muhammad Salman
Yah, kodenya ada di github, tapi kurasa tidak masuk akal untuk hanya menyalin-menempelnya di sini.
Nicolás Quesada
2
Jika cocok dengan jawaban, terutama jika Anda salah satu penulisnya, saya tidak berpikir ada yang salah dengan memposting solusi kompetitif yang dikaitkan dengan benar yang telah diterbitkan di tempat lain.
Dennis
@ NicolásQuesada Jawaban di situs ini harus lengkap jika memungkinkan, artinya kita tidak harus pergi ke situs lain untuk melihat jawaban / kode Anda.
mbomb007