El desafío es escribir codegolf para el hafniano de una matriz . El Hafnian de un 2n-by- 2nmatriz simétrica Ase define como:

Aquí S _2n representa el conjunto de todas las permutaciones de los enteros de 1a 2n, es decir [1, 2n].

El enlace de wikipedia habla sobre las matrices de adyacencia, pero su código debería funcionar para cualquier matriz de entrada simétrica de valor real.

Para aquellos interesados ​​en las aplicaciones de Hafnian, el enlace mathoverflow discute un poco más.

Su código puede recibir información de la forma que desee y proporcionar resultados en cualquier formato razonable, pero incluya en su respuesta un ejemplo completo que incluya instrucciones claras sobre cómo suministrar información a su código.

La matriz de entrada siempre es cuadrada y tendrá como máximo 16 por 16. No hay necesidad de poder manejar la matriz vacía o las matrices de dimensión impar.

Implementación de referencia

Aquí hay un ejemplo de código python del Sr. 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([[0, 4.5], [4.5, 0]]))
4.5
print(hafnian([[0, 4.7, 4.6, 4.5], [4.7, 0, 2.1, 0.4], [4.6, 2.1, 0, 1.2], [4.5, 0.4, 1.2, 0]])
16.93
print(hafnian([[1.3, 4.1, 1.2, 0.0, 0.9, 4.4], [4.1, 4.2, 2.7, 1.2, 0.4, 1.7], [1.2, 2.7, 4.9, 4.7, 4.0, 3.7], [0.0, 1.2, 4.7, 2.2, 3.3, 1.8], [0.9, 0.4, 4.0, 3.3, 0.5, 4.4], [4.4, 1.7, 3.7, 1.8, 4.4, 3.2]])
262.458

La página wiki ahora (2 de marzo de 2018) ha sido actualizada por ShreevatsaR para incluir una forma diferente de calcular el Hafnian. Sería muy interesante ver esto golfizado.

R , 150 142 127 119 bytes

function(A,N=nrow(A),k=1:(N/2)*2)sum(apply(gtools::permutations(N,N),1,function(r)prod(A[cbind(r[k-1],r[k])])))/prod(k)

Pruébalo en línea!

Utiliza el mismo truco que descubrí jugando golf en esta respuesta para indexar la matriz P, ¡y @Vlo sugirió un enfoque para eliminar por completo el forbucle para -6 bytes!

Para crear un nuevo caso de prueba, puede hacerlo matrix(c(values,separated,by,commas,going,across,rows),nrow=2n,ncol=2n,byrow=T).

Explicación: (el código es el mismo; utiliza un bucle en applylugar de un forbucle, pero la lógica es idéntica).

function(A){
N <- nrow(A)                   #N = 2*n
k <- 1:(N/2) * 2               #k = c(2,4,...,N) -- aka 2*j in the formula
P <- gtools::permutations(N,N) #all N-length permutations of 1:N
for(i in 1:nrow(P))
 F <- F + prod(A[cbind(P[i,k-1],P[i,k])]) # takes the product of all A_sigma(2j-1)sigma(2j) for fixed i and adds it to F (initialized to 0)
F / prod(k)                    #return value; prod(k) == n! * 2^n
}

Pyth , 24 bytes

sm*Fmc@@Qhkek2d{mScd2.pU

Pruébalo aquí!

Versión anterior, 35 bytes

*c1**FK/lQ2^2Ksm*Fm@@Q@[email protected]

Pruébalo aquí!

Stax , 23 22 19 17 bytes

ü;Y╙◘▌Φq↓ê²╧▐å↑┌C

Ejecútelo y depúrelo en línea

La representación ascii correspondiente del mismo programa es esta.

%r|TF2/{xsE@i^H/m:*+

El programa sufre algún error de redondeo de coma flotante. En particular, informa en 33673.5000000011lugar de 33673.5. Pero creo que la precisión es aceptable dado que este programa funciona con valores de coma flotante. También es muy lento, toma casi un minuto para las entradas de ejemplo en esta máquina.

%                             get size of matrix
 r|T                          get all permutations of [0 ... size-1]
    F                         for each, execute the rest of the program
     2/                       get consecutive pairs
       {        m             map each pair... 
        xsE@                      the matrix element at that location
            i^H/                  divided by 2*(i+1) where i=iteration index
                 :*           product of array
                   +          add to running total

05AB1E , 21 bytes

ā<œε2ô{}Ùεε`Isèsè;]PO

Pruébalo en línea!

Versión anterior, 32 bytes

āœvIg;©Lε·UIyX<èèyXèè}P}Oθ®!/®o/

Pruébalo en línea!

¿Cómo funciona?

āœvIg;©Lε·UIyX<èèyXèè}P}Oθ®!/®o/ – Full program. Argument: A matrix M.
ā                                – The range [1 ... len(M)].
 œ                               – Permutations.
  v                    }         – Iterate over the above with a variable y.
   Ig;©                          – Push len(M) / 2 and also store it in register c.
       Lε            }           – For each integer in the range [1 ... ^]:
         ·U                      – Double it and store it in a variable X.
            yX<                  – Push the element of y at index X-1.
           I   è                 – And index with the result into M.
                yXè              – Push the element of y at index X.
                   è             – And index with the result into ^^.
                      P          – Take the product of the resulting list.
                        O        – Sum the result of the mapping.
                         θ       – And take the last element*.
                          ®!     – Take the factorial of the last item in register c.
                             ®o  – Raise 2 to the power of the last item in register c.
                            /  / – And divide the sum of the mapping accordingly.

* – Yeah, this is needed because I mess up the stack when pushing so many values in the loop and not popping correctly ;P

Jalea , 19 bytes

LŒ!s€2Ṣ€QḅL_LịFHP€S

Pruébalo en línea!

Versión alternativa, 15 bytes, desafío de fecha posterior

LŒ!s€2Ṣ€QœịHP€S

Jelly finalmente obtuvo la indexación de matriz n-dimensional.

Pruébalo en línea!

Cómo funciona

LŒ!s€2Ṣ€QœiHP€S  Main link. Argument: M (matrix / 2D array)

L                Take the length, yielding 2n.
 Œ!              Generate all permutations of [1, ..., 2n].
   s€2           Split each permutation into pairs.
      Ṣ€         Sort the pair arrays.
        Q        Unique; deduplicate the array of pair arrays.
                 This avoids dividing by n! at the end.
           H     Halve; yield M, with all of its elements divided by 2.
                 This avoids dividing by 2**n at the end.
         œị      At-index (n-dimensional); take each pair of indices [i, j] and
                 yield M[i][j].
            P€   Take the product the results corresponding the same permutation.
              S  Take the sum of the products.

La versión de 19 bytes funciona de manera similar; solo tiene que implementarse œị.

...ḅL_LịFH...    Return value: Array of arrays of index pairs. Argument: M

    L            Length; yield 2n.
   ḅ             Convert each pair of indices [i, j] from base 2n to integer,
                 yielding ((2n)i + j).
     _L          Subtract 2n, yielding ((2n)(i - 1) + j).
                 This is necessary because indexing is 1-based in Jelly, so the
                 index pair [1, 1] must map to index 1.
        F        Yield M, flattened.
       ị         Take the indices to the left and get the element at these indices
                 from the array to the right.
         H       Halve; divide all retrieved elements by 2.

C (gcc) , 288 285 282 293 292 272 271 bytes

• Ahorró tres bytes jugando con dos incrementos posteriores y para la colocación de bucles.
• Guardado tres bytes por jugando con anothter de incremento posterior, moviéndose ambas inicializaciones variables antes de la rama - golfed if(...)...k=0...else...,j=0...a if(k=j=0,...)...else...- y realizado un cambio en el índice.
• Once bytes requeridos por floatmatrices de soporte .
• Ahorré un byte gracias al Sr. Xcoder ; golf 2*j+++1a j-~j++.
• Ahorró veinte bytes eliminando una intdeclaración de tipo variable superflua y no utilizando una función factorial, sino calculando el valor factorial utilizando un bucle for ya existente.
• Salvó un byte jugando S=S/F/(1<<n);al golf S/=F*(1<<n);.

float S,p,F;j,i;s(A,n,P,l,o,k)float*A;int*P;{if(k=j=0,o-l)for(;k<l;s(A,n,P,l,o+1))P[o]=k++;else{for(p=-l;j<l;j++)for(i=0;i<l;)p+=P[j]==P[i++];if(!p){for(F=p=1,j=0;j<n;F*=j)p*=A[P[2*j]*2*n+P[j-~j++]];S+=p;}}}float h(A,n)float*A;{int P[j=2*n];S=0;s(A,n,P,j,0);S/=F*(1<<n);}

Pruébalo en línea!

Explicación

float S,p,F;                    // global float variables: total sum, temporary, factorial
j,i;                            // global integer variables: indices
s(A,n,P,l,o,k)float*A;int*P;{   // recursively look at every permutation in S_n
 if(k=j=0,o-l)                  // initialize k and j, check if o != l (possible  permutation not yet fully generated)
  for(;k<l;s(A,n,P,l,o+1))      // loop through possible values for current possible  permuation position
   P[o]=k++;                    // set possible  permutation, recursively call (golfed into the for loop)
 else{for(p=-l;j<l;j++)         // there exists a possible permutation fully generated
  for(i=0;i<l;)                 // test if the possible permutation is a bijection
   p+=P[j]==P[i++];             // check for unique elements
  if(!p){                       // indeed, it is a permutation
   for(F=p=1,j=0;j<n;F*=j)      // Hafnian product loop and calculate the factorial (over and over to save bytes)
    p*=A[P[2*j]*2*n+P[j-~j++]]; // Hafnian product
   S+=p;}}}                     // add to sum
float h(A,n)float*A;{           // Hafnian function
 int P[j=2*n];S=0;              // allocate permutation memory, initialize sum
 s(A,n,P,j,0);                  // calculate Hafnian sum
 S/=F*(1<<n);}                  // calculate Hafnian

Pruébalo en línea!

En el núcleo del programa se encuentra el siguiente generador de permutación que se repite S_n. Todo el cálculo hafniano se basa simplemente en él, y se juega más.

j,i,p;Sn(A,l,o,k)int*A;{          // compute every element in S_n
 if(o-l)                          // o!=l, the permutation has not fully been generated
  for(k=0;k<l;k++)                // loop through the integers [0, n)
   A[o]=k,Sn(A,l,o+1);            // modify permutation, call recursively
 else{                            // possible permutation has been generated
  for(p=-l,j=0;j<l;j++)           // look at the entire possible permutation
   for(i=0;i<l;i++)p+=A[j]==A[i]; // check that all elements appear uniquely
  if(!p)                          // no duplicat elements, it is indeed a permutation
   for(printf("["),j=0;j<l        // print
   ||printf("]\n")*0;)            //  the
    printf("%d, ",A[j++]);}}      //   permutation
main(){int l=4,A[l];Sn(A,l,0);}   // all permutations in S_4

Pruébalo en línea!

R , 84 78 bytes

h=function(m)"if"(n<-nrow(m),{for(j in 2:n)F=F+m[1,j]*h(m[v<--c(1,j),v]);F},1)

Pruébalo en línea!

Editar: Gracias a Vlo por -6 bytes.

Parece que todos aquí están implementando el algoritmo de referencia estándar con permutaciones, pero traté de aprovechar el conocimiento de la comunidad adquirido en el desafío relacionado , que es básicamente la misma tarea dirigida al código más rápido en lugar del golf.

Resulta que para un lenguaje que es bueno para cortar matrices (como R), el algoritmo recursivo: hafnian(m) = sum(m[i,j] * hafnian(m[-rows and columns at i,j])no solo es más rápido, sino también bastante golfoso. Aquí está el código sin golf:

hafnian<-function(m)
{
    n=nrow(m)
    #Exits one step earlier than golfed version
    if(n == 2) return(m[1,2])
    h = 0
    for(j in 2:n) {
        if(m[1,j] == 0) next
        h = h + m[1,j] * hafnian(m[c(-1,-j),c(-1,-j)])
    }
    h
}

Jalea , 29 bytes

LHµ2*×!
LŒ!s€2;@€€Wị@/€€P€S÷Ç

Pruébalo en línea!

Creo que la ;@€€Wị@/€€P€parte probablemente se puede reducir. Necesito volver más tarde para verificar y agregar una explicación.

Haskell , 136 bytes

-14 bytes gracias a los ovs.

import Data.List
f m|n<-length m`div`2=sum[product[m!!(s!!(2*j-2)-1)!!(s!!(2*j-1)-1)/realToFrac(2*j)|j<-[1..n]]|s<-permutations[1..n*2]]

Pruébalo en línea!

Ugh ...

MATL , 29 24 22 bytes

Zy:Y@!"G@2eZ{)tn:E/pvs

Pruébalo en línea! O verifique todos los casos de prueba: 1 , 2 , 3 .

Cómo funciona

Zy       % Size of (implicit) input: pushes [2*n 2*n], where the
         % input is a 2*n × 2*n matrix. 
:        % Range: gives row vector [1 2 ... 2*n]
Y@       % All permutation of that vector as rows of a matrix
!"       % For each permutation 
  G      %   Push input matrix
  @      %   Push current permutation
  2e     %   Reshape as a 2-row array
  Z{     %   Split rows into a cell array of size 2
  )      %   Reference indexing. With a cell array as index this
         %   applies element-wise indexing (similar to sub2ind).
         %   Gives a row vector with the n matrix entries selected
         %   by the current permutation
  t      %   Duplicate
  n:     %   Number of elements, range: this gives [1 2 ... n]
  E      %   Double, element-wise: gives [2 4 ... 2*n]
  /      %   Divide, element-wise
  p      %   Product
  vs     %   Vertically concatenate and sum
         % End (implicit). Display (implicit)

Wolfram Language (Mathematica) , 85 bytes

#[[s[[j]],s[[n+j]]]]~Product~{j,n}~Sum~{s,n=Tr[1^#]/2;Permutations@Range[2n]}/n!/2^n&

Pruébalo en línea!

De coco , 165 145 128 127 bytes

-1 byte gracias al Sr. Xcoder

m->sum(reduce((o,j)->o*m[p[2*j]][p[j-~j]]/-~j/2,range(len(m)//2),1)for p in permutations(range(len(m))))
from itertools import*

Pruébalo en línea!

Perl 6, 86 bytes

{my \n=$^m/2;^$m .permutations.map({[*] .map(->\a,\b{$m[a][b]})}).sum/(2**n*[*] 1..n)}