diff --git a/smop/Makefile b/smop/Makefile
index e1c02161..5c2b8070 100644
--- a/smop/Makefile
+++ b/smop/Makefile
@@ -10,14 +10,9 @@ timer_py: a.py
 clean:
 	rm -f a.* *.pyc
 
-a.so: a.c
-	gcc -I /usr/include/python2.7 -L /usr/lib/python2.7 -O2 -shared -o a.so a.c
-
 a.py: solver.m r8_random.m 
 	python main.py $^
 
-a.c: a.py runtime.py
-	cython a.py
 check:
 	python runtime.py $(FLAGS)
 	python test_runtime.py $(FLAGS)
@@ -26,3 +21,10 @@ check:
 	python test_lexer.py $(FLAGS)
 	python main.py solver.m r8_random.m && python go.py | tail
 	#python main.py fastsolver.m && python go.py
+
+%.c: %.pyx
+	cython -a $^
+
+%.so: %.c
+	gcc -I /usr/include/python2.7 -O2 -shared -o $@ $^
+
diff --git a/smop/go.py b/smop/go.py
index 3e45d12b..11b3c497 100644
--- a/smop/go.py
+++ b/smop/go.py
@@ -1,5 +1,5 @@
 import numpy,time
-import a
+import solver as a
 from runtime import *
 
 def main():
diff --git a/smop/solver.pyx b/smop/solver.pyx
new file mode 100644
index 00000000..a111283d
--- /dev/null
+++ b/smop/solver.pyx
@@ -0,0 +1,55 @@
+from __future__ import division
+import numpy as np
+cimport numpy as np
+from runtime import *
+def solver_(np.ndarray ai,np.ndarray af,int w,int nargout=1):
+    cdef int nBlocks,m,n,i,j,r,bid,ni,nj,ti,tj,d,dn
+    cdef np.ndarray a,I,J,mv
+    rand_(1,2,3)
+    nBlocks=max_(ai[:])
+    m,n=size_(ai,nargout=2)
+    I=matlabarray([0,1,0,- 1])
+    J=matlabarray([1,0,- 1,0])
+    a=copy_(ai)
+    mv=matlabarray([])
+    while not isequal_(af,a):
+
+        bid=ceil_(rand_() * nBlocks)
+        i,j=find_(a == bid,nargout=2)
+        r=ceil_(rand_() * 4)
+        ni=i + I[r]
+        nj=j + J[r]
+        if (ni < 1) or (ni > m) or (nj < 1) or (nj > n):
+            continue
+        if a[ni,nj] > 0:
+            continue
+        ti,tj=find_(af == bid,nargout=2)
+        d=(ti - i) ** 2 + (tj - j) ** 2
+        dn=(ti - ni) ** 2 + (tj - nj) ** 2
+        if (d < dn) and (rand_() > 0.05):
+            continue
+        a[ni,nj]=bid
+        a[i,j]=0
+        mv[mv.shape[0] + 1,[1,2]]=[bid,r]
+
+    return mv
+def rand_(*args,**kwargs):
+    varargin = cellarray(args)
+    nargin = len(args)+0
+
+    global s1,s2,s3
+    if nargin != 0:
+        r=0
+        s1=varargin[1]
+        s2=varargin[2]
+        s3=varargin[3]
+    else:
+        r,s1,s2,s3=r8_random_(s1,s2,s3,nargout=4)
+    return r
+def r8_random_(double s1,double s2,double s3,nargout=1):
+    cdef double r
+    s1=171 * s1 % 30269
+    s2=172 * s2 % 30307
+    s3=170 * s3 % 30323
+    r=(s1 / 30269.0 + s2 / 30307.0 + s3 / 30323.0) % 1.0
+    return r,s1,s2,s3