--- a/pytouhou/utils/matrix.pyx
+++ b/pytouhou/utils/matrix.pyx
@@ -17,120 +17,123 @@ from libc.stdlib cimport malloc, free
 from libc.string cimport memcpy
 
 
-cdef float[16] identity
-identity[:] = [1, 0, 0, 0,
-               0, 1, 0, 0,
-               0, 0, 1, 0,
-               0, 0, 0, 1]
+cdef Matrix identity
+identity = Matrix(1, 0, 0, 0,
+                  0, 1, 0, 0,
+                  0, 0, 1, 0,
+                  0, 0, 0, 1)
 
 
-cdef class Matrix:
-    def __init__(self, data=None):
-        if data is not None:
-            for i in xrange(16):
-                self.data[i] = data[i]
-        else:
-            memcpy(self.data, identity, 16 * sizeof(float))
+cdef Matrix *new_matrix(Matrix *data) nogil:
+    mat = <Matrix*> malloc(sizeof(Matrix))
+    memcpy(mat, data, sizeof(Matrix))
+    return mat
+
+
+cdef Matrix *new_identity() nogil:
+    return new_matrix(&identity)
 
 
-    def __mul__(Matrix self, Matrix other):
-        cdef float *d1, *d2, *d3
+cdef void mul(Matrix *mat1, Matrix *mat2) nogil:
+    cdef float *d3
 
-        out = Matrix()
-        d1 = self.data
-        d2 = other.data
-        d3 = out.data
-        for i in xrange(4):
-            for j in xrange(4):
-                d3[4*i+j] = 0
-                for k in xrange(4):
-                    d3[4*i+j] += d1[4*i+k] * d2[4*k+j]
-        return out
+    out = <Matrix*> malloc(sizeof(Matrix))
+    d1 = <float*>mat1
+    d2 = <float*>mat2
+    d3 = <float*>out
+    for i in xrange(4):
+        for j in xrange(4):
+            d3[4*i+j] = 0
+            for k in xrange(4):
+                d3[4*i+j] += d1[4*i+k] * d2[4*k+j]
+    memcpy(mat1, out, sizeof(Matrix))
+    free(out)
 
 
-    cdef void flip(self) nogil:
-        cdef float *data
-
-        data = self.data
-        for i in xrange(4):
-            data[i] = -data[i]
+cdef void flip(Matrix *mat) nogil:
+    data = <float*>mat
+    for i in xrange(4):
+        data[i] = -data[i]
 
 
-    cdef void scale(self, float x, float y, float z) nogil:
-        cdef float *data, coordinate[3]
+cdef void scale(Matrix *mat, float x, float y, float z) nogil:
+    cdef float coordinate[3]
+
+    data = <float*>mat
+    coordinate[0] = x
+    coordinate[1] = y
+    coordinate[2] = z
 
-        data = self.data
-        coordinate[0] = x
-        coordinate[1] = y
-        coordinate[2] = z
+    for i in xrange(3):
+        for j in xrange(4):
+            data[4*i+j] *= coordinate[i]
+
 
-        for i in xrange(3):
-            for j in xrange(4):
-                data[4*i+j] *= coordinate[i]
+cdef void scale2d(Matrix *mat, float x, float y) nogil:
+    data = <float*>mat
+    for i in xrange(4):
+        data[  i] *= x
+        data[4+i] *= y
 
 
-    cdef void scale2d(self, float x, float y) nogil:
-        cdef float *data
+cdef void translate(Matrix *mat, float x, float y, float z) nogil:
+    cdef float offset[3], item[3]
+
+    offset[0] = x
+    offset[1] = y
+    offset[2] = z
 
-        data = self.data
-        for i in xrange(4):
-            data[  i] *= x
-            data[4+i] *= y
+    data = <float*>mat
+    for i in xrange(3):
+        item[i] = data[12+i] * offset[i]
+
+    for i in xrange(3):
+        for j in xrange(4):
+            data[4*i+j] += item[i]
 
 
-    cdef void translate(self, float x, float y, float z) nogil:
-        cdef float *data, coordinate[3], item[3]
-
-        data = self.data
-        coordinate[0] = x
-        coordinate[1] = y
-        coordinate[2] = z
-        for i in xrange(3):
-            item[i] = data[12+i] * coordinate[i]
-
-        for i in xrange(3):
-            for j in xrange(4):
-                data[4*i+j] += item[i]
+cdef void translate2d(Matrix *mat, float x, float y) nogil:
+    translate(mat, x, y, 0)
 
 
-    cdef void rotate_x(self, float angle) nogil:
-        cdef float cos_a, sin_a
-        cdef float lines[8], *data
+cdef void rotate_x(Matrix *mat, float angle) nogil:
+    cdef float cos_a, sin_a
+    cdef float lines[8]
 
-        data = self.data
-        cos_a = cos(angle)
-        sin_a = sin(angle)
-        for i in xrange(8):
-            lines[i] = data[i+4]
-        for i in xrange(4):
-            data[4+i] = cos_a * lines[i] - sin_a * lines[4+i]
-            data[8+i] = sin_a * lines[i] + cos_a * lines[4+i]
+    data = <float*>mat
+    cos_a = cos(angle)
+    sin_a = sin(angle)
+    for i in xrange(8):
+        lines[i] = data[i+4]
+    for i in xrange(4):
+        data[4+i] = cos_a * lines[i] - sin_a * lines[4+i]
+        data[8+i] = sin_a * lines[i] + cos_a * lines[4+i]
 
 
-    cdef void rotate_y(self, float angle) nogil:
-        cdef float cos_a, sin_a
-        cdef float lines[8], *data
+cdef void rotate_y(Matrix *mat, float angle) nogil:
+    cdef float cos_a, sin_a
+    cdef float lines[8]
 
-        data = self.data
-        cos_a = cos(angle)
-        sin_a = sin(angle)
-        for i in xrange(4):
-            lines[i] = data[i]
-            lines[i+4] = data[i+8]
-        for i in xrange(4):
-            data[  i] =  cos_a * lines[i] + sin_a * lines[4+i]
-            data[8+i] = -sin_a * lines[i] + cos_a * lines[4+i]
+    data = <float*>mat
+    cos_a = cos(angle)
+    sin_a = sin(angle)
+    for i in xrange(4):
+        lines[i] = data[i]
+        lines[i+4] = data[i+8]
+    for i in xrange(4):
+        data[  i] =  cos_a * lines[i] + sin_a * lines[4+i]
+        data[8+i] = -sin_a * lines[i] + cos_a * lines[4+i]
 
 
-    cdef void rotate_z(self, float angle) nogil:
-        cdef float cos_a, sin_a
-        cdef float lines[8], *data
+cdef void rotate_z(Matrix *mat, float angle) nogil:
+    cdef float cos_a, sin_a
+    cdef float lines[8]
 
-        data = self.data
-        cos_a = cos(angle)
-        sin_a = sin(angle)
-        for i in xrange(8):
-            lines[i] = data[i]
-        for i in xrange(4):
-            data[  i] = cos_a * lines[i] - sin_a * lines[4+i]
-            data[4+i] = sin_a * lines[i] + cos_a * lines[4+i]
+    data = <float*>mat
+    cos_a = cos(angle)
+    sin_a = sin(angle)
+    for i in xrange(8):
+        lines[i] = data[i]
+    for i in xrange(4):
+        data[  i] = cos_a * lines[i] - sin_a * lines[4+i]
+        data[4+i] = sin_a * lines[i] + cos_a * lines[4+i]