diff --git a/makefile b/makefile index 634105c..207db8f 100644 --- a/makefile +++ b/makefile @@ -7,7 +7,7 @@ setup = setup.py default : $(target) test # Compile -$(target): $(srcdir)/permanent.c +$(target): $(srcdir)/* python $(setup) build_ext --inplace test : diff --git a/run-tests.py b/run-tests.py index 3eef926..4099c8d 100644 --- a/run-tests.py +++ b/run-tests.py @@ -4,6 +4,7 @@ import multiprocessing as mp import numpy as np import lib import time +from matplotlib import pyplot as plt def perm_ryser(a): ''' the permanent calculated using the ryser formula. much faster than the naive approach ''' @@ -16,23 +17,30 @@ def perm_ryser(a): terms=map(get_term, indeces) return np.sum(terms)*((-1)**n) -dimension=9 -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 +maxtime=1 +dimensions=range(1,11) -#print lib.permanent(submatrix), perm_ryser(submatrix) +for (function, label) in zip((lib.permanent, perm_ryser), ("C", "Python")): + counts=[] + for dimension in dimensions: + print dimension + 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 -#sys.exit(0) -#t=time.clock() -#for i in range(1000): - #perm_ryser(submatrix) -#t1=time.clock()-t + t=time.clock() + n=0 + while time.clock()-t < maxtime: + for i in range(5): + function(submatrix) + n+=5 + counts.append(n) -t=time.clock() -for i in range(1000): - lib.permanent(submatrix) -t2=time.clock()-t -print t2 + plt.plot(dimensions, counts, '.-', label=label) -#print t1/t2 +plt.ylabel('Number of permanents per second') +plt.xlabel('Dimension') +plt.xlim(min(dimensions), max(dimensions)) +plt.legend() +plt.semilogy() +plt.savefig('out.pdf') diff --git a/src/bithacks.h b/src/bithacks.h index 84d28c4..ee87c03 100644 --- a/src/bithacks.h +++ b/src/bithacks.h @@ -8,7 +8,7 @@ int countbits(unsigned int n) 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. + q = (q & 0x00000000FFFFFFFF) + ((q & 0xFFFFFFFF00000000) >> 32); // This last & isn't strictly qecessary. return q; } diff --git a/src/npy_util.h b/src/npy_util.h index fff6c30..4857cc5 100644 --- a/src/npy_util.h +++ b/src/npy_util.h @@ -40,8 +40,8 @@ void complex_inc(npy_complex128 *a, npy_complex128 b) { // 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; + double c = a->real; + a->real = a->real*b.real-a->imag*b.imag; + a->imag = c*b.imag+a->imag*b.real; } diff --git a/src/permanent.c b/src/permanent.c index 7c40319..782bce1 100644 --- a/src/permanent.c +++ b/src/permanent.c @@ -4,7 +4,7 @@ #include #include "npy_util.h" #include "bithacks.h" -#include "permanents.h" +#include "ryser.h" // Forward function declaration static PyObject *permanent(PyObject *self, PyObject *args); @@ -28,6 +28,6 @@ static PyObject *permanent(PyObject *self, PyObject *args) { if (!PyArg_ParseTuple(args, "O!", &PyArray_Type, &submatrix)) {return NULL;} // Compute the permanent - npy_complex128 p = perm_ryser(submatrix); + npy_complex128 p = ryser(submatrix); return PyComplex_FromDoubles(p.real, p.imag); } diff --git a/src/permanents.h b/src/permanents.h deleted file mode 100644 index 153ce14..0000000 --- a/src/permanents.h +++ /dev/null @@ -1,22 +0,0 @@ -// Ryser's algorithm -static npy_complex128 perm_ryser(PyArrayObject *submatrix) { - int n = (int) PyArray_DIM(submatrix, 0); - npy_complex128 rowsum, rowsumprod; - npy_complex128 perm = complex_zero; - int exp = 1 << n; - int i, y, z; - for (i=0; i