Still broken. Added a showdown

run-tests.py

@@ -5,12 +5,32 @@ import numpy as np
 import lib
 import time
 
-dimension=2
+def perm_ryser(a):
+    ''' the permanent calculated using the ryser formula. much faster than the naive approach '''
+    n,n2=a.shape
+    z=np.arange(n)
+    irange=xrange(2**n)
+    get_index=lambda i: (i & (1 << z)) != 0 
+    get_term=lambda index: ((-1)**np.sum(index))*np.prod(np.sum(a[index,:], 0))
+    indeces=map(get_index, irange)
+    terms=map(get_term, indeces) 
+    return np.sum(terms)*((-1)**n)
+
+dimension=8
 real=np.ones((dimension, dimension))
 imag=np.ones((dimension, dimension))*0
 submatrix=real+1j*imag
 submatrix[0,0]*=2
 print submatrix
 
-p=lib.permanent(submatrix)
-print p
+t=time.clock()
+for i in range(10000):
+    lib.permanent(submatrix)
+print time.clock()-t
+
+#t=time.clock()
+#for i in range(4000):
+    #perm_ryser(submatrix)
+#print time.clock()-t
+
+

src/permanent.c

@@ -12,7 +12,12 @@
 // Complex numbers
 static const npy_complex128 complex_one = {.real=1, .imag=0};
 static const npy_complex128 complex_zero = {.real=0, .imag=0};
-static void complex_inc(npy_complex128 a, npy_complex128 b) { a.real+=b.real; a.imag+=b.imag;}
+static npy_complex128 complex_add(npy_complex128 a, npy_complex128 b) { 
+    npy_complex128 x;
+    x.real = a.real+b.real;
+    x.imag = a.imag+b.imag;
+    return x;
+}
 static npy_complex128 complex_prod(npy_complex128 a, npy_complex128 b) { 
     npy_complex128 x;
     x.real = a.real*b.real+a.imag*b.imag;
@@ -31,6 +36,7 @@ PyMODINIT_FUNC initpermanent(void) {                            // Module initia
   import_array();
 }
 
+
 // Ryser's algorithm 
 static npy_complex128 perm_ryser(PyArrayObject *submatrix) {
     int n = (int) PyArray_DIM(submatrix, 0);
@@ -39,27 +45,23 @@ static npy_complex128 perm_ryser(PyArrayObject *submatrix) {
     npy_complex128 perm = complex_zero;
     npy_complex128 sum = complex_zero;
     int exp = 1 << n; 
-    printf("2^%d = %d\n", n, exp);
 
     // Iterate over exponentially many index strings
     for (i=0; i<exp; ++i) {
         prod = complex_one;
-        printf("prod %.3f %.3f \n", prod.real, prod.imag);
         for (y=0; y<n; ++y) {               // Rows
             sum = complex_zero;
             for (z=0; z<n; ++z) {           // Columns
-                printf("element %.3f %.3f \n", SM(z,y).real, SM(z,y).imag);
-                // Below is a suspicious line
-                if ((i & (1 << z)) > 0) { complex_inc(sum, SM(z,y)); }
-                printf("sum %.3f %.3f \n", sum.real, sum.imag);
+                if ((i & (1 << z)) > 0) { 
+                    sum = complex_add(sum, SM(z,y)); 
+                }
             }
             prod = complex_prod(prod, sum);
         }
-        printf("prod %.3f %.3f \n", prod.real, prod.imag);
-        if (i%2 == 0) {prod.real*=-1; prod.imag*=-1;}
-        complex_inc(perm, prod);
+        /*if (i%2 == 0) {prod.real*=-1; prod.imag*=-1;}*/
+        perm = complex_add(perm, prod);
     }
-    if (i%2 == 0) {perm.real*=-1; perm.imag*=-1;}
+    /*if (i%2 == 0) {perm.real*=-1; perm.imag*=-1;}*/
     return perm;
 }