Adding the real permanent now

src/permanent.c

@@ -6,6 +6,10 @@
 #include <math.h>
 #include <complex.h> 
 
+// Globals 
+PyArrayObject *submatrix;
+int size;
+
 // Boilerplate: Forward function declaration.
 static PyObject *permanent(PyObject *self, PyObject *args); 
 
@@ -27,32 +31,36 @@ PyMODINIT_FUNC initpermanent(void) {
                     (x1) * PyArray_STRIDES(submatrix)[1])))
 #define SM_shape(x0) (int) PyArray_DIM(submatrix, x0)
 
+// 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) {
+  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_complex128 q = SM(0,0);
+      p.real += q.real;
+      p.imag += q.imag;
+    }
+  }
+  return p;
+}
+
+// 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;}
-  int d = (int) PyArray_DIM(submatrix, 0);
-  if ((int) PyArray_DIM(submatrix, 1) != d) {return NULL;}
-
-  // This number ends up being the permanent
-  npy_complex64 p;
-  p.real=0; p.imag=0;
-
-  int i, j;
-  for (i = 0; i < d; ++i) {
-    for (j = 0; j<d; ++j) {
-      npy_complex128 q = SM(0,0);
-      p.real += q.real;
-      p.imag += q.imag;
-    }
-  }
+  size = (int) PyArray_DIM(submatrix, 0);
+  if ((int) PyArray_DIM(submatrix, 1) != size) {return NULL;}
 
-  // Convert to a python complex number and return
-  PyObject *output=PyComplex_FromDoubles(p.real,p.imag);
-  return output;
+  // Get the permanent, convert to a python object, and return
+  npy_complex64 p = perm_ryser();
+  return PyComplex_FromDoubles(p.real,p.imag);
 }