diff --git a/run-tests.py b/run-tests.py index 634b0b6..9513f17 100644 --- a/run-tests.py +++ b/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) diff --git a/src/bithacks.h b/src/bithacks.h new file mode 100644 index 0000000..84d28c4 --- /dev/null +++ b/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; } + diff --git a/src/npy_util.h b/src/npy_util.h new file mode 100644 index 0000000..fff6c30 --- /dev/null +++ b/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; +} + diff --git a/src/old.c b/src/old.c deleted file mode 100644 index da69ce7..0000000 --- a/src/old.c +++ /dev/null @@ -1,90 +0,0 @@ -/* Computes the permanent, given a numpy array */ - -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION -#include -#include -#include -#include - -// 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 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); -} diff --git a/src/oldold.c b/src/oldold.c deleted file mode 100644 index ae162ab..0000000 --- a/src/oldold.c +++ /dev/null @@ -1,68 +0,0 @@ -/* Computes the permanent, given a numpy array */ - -#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION -#include -#include -#include -#include - -// 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); -} diff --git a/src/permanent.c b/src/permanent.c index e97452d..39225c1 100644 --- a/src/permanent.c +++ b/src/permanent.c @@ -2,51 +2,20 @@ #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION #include #include +#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