Still not quite there...

run-tests.py

@@ -3,10 +3,12 @@ import time
 import multiprocessing as mp
 import numpy as np
 import lib
+import time
 
-dimension=8
-real=np.random.uniform(0,1,(dimension, dimension))
-imag=np.random.uniform(0,1,(dimension, dimension))
+dimension=2
+real=np.ones((dimension, dimension))
+imag=np.ones((dimension, dimension))
 submatrix=real+1j*imag
+print submatrix
 
-print "OUTPUT:", lib.permanent(submatrix)
+p=lib.permanent(submatrix)

src/permanent.c

@@ -7,8 +7,6 @@
 #include <complex.h> 
 
 // Globals 
-PyArrayObject *submatrix;
-int size;
 
 // Boilerplate: Forward function declaration.
 static PyObject *permanent(PyObject *self, PyObject *args); 
@@ -34,13 +32,16 @@ PyMODINIT_FUNC initpermanent(void) {
 // Ryser's algorithm takes exponential time
 // The advantage over naive perm() only kicks in at around 6x6 matrices 
 // TODO: should take matrix as arg really, get rid of consts
-static npy_complex64 perm_ryser(void) {
+static npy_complex64 perm_ryser(PyArrayObject *submatrix) {
+  int size = (int) PyArray_DIM(submatrix, 0);
   npy_complex64 p;
   p.real=0; p.imag=0; 
   int i, j;
   for (i = 0; i < size; ++i) {
     for (j = 0; j < size; ++j) {
-      npy_complex64 q = SM(0,0);
+      npy_complex64 q = SM(i,j);
+      printf("real: %f\n", q.real);
+      printf("imag: %f\n", q.imag);
       p.real += q.real;
       p.imag += q.imag;
     }
@@ -48,43 +49,20 @@ static npy_complex64 perm_ryser(void) {
   return p;
 }
 
-void TODO(void) {
-    int n = size;
-    int i = 0; int z = 0; int y = 0;
-    npy_complex64 perm; perm.real=0; perm.imag=0; 
-    npy_complex64 prod; 
-    npy_complex64 sum; sum.real=0; sum.imag=0; 
-    npy_complex64 element;
-    int exp=pow(2.0, n);
-
-    // Iterate over exponentially many index strings
-    for (i=0; i<exp; ++i) {
-        prod.real = 1; prod.imag=0;
-        for (y=0; y<n; ++y) {               // Rows
-            sum.real = 0; sum.imag = 0;
-            for (z=0; z<n; ++z) {           // Columns
-                if ((i & (1 << z)) > 0) { sum.real += SM(z,y).real;  sum.imag += SM(z,y).imag; }
-            }
-            prod = prod*sum;
-        }
-        perm += parity(i) * prod;
-    }
-    return (pow(-1,n))*perm;
-}
-
 // This is basically a wrapper which chooses the optimal permanent function
 static PyObject *permanent(PyObject *self, PyObject *args) {
   // Parse input
+  PyArrayObject *submatrix;
   if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &submatrix)) {
     return NULL;
   }
 
   // Check for stupid mistakes
   if ((int) PyArray_NDIM(submatrix) != 2) {return NULL;}
-  size = (int) PyArray_DIM(submatrix, 0);
+  int size = (int) PyArray_DIM(submatrix, 0);
   if ((int) PyArray_DIM(submatrix, 1) != size) {return NULL;}
 
   // Get the permanent, convert to a python object, and return
-  npy_complex64 p = perm_ryser();
-  return PyComplex_FromDoubles(p.real,p.imag);
+  npy_complex64 p = perm_ryser(submatrix);
+  return PyComplex_FromDoubles(p.real, p.imag);
 }