diff --git a/Programming/NLP/processing.py b/Programming/NLP/processing.py
index ae4efa32..2342ebd2 100755
--- a/Programming/NLP/processing.py
+++ b/Programming/NLP/processing.py
@@ -1,9 +1,9 @@
 #!/usr/bin/env python3
 """Load the SemEval corpus: sentences, chunkings and gold alignments."""
 
-import sys
 import logging
 import re
+import sys
 from pathlib import Path
 
 import numpy as np
diff --git a/Programming/PORR/code/linear_algebra_utils.py b/Programming/PORR/code/linear_algebra_utils.py
old mode 100644
new mode 100755
index 0d907e00..2108b0bf
--- a/Programming/PORR/code/linear_algebra_utils.py
+++ b/Programming/PORR/code/linear_algebra_utils.py
@@ -1,121 +1,154 @@
+#!/usr/bin/env python3
+"""Four linear-algebra backends with the same interface.
+
+Sequential, threads, processes and dask distributed arrays -- comparing them
+is what this course is about.
+"""
+
 import cmath
 import itertools
 import logging
 import math
 import multiprocessing
 import operator
-from abc import ABC, abstractmethod
+from abc import abstractmethod
 from concurrent.futures import ThreadPoolExecutor
 from functools import partial
 from multiprocessing import Pool
 
 import dask.array as da
+import numpy as np
 from time_measurement import time_accumulator, time_measurement
 
+# These helpers are called with numpy arrays in some backends and plain
+# Python lists in others, which is the point of the comparison this course
+# is about, so the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
 logger = logging.getLogger(__name__)
 
 
 class LinearAlgebraUtils:
-    @staticmethod
-    @abstractmethod
-    def dot_product(v1, v2) -> None:
-        pass
+    """The interface every linear-algebra backend implements."""
 
     @staticmethod
     @abstractmethod
-    def matrix_vector_multiply(A, x) -> None:
-        pass
+    def dot_product(v1: Vector, v2: Vector) -> None:
+        """Return the dot product of two vectors."""
 
     @staticmethod
     @abstractmethod
-    def vector_norm(v) -> None:
-        pass
+    def matrix_vector_multiply(matrix: Matrix, x: Vector) -> None:
+        """Return the matrix-vector product."""
 
     @staticmethod
     @abstractmethod
-    def vector_scalar_divide(x, scalar) -> None:
-        pass
+    def vector_norm(v: Vector) -> None:
+        """Return the Euclidean norm of a vector."""
 
     @staticmethod
     @abstractmethod
-    def matrix_scalar_multiply(A, w) -> None:
-        pass
+    def vector_scalar_divide(x: Vector, scalar: float) -> None:
+        """Divide every element of a vector by a scalar."""
 
     @staticmethod
     @abstractmethod
-    def vector_vector_subtraction(v1, v2) -> None:
-        pass
+    def matrix_scalar_multiply(matrix: Matrix, w: float) -> None:
+        """Multiply every element of a matrix by a scalar."""
 
     @staticmethod
     @abstractmethod
-    def vector_vector_addition(v1, v2) -> None:
-        pass
+    def vector_vector_subtraction(v1: Vector, v2: Vector) -> None:
+        """Subtract one vector from another, elementwise."""
 
     @staticmethod
     @abstractmethod
-    def scalar_vector_multiply(omega, vector) -> None:
-        pass
+    def vector_vector_addition(v1: Vector, v2: Vector) -> None:
+        """Add two vectors elementwise."""
 
     @staticmethod
     @abstractmethod
-    def matrix_norm(A) -> None:
-        pass
+    def scalar_vector_multiply(omega: float, vector: Vector) -> None:
+        """Multiply every element of a vector by a scalar."""
 
     @staticmethod
     @abstractmethod
-    def matrix_matrix_subtraction(A, B) -> None:
-        pass
+    def matrix_norm(matrix: Matrix) -> None:
+        """Return the Frobenius norm of a matrix."""
+
+    @staticmethod
+    @abstractmethod
+    def matrix_matrix_subtraction(matrix: Matrix, other_matrix: Matrix) -> None:
+        """Subtract one matrix from another, elementwise."""
 
 
 class SequentialLinearAlgebraUtils:
-    @staticmethod
-    def dot_product(v1, v2):
-        return sum(x * y for x, y in zip(v1, v2))
+    """Linear algebra in plain Python, one core."""
 
     @staticmethod
-    def matrix_vector_multiply(A, x):
-        return [SequentialLinearAlgebraUtils.dot_product(row, x) for row in A]
+    def dot_product(v1: Vector, v2: Vector) -> float:
+        """Return the dot product of two vectors."""
+        return sum(x * y for x, y in zip(v1, v2, strict=False))
 
     @staticmethod
-    def vector_norm(v):
+    def matrix_vector_multiply(matrix: Matrix, x: Vector) -> Vector:
+        """Return the matrix-vector product."""
+        return [SequentialLinearAlgebraUtils.dot_product(row, x) for row in matrix]
+
+    @staticmethod
+    def vector_norm(v: Vector) -> float:
+        """Return the Euclidean norm of a vector."""
         x_values = (x * x for x in v)
         x_values_sum = sum(x_values)
         return cmath.sqrt(x_values_sum).real
 
     @staticmethod
-    def vector_scalar_divide(x, scalar):
+    def vector_scalar_divide(x: Vector, scalar: float) -> Vector:
+        """Divide every element of a vector by a scalar."""
         return [xi / scalar for xi in x]
 
     @staticmethod
-    def matrix_scalar_multiply(A, w):
-        return A * w
+    def matrix_scalar_multiply(matrix: Matrix, w: float) -> Matrix:
+        """Multiply every element of a matrix by a scalar."""
+        return matrix * w
 
     @staticmethod
-    def vector_vector_subtraction(v1, v2):
-        return [x - y for x, y in zip(v1, v2)]
+    def vector_vector_subtraction(v1: Vector, v2: Vector) -> Vector:
+        """Subtract one vector from another, elementwise."""
+        return [x - y for x, y in zip(v1, v2, strict=False)]
 
     @staticmethod
-    def vector_vector_addition(v1, v2):
-        return [x + y for x, y in zip(v1, v2)]
+    def vector_vector_addition(v1: Vector, v2: Vector) -> Vector:
+        """Add two vectors elementwise."""
+        return [x + y for x, y in zip(v1, v2, strict=False)]
 
     @staticmethod
-    def scalar_vector_multiply(omega, vector):
+    def scalar_vector_multiply(omega: float, vector: Vector) -> Vector:
+        """Multiply every element of a vector by a scalar."""
         return [omega * x for x in vector]
 
     @staticmethod
-    def matrix_norm(A):
-        return math.sqrt(sum(sum(element**2 for element in row) for row in A))
+    def matrix_norm(matrix: Matrix) -> float:
+        """Return the Frobenius norm of a matrix."""
+        return math.sqrt(sum(sum(element**2 for element in row) for row in matrix))
 
     @staticmethod
-    def matrix_matrix_subtraction(A, B):
-        return [[A[i][j] - B[i][j] for j in range(len(A[0]))] for i in range(len(A))]
+    def matrix_matrix_subtraction(matrix: Matrix, other_matrix: Matrix) -> Matrix:
+        """Subtract one matrix from another, elementwise."""
+        return [
+            [matrix[i][j] - other_matrix[i][j] for j in range(len(matrix[0]))]
+            for i in range(len(matrix))
+        ]
 
 
 class ThreadsLinearAlgebraUtils:
+    """Linear algebra split across threads."""
+
     NUM_THREADS = multiprocessing.cpu_count()
 
     @staticmethod
-    def get_chunk_size(data):
+    def get_chunk_size(data: Vector) -> tuple[int, int, int]:
+        """Work out the chunk size, thread count and remainder."""
         num_elements = len(data)
         num_threads = min(ThreadsLinearAlgebraUtils.NUM_THREADS, num_elements)
         chunk_size = num_elements // num_threads
@@ -123,7 +156,8 @@ class ThreadsLinearAlgebraUtils:
         return chunk_size, num_threads, remainder
 
     @staticmethod
-    def divide_vectors_to_chunks(v1, v2):
+    def divide_vectors_to_chunks(v1: Vector, v2: Vector) -> list[tuple[Vector, Vector]]:
+        """Split two vectors into matching chunks."""
         chunk_size, num_threads, remainder = ThreadsLinearAlgebraUtils.get_chunk_size(
             v1
         )
@@ -138,7 +172,8 @@ class ThreadsLinearAlgebraUtils:
         return chunks
 
     @staticmethod
-    def divide_vector_or_matrix_to_chunks(v):
+    def divide_vector_or_matrix_to_chunks(v: Vector) -> list[Vector]:
+        """Split a vector or matrix into chunks."""
         chunk_size, num_threads, remainder = ThreadsLinearAlgebraUtils.get_chunk_size(v)
 
         chunks = []
@@ -152,8 +187,9 @@ class ThreadsLinearAlgebraUtils:
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_vector_multiply(A, x):
-        chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(A)
+    def matrix_vector_multiply(matrix: Matrix, x: Vector) -> Vector:
+        """Return the matrix-vector product."""
+        chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(matrix)
         with ThreadPoolExecutor(
             max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS
         ) as executor:
@@ -163,10 +199,12 @@ class ThreadsLinearAlgebraUtils:
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_norm(v):
+    def vector_norm(v: Vector) -> float:
+        """Return the Euclidean norm of a vector."""
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(v)
 
-        def partial_norm(chunk):
+        def partial_norm(chunk: Vector) -> float:
+            """Return the sum of squares of one chunk, for a parallel norm."""
             return sum(x * x for x in chunk)
 
         with ThreadPoolExecutor(
@@ -178,7 +216,8 @@ class ThreadsLinearAlgebraUtils:
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_vector_subtraction(v1, v2):
+    def vector_vector_subtraction(v1: Vector, v2: Vector) -> Vector:
+        """Subtract one vector from another, elementwise."""
         chunks = ThreadsLinearAlgebraUtils.divide_vectors_to_chunks(v1, v2)
         with ThreadPoolExecutor(
             max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS
@@ -193,7 +232,8 @@ class ThreadsLinearAlgebraUtils:
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_vector_addition(v1, v2):
+    def vector_vector_addition(v1: Vector, v2: Vector) -> Vector:
+        """Add two vectors elementwise."""
         chunks = ThreadsLinearAlgebraUtils.divide_vectors_to_chunks(v1, v2)
         with ThreadPoolExecutor(
             max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS
@@ -206,7 +246,8 @@ class ThreadsLinearAlgebraUtils:
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def scalar_vector_multiply(omega, vector):
+    def scalar_vector_multiply(omega: float, vector: Vector) -> Vector:
+        """Multiply every element of a vector by a scalar."""
         chunks = ThreadsLinearAlgebraUtils.divide_vector_or_matrix_to_chunks(vector)
         with ThreadPoolExecutor(
             max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS
@@ -222,77 +263,90 @@ class ThreadsLinearAlgebraUtils:
 
 
 @time_measurement(time_accumulator)
-def process_row(params):
-    A, k, i = params
-    factor = A[i][k] / A[k][k]
-    return [A[i][j] - factor * A[k][j] for j in range(len(A[0]))]
+def process_row(params: tuple[object, ...]) -> Vector:
+    """Eliminate one row against the pivot row."""
+    matrix, k, i = params
+    factor = matrix[i][k] / matrix[k][k]
+    return [matrix[i][j] - factor * matrix[k][j] for j in range(len(matrix[0]))]
 
 
 @time_measurement(time_accumulator)
-def divide_by_scalar(pair):
+def divide_by_scalar(pair: tuple[Vector, float]) -> Vector:
+    """Divide one chunk of a vector by its scalar."""
     xi, scalar = pair
     return xi / scalar
 
 
 @time_measurement(time_accumulator)
-def multiply_by_scalar(pair):
+def multiply_by_scalar(pair: tuple[Vector, float]) -> Vector:
+    """Multiply one chunk of a vector by its scalar."""
     element, scalar = pair
     return element * scalar
 
 
 class ProcessLinearAlgebraUtils:
+    """Linear algebra split across a process pool."""
+
     @staticmethod
     @time_measurement(time_accumulator)
-    def dot_product(v1, v2):
+    def dot_product(v1: Vector, v2: Vector) -> float:
+        """Return the dot product of two vectors."""
         with Pool() as pool:
             result = pool.starmap(
-                ProcessLinearAlgebraUtils.multiply_elements, zip(v1, v2)
+                ProcessLinearAlgebraUtils.multiply_elements, zip(v1, v2, strict=False)
             )
         return sum(result)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def multiply_elements(x, y):
+    def multiply_elements(x: Vector, y: float) -> float:
+        """Return the product of two numbers."""
         return x * y
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_vector_multiply_row(params):
+    def matrix_vector_multiply_row(params: tuple[object, ...]) -> Vector:
+        """Return one row of a matrix-vector product."""
         row, vector = params
         return SequentialLinearAlgebraUtils.dot_product(row, vector)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_vector_multiply(A, x):
+    def matrix_vector_multiply(matrix: Matrix, x: Vector) -> Vector:
+        """Return the matrix-vector product."""
         with Pool() as pool:
             result = pool.map(
                 ProcessLinearAlgebraUtils.matrix_vector_multiply_row,
-                [(row, x) for row in A],
+                [(row, x) for row in matrix],
             )
         return list(result)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_norm(v):
+    def vector_norm(v: Vector) -> float:
+        """Return the Euclidean norm of a vector."""
         with Pool() as pool:
             squared = pool.map(ProcessLinearAlgebraUtils.square, v)
         return math.sqrt(sum(squared))
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def square(x):
+    def square(x: Vector) -> float:
+        """Return the square of a number."""
         return x * x
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_scalar_divide(x, scalar):
+    def vector_scalar_divide(x: Vector, scalar: float) -> Vector:
+        """Divide every element of a vector by a scalar."""
         with Pool() as pool:
             result = pool.map(divide_by_scalar, [(xi, scalar) for xi in x])
         return list(result)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def divide_vector_by_scalar(x, scalar):
+    def divide_vector_by_scalar(x: Vector, scalar: float) -> Vector:
+        """Divide every element of a vector by a scalar."""
         with Pool() as pool:
             result = pool.map(
                 ProcessLinearAlgebraUtils.vector_scalar_divide,
@@ -302,49 +356,54 @@ class ProcessLinearAlgebraUtils:
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_scalar_multiply_row(params):
+    def matrix_scalar_multiply_row(params: tuple[object, ...]) -> Vector:
+        """Multiply one row of a matrix by a scalar."""
         row, w = params
         return [w * element for element in row]
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_scalar_multiply(A, w):
+    def matrix_scalar_multiply(matrix: Matrix, w: float) -> Matrix:
+        """Multiply every element of a matrix by a scalar."""
         with Pool() as pool:
-            result = pool.map(
+            return pool.map(
                 ProcessLinearAlgebraUtils.matrix_scalar_multiply_row,
-                [(row, w) for row in A],
+                [(row, w) for row in matrix],
             )
-        return result
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_vector_operation(params):
+    def vector_vector_operation(params: tuple[object, ...]) -> Vector:
+        """Apply one elementwise operation to a pair of chunks."""
         v1, v2, op = params
         return op(v1, v2)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_vector_subtraction(v1, v2):
+    def vector_vector_subtraction(v1: Vector, v2: Vector) -> Vector:
+        """Subtract one vector from another, elementwise."""
         with Pool() as pool:
             result = pool.map(
                 ProcessLinearAlgebraUtils.vector_vector_operation,
-                zip(v1, v2, itertools.repeat(operator.sub)),
+                zip(v1, v2, itertools.repeat(operator.sub), strict=False),
             )
         return list(result)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_vector_addition(v1, v2):
+    def vector_vector_addition(v1: Vector, v2: Vector) -> Vector:
+        """Add two vectors elementwise."""
         with Pool() as pool:
             result = pool.map(
                 ProcessLinearAlgebraUtils.vector_vector_operation,
-                zip(v1, v2, itertools.repeat(operator.add)),
+                zip(v1, v2, itertools.repeat(operator.add), strict=False),
             )
         return list(result)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def scalar_vector_multiply(omega, vector):
+    def scalar_vector_multiply(omega: float, vector: Vector) -> Vector:
+        """Multiply every element of a vector by a scalar."""
         with Pool() as pool:
             result = pool.map(
                 multiply_by_scalar, [(element, omega) for element in vector]
@@ -353,120 +412,142 @@ class ProcessLinearAlgebraUtils:
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def sum_of_squares(row):
+    def sum_of_squares(row: Vector) -> float:
+        """Return the sum of the squares of a row."""
         return sum(x**2 for x in row)
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_norm(A):
+    def matrix_norm(matrix: Matrix) -> float:
+        """Return the Frobenius norm of a matrix."""
         with Pool() as pool:
-            row_sums = pool.map(ProcessLinearAlgebraUtils.sum_of_squares, A)
+            row_sums = pool.map(ProcessLinearAlgebraUtils.sum_of_squares, matrix)
         return math.sqrt(sum(row_sums))
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def subtract_rows(row_from_A, row_from_B):
-        return [a - b for a, b in zip(row_from_A, row_from_B)]
+    def subtract_rows(row_from_matrix: Vector, row_from_other: Vector) -> Vector:
+        """Subtract one row from another, elementwise."""
+        return [a - b for a, b in zip(row_from_matrix, row_from_other, strict=False)]
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_matrix_subtraction(A, B):
+    def matrix_matrix_subtraction(matrix: Matrix, other_matrix: Matrix) -> Matrix:
+        """Subtract one matrix from another, elementwise."""
         with Pool() as pool:
-            result = pool.starmap(ProcessLinearAlgebraUtils.subtract_rows, zip(A, B))
-        return result
+            return pool.starmap(
+                ProcessLinearAlgebraUtils.subtract_rows,
+                zip(matrix, other_matrix, strict=False),
+            )
 
 
-class DistributedArraysLinearAlgebraUtils(ABC):
+class DistributedArraysLinearAlgebraUtils:
+    """Linear algebra over dask distributed arrays."""
+
     @staticmethod
     @time_measurement(time_accumulator)
-    def dot_product(v1, v2):
+    def dot_product(v1: Vector, v2: Vector) -> float:
+        """Return the dot product of two vectors."""
         dv1 = da.from_array(v1, chunks="auto")
         dv2 = da.from_array(v2, chunks="auto")
         return da.dot(dv1, dv2).compute()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_vector_multiply(A, x):
-        dA = da.from_array(A, chunks="auto")
+    def matrix_vector_multiply(matrix: Matrix, x: Vector) -> Vector:
+        """Return the matrix-vector product."""
+        dask_matrix = da.from_array(matrix, chunks="auto")
         dx = da.from_array(x, chunks="auto")
-        return da.dot(dA, dx).compute().tolist()
+        return da.dot(dask_matrix, dx).compute().tolist()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_norm(v):
+    def vector_norm(v: Vector) -> float:
+        """Return the Euclidean norm of a vector."""
         dv = da.from_array(v, chunks="auto")
         return da.linalg.norm(dv).compute()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_scalar_divide(x, scalar):
+    def vector_scalar_divide(x: Vector, scalar: float) -> Vector:
+        """Divide every element of a vector by a scalar."""
         dx = da.from_array(x, chunks="auto")
         return (dx / scalar).compute().tolist()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_scalar_multiply(A, w):
-        dA = da.from_array(A, chunks="auto")
-        return (dA * w).compute().tolist()
+    def matrix_scalar_multiply(matrix: Matrix, w: float) -> Matrix:
+        """Multiply every element of a matrix by a scalar."""
+        dask_matrix = da.from_array(matrix, chunks="auto")
+        return (dask_matrix * w).compute().tolist()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_vector_subtraction(v1, v2):
+    def vector_vector_subtraction(v1: Vector, v2: Vector) -> Vector:
+        """Subtract one vector from another, elementwise."""
         dv1 = da.from_array(v1, chunks="auto")
         dv2 = da.from_array(v2, chunks="auto")
         return (dv1 - dv2).compute().tolist()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def vector_vector_addition(v1, v2):
+    def vector_vector_addition(v1: Vector, v2: Vector) -> Vector:
+        """Add two vectors elementwise."""
         dv1 = da.from_array(v1, chunks="auto")
         dv2 = da.from_array(v2, chunks="auto")
         return (dv1 + dv2).compute().tolist()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def scalar_vector_multiply(omega, vector):
+    def scalar_vector_multiply(omega: float, vector: Vector) -> Vector:
+        """Multiply every element of a vector by a scalar."""
         dvector = da.from_array(vector, chunks="auto")
         return (omega * dvector).compute().tolist()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_norm(A):
-        dA = da.from_array(A, chunks="auto")
-        return da.linalg.norm(dA).compute()
+    def matrix_norm(matrix: Matrix) -> float:
+        """Return the Frobenius norm of a matrix."""
+        dask_matrix = da.from_array(matrix, chunks="auto")
+        return da.linalg.norm(dask_matrix).compute()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def matrix_matrix_subtraction(A, B):
-        dA = da.from_array(A, chunks="auto")
-        dB = da.from_array(B, chunks="auto")
-        return (dA - dB).compute().tolist()
+    def matrix_matrix_subtraction(matrix: Matrix, other_matrix: Matrix) -> Matrix:
+        """Subtract one matrix from another, elementwise."""
+        dask_matrix = da.from_array(matrix, chunks="auto")
+        dask_other = da.from_array(other_matrix, chunks="auto")
+        return (dask_matrix - dask_other).compute().tolist()
 
     @staticmethod
     @time_measurement(time_accumulator)
-    def gaussian_elimination(A, b):
+    def gaussian_elimination(matrix: Matrix, b: Vector) -> Vector | None:
+        """Solve the system by Gaussian elimination."""
         try:
-            dA = da.from_array(A, chunks="auto")
+            dask_matrix = da.from_array(matrix, chunks="auto")
             db = da.from_array(b, chunks="auto")
-            Ab = da.hstack([dA, db[:, None]])
-            Ab = Ab.persist()
+            augmented = da.hstack([dask_matrix, db[:, None]])
+            augmented = augmented.persist()
 
-            def elimination_step(Ab, k):
-                n = Ab.shape[0]
-                max_index = da.argmax(da.abs(Ab[k:, k])) + k
-                Ab[[k, max_index]] = Ab[[max_index, k]]
-                Ab = Ab.persist()
-                factor = Ab[k + 1 :, k] / Ab[k, k]
-                Ab[k + 1 :] -= factor[:, None] * Ab[k]
-                return Ab
+            def elimination_step(augmented: Matrix, k: int) -> Matrix:
+                """Run one Gaussian elimination step over the augmented matrix."""
+                max_index = da.argmax(da.abs(augmented[k:, k])) + k
+                augmented[[k, max_index]] = augmented[[max_index, k]]
+                augmented = augmented.persist()
+                factor = augmented[k + 1 :, k] / augmented[k, k]
+                augmented[k + 1 :] -= factor[:, None] * augmented[k]
+                return augmented
 
-            for k in range(A.shape[0]):
-                Ab = elimination_step(Ab, k)
+            for k in range(matrix.shape[0]):
+                augmented = elimination_step(augmented, k)
 
-            x = da.zeros(A.shape[0])
-            for i in range(A.shape[0] - 1, -1, -1):
-                x[i] = (Ab[i, -1] - da.dot(Ab[i, i + 1 : -1], x[i + 1 :])) / Ab[i, i]
+            x = da.zeros(matrix.shape[0])
+            for i in range(matrix.shape[0] - 1, -1, -1):
+                x[i] = (
+                    augmented[i, -1] - da.dot(augmented[i, i + 1 : -1], x[i + 1 :])
+                ) / augmented[i, i]
             return x.compute().tolist()
-        except Exception as e:
-            logger.info("Error during Gaussian elimination: %s", e)
+        except (ValueError, ArithmeticError) as exc:
+            # A singular or badly-scaled system; dask surfaces both here.
+            logger.info("Error during Gaussian elimination: %s", exc)
             return None
diff --git a/Programming/PORR/code/main.py b/Programming/PORR/code/main.py
old mode 100644
new mode 100755
index 10263d27..b5fc1002
--- a/Programming/PORR/code/main.py
+++ b/Programming/PORR/code/main.py
@@ -1,3 +1,4 @@
+#!/usr/bin/env python3
 """Entry point for the PORR benchmark suite.
 
 Runs ``tests.py``, which compares this project's Richardson implementations
diff --git a/Programming/PORR/code/matrix_generator.py b/Programming/PORR/code/matrix_generator.py
old mode 100644
new mode 100755
index 0d379652..10aa35b7
--- a/Programming/PORR/code/matrix_generator.py
+++ b/Programming/PORR/code/matrix_generator.py
@@ -1,53 +1,71 @@
-import os
+#!/usr/bin/env python3
+"""Build and cache the SPD, nemeth12 and poli3 test problems."""
+
+from pathlib import Path
 
 import numpy as np
 import scipy.io
 
+# Generator rather than the legacy global np.random state (NPY002).
+_RNG = np.random.default_rng()
+# These helpers are called with numpy arrays in some backends and plain
+# Python lists in others, which is the point of the comparison this course
+# is about, so the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
+
 
 class MatrixGenerator:
+    """Build or load the test problems this course measures against."""
+
     @staticmethod
     def generate_spd_matrix(n: int) -> np.ndarray:
-        A = np.random.rand(n, n)
-        return np.dot(A, A.T) + n * MatrixGenerator.generate_identity_matrix(
+        """Generate a symmetric positive-definite matrix."""
+        matrix = _RNG.random((n, n))
+        return np.dot(matrix, matrix.T) + n * MatrixGenerator.generate_identity_matrix(
             n
         )  # Adding n*I ensures positive definiteness
 
     @staticmethod
-    def generate_identity_matrix(size):
+    def generate_identity_matrix(size: int) -> Matrix:
+        """Return an identity matrix of the given size."""
         return np.eye(size)
 
     @staticmethod
-    def generate_alternating_vector(size):
+    def generate_alternating_vector(size: int) -> Vector:
+        """Return a vector alternating between 1 and 2."""
         return np.tile([1, 2], int(np.ceil(size / 2)))[:size]
 
     @staticmethod
-    def get_matrix_from_file(file_path, problem):
+    def get_matrix_from_file(file_path: str, problem: int) -> Matrix:
+        """Load a matrix from a MatrixMarket .mat file."""
         mat_contents = scipy.io.loadmat(file_path)
         problem_record = mat_contents["Problem"][0][0]
-        A = problem_record[problem]
-        if scipy.sparse.issparse(A):
-            A_dense = A.todense()
-        else:
-            A_dense = A
-        return np.array(A_dense)
+        matrix = problem_record[problem]
+        dense = matrix.todense() if scipy.sparse.issparse(matrix) else matrix
+        return np.array(dense)
 
     @staticmethod
-    def generate_matrix_and_vector(type, size=None):
-        if type == "spd":
+    def generate_matrix_and_vector(
+        kind: str, size: int | None = None
+    ) -> tuple[Matrix, Vector, float, float]:
+        """Return the matrix, vector and eigenvalue bounds."""
+        if kind == "spd":
             if size is None:
-                raise ValueError("Size must be provided for SPD matrix generation.")
+                msg = "Size must be provided for SPD matrix generation."
+                raise ValueError(msg)
             try:
                 matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file(
                     "spd_" + str(size) + ".npz"
                 )
             except FileNotFoundError:
                 matrix = MatrixGenerator.generate_spd_matrix(size)
-                vector = np.random.uniform(-1, 1, size)
+                vector = _RNG.uniform(-1, 1, size)
                 lambda_min, lambda_max = MatrixGenerator.calculate_eigenvalues(matrix)
                 MatrixGenerator.save_to_file(
                     matrix, vector, lambda_min, lambda_max, "spd_" + str(size) + ".npz"
                 )
-        elif type == "nemeth12":
+        elif kind == "nemeth12":
             try:
                 matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file(
                     "nemeth12.npz"
@@ -60,7 +78,7 @@ class MatrixGenerator:
                 MatrixGenerator.save_to_file(
                     matrix, vector, lambda_min, lambda_max, "nemeth12.npz"
                 )
-        elif type == "poli3":
+        elif kind == "poli3":
             try:
                 matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file(
                     "poli3.npz"
@@ -74,21 +92,28 @@ class MatrixGenerator:
                     matrix, vector, lambda_min, lambda_max, "poli3.npz"
                 )
         else:
-            raise ValueError(
-                "Invalid type specified. Choose 'spd', 'nemeth12', or 'poli3'."
-            )
+            msg = "Invalid type specified. Choose 'spd', 'nemeth12', or 'poli3'."
+            raise ValueError(msg)
 
         return matrix, vector, lambda_min, lambda_max
 
     @staticmethod
-    def calculate_eigenvalues(A):
-        eigenvalues = np.linalg.eigvals(A)
+    def calculate_eigenvalues(matrix: Matrix) -> tuple[float, float]:
+        """Return the smallest and largest eigenvalue."""
+        eigenvalues = np.linalg.eigvals(matrix)
         lambda_min = np.min(eigenvalues)
         lambda_max = np.max(eigenvalues)
         return lambda_min, lambda_max
 
     @staticmethod
-    def save_to_file(matrix, vector, lambda_min, lambda_max, file_path) -> None:
+    def save_to_file(
+        matrix: Matrix,
+        vector: Vector,
+        lambda_min: float,
+        lambda_max: float,
+        file_path: str,
+    ) -> None:
+        """Cache a generated problem so the next run can reuse it."""
         np.savez(
             file_path,
             matrix=matrix,
@@ -97,9 +122,12 @@ class MatrixGenerator:
             lambda_max=lambda_max,
         )
 
-    def load_from_file(file_path):
-        if not os.path.exists(file_path):
-            raise FileNotFoundError(f"The file {file_path} does not exist.")
+    @staticmethod
+    def load_from_file(file_path: str) -> tuple[Matrix, Vector, float, float]:
+        """Load a cached problem."""
+        if not Path(file_path).exists():
+            msg = f"The file {file_path} does not exist."
+            raise FileNotFoundError(msg)
         data = np.load(file_path)
         matrix = data["matrix"]
         vector = data["vector"]
diff --git a/Programming/PORR/code/processing_type.py b/Programming/PORR/code/processing_type.py
old mode 100644
new mode 100755
index 8bedd46d..ccaea421
--- a/Programming/PORR/code/processing_type.py
+++ b/Programming/PORR/code/processing_type.py
@@ -1,7 +1,12 @@
+#!/usr/bin/env python3
+"""The enum naming the four processing backends."""
+
 from enum import Enum, auto
 
 
 class ProcessingType(Enum):
+    """Which backend a run should use."""
+
     SEQUENTIAL = auto()
     THREADS = auto()
     PROCESSES = auto()
diff --git a/Programming/PORR/code/richardson_method.py b/Programming/PORR/code/richardson_method.py
old mode 100644
new mode 100755
index 7a22ee4e..f1c8c2f1
--- a/Programming/PORR/code/richardson_method.py
+++ b/Programming/PORR/code/richardson_method.py
@@ -1,79 +1,100 @@
+#!/usr/bin/env python3
+"""Richardson iteration, parameterised by linear-algebra backend."""
+
 import gc
 import logging
 import time
 
-import linear_algebra_utils as linAlg
+import linear_algebra_utils as lin_alg
+import numpy as np
 from processing_type import ProcessingType
+from richardson_problem import Problem, Settings
 from time_measurement import tests_time, time_accumulator, time_measurement
 
+# These helpers are called with numpy arrays in some backends and plain
+# Python lists in others, which is the point of the comparison this course
+# is about, so the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
 logger = logging.getLogger(__name__)
 
 
 class RichardsonMethod:
+    """Richardson iteration over a pluggable linear-algebra backend."""
+
     @time_measurement(time_accumulator)
     def __init__(
         self,
         method: ProcessingType,
-        A,
-        b,
-        lambda_min,
-        lambda_max,
-        max_iterations,
-        size: int,
-        x0=None,
-        tol=1e-5,
+        problem: Problem,
+        settings: Settings,
     ) -> None:
-        self.LinAlg = self.assign_LinAlgType(method)
-        self.A = A
-        self.b = b
-        self.x0 = x0 if x0 is not None else [0.0] * len(b)
-        self.max_iterations = max_iterations
-        self.tol = tol
-        self.lambda_min = lambda_min
-        self.lambda_max = lambda_max
+        """Set up the solver for one problem.
+
+        `size` used to be a parameter and was never read -- len(b) is the
+        size, and it is right here.
+        """
+        self.LinAlg = self.assign_lin_alg_type(method)
+        self.matrix = problem.matrix
+        self.b = problem.rhs
+        self.x0 = settings.x0 if settings.x0 is not None else [0.0] * len(problem.rhs)
+        self.max_iterations = settings.max_iterations
+        self.tol = settings.tol
+        self.lambda_min = problem.lambda_min
+        self.lambda_max = problem.lambda_max
         if self.lambda_min < 0:
-            raise ValueError("Matrix A is not positive semi-definite.")
+            msg = "Matrix A is not positive semi-definite."
+            raise ValueError(msg)
         self.omega = RichardsonMethod.calculate_omega(self.lambda_min, self.lambda_max)
 
     @staticmethod
-    def calculate_omega(lambda_min, lambda_max):
+    def calculate_omega(lambda_min: float, lambda_max: float) -> float:
+        """Return the optimal Richardson relaxation factor."""
         return 2 / (lambda_min + lambda_max)
 
     @staticmethod
-    def convergence_norm(LinAlgType, A, omega, I) -> bool:
-        wA = LinAlgType.matrix_scalar_multiply(A, omega)
-        IMinuswA = LinAlgType.matrix_matrix_subtraction(I, wA)
-        return LinAlgType.matrix_norm(IMinuswA)
+    def convergence_norm(
+        lin_alg_type: type, matrix: Matrix, omega: float, identity: Matrix
+    ) -> bool:
+        """Report whether the iteration matrix norm is under 1."""
+        scaled_matrix = lin_alg_type.matrix_scalar_multiply(matrix, omega)
+        identity_minus_scaled = lin_alg_type.matrix_matrix_subtraction(
+            identity, scaled_matrix
+        )
+        return lin_alg_type.matrix_norm(identity_minus_scaled)
 
     @staticmethod
-    def assign_LinAlgType(method):
+    def assign_lin_alg_type(method: ProcessingType) -> type:
+        """Return the linear-algebra backend for a processing type."""
         methods = {
-            ProcessingType.SEQUENTIAL: linAlg.SequentialLinearAlgebraUtils,
-            ProcessingType.THREADS: linAlg.ThreadsLinearAlgebraUtils,
-            ProcessingType.PROCESSES: linAlg.ProcessLinearAlgebraUtils,
-            ProcessingType.DISTRIBUTED_ARRAYS: linAlg.DistributedArraysLinearAlgebraUtils,
+            ProcessingType.SEQUENTIAL: lin_alg.SequentialLinearAlgebraUtils,
+            ProcessingType.THREADS: lin_alg.ThreadsLinearAlgebraUtils,
+            ProcessingType.PROCESSES: lin_alg.ProcessLinearAlgebraUtils,
+            ProcessingType.DISTRIBUTED_ARRAYS: (
+                lin_alg.DistributedArraysLinearAlgebraUtils
+            ),
         }
 
         try:
             return methods[method]
-        except KeyError:
-            raise ValueError(
-                "Unknown method, please use 'SEQUENTIAL', 'THREADS' or 'PROCESSES'."
-            )
+        except KeyError as exc:
+            msg = "Unknown method, please use 'SEQUENTIAL', 'THREADS' or 'PROCESSES'."
+            raise ValueError(msg) from exc
 
-    def solve(self):
+    def solve(self) -> tuple[Vector, str | None]:
+        """Run the Richardson iteration to convergence, or give up."""
         gc.disable()
         time_accumulator.total_time = 0
         start = time.perf_counter()
         x = self.x0[:]
 
-        for iteration in range(self.max_iterations):
-            Ax = self.LinAlg.matrix_vector_multiply(self.A, x)
-            residual = self.LinAlg.vector_vector_subtraction(self.b, Ax)
+        for _iteration in range(self.max_iterations):
+            product = self.LinAlg.matrix_vector_multiply(self.matrix, x)
+            residual = self.LinAlg.vector_vector_subtraction(self.b, product)
             x = self.LinAlg.vector_vector_addition(
                 x, self.LinAlg.scalar_vector_multiply(self.omega, residual)
             )
-            if linAlg.SequentialLinearAlgebraUtils.vector_norm(residual) < self.tol:
+            if lin_alg.SequentialLinearAlgebraUtils.vector_norm(residual) < self.tol:
                 break
 
         end = time.perf_counter()
@@ -81,11 +102,11 @@ class RichardsonMethod:
         gc.enable()
 
         match self.LinAlg:
-            case linAlg.SequentialLinearAlgebraUtils:
+            case lin_alg.SequentialLinearAlgebraUtils:
                 logger.info(
                     "Total: %ss, Tests time: %ss", total_time, tests_time.total_time
                 )
-            case linAlg.ThreadsLinearAlgebraUtils:
+            case lin_alg.ThreadsLinearAlgebraUtils:
                 sequential_time = total_time - time_accumulator.total_time
                 logger.info(
                     "Total: %ss, Seq: %ss, Parallel (threads): %ss, Tests time: %ss",
@@ -94,7 +115,7 @@ class RichardsonMethod:
                     time_accumulator.total_time,
                     tests_time.total_time,
                 )
-            case linAlg.ProcessLinearAlgebraUtils:
+            case lin_alg.ProcessLinearAlgebraUtils:
                 sequential_time = total_time - time_accumulator.total_time
                 logger.info(
                     "Total: %ss, Seq: %ss, Parallel (processes): %ss, Tests time: %ss",
@@ -103,10 +124,11 @@ class RichardsonMethod:
                     time_accumulator.total_time,
                     tests_time.total_time,
                 )
-            case linAlg.DistributedArraysLinearAlgebraUtils:
+            case lin_alg.DistributedArraysLinearAlgebraUtils:
                 sequential_time = total_time - time_accumulator.total_time
                 logger.info(
-                    "Total: %ss, Seq: %ss, Parallel (distributed arrays): %ss, Tests time: %ss",
+                    "Total: %ss, Seq: %ss, Parallel (distributed arrays): %ss, "
+                    "Tests time: %ss",
                     total_time,
                     sequential_time,
                     time_accumulator.total_time,
diff --git a/Programming/PORR/code/richardson_mpi.py b/Programming/PORR/code/richardson_mpi.py
old mode 100644
new mode 100755
index b69d476e..7ac9e919
--- a/Programming/PORR/code/richardson_mpi.py
+++ b/Programming/PORR/code/richardson_mpi.py
@@ -1,25 +1,42 @@
+#!/usr/bin/env python3
+"""Richardson iteration split across MPI ranks."""
+
 import logging
 import time
 
 import numpy as np
 from mpi4py import MPI
+from richardson_problem import Problem, Settings
 
+# These helpers are called with numpy arrays in some backends and plain
+# Python lists in others, which is the point of the comparison this course
+# is about, so the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
 logger = logging.getLogger(__name__)
 
 
-def richardson_parallel(A, b, lambda_min, lambda_max, tol=1e-5, max_iter=10000):
+def richardson_parallel(
+    problem: Problem, settings: Settings | None = None
+) -> tuple[Vector, float]:
+    """Run Richardson across MPI ranks."""
+    if settings is None:
+        settings = Settings(max_iterations=1000)
+    matrix = problem.matrix
+    b = problem.rhs
+    lambda_min = problem.lambda_min
+    lambda_max = problem.lambda_max
+    tol = settings.tol
+    max_iter = settings.max_iterations
     comm = MPI.COMM_WORLD
     rank = comm.Get_rank()
     size = comm.Get_size()
 
     # Rozmiar macierzy A
-    n = A.shape[0]
+    n = matrix.shape[0]
 
     # Obliczanie wartości własnych tylko na jednym procesie
-    if rank == 0:
-        omega = 2 / (lambda_min + lambda_max)
-    else:
-        omega = None
+    omega = 2 / (lambda_min + lambda_max) if rank == 0 else None
 
     # Rozgłoszenie omega do wszystkich procesów
     omega = comm.bcast(omega, root=0)
@@ -33,7 +50,7 @@ def richardson_parallel(A, b, lambda_min, lambda_max, tol=1e-5, max_iter=10000):
     end_row = start_row + local_rows if rank != size - 1 else n
 
     # Przydzielenie lokalnych porcji A i b
-    local_A = A[start_row:end_row, :]
+    local_matrix = matrix[start_row:end_row, :]
     local_b = b[start_row:end_row]
 
     # Lokalny wektor residuum
@@ -44,9 +61,9 @@ def richardson_parallel(A, b, lambda_min, lambda_max, tol=1e-5, max_iter=10000):
 
     start_time = time.time()
 
-    for i in range(max_iter):
+    for _i in range(max_iter):
         # Oblicz lokalny residuum r = b - A @ x
-        local_r[:] = local_b - np.dot(local_A, x)
+        local_r[:] = local_b - np.dot(local_matrix, x)
 
         # Tworzymy tymczasowy wektor o pełnym rozmiarze i kopiujemy lokalne dane
         temp_r = np.zeros_like(b, dtype=np.float64)
@@ -69,8 +86,11 @@ def richardson_parallel(A, b, lambda_min, lambda_max, tol=1e-5, max_iter=10000):
     return x, execution_time
 
 
-def check_solution(A, b, x_approx, tolerance=8e-3):
-    x_true = np.linalg.solve(A, b)
+def check_solution(
+    matrix: Matrix, b: Vector, x_approx: Vector, tolerance: float = 0.008
+) -> tuple[bool, float]:
+    """Report whether an approximate solution is within tolerance."""
+    x_true = np.linalg.solve(matrix, b)
     error = np.linalg.norm(x_true - x_approx)
     return error < tolerance, error
 
@@ -85,25 +105,25 @@ if __name__ == "__main__":
     sizes = [2, 5, 10, 50, 80, 100, 300, 500, 750, 1000, 5000, 10000]
     for i in sizes:
         if rank == 0:
-            A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
-                "spd", size=i
+            matrix, b, lambda_min, lambda_max = (
+                MatrixGenerator.generate_matrix_and_vector("spd", size=i)
             )
         else:
-            A = None
+            matrix = None
             b = None
             lambda_min = None
             lambda_max = None
 
-        A = comm.bcast(A, root=0)
+        matrix = comm.bcast(matrix, root=0)
         b = comm.bcast(b, root=0)
 
         # Rozwiązanie przy użyciu zrównoleglonej metody Richardsona
-        x, time_taken = richardson_parallel(A, b, lambda_min, lambda_max)
+        x, time_taken = richardson_parallel(Problem(matrix, b, lambda_min, lambda_max))
 
         # Sprawdzanie poprawności rozwiązania (na procesie 0)
         if rank == 0:
             logger.info("Spd matrix with size %s", i)
-            is_correct, error = check_solution(A, b, x)
+            is_correct, error = check_solution(matrix, b, x)
             logger.info("Czas wykonania [s]: %s", time_taken)
             logger.info(
                 "Czy rozwiązanie jest poprawne: %s", "Tak" if is_correct else "Nie"
@@ -111,49 +131,49 @@ if __name__ == "__main__":
             logger.info("Błąd rozwiązania: %s", error)
 
     if rank == 0:
-        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
+        matrix, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
             "nemeth12"
         )
     else:
-        A = None
+        matrix = None
         b = None
         lambda_min = None
         lambda_max = None
 
-    A = comm.bcast(A, root=0)
+    matrix = comm.bcast(matrix, root=0)
     b = comm.bcast(b, root=0)
 
     # Rozwiązanie przy użyciu zrównoleglonej metody Richardsona
-    x, time_taken = richardson_parallel(A, b, lambda_min, lambda_max)
+    x, time_taken = richardson_parallel(Problem(matrix, b, lambda_min, lambda_max))
 
     # Sprawdzanie poprawności rozwiązania (na procesie 0)
     if rank == 0:
         logger.info("Nemeth12 matrix")
-        is_correct, error = check_solution(A, b, x)
+        is_correct, error = check_solution(matrix, b, x)
         logger.info("Czas wykonania [s]: %s", time_taken)
         logger.info("Czy rozwiązanie jest poprawne: %s", "Tak" if is_correct else "Nie")
         logger.info("Błąd rozwiązania: %s", error)
 
     if rank == 0:
-        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
+        matrix, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
             "poli3"
         )
     else:
-        A = None
+        matrix = None
         b = None
         lambda_min = None
         lambda_max = None
 
-    A = comm.bcast(A, root=0)
+    matrix = comm.bcast(matrix, root=0)
     b = comm.bcast(b, root=0)
 
     # Rozwiązanie przy użyciu zrównoleglonej metody Richardsona
-    x, time_taken = richardson_parallel(A, b, lambda_min, lambda_max)
+    x, time_taken = richardson_parallel(Problem(matrix, b, lambda_min, lambda_max))
 
     # Sprawdzanie poprawności rozwiązania (na procesie 0)
     if rank == 0:
         logger.info("Poli3 matrix")
-        is_correct, error = check_solution(A, b, x)
+        is_correct, error = check_solution(matrix, b, x)
         logger.info("Czas wykonania [s]: %s", time_taken)
         logger.info("Czy rozwiązanie jest poprawne: %s", "Tak" if is_correct else "Nie")
         logger.info("Błąd rozwiązania: %s", error)
diff --git a/Programming/PORR/code/richardson_problem.py b/Programming/PORR/code/richardson_problem.py
new file mode 100755
index 00000000..4edfc9d3
--- /dev/null
+++ b/Programming/PORR/code/richardson_problem.py
@@ -0,0 +1,34 @@
+#!/usr/bin/env python3
+"""The system to solve, and the settings to solve it with.
+
+Every Richardson backend in this project took the same seven or nine loose
+arguments -- matrix, right-hand side, both eigenvalue bounds, iteration cap,
+starting guess and tolerance -- which is what put all of them over ruff's
+argument limit. They travel together, so they are two dataclasses.
+"""
+
+from dataclasses import dataclass
+
+import numpy as np
+
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
+
+
+@dataclass(frozen=True)
+class Problem:
+    """One linear system, with the eigenvalue bounds Richardson needs."""
+
+    matrix: Matrix
+    rhs: Vector
+    lambda_min: float
+    lambda_max: float
+
+
+@dataclass(frozen=True)
+class Settings:
+    """How hard to try before giving up."""
+
+    max_iterations: int
+    x0: Vector | None = None
+    tol: float = 1e-05
diff --git a/Programming/PORR/code/tests.py b/Programming/PORR/code/tests.py
old mode 100644
new mode 100755
index 8f049ab0..5cf41117
--- a/Programming/PORR/code/tests.py
+++ b/Programming/PORR/code/tests.py
@@ -1,3 +1,6 @@
+#!/usr/bin/env python3
+"""Compare every Richardson backend against numpy's direct solver."""
+
 import logging
 
 import numpy as np
@@ -5,34 +8,48 @@ import pytest
 from matrix_generator import MatrixGenerator
 from processing_type import ProcessingType
 from richardson_method import RichardsonMethod
-from threads_indep import RichardsonMethodThreads
+from richardson_problem import Problem, Settings
+from threads_indep import richardson_method_threads
 from time_measurement import tests_time, time_measurement
 
+# These helpers are called with numpy arrays in some backends and plain
+# Python lists in others, which is the point of the comparison this course
+# is about, so the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
 logger = logging.getLogger(__name__)
 
+# nemeth12 and poli3 come from a file, so their size is fixed; the suite
+# only runs them under the first n to avoid repeating identical work.
+FIXED_SIZE_N = 2
 
-def calculate_norm_numpy(I, w, A):
-    difference = I - w * A
+
+def calculate_norm_numpy(identity: Matrix, w: float, matrix: Matrix) -> float:
+    """Return the norm of I - wA, computed with numpy."""
+    difference = identity - w * matrix
     return np.linalg.norm(difference)
 
 
-def calculate_eigenvalues(A):
-    eigenvalues = np.linalg.eigvals(A)
+def calculate_eigenvalues(matrix: Matrix) -> tuple[float, float]:
+    """Return the smallest and largest eigenvalue."""
+    eigenvalues = np.linalg.eigvals(matrix)
     lambda_min = np.min(eigenvalues)
     lambda_max = np.max(eigenvalues)
     return lambda_min, lambda_max
 
 
-def calcualte_norm_from_matrix_numpy(A, n):
-    lambda_min, lambda_max = calculate_eigenvalues(A)
+def calcualte_norm_from_matrix_numpy(matrix: Matrix, n: int) -> float:
+    """Return the iteration-matrix norm for a matrix."""
+    lambda_min, lambda_max = calculate_eigenvalues(matrix)
     omega = 2 / (lambda_min + lambda_max)
-    I = np.eye(n)
-    return calculate_norm_numpy(I, omega, A)
+    identity = np.eye(n)
+    return calculate_norm_numpy(identity, omega, matrix)
 
 
 @time_measurement(tests_time)
-def solution_lib(A, b):
-    return np.linalg.solve(A, b)
+def solution_lib(matrix: Matrix, b: Vector) -> Vector:
+    """Solve the system with numpy, as the reference answer."""
+    return np.linalg.solve(matrix, b)
 
 
 @pytest.mark.parametrize("n", [2, 5, 10, 50, 100, 300, 500, 750, 1000, 5000, 10000])
@@ -47,86 +64,99 @@ def solution_lib(A, b):
 )
 @pytest.mark.parametrize("matrix_type", ["spd", "nemeth12", "poli3"])
 def test_richardson_vs_cg(
-    n: int, processing_type: ProcessingType, matrix_type: str, capsys
+    n: int,
+    processing_type: ProcessingType,
+    matrix_type: str,
+    capsys: pytest.CaptureFixture[str],
 ) -> None:
+    """Every backend must agree with numpy's direct solver, at every size."""
     logger.info("matrix type:  %s", matrix_type)
     logger.info("matrix size:  %s", n if matrix_type == "spd" else "fixed")
     tolerance = 8e-3
     max_iterations = 100
-    if matrix_type in ["nemeth12", "poli3"] and n != 2:
+    if matrix_type in ("nemeth12", "poli3") and n != FIXED_SIZE_N:
         pytest.skip(
             "Fixed matrix size for nemeth12 and poli3, skipping redundant runs."
         )
 
     if matrix_type == "spd":
-        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
+        matrix, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
             "spd", size=n
         )
     elif matrix_type == "poli3":
-        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
+        matrix, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
             "poli3"
         )
     elif matrix_type == "nemeth12":
-        A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
+        matrix, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector(
             "nemeth12"
         )
     else:
-        raise ValueError(
-            "Invalid matrix type specified. Choose 'spd', 'poli3', or 'nemeth12'."
-        )
+        msg = "Invalid matrix type specified. Choose 'spd', 'poli3', or 'nemeth12'."
+        raise ValueError(msg)
 
     solution_richardson, info_richardson = None, None
 
     if processing_type != ProcessingType.THREADS:
         richardson_solver = RichardsonMethod(
             processing_type,
-            A,
-            b,
-            lambda_min,
-            lambda_max,
-            max_iterations,
-            size=A.shape[0],
-            tol=1e-7,
+            Problem(matrix, b, lambda_min, lambda_max),
+            Settings(max_iterations, tol=1e-7),
         )
         with capsys.disabled():
             solution_richardson, info_richardson = richardson_solver.solve()
     else:
         with capsys.disabled():
-            solution_richardson, info_richardson = RichardsonMethodThreads(
-                A, b, lambda_min, lambda_max, max_iterations, tol=1e-7
+            solution_richardson, info_richardson = richardson_method_threads(
+                Problem(matrix, b, lambda_min, lambda_max),
+                Settings(max_iterations, tol=1e-7),
             )
 
-    # Przechwytywanie wyjścia po solve
+            # Przechwytywanie wyjścia po solve
     captured = capsys.readouterr()
     logger.info("Captured output: %s", captured.out)
 
-    solution = solution_lib(A, b)
+    solution = solution_lib(matrix, b)
 
-    assert_converged(solution_richardson, info_richardson, solution, tolerance, A, n)
+    assert_converged(
+        (solution_richardson, info_richardson),
+        solution,
+        tolerance,
+        (matrix, n),
+    )
 
 
 def assert_converged(
-    solution_richardson, info_richardson, solution, tolerance, A, n
+    result: tuple[Vector, str | None],
+    solution: Vector,
+    tolerance: float,
+    problem: tuple[Matrix, int],
 ) -> None:
+    """Raise unless Richardson agrees with numpy within tolerance.
+
+    ``result`` is what a backend returned, ``problem`` is (matrix, n).
+    """
+    solution_richardson, info_richardson = result
+    matrix, n = problem
     if info_richardson == "Richardson method for those values will NOT converge":
-        numpy_norm = calcualte_norm_from_matrix_numpy(A, n)
+        numpy_norm = calcualte_norm_from_matrix_numpy(matrix, n)
         logger.info(
             "Numpy norm:  %s  Richardson norm:  %s", numpy_norm, solution_richardson
         )
-        assert False, "Richardson did not converge"
-    else:
-        difference = np.linalg.norm(solution_richardson - solution)
-        logger.info("Difference between Richardson and numpy solutions: %s", difference)
-        if difference < tolerance:
-            logger.info(
-                "Both Richardson and numpy method converged and calculated correct values."
-            )
-        else:
-            logger.info("Solution numpy:\n %s", solution)
-            logger.info("Solution Richardson:\n %s", solution_richardson)
-        assert difference < tolerance, (
-            f"The solutions are different! Difference: {difference:.8f}"
+        msg = "Richardson did not converge"
+        raise AssertionError(msg)
+    difference = np.linalg.norm(solution_richardson - solution)
+    logger.info("Difference between Richardson and numpy solutions: %s", difference)
+    if difference < tolerance:
+        logger.info(
+            "Both Richardson and numpy method converged and calculated correct values."
         )
+    else:
+        logger.info("Solution numpy:\n %s", solution)
+        logger.info("Solution Richardson:\n %s", solution_richardson)
+    if difference >= tolerance:
+        msg = f"The solutions are different! Difference: {difference:.8f}"
+        raise AssertionError(msg)
 
 
 if __name__ == "__main__":
diff --git a/Programming/PORR/code/threads_indep.py b/Programming/PORR/code/threads_indep.py
old mode 100644
new mode 100755
index 7a7421cd..139f78f7
--- a/Programming/PORR/code/threads_indep.py
+++ b/Programming/PORR/code/threads_indep.py
@@ -1,9 +1,14 @@
+#!/usr/bin/env python3
+"""Threaded Richardson with njit-compiled kernels."""
+
 import gc
 import logging
 import time
 
-import linear_algebra_utils as linAlg
+import linear_algebra_utils as lin_alg
+import numpy as np
 from numba import njit, prange
+from richardson_problem import Problem, Settings
 from time_measurement import (
     longest_threads_time_accumulator,
     tests_time,
@@ -12,60 +17,89 @@ from time_measurement import (
 
 logger = logging.getLogger(__name__)
 
+# Called with numpy arrays here and with plain lists in the sequential
+# backend; the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
+
 
 @njit(parallel=True)
-def numba_matrix_vector_multiply(A, input_x, Ax) -> None:
-    for i in prange(len(A)):
+def numba_matrix_vector_multiply(
+    matrix: Matrix, input_x: Vector, product: Vector
+) -> None:
+    """Matrix-vector product, njit-compiled."""
+    for i in prange(len(matrix)):
         acc = 0.0
         for j in range(len(input_x)):
-            acc += A[i][j] * input_x[j]
-        Ax[i] = acc
+            acc += matrix[i][j] * input_x[j]
+        product[i] = acc
 
 
 @njit(parallel=True)
-def numba_vector_vector_subtraction(b, Ax, residual) -> None:
+def numba_vector_vector_subtraction(
+    b: Vector, product: Vector, residual: Vector
+) -> None:
+    """Elementwise subtraction, njit-compiled."""
     for i in prange(len(b)):
-        residual[i] = b[i] - Ax[i]
+        residual[i] = b[i] - product[i]
 
 
 @njit(nopython=True)
-def numba_scalar_vector_multiply(omega, vector, result) -> None:
+def numba_scalar_vector_multiply(omega: float, vector: Vector, result: Vector) -> None:
+    """Scalar-vector product, njit-compiled."""
     omega_real = omega.real
     for i in range(len(vector)):
         result[i] = omega_real * vector[i]
 
 
 @njit(parallel=True)
-def numba_vector_vector_addition(input_x, vector, output_x) -> None:
+def numba_vector_vector_addition(
+    input_x: Vector, vector: Vector, output_x: Vector
+) -> None:
+    """Elementwise addition, njit-compiled."""
     for i in prange(len(input_x)):
         output_x[i] = input_x[i] + vector[i]
 
-
-# Funkcje z dekoratorem
-@time_measurement_longest(longest_threads_time_accumulator)
-def matrix_vector_multiply(A, input_x, Ax) -> None:
-    numba_matrix_vector_multiply(A, input_x, Ax)
+        # Funkcje z dekoratorem
 
 
 @time_measurement_longest(longest_threads_time_accumulator)
-def vector_vector_subtraction(b, Ax, residual) -> None:
-    numba_vector_vector_subtraction(b, Ax, residual)
+def matrix_vector_multiply(matrix: Matrix, input_x: Vector, product: Vector) -> None:
+    """Return the matrix-vector product."""
+    numba_matrix_vector_multiply(matrix, input_x, product)
 
 
 @time_measurement_longest(longest_threads_time_accumulator)
-def scalar_vector_multiply(omega, vector, result) -> None:
+def vector_vector_subtraction(b: Vector, product: Vector, residual: Vector) -> None:
+    """Subtract one vector from another, elementwise."""
+    numba_vector_vector_subtraction(b, product, residual)
+
+
+@time_measurement_longest(longest_threads_time_accumulator)
+def scalar_vector_multiply(omega: float, vector: Vector, result: Vector) -> None:
+    """Multiply every element of a vector by a scalar."""
     numba_scalar_vector_multiply(omega, vector, result)
 
 
 @time_measurement_longest(longest_threads_time_accumulator)
-def vector_vector_addition(input_x, vector, output_x) -> None:
+def vector_vector_addition(input_x: Vector, vector: Vector, output_x: Vector) -> None:
+    """Add two vectors elementwise."""
     numba_vector_vector_addition(input_x, vector, output_x)
 
+    # Metoda Richardson z obsługą wątków
 
-# Metoda Richardson z obsługą wątków
-def RichardsonMethodThreads(
-    A, b, lambda_min, lambda_max, max_iterations, x0=None, tol=1e-5
-):
+
+def richardson_method_threads(
+    problem: Problem, settings: Settings
+) -> tuple[Vector, str | None]:
+    """Run Richardson with one thread per chunk."""
+    matrix = problem.matrix
+    b = problem.rhs
+    lambda_min = problem.lambda_min
+    lambda_max = problem.lambda_max
+    max_iterations = settings.max_iterations
+    x0 = settings.x0
+    tol = settings.tol
     longest_threads_time_accumulator.hard_reset()
 
     gc.disable()
@@ -77,13 +111,13 @@ def RichardsonMethodThreads(
     omega = 2 / (lambda_min + lambda_max)
     n = len(b)
 
-    for iteration in range(max_iterations):
-        Ax = [0.0] * n
-        matrix_vector_multiply(A, x, Ax)
+    for _iteration in range(max_iterations):
+        product = [0.0] * n
+        matrix_vector_multiply(matrix, x, product)
         longest_threads_time_accumulator.save_lap_and_reset()
 
         residual = [0.0] * n
-        vector_vector_subtraction(b, Ax, residual)
+        vector_vector_subtraction(b, product, residual)
         longest_threads_time_accumulator.save_lap_and_reset()
 
         change_vector = [0.0] * n
@@ -96,7 +130,7 @@ def RichardsonMethodThreads(
 
         x = _x[:]
 
-        if linAlg.SequentialLinearAlgebraUtils.vector_norm(residual) < tol:
+        if lin_alg.SequentialLinearAlgebraUtils.vector_norm(residual) < tol:
             break
 
     end_time = time.perf_counter()
diff --git a/Programming/PORR/code/threads_indep_old.py b/Programming/PORR/code/threads_indep_old.py
old mode 100644
new mode 100755
index a97c29d3..c68a662f
--- a/Programming/PORR/code/threads_indep_old.py
+++ b/Programming/PORR/code/threads_indep_old.py
@@ -1,10 +1,16 @@
+#!/usr/bin/env python3
+"""Threaded Richardson, the earlier index-slicing version."""
+
 import gc
 import logging
 import multiprocessing
 import time
+from collections.abc import Callable
 from concurrent.futures import ThreadPoolExecutor
 
-import linear_algebra_utils as linAlg
+import linear_algebra_utils as lin_alg
+import numpy as np
+from richardson_problem import Problem, Settings
 from time_measurement import (
     longest_threads_time_accumulator,
     tests_time,
@@ -13,30 +19,89 @@ from time_measurement import (
 
 logger = logging.getLogger(__name__)
 
-
-@time_measurement_longest(longest_threads_time_accumulator)
-def matrix_vector_multiply(A, input_x, start, end, Ax) -> None:
-    Ax[start:end] = [sum(x * y for x, y in zip(row, input_x)) for row in A[start:end]]
+# Called with numpy arrays here and with plain lists in the sequential
+# backend; the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
 
 
 @time_measurement_longest(longest_threads_time_accumulator)
-def vector_vector_subtraction(b, Ax, start, end, residual) -> None:
-    residual[start:end] = [x - y for x, y in zip(b[start:end], Ax[start:end])]
+def matrix_vector_multiply(
+    matrix: Matrix, input_x: Vector, start: int, end: int, product: Vector
+) -> None:
+    """Return the matrix-vector product."""
+    product[start:end] = [
+        sum(x * y for x, y in zip(row, input_x, strict=False))
+        for row in matrix[start:end]
+    ]
 
 
 @time_measurement_longest(longest_threads_time_accumulator)
-def scalar_vector_multiply(omega, vector, start, end, result) -> None:
+def vector_vector_subtraction(
+    b: Vector, product: Vector, start: int, end: int, residual: Vector
+) -> None:
+    """Subtract one vector from another, elementwise."""
+    residual[start:end] = [
+        x - y for x, y in zip(b[start:end], product[start:end], strict=False)
+    ]
+
+
+@time_measurement_longest(longest_threads_time_accumulator)
+def scalar_vector_multiply(
+    omega: float, vector: Vector, start: int, end: int, result: Vector
+) -> None:
+    """Multiply every element of a vector by a scalar."""
     result[start:end] = [omega * x for x in vector[start:end]]
 
 
 @time_measurement_longest(longest_threads_time_accumulator)
-def vector_vector_addition(input_x, vector, start, end, output_x) -> None:
-    output_x[start:end] = [x + y for x, y in zip(input_x[start:end], vector[start:end])]
+def vector_vector_addition(
+    input_x: Vector, vector: Vector, start: int, end: int, output_x: Vector
+) -> None:
+    """Add two vectors elementwise."""
+    output_x[start:end] = [
+        x + y for x, y in zip(input_x[start:end], vector[start:end], strict=False)
+    ]
 
 
-def RichardsonMethodThreads(
-    A, b, lambda_min, lambda_max, max_iterations, x0=None, tol=1e-5
-):
+def _run_chunked(
+    executor: ThreadPoolExecutor,
+    layout: tuple[int, int, int],
+    func: Callable[..., None],
+    args: tuple[object, ...],
+) -> None:
+    """Split a row range into per-thread chunks, submit them, and wait.
+
+    ``layout`` is (num_threads, n, chunk_size). The four stages of one
+    Richardson iteration each did this by hand, which is what made the
+    iteration too complex and too long for ruff.
+    """
+    num_threads, n, chunk_size = layout
+    futures = []
+    for i in range(num_threads):
+        start = i * chunk_size
+        end = n if i == num_threads - 1 else (i + 1) * chunk_size
+        futures.append(executor.submit(func, *args[:-1], start, end, args[-1]))
+    for future in futures:
+        future.result()
+    longest_threads_time_accumulator.save_lap_and_reset()
+
+
+def richardson_method_threads(
+    problem: Problem, settings: Settings
+) -> tuple[Vector, str | None]:
+    """Run Richardson with one thread per chunk.
+
+    The unpacking below is deliberate: the body predates Problem/Settings and
+    reads these as locals throughout.
+    """
+    matrix = problem.matrix
+    b = problem.rhs
+    lambda_min = problem.lambda_min
+    lambda_max = problem.lambda_max
+    max_iterations = settings.max_iterations
+    x0 = settings.x0
+    tol = settings.tol
     longest_threads_time_accumulator.hard_reset()
 
     gc.disable()
@@ -48,89 +113,36 @@ def RichardsonMethodThreads(
     omega = 2 / (lambda_min + lambda_max)
     num_threads = multiprocessing.cpu_count()
     chunk_size = n // num_threads
+    layout = (num_threads, n, chunk_size)
 
     with ThreadPoolExecutor(
         max_workers=num_threads
     ) as executor:  # wątki są tworzone raz i nie są niszczone
-        for iteration in range(max_iterations):
-            Ax = [0] * len(
-                x
-            )  # tutaj zostanie przypisany wynik z mnożenia macierzy A z wektorem x
-            futures = []
-
-            for i in range(num_threads):
-                start = i * chunk_size
-                end = n if i == num_threads - 1 else (i + 1) * chunk_size
-                futures.append(
-                    executor.submit(matrix_vector_multiply, A, x, start, end, Ax)
-                )
-
-            for future in futures:
-                future.result()
-
-            longest_threads_time_accumulator.save_lap_and_reset()
-            residual = [0] * len(
-                b
-            )  # tutaj zostanie przypisany wynik z vector_vector_subtraction
-            futures = []
-
-            for i in range(num_threads):
-                start = i * chunk_size
-                end = n if i == num_threads - 1 else (i + 1) * chunk_size
-                futures.append(
-                    executor.submit(
-                        vector_vector_subtraction, b, Ax, start, end, residual
-                    )
-                )
-
-            for future in futures:
-                future.result()
-
-            longest_threads_time_accumulator.save_lap_and_reset()
-            change_vector = [0] * len(
-                residual
-            )  # zostanie tu przypisany wynik scalar_vector_multiply po pracy wątków
-            futures = []
-
-            for i in range(num_threads):
-                start = i * chunk_size
-                end = n if i == num_threads - 1 else (i + 1) * chunk_size
-                futures.append(
-                    executor.submit(
-                        scalar_vector_multiply,
-                        omega,
-                        residual,
-                        start,
-                        end,
-                        change_vector,
-                    )
-                )
-
-            for future in futures:
-                future.result()
-
-            longest_threads_time_accumulator.save_lap_and_reset()
-            _x = x[
-                :
-            ]  # do _x zostanie przez wątki przypisany wynik pracy w danej iteracji
-            futures = []
-
-            for i in range(num_threads):
-                start = i * chunk_size
-                end = n if i == num_threads - 1 else (i + 1) * chunk_size
-                futures.append(
-                    executor.submit(
-                        vector_vector_addition, x, change_vector, start, end, _x
-                    )
-                )
-
-            for future in futures:
-                future.result()
-
-            longest_threads_time_accumulator.save_lap_and_reset()
+        for _iteration in range(max_iterations):
+            # holds the result of multiplying A by x
+            product = [0] * len(x)
+            _run_chunked(executor, layout, matrix_vector_multiply, (matrix, x, product))
+            # holds the result of b - Ax
+            residual = [0] * len(b)
+            _run_chunked(
+                executor, layout, vector_vector_subtraction, (b, product, residual)
+            )
+            # holds omega * residual
+            change_vector = [0] * len(residual)
+            _run_chunked(
+                executor,
+                layout,
+                scalar_vector_multiply,
+                (omega, residual, change_vector),
+            )
+            # the threads write this iteration's result into _x
+            _x = x[:]
+            _run_chunked(
+                executor, layout, vector_vector_addition, (x, change_vector, _x)
+            )
             x = _x[:]
 
-            if linAlg.SequentialLinearAlgebraUtils.vector_norm(residual) < tol:
+            if lin_alg.SequentialLinearAlgebraUtils.vector_norm(residual) < tol:
                 break
 
     end_time = time.perf_counter()
diff --git a/Programming/PORR/code/threads_too_simple.py b/Programming/PORR/code/threads_too_simple.py
old mode 100644
new mode 100755
index dc20eaff..fbdd0835
--- a/Programming/PORR/code/threads_too_simple.py
+++ b/Programming/PORR/code/threads_too_simple.py
@@ -1,3 +1,6 @@
+#!/usr/bin/env python3
+"""Threaded Richardson, the first attempt (one thread per slice)."""
+
 import gc
 import logging
 import multiprocessing
@@ -6,6 +9,7 @@ import threading
 import time
 
 import numpy as np
+from richardson_problem import Problem, Settings
 from time_measurement import (
     longest_time_accumulator,
     tests_time,
@@ -14,17 +18,40 @@ from time_measurement import (
 
 logger = logging.getLogger(__name__)
 
+# Called with numpy arrays here and with plain lists in the sequential
+# backend; the aliases name both.
+Matrix = np.ndarray | list[list[float]]
+Vector = np.ndarray | list[float]
+
 
 @time_measurement_longest(longest_time_accumulator)
-def RichardsonThread(A, b, x, _x, omega, start, end) -> None:
+def richardson_thread(
+    matrix: Matrix,
+    b: Vector,
+    vectors: tuple[Vector, Vector],
+    omega: float,
+    span: tuple[int, int],
+) -> None:
+    """Run one thread's slice of a Richardson iteration.
+
+    ``vectors`` is (x, next_x) and ``span`` is the half-open row range this
+    thread owns.
+    """
+    x, _x = vectors
+    start, end = span
     for i in range(start, end):
-        sigma = np.dot(A[i, :], _x) - A[i, i] * _x[i]
-        x[i] = (1 - omega) * _x[i] + omega * (b[i] - sigma) / A[i, i]
+        sigma = np.dot(matrix[i, :], _x) - matrix[i, i] * _x[i]
+        x[i] = (1 - omega) * _x[i] + omega * (b[i] - sigma) / matrix[i, i]
 
 
-def RichardsonMethodThreads(
-    A, b, lambda_min, lambda_max, max_iterations, x0=None, tol=1e-5
-):
+def richardson_method_threads(
+    problem: Problem, settings: Settings
+) -> tuple[Vector, str | None]:
+    """Run Richardson with one thread per chunk."""
+    matrix = problem.matrix
+    b = problem.rhs
+    max_iterations = settings.max_iterations
+    x0 = settings.x0
     longest_time_accumulator.total_time = 0
     longest_time_accumulator.start = sys.float_info.max
     longest_time_accumulator.end = 0
@@ -44,21 +71,18 @@ def RichardsonMethodThreads(
     for _ in range(max_iterations):
         _x = x[:]
         for i in range(num_threads):
-            start = (
-                i * chunk_size
-            )  # start jest indeksem w A. Wątki otrzymują kolejny punkt startowy będący wielokrotnością rozmiaru porcji na wątek
+            # start is an index into A: each thread gets the next multiple of
+            # the per-thread chunk size as its starting point.
+            start = i * chunk_size
             end = n if i == num_threads - 1 else (i + 1) * chunk_size
             thread = threading.Thread(
-                target=RichardsonThread, args=(A, b, x, _x, omega, start, end)
+                target=richardson_thread, args=(matrix, b, (x, _x), omega, (start, end))
             )
             threads.append(thread)
             thread.start()
 
         for thread in threads:
             thread.join()
-
-        #         break
-
     end_time = time.perf_counter()
     gc.enable()
     total_time = end_time - start_time
diff --git a/Programming/PORR/code/time_measurement.py b/Programming/PORR/code/time_measurement.py
old mode 100644
new mode 100755
index 43fc8b5a..997daab4
--- a/Programming/PORR/code/time_measurement.py
+++ b/Programming/PORR/code/time_measurement.py
@@ -1,27 +1,39 @@
+#!/usr/bin/env python3
+"""Decorators that accumulate how long the measured calls take."""
+
 import sys
 import time
-from functools import wraps
+from collections.abc import Callable
 
 
 class TimeAccumulator:
+    """Total time spent inside the decorated functions."""
+
     def __init__(self) -> None:
+        """Set up the object."""
         self.total_time = 0
 
 
 class ComplexTimeAcumulator:
+    """Per-lap and longest times for the decorated functions."""
+
     def __init__(self) -> None:
+        """Set up the object."""
         self.hard_reset()
 
     def hard_reset(self) -> None:
+        """Clear every recorded time."""
         self.total_time = 0
         self.reset()
 
     def reset(self) -> None:
+        """Clear the current lap."""
         self.lap_time = 0
         self.start = sys.float_info.max
         self.end = 0
 
     def save_lap_and_reset(self) -> None:
+        """Store the current lap and start a new one."""
         self.total_time += self.lap_time
         self.reset()
 
@@ -32,10 +44,14 @@ tests_time = TimeAccumulator()
 longest_threads_time_accumulator = ComplexTimeAcumulator()
 
 
-def time_measurement(accumulator: TimeAccumulator):
-    def decorator(func):
-        @wraps(func)
-        def inner(*args, **kwargs):
+def time_measurement(accumulator: TimeAccumulator) -> Callable[..., object]:
+    """Decorate a function so its runtime accumulates."""
+
+    def decorator(func: Callable[..., object]) -> Callable[..., object]:
+        """Wrap the function so its runtime is recorded."""
+
+        def inner(*args: object, **kwargs: object) -> object:
+            """Call the wrapped function and record how long it took."""
             start = time.perf_counter()
             result = func(*args, **kwargs)
             end = time.perf_counter()
@@ -47,10 +63,16 @@ def time_measurement(accumulator: TimeAccumulator):
     return decorator
 
 
-def time_measurement_longest(accumulator: ComplexTimeAcumulator):
-    def decorator(func):
-        @wraps(func)
-        def inner(*args, **kwargs):
+def time_measurement_longest(
+    accumulator: ComplexTimeAcumulator,
+) -> Callable[..., object]:
+    """Decorate a function so its longest runtime is kept."""
+
+    def decorator(func: Callable[..., object]) -> Callable[..., object]:
+        """Wrap the function so its runtime is recorded."""
+
+        def inner(*args: object, **kwargs: object) -> object:
+            """Call the wrapped function and record how long it took."""
             start = time.perf_counter()
             result = func(*args, **kwargs)
             end = time.perf_counter()