diff -r 303c51137cb2 -r fea8e9ae6f29 geometry.py
--- a/geometry.py	Thu Dec 21 12:27:16 2017 +0200
+++ b/geometry.py	Fri Dec 22 15:13:43 2017 +0200
@@ -205,28 +205,107 @@
 def is_real(number):
     return isinstance(number, int) or isinstance(number, float)
 
+def matrix_indices():
+    '''
+        Yields index pairs for matrix iteration
+    '''
+    from itertools import product
+    yield from product(range(3), range(3))
+
 class Matrix3x3:
     '''
         A 3×3 matrix forming the top-left portion of a full 4×4 transformation
         matrix.
     '''
-    def __init__(self, values):
-        assert(all(is_real(x) for x in values))
-        assert len(values) == 9
-        self.values = values
+    class Row:
+        def __init__(self, matrix_ref, i1, i2, i3):
+            self.matrix_ref = matrix_ref
+            self.matrix_indices = i1, i2, i3
+        def __getitem__(self, row_index):
+            return self.matrix_ref[self.matrix_indices[row_index]]
+        def __setitem__(self, row_index, value):
+            self.matrix_ref[self.matrix_indices[row_index]] = value
+    def __init__(self, values = None):
+        if values is not None:
+            assert(all(is_real(x) for x in values))
+            assert len(values) == 9
+            self.values = values
+        else:
+            self.values = [1, 0, 0, 0, 1, 0, 0, 0, 1]
     def __repr__(self):
-        return str.format('Matrix3x3({!r})', self.values)
+        if self != Matrix3x3():
+            return str.format('Matrix3x3({!r})', self.values)
+        else:
+            return 'Matrix3x3()'
     def __getitem__(self, index):
-        return self.values[index]
+        if isinstance(index, tuple):
+            assert all(k in [0, 1, 2] for k in index)
+            return self.values[index[0] * 3 + index[1]]
+        else:
+            return Matrix3x3.Row(self, index * 3, index * 3 + 1, index * 3 + 2)
     def determinant(self):
-        v = self.values
         return sum([
-            +(v[0] * v[4] * v[8]),
-            +(v[1] * v[5] * v[6]),
-            +(v[2] * v[3] * v[7]),
-            -(v[2] * v[4] * v[6]),
-            -(v[1] * v[3] * v[8]),
-            -(v[0] * v[5] * v[7]),
+            +(self[0, 0] * self[1, 1] * self[2, 2]),
+            +(self[0, 1] * self[1, 2] * self[2, 0]),
+            +(self[0, 2] * self[1, 0] * self[2, 1]),
+            -(self[0, 2] * self[1, 1] * self[2, 0]),
+            -(self[0, 1] * self[1, 0] * self[2, 2]),
+            -(self[0, 0] * self[1, 2] * self[2, 1]),
+        ])
+    @property
+    def scaling_vector(self):
+        '''
+            Extracts scaling factors from this transformation matrix.
+        '''
+        from math import sqrt
+        return Vector(
+            x = sqrt(self[0, 0] ** 2 + self[1, 0] ** 2 + self[2, 0] ** 2),
+            y = sqrt(self[0, 1] ** 2 + self[1, 1] ** 2 + self[2, 1] ** 2),
+            z = sqrt(self[0, 2] ** 2 + self[1, 2] ** 2 + self[2, 2] ** 2),
+        )
+    @property
+    def rotation_component(self):
+        '''
+            Extracts rotation from this matrix.
+        '''
+        return self.split()[0]
+    def split(self):
+        '''
+            Extracts the rotation matrix and scaling vector.
+        '''
+        vec = self.scaling_vector
+        return Matrix3x3([
+            self[i, j] / vec.coordinates[j]
+            for i, j in matrix_indices()
+        ]), vec
+    @property
+    def contains_scaling(self):
+        '''
+            Returns whether this matrix contains scaling factors.
+        '''
+        vec = self.scaling_vector
+        return abs((vec.x * vec.y * vec.z) - 1) >= 1e-10
+    @property
+    def contains_rotation(self):
+        '''
+            Returns whether this matrix contains rotation factors.
+        '''
+        return self.rotation_component != Matrix3x3()
+    def __eq__(self, other):
+        '''
+            Returns whether two matrices are equivalent.
+        '''
+        return all(
+            abs(self[(i, j)] - other[(i, j)]) < 1e-10
+            for i, j in matrix_indices()
+        )
+    def __matmul__(self, other):
+        '''
+            Computes the matrix multiplication self × other.
+        '''
+        return Matrix3x3([
+            sum(self[i, k] * other[k, j] for k in range(3))
+            for i, j in matrix_indices()
         ])
 
 def complete_matrix(matrix, anchor):
@@ -235,9 +314,9 @@
         transformation matrix.
     '''
     return [
-        matrix[0], matrix[1], matrix[2], anchor.x,
-        matrix[3], matrix[4], matrix[5], anchor.y,
-        matrix[6], matrix[7], matrix[8], anchor.z,
+        matrix[0, 0], matrix[0, 1], matrix[0, 2], anchor.x,
+        matrix[1, 0], matrix[1, 1], matrix[1, 2], anchor.y,
+        matrix[2, 0], matrix[2, 1], matrix[2, 2], anchor.z,
         0, 0, 0, 1,
     ]
 
@@ -245,16 +324,17 @@
     '''
         Transforms a vertex by a 4×4 transformation matrix.
     '''
-    u = transformation_matrix[0] * vertex.x \
-        + transformation_matrix[1] * vertex.y \
-        + transformation_matrix[2] * vertex.z \
-        + transformation_matrix[3]
-    v = transformation_matrix[4] * vertex.x \
-        + transformation_matrix[5] * vertex.y \
-        + transformation_matrix[6] * vertex.z \
-        + transformation_matrix[7]
-    w = transformation_matrix[8] * vertex.x \
-        + transformation_matrix[9] * vertex.y \
-        + transformation_matrix[10] * vertex.z \
-        + transformation_matrix[11]
-    return Vertex(u, v, w)
+    return Vertex(
+        x = transformation_matrix[0] * vertex.x \
+            + transformation_matrix[1] * vertex.y \
+            + transformation_matrix[2] * vertex.z \
+            + transformation_matrix[3],
+        y = transformation_matrix[4] * vertex.x \
+            + transformation_matrix[5] * vertex.y \
+            + transformation_matrix[6] * vertex.z \
+            + transformation_matrix[7],
+        z = transformation_matrix[8] * vertex.x \
+            + transformation_matrix[9] * vertex.y \
+            + transformation_matrix[10] * vertex.z \
+            + transformation_matrix[11],
+    )