Had a tidy up, all works well. TODO: hard-coded versions

run-tests.py

@@ -1,4 +1,4 @@
-import os
+import os, sys
 import time
 import multiprocessing as mp
 import numpy as np
@@ -16,26 +16,14 @@ def perm_ryser(a):
     terms=map(get_term, indeces) 
     return np.sum(terms)*((-1)**n)
 
-def explain_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 
-    for q in irange:
-        print get_index(q)
-    #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=5
 real=np.random.uniform(-1, 1, dimension*dimension).reshape((dimension, dimension))
 imag=np.random.uniform(-1, 1, dimension*dimension).reshape((dimension, dimension))
 submatrix=real+1j*imag
 
+print lib.permanent(submatrix), perm_ryser(submatrix)
+
+sys.exit(0)
 t=time.clock()
 for i in range(1000):
     perm_ryser(submatrix)

src/bithacks.h

@@ -0,0 +1,16 @@
+
+// Count the number of set bits in a binary string
+int countbits(unsigned int n) 
+{
+    int q=n;
+    q = (q & 0x5555555555555555) + ((q & 0xAAAAAAAAAAAAAAAA) >> 1);
+    q = (q & 0x3333333333333333) + ((q & 0xCCCCCCCCCCCCCCCC) >> 2);
+    q = (q & 0x0F0F0F0F0F0F0F0F) + ((q & 0xF0F0F0F0F0F0F0F0) >> 4);
+    q = (q & 0x00FF00FF00FF00FF) + ((q & 0xFF00FF00FF00FF00) >> 8);
+    q = (q & 0x0000FFFF0000FFFF) + ((q & 0xFFFF0000FFFF0000) >> 16);
+    q = (q & 0x00000000FFFFFFFF) + ((q & 0xFFFFFFFF00000000) >> 32); // This last & isq't strictly qecessary.
+    return q;
+}
+
+int bitparity (unsigned int n) { return 1 - (countbits(n) & 1)*2; }
+

src/npy_util.h

@@ -0,0 +1,47 @@
+// Array access macros.
+#define SM(x0, x1) (*(npy_complex128*)((PyArray_DATA(submatrix) + \
+                    (x0) * PyArray_STRIDES(submatrix)[0] +  \
+                    (x1) * PyArray_STRIDES(submatrix)[1])))
+#define SM_shape(x0) (int) PyArray_DIM(submatrix, x0)
+
+// Complex numbers
+static const npy_complex128 complex_one = {.real=1, .imag=0};
+static const npy_complex128 complex_zero = {.real=0, .imag=0};
+
+// Add two numbers
+npy_complex128 complex_add(npy_complex128 a, npy_complex128 b) { 
+    npy_complex128 x;
+    x.real = a.real+b.real;
+    x.imag = a.imag+b.real;
+    return x;
+}
+
+// Product of two numbers
+npy_complex128 complex_prod(npy_complex128 a, npy_complex128 b) { 
+    npy_complex128 x;
+    x.real = a.real*b.real - a.imag*b.imag;
+    x.imag = a.imag*b.real + a.real*b.imag;
+    return x;
+}
+
+// Product of complex and float
+npy_complex128 complex_float_prod(npy_complex128 a, float b) { 
+    npy_complex128 x;
+    x.real = a.real*b;
+    x.imag = a.imag*b;
+    return x;
+}
+
+// Increment a number
+void complex_inc(npy_complex128 *a, npy_complex128 b) { 
+    a->real += b.real;
+    a->imag += b.imag;
+}
+
+// Multipy a number by another one
+void complex_multiply(npy_complex128 *a, npy_complex128 b) { 
+    npy_complex128 c = {.real=a->real, .imag=a->imag};
+    a->real  = c.real*b.real-c.imag*b.imag;
+    a->imag  = c.real*b.imag+c.imag*b.real;
+}
+

src/old.c

@@ -1,90 +0,0 @@
-/* Computes the permanent, given a numpy array */
-
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
-#include <Python.h>
-#include <numpy/arrayobject.h>
-#include <math.h>
-#include <complex.h> 
-
-// Globals 
-PyArrayObject *submatrix;
-int size;
-
-// Boilerplate: Forward function declaration.
-static PyObject *permanent(PyObject *self, PyObject *args); 
-
-// Boilerplate: method list.
-static PyMethodDef methods[] = {
-  { "permanent", permanent, METH_VARARGS, "Compute the permanent"},
-  { NULL, NULL, 0, NULL } /* Sentinel */
-};
-
-// Boilerplate: Module initialization.
-PyMODINIT_FUNC initpermanent(void) {
-  (void) Py_InitModule("permanent", methods);
-  import_array();
-}
-
-// Array access macros.
-#define SM(x0, x1) (*(npy_complex64*)((PyArray_DATA(submatrix) + \
-                    (x0) * PyArray_STRIDES(submatrix)[0] +  \
-                    (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_complex64 q = SM(0,0);
-      p.real += q.real;
-      p.imag += q.imag;
-    }
-  }
-  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 = c_prod(prod, sum);
-        }
-        perm += parity(i) * prod;
-    }
-    return c_prod((pow(-1,n)), perm);
-}
-
-// This is basically a wrapper which chooses the optimal permanent function
-static PyObject *permanent(PyObject *self, PyObject *args) {
-  // Parse input
-  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);
-  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);
-}

src/oldold.c

@@ -1,68 +0,0 @@
-/* Computes the permanent, given a numpy array */
-
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
-#include <Python.h>
-#include <numpy/arrayobject.h>
-#include <math.h>
-#include <complex.h> 
-
-// Globals 
-
-// Boilerplate: Forward function declaration.
-static PyObject *permanent(PyObject *self, PyObject *args); 
-
-// Boilerplate: method list.
-static PyMethodDef methods[] = {
-  { "permanent", permanent, METH_VARARGS, "Compute the permanent"},
-  { NULL, NULL, 0, NULL } /* Sentinel */
-};
-
-// Boilerplate: Module initialization.
-PyMODINIT_FUNC initpermanent(void) {
-  (void) Py_InitModule("permanent", methods);
-  import_array();
-}
-
-// Array access macros.
-#define SM(x0, x1) (*(npy_complex64*)((PyArray_DATA(submatrix) + \
-                    (x0) * PyArray_STRIDES(submatrix)[0] +  \
-                    (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(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(i,j);
-      printf("real: %f\n", q.real);
-      printf("imag: %f\n", q.imag);
-      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 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(submatrix);
-  return PyComplex_FromDoubles(p.real, p.imag);
-}

src/permanent.c

@@ -2,51 +2,20 @@
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #include <Python.h>
 #include <numpy/arrayobject.h>
+#include "bithacks.h"
+#include "npy_util.h"
 
-// Array access macros.
-#define SM(x0, x1) (*(npy_complex128*)((PyArray_DATA(submatrix) + \
-                    (x0) * PyArray_STRIDES(submatrix)[0] +  \
-                    (x1) * PyArray_STRIDES(submatrix)[1])))
-#define SM_shape(x0) (int) PyArray_DIM(submatrix, x0)
+// Forward function declaration 
+static PyObject *permanent(PyObject *self, PyObject *args);    
 
-int countbits(unsigned int n) 
-{
-    int q=n;
-    q = (q & 0x5555555555555555) + ((q & 0xAAAAAAAAAAAAAAAA) >> 1);
-    q = (q & 0x3333333333333333) + ((q & 0xCCCCCCCCCCCCCCCC) >> 2);
-    q = (q & 0x0F0F0F0F0F0F0F0F) + ((q & 0xF0F0F0F0F0F0F0F0) >> 4);
-    q = (q & 0x00FF00FF00FF00FF) + ((q & 0xFF00FF00FF00FF00) >> 8);
-    q = (q & 0x0000FFFF0000FFFF) + ((q & 0xFFFF0000FFFF0000) >> 16);
-    q = (q & 0x00000000FFFFFFFF) + ((q & 0xFFFFFFFF00000000) >> 32); // This last & isq't strictly qecessary.
-    return q;
-}
-
-int bitparity (unsigned int n) { return 1 - (countbits(n) & 1)*2; }
-
-
-// Complex numbers
-static const npy_complex128 complex_one = {.real=1, .imag=0};
-static const npy_complex128 complex_zero = {.real=0, .imag=0};
-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;
-    x.imag = a.imag*b.real + a.real*b.imag;
-    return x;
-}
-
-// Boilerplate
-static PyObject *permanent(PyObject *self, PyObject *args);     // Forward function declaration
-static PyMethodDef methods[] = {                                // Method list
+// Method list
+static PyMethodDef methods[] = {                                
   { "permanent", permanent, METH_VARARGS, "Compute the permanent"},
-  { NULL, NULL, 0, NULL } /* Sentinel */
+  { NULL, NULL, 0, NULL } // Sentinel
 };
-PyMODINIT_FUNC initpermanent(void) {                            // Module initialization
+
+// Module initialization
+PyMODINIT_FUNC initpermanent(void) {                            
   (void) Py_InitModule("permanent", methods);
   import_array();
 }
@@ -58,21 +27,18 @@ static npy_complex128 perm_ryser(PyArrayObject *submatrix) {
     npy_complex128 perm = complex_zero;
     int exp = 1 << n; 
     int i, y, z;
-
-    // Iterate over exponentially many index strings
     for (i=0; i<exp; ++i) {
         rowsumprod = complex_one;
-        for (y=0; y<n; ++y) {               // Rows
+        for (y=0; y<n; ++y) {               
             rowsum = complex_zero;
-            for (z=0; z<n; ++z) {           // Columns
-                if ((i & (1 << z)) != 0) { rowsum = complex_add(rowsum, SM(z,y)); }
+            for (z=0; z<n; ++z) { 
+                if ((i & (1 << z)) != 0) { complex_inc(&rowsum, SM(z, y)); }
             }
-            rowsumprod = complex_prod(rowsumprod, rowsum);
+            complex_multiply(&rowsumprod, rowsum);
         }
-        int sign = bitparity(i); 
-        perm.real+=sign*rowsumprod.real; perm.imag+=sign*rowsumprod.imag;
+        complex_inc(&perm, complex_float_prod(rowsumprod, bitparity(i)));
     }
-    if (n%2 == 1) {perm.real*=-1; perm.imag*=-1;}
+    if (n%2 == 1) {perm=complex_float_prod(perm, -1);}
     return perm;
 }