diff --git a/pyRCV2/constraints.py b/pyRCV2/constraints.py
new file mode 100644
index 0000000..3c58823
--- /dev/null
+++ b/pyRCV2/constraints.py
@@ -0,0 +1,144 @@
+#   pyRCV2: Preferential vote counting
+#   Copyright © 2020–2021  Lee Yingtong Li (RunasSudo)
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU Affero General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU Affero General Public License for more details.
+#
+#   You should have received a copy of the GNU Affero General Public License
+#   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+from pyRCV2.method.base_stv import STVException
+
+class ConstraintMatrix:
+	def __init__(self, dimensions):
+		self.dimensions = dimensions
+		
+		num_cells = 1
+		for d in dimensions:
+			num_cells *= (d + 1)
+		
+		# FIXME: This is not the most efficient packing, as we require margins only along one dimension at a time
+		self.matrix = [ConstraintMatrixCell() for _ in range(num_cells)]
+	
+	# Pass -1 to get margins
+	def index_of(self, address):
+		index = 0
+		for i, d in enumerate(self.dimensions):
+			index *= (d + 1)
+			index += (address[i] + 1)
+		return index
+	
+	def get(self, address):
+		return self.matrix[self.index_of(address)]
+
+class ConstraintMatrixCell:
+	def __init__(self):
+		self.elected = 0
+		self.min = 0
+		self.max = 0
+		self.cands = 0
+	
+	def __repr__(self):
+		return '<ConstraintMatrixCell: E={}, Max={}, Min={}, C={}>'.format(self.elected, self.min, self.max, self.cands)
+
+def ndrange(dimensions):
+	# n-dimensional range function
+	if len(dimensions) == 0:
+		yield []
+		return
+	
+	if len(dimensions) == 1:
+		yield from ([n] for n in range(dimensions[0]))
+		return
+	
+	yield from ([n] + j for n in range(dimensions[0]) for j in ndrange(dimensions[1:]))
+
+def init_matrix(counter):
+	counter._constraint_matrix = ConstraintMatrix([len(category) for category_name, category in counter.election.constraints.items()])
+
+def update_matrix(counter):
+	# Update matrix with elected numbers
+	...
+
+def step_matrix(counter):
+	cm = counter._constraint_matrix
+	
+	# Rule 1
+	for cell in cm.matrix:
+		if cell.elected > cell.min:
+			cell.min = cell.elected
+			return False
+		if cell.cands < cell.max:
+			cell.max = cell.cands
+			return False
+		if cell.min > cell.max:
+			raise STVException('No result conformant with the constraints is possible')
+	
+	# Rule 2/3
+	for address in ndrange(cm.dimensions):
+		cell = cm.get(address)
+		
+		for dimension in range(len(cm.dimensions)):
+			tot_min_others = 0
+			tot_max_others = 0
+			
+			addr2 = [x for x in address]
+			for i in range(cm.dimensions[dimension]):
+				if i == address[dimension]:
+					continue
+				addr2[dimension] = i
+				
+				cell2 = cm.get(addr2)
+				tot_min_others += cell2.min
+				tot_max_others += cell2.max
+			
+			addr2[dimension] = -1
+			cell_dimension = cm.get(addr2)
+			min_dimension = cell_dimension.min
+			max_dimension = cell_dimension.max
+			
+			# This many must be elected from this cell at least
+			this_cell_min = min_dimension - tot_max_others
+			this_cell_max = max_dimension - tot_min_others
+			
+			# Rule 2
+			if this_cell_min > cell.min:
+				cell.min = this_cell_min
+				return False
+			
+			# Rule 3
+			if this_cell_max < cell.max:
+				cell.max = this_cell_max
+				return False
+	
+	# Rule 4/5
+	for dimension in range(len(cm.dimensions)):
+		for addr1 in ndrange(cm.dimensions[:dimension]):
+			for addr2 in ndrange(cm.dimensions[dimension+1:]):
+				tot_min = 0
+				tot_max = 0
+				
+				for addr_d in range(cm.dimensions[dimension]):
+					address = addr1 + [addr_d] + addr2
+					tot_min += cm.get(address).min
+					tot_max += cm.get(address).max
+				
+				address = addr1 + [-1] + addr2
+				cell = cm.get(address)
+				
+				# Rule 4
+				if cell.min < tot_min:
+					cell.min = tot_min
+				
+				# Rule 5
+				if cell.max > tot_max:
+					cell.max = tot_max
+	
+	return True
diff --git a/tests/test_constraints_otten.py b/tests/test_constraints_otten.py
new file mode 100644
index 0000000..cf65ec7
--- /dev/null
+++ b/tests/test_constraints_otten.py
@@ -0,0 +1,132 @@
+#   pyRCV2: Preferential vote counting
+#   Copyright © 2020–2021  Lee Yingtong Li (RunasSudo)
+#
+#   This program is free software: you can redistribute it and/or modify
+#   it under the terms of the GNU Affero General Public License as published by
+#   the Free Software Foundation, either version 3 of the License, or
+#   (at your option) any later version.
+#
+#   This program is distributed in the hope that it will be useful,
+#   but WITHOUT ANY WARRANTY; without even the implied warranty of
+#   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#   GNU Affero General Public License for more details.
+#
+#   You should have received a copy of the GNU Affero General Public License
+#   along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+from pyRCV2 import constraints
+
+class CounterStub:
+	pass
+
+def test_constraints_otten():
+	counter = CounterStub()
+	counter._constraint_matrix = constraints.ConstraintMatrix([2, 3])
+	cm = counter._constraint_matrix
+	
+	# Fill in details
+	c = cm.get([0, 0]); c.max = 7; c.cands = 4
+	c = cm.get([0, 1]); c.max = 6; c.cands = 11
+	c = cm.get([0, 2]); c.max = 1; c.cands = 2
+	c = cm.get([0, -1]); c.min = 7; c.max = 7; c.cands = 17
+	
+	c = cm.get([1, 0]); c.max = 7; c.cands = 7
+	c = cm.get([1, 1]); c.max = 6; c.cands = 3
+	c = cm.get([1, 2]); c.max = 1; c.cands = 1
+	c = cm.get([1, -1]); c.min = 7; c.max = 7; c.cands = 11
+	
+	c = cm.get([-1, 0]); c.min = 7; c.max = 7; c.cands = 11
+	c = cm.get([-1, 1]); c.min = 6; c.max = 6; c.cands = 14
+	c = cm.get([-1, 2]); c.min = 1; c.max = 1; c.cands = 3
+	c = cm.get([-1, -1]); c.min = 14; c.max = 14; c.cands = 28
+	
+	assert cm.get([0, 0]).max != 4
+	assert cm.get([1, 1]).max != 3
+	assert cm.get([0, 1]).min != 3
+	assert cm.get([1, 0]).min != 3
+	
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == True
+	
+	assert cm.get([0, 0]).max == 4
+	assert cm.get([1, 1]).max == 3
+	assert cm.get([0, 1]).min == 3
+	assert cm.get([1, 0]).min == 3
+	
+	# Election of Welsh Man
+	cm.get([0, 2]).elected += 1
+	cm.get([-1, 2]).elected += 1
+	cm.get([0, -1]).elected += 1
+	
+	assert cm.get([0, 2]).min != 1
+	assert cm.get([0, 0]).max != 3
+	assert cm.get([1, 0]).min != 4
+	assert cm.get([1, 2]).max != 0
+	
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == True
+	
+	assert cm.get([0, 2]).min == 1
+	assert cm.get([0, 0]).max == 3
+	assert cm.get([1, 0]).min == 4
+	assert cm.get([1, 2]).max == 0
+	
+	# Welsh Man and Welsh Woman are doomed
+	assert cm.get([0, 2]).elected == cm.get([0, 2]).max
+	assert cm.get([1, 2]).elected == cm.get([1, 2]).max
+	cm.get([0, 2]).cands = 1
+	cm.get([1, 2]).cands = 0
+	
+	# Election of 2 English Men
+	cm.get([0, 0]).elected += 2
+	cm.get([-1, 0]).elected += 2
+	cm.get([0, -1]).elected += 2
+	
+	# Election of 2 English Women
+	cm.get([1, 0]).elected += 2
+	cm.get([-1, 0]).elected += 2
+	cm.get([1, -1]).elected += 2
+	
+	# Exclusion of Scottish Woman
+	cm.get([1, 1]).cands -= 1
+	cm.get([-1, 1]).cands -= 1
+	cm.get([1, -1]).cands -= 1
+	
+	assert cm.get([0, 0]).min != 2
+	assert cm.get([1, 1]).max != 2
+	assert cm.get([0, 1]).min != 4
+	assert cm.get([1, 0]).min != 5
+	assert cm.get([0, 0]).max != 2
+	assert cm.get([0, 1]).max != 4
+	assert cm.get([1, 1]).min != 2
+	assert cm.get([1, 0]).max != 5
+	
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == False
+	assert constraints.step_matrix(counter) == True
+	
+	assert cm.get([0, 0]).min == 2
+	assert cm.get([1, 1]).max == 2
+	assert cm.get([0, 1]).min == 4
+	assert cm.get([1, 0]).min == 5
+	assert cm.get([0, 0]).max == 2
+	assert cm.get([0, 1]).max == 4
+	assert cm.get([1, 1]).min == 2
+	assert cm.get([1, 0]).max == 5
+	
+	# English Men doomed
+	assert cm.get([0, 0]).elected == cm.get([0, 0]).max
+	# Scottish Women guarded
+	assert cm.get([1, 1]).cands == cm.get([1, 1]).min