diff --git a/.devcontainer/devcontainer.json b/.devcontainer/devcontainer.json
index 0cccda3..97a4b68 100644
--- a/.devcontainer/devcontainer.json
+++ b/.devcontainer/devcontainer.json
@@ -19,5 +19,5 @@
         }
     },
     "remoteUser": "vscode",
-    "postCreateCommand": "python -m pip install -r ./requirements.txt && python -m pip install -e ./[all] && pre-commit install && pre-commit run -a"
+    "postCreateCommand": "python -m pip install -e ./[all] && pre-commit install && pre-commit run -a"
 }
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 4fa1193..877ab0d 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -16,7 +16,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: "3.12"
+          python-version: "3.13"
           cache: "pip"
       - name: Install dependencies
         run: |
@@ -29,7 +29,7 @@ jobs:
     needs: pre-commit-ci
     strategy:
       matrix:
-        python-version: ["3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"]
     steps:
       - uses: actions/checkout@v4
       - name: Set up Python
@@ -37,17 +37,15 @@ jobs:
         with:
           python-version: ${{ matrix.python-version }}
           cache: "pip"
-      - name: Install dependencies
+      - name: Build
         run: |
           python -m pip install --upgrade pip
-          python -m pip install -r requirements.txt
+          python -m pip install -e .
       - name: Analyze code with pylint
         run: |
           python -m pip install pylint
           pylint $(git ls-files '*.py')
         continue-on-error: true
-      - name: Build
-        run: python -m pip install -e .
       - name: Test and coverage
         run: |
           python -m pip install pytest pytest-cov
@@ -69,7 +67,7 @@ jobs:
       - name: Set up Python
         uses: actions/setup-python@v5
         with:
-          python-version: "3.12"
+          python-version: "3.13"
           cache: "pip"
       - name: Install dependencies and tqec package
         run: |
diff --git a/.github/workflows/exotic-oses.yml b/.github/workflows/exotic-oses.yml
index 0ec70df..86bf02d 100644
--- a/.github/workflows/exotic-oses.yml
+++ b/.github/workflows/exotic-oses.yml
@@ -21,21 +21,19 @@ jobs:
         uses: actions/setup-python@v5
         with:
           # Use smallest supported Python version as the common denominator.
-          # In theory, if everything succeed in 3.9, everything should succeed in
-          # 3.10, 3.11 and follow-up version.
-          python-version: "3.9"
+          # In theory, if everything succeed in 3.10, everything should succeed in
+          # 3.11 and follow-up versions.
+          python-version: "3.10"
           cache: "pip"
-      - name: Install dependencies
+      - name: Build
         run: |
           python -m pip install --upgrade pip
-          python -m pip install -r requirements.txt
+          python -m pip install '.[all]'
       - name: Analyze code with pylint
         run: |
           python -m pip install pylint
           pylint $(git ls-files '*.py')
         continue-on-error: true
-      - name: Build
-        run: python -m pip install '.[all]'
       - name: Test and coverage
         run: |
           python -m pip install pytest pytest-cov
diff --git a/.github/workflows/gh-pages.yml b/.github/workflows/gh-pages.yml
index c4b9384..751b4cd 100644
--- a/.github/workflows/gh-pages.yml
+++ b/.github/workflows/gh-pages.yml
@@ -26,7 +26,6 @@ jobs:
           sudo apt-get update
           sudo apt-get install -y --no-install-recommends pandoc
           python -m pip install --upgrade pip
-          python -m pip install -r requirements.txt
           python -m pip install '.[all]'
       - name: Generate the documentation
         run: |
diff --git a/README.md b/README.md
index 72dfeb6..2d9ec3e 100644
--- a/README.md
+++ b/README.md
@@ -16,8 +16,6 @@ Currently, `tqecd` needs to be installed from source using
 python -m pip install git+https://github.com/tqec/tqecd.git
 ```
 
-The `tqecd` package has some dependencies that might be harder to install than a simple `pip install`. If you have any issues with the simple installation method above, please look at the [full installation page](https://tqec.github.io/tqecd/user_guide/installation.html).
-
 ## Basic usage
 
 ```py
diff --git a/docs/user_guide/installation.rst b/docs/user_guide/installation.rst
index 63a5faf..d7a6b68 100644
--- a/docs/user_guide/installation.rst
+++ b/docs/user_guide/installation.rst
@@ -10,21 +10,6 @@ Python version
 The ``tqecd`` package only supports Python 3.10 and onward. If you have Python 3.9 or below,
 please update your Python installation.
 
-Additional toolchains
-~~~~~~~~~~~~~~~~~~~~~
-
-Some of the dependencies of ``tqecd`` are implemented using compiled languages. This is for
-example the case of the `pycryptosat <https://pypi.org/project/pycryptosat/>`_ dependency.
-Pre-compiled Python packages that should be compatible with any GNU/Linux are provided
-by the author, but no pre-compiled package exist for Windows or MacOS.
-
-This means that, if you try to install the ``tqecd`` package on Windows or MacOS, a working
-C++ toolchain should also be installed on your system.
-
-Here is a list of potential issues you might encounter and how to solve them:
-
-- `Failed building wheel for pycryptosat <https://github.com/tqec/tqec/issues/311>`_
-
 Installation procedure
 ----------------------
 
diff --git a/pyproject.toml b/pyproject.toml
index 0dd9b2a..312d148 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "tqecd"
-version = "0.1.3"
+version = "0.1.4"
 authors = [
     { name = "TQEC community", email = "tqec-design-automation@googlegroups.com" },
 ]
@@ -12,7 +12,7 @@ description = "Automatically find detectors in a topologically quantum error cor
 readme = "README.md"
 license = { file = "LICENSE" }
 keywords = ["topological quantum error correction", "qec", "detectors"]
-requires-python = ">= 3.9"
+requires-python = ">= 3.10"
 classifiers = [
     "Development Status :: 4 - Beta",
     "Intended Audience :: Developers",
@@ -28,12 +28,19 @@ classifiers = [
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
+    "Programming Language :: Python :: 3.14",
     "Topic :: Scientific/Engineering",
     "Topic :: Software Development :: Libraries",
     "Topic :: Utilities",
     "Typing :: Typed",
 ]
-dynamic = ["dependencies"]
+dependencies = [
+    "numpy>=1.22,<2.3",                      # Upper bound: https://github.com/tqec/tqec/pull/659
+    "numpy>=1.24; python_version >= '3.11'",
+    "numpy>=1.26; python_version >= '3.12'",
+    "numpy>=2.1; python_version >= '3.13'",
+    "stim>=1.15.0",
+]
 
 [project.urls]
 Documentation = "https://tqec.github.io/tqecd/"
@@ -64,10 +71,6 @@ exclude = ["*test.py"]
 [tool.setuptools.package-data]
 tqecd = ["py.typed"]
 
-[tool.setuptools.dynamic]
-dependencies = { file = "requirements.txt" }
-
-
 [tool.mypy]
 # See https://numpy.org/devdocs/reference/typing.html
 plugins = "numpy.typing.mypy_plugin"
@@ -106,7 +109,7 @@ warn_return_any = true
 
 # See https://mypy.readthedocs.io/en/stable/config_file.html#using-a-pyproject-toml
 [[tool.mypy.overrides]]
-module = ["stim", "pysat", "pysat.solvers"]
+module = ["stim"]
 ignore_missing_imports = true
 
 [tool.coverage.run]
diff --git a/requirements.txt b/requirements.txt
deleted file mode 100644
index b51c8ed..0000000
--- a/requirements.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-numpy
-stim
-pycryptosat!=5.11.23
-python-sat[cryptosat]
diff --git a/setup.py b/setup.py
deleted file mode 100644
index 6068493..0000000
--- a/setup.py
+++ /dev/null
@@ -1,3 +0,0 @@
-from setuptools import setup
-
-setup()
diff --git a/src/tqecd/cover.py b/src/tqecd/cover.py
new file mode 100644
index 0000000..95c2a07
--- /dev/null
+++ b/src/tqecd/cover.py
@@ -0,0 +1,150 @@
+"""Provides utility functions to find "covers" of Pauli strings."""
+
+from __future__ import annotations
+
+from collections.abc import Callable, Iterable
+from typing import Any
+
+from tqecd.pauli import PauliString
+
+
+def _find_cover(
+    target: PauliString, sources: list[PauliString], on_qubits: frozenset[int]
+) -> list[int] | None:
+    """Try to find a set of boundary stabilizers from ``sources`` that generate
+    target on qubits ``on_qubits``.
+
+    If multiple valid covers exist, the covers involving the lowest number of
+    :class:`~tqecd.pauli.PauliString` instances from ``sources`` are listed, and
+    a random cover is picked from that list.
+
+    Args:
+        target: the stabilizers to cover with stabilizers from ``sources``.
+        sources: stabilizers that can be used to cover `target`.
+        on_qubits: qubits to consider when trying to cover ``target`` with
+            ``sources``.
+
+    Returns:
+        Either a list of indices over ``sources`` that, when combined, cover
+        exactly the provided ``target`` on all the qubits provided in
+        ``on_qubits``, or ``None`` if such a list could not be found.
+    """
+    return _solve_linear_system(
+        _construct_basis({}, sources, lambda s: s.to_int(on_qubits)),
+        target.to_int(on_qubits),
+    )
+
+
+def find_exact_cover(
+    target: PauliString, sources: list[PauliString]
+) -> list[int] | None:
+    """Try to find a set of pauli strings from ``sources`` that generate exactly
+    ``target``.
+
+    The Pauli strings returned (via indices over the provided ``sources``), once
+    multiplied together, should be exactly equal to ``target``. In particular, the
+    following post-condition should hold:
+
+    .. code-block:: python
+
+        target = None     # to replace
+        sources = [None]  # to replace
+        cover_indices = find_exact_cover(target, sources)
+        resulting_pauli_string = PauliString({})
+        for i in cover_indices:
+            resulting_pauli_string = resulting_pauli_string * sources[i]
+        assert resulting_pauli_string == target, "Should hold"
+
+    Args:
+        target: the stabilizers to cover with stabilizers from ``sources``.
+        sources: stabilizers that can be used to cover ``target``.
+
+    Returns:
+        Either a list of indices over ``sources`` that, when combined, cover
+        exactly the provided ``target``, or ``None`` if such a list could not be
+        found.
+    """
+    # If target is the identity, pick no Pauli string.
+    # Note: we might want to disallow an empty return in the future.
+    if target.non_trivial_pauli_count == 0:
+        return []
+
+    # Else, if there are no sources, we cannot find a solution.
+    if not sources:
+        return None
+
+    # We want an exact (i.e., equality) cover on all qubits, to be sure that
+    # the post-condition in the docstring holds. For that, it is sufficient to
+    # only consider all the qubits where either `target` or at least one of the
+    # items of `sources` acts non-trivially (i.e., something else than the identity).
+    involved_qubits = frozenset(target.qubits)
+    for source in sources:
+        involved_qubits |= frozenset(source.qubits)
+
+    return _find_cover(target, sources, involved_qubits)
+
+
+def find_commuting_cover_on_target_qubits(
+    target: PauliString, sources: list[PauliString]
+) -> list[int] | None:
+    """Try to find a set of boundary stabilizers from ``sources`` that generate a
+    superset of ``target``.
+
+    This function try to find a set of Pauli strings from ``sources`` that
+    includes ``target`` (i.e., on every qubit where ``target`` is non-trivial,
+    the product of each of the returned Pauli strings should commute with
+    ``target``).
+
+    The differences with :func:`find_cover` are:
+
+    1. this function does not restrict the output of the product of each of
+       the returned Pauli string on qubits where ``target`` acts trivially (i.e.
+       "I"). So in practice, on qubits where ``target[qubit] == "I"``, the value
+       of the returned Pauli string can be anything.
+    2. this function does not restrict the output of the product of each of
+       the returned Pauli string to exactly match ``target`` on qubits where
+       it acts non-trivially, but rather requires the output to commute with
+       ``target`` on those qubits.
+
+    Args:
+        target: the stabilizers to cover with stabilizers from ``sources``.
+        sources: stabilizers that can be used to cover ``target``.
+
+    Returns:
+        Either a list of a stabilizers that, when combined, commute with
+        the provided ``target``, or ``None`` if such a list could not be found.
+    """
+    if not sources:
+        return None
+    return _find_cover(target, sources, frozenset(target.qubits))
+
+
+def _solve_linear_system(
+    basis: dict[int, tuple[int, int]], x: int, update_basis: bool = True
+) -> list[int] | None:
+    """Gaussian elimination over GF(2)."""
+    mask = 1 << len(basis)
+    while x:
+        highest_bit = x.bit_length() - 1
+        if highest_bit not in basis:
+            if update_basis:
+                basis[highest_bit] = (x, mask)
+            return None
+        pivot, pivot_mask = basis[highest_bit]
+        x ^= pivot
+        mask ^= pivot_mask
+    return _int_to_bit_indices(mask)[:-1]
+
+
+def _int_to_bit_indices(x: int) -> list[int]:
+    """Convert an integer to a list of indices where the bits are set."""
+    return [i for i in range(x.bit_length()) if (x >> i) & 1]
+
+
+def _construct_basis(
+    basis: dict[int, tuple[int, int]], items: Iterable[Any], func: Callable[[Any], int]
+) -> dict[int, tuple[int, int]]:
+    """Construct a linear basis from the given items using the provided function."""
+    for item in items:
+        _solve_linear_system(basis, func(item))
+    return basis
diff --git a/src/tqecd/match_utils/cover_test.py b/src/tqecd/cover_test.py
similarity index 64%
rename from src/tqecd/match_utils/cover_test.py
rename to src/tqecd/cover_test.py
index f2dcdf2..742370d 100644
--- a/src/tqecd/match_utils/cover_test.py
+++ b/src/tqecd/cover_test.py
@@ -1,15 +1,11 @@
 from __future__ import annotations
 
-import itertools
-from typing import Iterator
-
 import pytest
 import stim
 
-from tqecd.match_utils.cover import (
-    _all_pauli_string_combination_results,  # pyright: ignore[reportPrivateUsage]
-    find_commuting_cover_on_target_qubits_sat,
-    find_exact_cover_sat,
+from tqecd.cover import (
+    find_commuting_cover_on_target_qubits,
+    find_exact_cover,
 )
 from tqecd.pauli import PauliString, pauli_product
 
@@ -18,48 +14,6 @@ def _pss(pauli_string: str) -> PauliString:
     return PauliString.from_stim_pauli_string(stim.PauliString(pauli_string))
 
 
-def _all_pauli_string_combination_results_alternative_implementation(
-    pauli_string_list: list[PauliString],
-) -> Iterator[PauliString]:
-    yield from (
-        pauli_product([p for p, c in zip(pauli_string_list, choices) if c])
-        for choices in itertools.product([True, False], repeat=len(pauli_string_list))
-    )
-
-
-def test_all_pauli_string_combination_results_randomized() -> None:
-    paulis = [
-        PauliString.from_stim_pauli_string(stim.PauliString.random(10))
-        for _ in range(5)
-    ]
-
-    from_function = {
-        pauli for _, pauli in _all_pauli_string_combination_results(paulis)
-    }
-    from_alternative_implementation = set(
-        _all_pauli_string_combination_results_alternative_implementation(paulis)
-    )
-    assert from_function == from_alternative_implementation
-
-
-def test_all_pauli_string_combination_results() -> None:
-    paulis = [PauliString({i: "Z"}) for i in range(3)]
-    all_combinations = list(
-        pauli for _, pauli in _all_pauli_string_combination_results(paulis)
-    )
-    assert len(all_combinations) == 2 ** len(paulis)
-    assert set(all_combinations) == {
-        PauliString({}),
-        PauliString({0: "Z"}),
-        PauliString({1: "Z"}),
-        PauliString({2: "Z"}),
-        PauliString({0: "Z", 1: "Z"}),
-        PauliString({0: "Z", 2: "Z"}),
-        PauliString({1: "Z", 2: "Z"}),
-        PauliString({0: "Z", 1: "Z", 2: "Z"}),
-    }
-
-
 @pytest.mark.parametrize(
     "target,sources,expected_result",
     [
@@ -117,7 +71,7 @@ def test_all_pauli_string_combination_results() -> None:
 def test_exact_match(
     target: PauliString, sources: list[PauliString], expected_result: list[int] | None
 ) -> None:
-    obtained_result = find_exact_cover_sat(target, sources)
+    obtained_result = find_exact_cover(target, sources)
     # We expect the results to either both be None, or both be a list.
     assert (obtained_result is None) == (expected_result is None)
     # If they are both a list, compare them and check the post-condition documented.
@@ -154,7 +108,7 @@ def test_exact_match(
 def test_commuting_match(
     target: PauliString, sources: list[PauliString], expected_result: list[int] | None
 ) -> None:
-    obtained_result = find_commuting_cover_on_target_qubits_sat(target, sources)
+    obtained_result = find_commuting_cover_on_target_qubits(target, sources)
     # We expect the results to either both be None, or both be a list.
     assert (obtained_result is None) == (expected_result is None)
     # If they are both a list, compare them and check the post-condition documented.
diff --git a/src/tqecd/flow.py b/src/tqecd/flow.py
index a06e215..e8c5a46 100644
--- a/src/tqecd/flow.py
+++ b/src/tqecd/flow.py
@@ -4,9 +4,9 @@
 from dataclasses import dataclass
 
 from tqecd.boundary import BoundaryStabilizer
+from tqecd.cover import find_commuting_cover_on_target_qubits
 from tqecd.exceptions import TQECDException
 from tqecd.fragment import Fragment, FragmentLoop
-from tqecd.match_utils.cover import find_commuting_cover_on_target_qubits_sat
 from tqecd.measurement import get_relative_measurement_index
 from tqecd.pauli import PauliString, pauli_product
 
@@ -73,7 +73,7 @@ def _try_merge_anticommuting_flows_inplace(flows: list[BoundaryStabilizer]) -> N
         flows[fi].before_collapse for fi in anti_commuting_index_to_flows_index
     ]
     indices_of_anti_commuting_stabilizers_to_merge = (
-        find_commuting_cover_on_target_qubits_sat(
+        find_commuting_cover_on_target_qubits(
             collapsing_pauli, anticommuting_stabilizers
         )
     )
@@ -105,7 +105,7 @@ def _try_merge_anticommuting_flows_inplace(flows: list[BoundaryStabilizer]) -> N
         flows.append(new_commuting_stabilizer)
         # Update for loop condition
         indices_of_anti_commuting_stabilizers_to_merge = (
-            find_commuting_cover_on_target_qubits_sat(
+            find_commuting_cover_on_target_qubits(
                 collapsing_pauli, anticommuting_stabilizers
             )
         )
diff --git a/src/tqecd/match.py b/src/tqecd/match.py
index b6db721..753cf41 100644
--- a/src/tqecd/match.py
+++ b/src/tqecd/match.py
@@ -8,9 +8,9 @@
 import stim
 
 from tqecd.boundary import BoundaryStabilizer
+from tqecd.cover import find_exact_cover
 from tqecd.exceptions import TQECDException
 from tqecd.flow import FragmentFlows, FragmentLoopFlows
-from tqecd.match_utils.cover import find_exact_cover_sat
 from tqecd.measurement import RelativeMeasurementLocation
 
 
@@ -457,7 +457,7 @@ def _match_boundary_stabilizers_by_disjoint_cover(
         # Try to find the cover. If unsuccessful, skip `target` and continue the loop.
         # Else, a detector is matched, so insert it into the returned list and mark
         # `target` as "to remove".
-        cover_indices = find_exact_cover_sat(
+        cover_indices = find_exact_cover(
             target.after_collapse, [cs.after_collapse for cs in covering_stabilizers]
         )
         if cover_indices is None:
diff --git a/src/tqecd/match_utils/__init__.py b/src/tqecd/match_utils/__init__.py
deleted file mode 100644
index e2a7e00..0000000
--- a/src/tqecd/match_utils/__init__.py
+++ /dev/null
@@ -1,27 +0,0 @@
-"""Provides a few utility functions to match flows.
-
-This module provides functions to encode the problem of finding matching flows
-as a SAT problem, as well as a function to solve that SAT problem.
-
-The three functions re-exported by this module all aims at finding a list of
-Pauli strings from a given list that "cover" a target Pauli string:
-
-- :func:`~.cover.find_cover` performs a brute-force search to find Pauli strings
-  from a provided candidate list that, once multiplied together, exactly match
-  the provided target Pauli string,
-- :func:`~.cover.find_exact_cover_sat` performs the same task as
-  :func:`~.cover.find_cover` but encodes the problem as a SAT problem and uses a
-  SAT solver to find potential solutions instead of performing a brute-force
-  search.
-- :func:`~.cover.find_commuting_cover_on_target_qubits_sat` uses the same
-  strategy as :func:`~.cover.find_exact_cover_sat` by using a SAT solver, but
-  solves a slightly different problem that consists in finding Pauli strings
-  that, once multiplied together, commute with a provided target Pauli string
-  (instead of being exactly equal to it).
-"""
-
-from .cover import (
-    find_commuting_cover_on_target_qubits_sat as find_commuting_cover_on_target_qubits_sat,
-)
-from .cover import find_cover as find_cover
-from .cover import find_exact_cover_sat as find_exact_cover_sat
diff --git a/src/tqecd/match_utils/cover.py b/src/tqecd/match_utils/cover.py
deleted file mode 100644
index 6a59a15..0000000
--- a/src/tqecd/match_utils/cover.py
+++ /dev/null
@@ -1,319 +0,0 @@
-"""Provides utility functions to find "covers" of Pauli strings."""
-
-from __future__ import annotations
-
-import time
-from typing import Iterator
-
-import pysat.solvers
-
-from tqecd.boundary import BoundaryStabilizer, manhattan_distance
-from tqecd.match_utils.sat import (
-    encode_pauli_string_commuting_cover_sat_problem_in_solver,
-    encode_pauli_string_exact_cover_sat_problem_in_solver,
-)
-from tqecd.pauli import PauliString
-
-
-def _all_pauli_string_combination_results(
-    pauli_strings: list[PauliString],
-) -> Iterator[tuple[list[bool], PauliString]]:
-    """Iterate over all the possible Pauli string products.
-
-    This function iterates over all the ``2**len(pauli_strings)`` products
-    that may be generated by either picking or not picking any of the
-    provided Pauli strings.
-
-    It is efficient in the sense that it performs the minimum number of
-    Pauli products required to do so.
-
-    Args:
-        pauli_strings: a list of Pauli strings that will be considered
-            in the returned products.
-
-    Yields:
-        all the ``2**len(pauli_strings)`` possible Pauli products, one by one.
-    """
-    yield from _all_pauli_string_combination_results_impl(
-        pauli_strings, PauliString({})
-    )
-
-
-def _all_pauli_string_combination_results_impl(
-    pauli_strings: list[PauliString], current_pauli_string: PauliString
-) -> Iterator[tuple[list[bool], PauliString]]:
-    """Iterate over all the possible Pauli string products.
-
-    This function iterates over all the ``2**len(pauli_strings)`` products
-    that may be generated by either picking or not picking any of the
-    provided Pauli strings.
-
-    It is efficient in the sense that it performs the minimum number of
-    Pauli products required to do so.
-
-    Args:
-        pauli_strings: a list of Pauli strings that will be considered
-            in the returned products.
-
-    Yields:
-        tuples containing a list of boolean of the exact same size as the
-        provided ``pauli_strings`` list and representing the Pauli strings
-        that have been considered in the product returned as the second
-        entry of the tuple.
-    """
-    if len(pauli_strings) == 0:
-        yield [], current_pauli_string
-        return
-    yield from (
-        ([True] + choices, res_pauli_string)
-        for choices, res_pauli_string in _all_pauli_string_combination_results_impl(
-            pauli_strings[1:],
-            current_pauli_string * pauli_strings[0],
-        )
-    )
-    yield from (
-        ([False] + choices, res_pauli_string)
-        for choices, res_pauli_string in _all_pauli_string_combination_results_impl(
-            pauli_strings[1:], current_pauli_string
-        )
-    )
-
-
-def find_cover(
-    target: BoundaryStabilizer,
-    sources: list[BoundaryStabilizer],
-    qubit_coordinates: dict[int, tuple[float, ...]],
-    maximum_qubit_distance: int = 5,
-) -> list[BoundaryStabilizer] | None:
-    """Try to cover the provided ``target`` stabilizer with stabilizers from
-    ``sources``.
-
-    This function is currently performing a bruteforce search: it tries all
-    combinations of stabilizers from ``sources`` and check if one matches with
-    the provided ``target``.
-    This approach blows up very quickly, as it may test up to ``2**len(sources)``
-    different combinations before being able to tell that no match has been found.
-
-    Note that **all** the combinations are tested, even the only combination that
-    consist in "picking nothing from the list", resulting in an empty PauliString.
-    To avoid any surprise, and because we do not expect an empty ``target`` to make
-    sense here, this function will raise on such a case.
-
-    Args:
-        target: the boundary stabilizer to try to match with the provided
-            ``sources``.
-        sources: the boundary stabilizers that this function will try to combine to
-            find a stabilizer involving ``target``.
-        qubit_coordinates: a mapping from qubit indices to coordinates. Used to
-            annotate the matched detectors with the coordinates from the qubits
-            involved in the measurement forming the detector.
-        maximum_qubit_distance: radius (in number of qubits) to consider when
-            searching for covering boundary stabilizers. Any boundary stabilizer
-            in ``sources`` that has coordinates outside of that radius from
-            ``target`` will not be considered, reducing the overall complexity
-            of this function. The radius is computed with the Manhattan distance.
-
-    Raises:
-        TQECDException: if ``target`` is the empty PauliString.
-        TQECDException: if any of the provided instances in ``target`` or
-            ``sources`` has anti-commuting stabilizers.
-
-    Returns:
-        a matching set of boundary stabilizers or None if no matching stabilizers
-        could be found.
-    """
-
-    sources = [
-        s
-        for s in sources
-        if manhattan_distance(target, s, qubit_coordinates) <= maximum_qubit_distance
-    ]
-
-    after_collapse_sources = [s.after_collapse for s in sources]
-    for (
-        picked_stabilizers,
-        resulting_pauli_string,
-    ) in _all_pauli_string_combination_results(after_collapse_sources):
-        if target.after_collapse == resulting_pauli_string:
-            return [
-                boundary_stabilizer
-                for i, boundary_stabilizer in enumerate(sources)
-                if picked_stabilizers[i]
-            ]
-    return None
-
-
-def _smallest_solution_shortcircuit(
-    solutions: Iterator[list[int]], lower_length_bound: int = 0, timeout: float = 0.1
-) -> list[int] | None:
-    """Iterate over the provided ``solutions`` iterator to find the smallest
-    possible solution within the provided ``timeout``.
-
-    Args:
-        solutions: iterator yielding lists of integers representing the indices of
-            Pauli strings that should be used to cover a specific target.
-        lower_length_bound: The minimum expected length. This parameter is used to
-            quit early if a solution with this length is found. Default to exploring
-            all the solutions to find the best one.
-        timeout: maximum time in seconds that can be spent finding the smallest
-            solution. The provided ``solutions`` iterator might yield millions of
-            inputs, so this parameter is here to stop the search and return the best
-            found solution when more than ``timeout`` seconds have been spent within
-            this function.
-
-    Returns:
-        the smallest solution found, or ``None`` if no solution was found.
-    """
-    start_time: float = time.monotonic()
-    smallest_solution = next(solutions, None)
-
-    if smallest_solution is None:
-        return None
-    if len(smallest_solution) == lower_length_bound:
-        return smallest_solution
-
-    for solution in solutions:
-        if len(solution) == lower_length_bound:
-            return solution
-        smallest_solution = min((smallest_solution, solution), key=len)
-        if time.monotonic() - start_time > timeout:
-            return smallest_solution
-    return smallest_solution
-
-
-def _find_cover_sat(
-    target: PauliString, sources: list[PauliString], on_qubits: frozenset[int]
-) -> list[int] | None:
-    """Try to find a set of boundary stabilizers from ``sources`` that generate
-    target on qubits ``on_qubits``.
-
-    If multiple valid covers exist, the covers involving the lowest number of
-    :class:`~tqecd.pauli.PauliString` instances from ``sources`` are listed, and
-    a random cover is picked from that list.
-
-    Args:
-        target: the stabilizers to cover with stabilizers from ``sources``.
-        sources: stabilizers that can be used to cover `target`.
-        on_qubits: qubits to consider when trying to cover ``target`` with
-            ``sources``.
-
-    Returns:
-        Either a list of indices over ``sources`` that, when combined, cover
-        exactly the provided ``target`` on all the qubits provided in
-        ``on_qubits``, or ``None`` if such a list could not be found.
-    """
-    with pysat.solvers.CryptoMinisat() as solver:
-        encode_pauli_string_exact_cover_sat_problem_in_solver(
-            solver, target, sources, on_qubits
-        )
-        return _smallest_solution_shortcircuit(
-            (
-                [i - 1 for i in satisfying_proof if i > 0]
-                for satisfying_proof in solver.enum_models()
-            ),
-            2,
-        )
-
-
-def find_exact_cover_sat(
-    target: PauliString, sources: list[PauliString]
-) -> list[int] | None:
-    """Try to find a set of pauli strings from ``sources`` that generate exactly
-    ``target``.
-
-    The Pauli strings returned (via indices over the provided ``sources``), once
-    multiplied together, should be exactly equal to ``target``. In particular, the
-    following post-condition should hold:
-
-    If multiple valid covers exist, the covers involving the lowest number of
-    :class:`~tqecd.pauli.PauliString` instances from ``sources`` are listed, and
-    a random cover is picked from that list.
-
-    .. code-block:: python
-
-        target = None     # to replace
-        sources = [None]  # to replace
-        cover_indices = find_exact_cover_sat(target, sources)
-        resulting_pauli_string = PauliString({})
-        for i in cover_indices:
-            resulting_pauli_string = resulting_pauli_string * sources[i]
-        assert resulting_pauli_string == target, "Should hold"
-
-    Args:
-        target: the stabilizers to cover with stabilizers from ``sources``.
-        sources: stabilizers that can be used to cover ``target``.
-
-    Returns:
-        Either a list of indices over ``sources`` that, when combined, cover
-        exactly the provided ``target``, or ``None`` if such a list could not be
-        found.
-    """
-    # If target is the identity, we do not have to call a SAT solver to find
-    # a solution: pick no Pauli string.
-    # Note: we might want to disallow an empty return in the future.
-    if target.non_trivial_pauli_count == 0:
-        return []
-
-    # Else, if there are no sources, we cannot find a solution.
-    if not sources:
-        return None
-
-    # We want an exact (i.e., equality) cover on all qubits, to be sure that
-    # the post-condition in the docstring holds. For that, it is sufficient to
-    # only consider all the qubits where either `target` or at least one of the
-    # items of `sources` acts non-trivially (i.e., something else than the identity).
-    involved_qubits = frozenset(target.qubits)
-    for source in sources:
-        involved_qubits |= frozenset(source.qubits)
-
-    return _find_cover_sat(target, sources, involved_qubits)
-
-
-def find_commuting_cover_on_target_qubits_sat(
-    target: PauliString, sources: list[PauliString]
-) -> list[int] | None:
-    """Try to find a set of boundary stabilizers from ``sources`` that generate a
-    superset of ``target``.
-
-    This function try to find a set of Pauli strings from ``sources`` that
-    includes ``target`` (i.e., on every qubit where ``target`` is non-trivial,
-    the product of each of the returned Pauli strings should commute with
-    ``target``).
-
-    If multiple valid covers exist, the covers involving the lowest number of
-    :class:`~tqecd.pauli.PauliString` instances from ``sources`` are listed, and
-    a random cover is picked from that list.
-
-    The differences with :func:`find_cover_sat` are:
-
-    1. this function does not restrict the output of the product of each of
-       the returned Pauli string on qubits where ``target`` acts trivially (i.e.
-       "I"). So in practice, on qubits where ``target[qubit] == "I"``, the value
-       of the returned Pauli string can be anything.
-    2. this function does not restrict the output of the product of each of
-       the returned Pauli string to exactly match ``target`` on qubits where
-       it acts non-trivially, but rather requires the output to commute with
-       ``target`` on those qubits.
-
-    Args:
-        target: the stabilizers to cover with stabilizers from ``sources``.
-        sources: stabilizers that can be used to cover ``target``.
-
-    Returns:
-        Either a list of a stabilizers that, when combined, commute with
-        the provided ``target``, or ``None`` if such a list could not be found.
-    """
-    if not sources:
-        return None
-
-    with pysat.solvers.CryptoMinisat() as solver:
-        encode_pauli_string_commuting_cover_sat_problem_in_solver(
-            solver, target, sources, frozenset(target.qubits)
-        )
-        return _smallest_solution_shortcircuit(
-            (
-                [i - 1 for i in satisfying_proof if i > 0]
-                for satisfying_proof in solver.enum_models()
-            ),
-            2,
-        )
diff --git a/src/tqecd/match_utils/sat.py b/src/tqecd/match_utils/sat.py
deleted file mode 100644
index ba03c7e..0000000
--- a/src/tqecd/match_utils/sat.py
+++ /dev/null
@@ -1,238 +0,0 @@
-"""Semi-private module providing functions to construct SAT problems.
-
-This module is semi-private in the sense that even though the two functions it
-implements are not private, they are quite low-level and are not expected to be
-used by external users of the library.
-
-The module defines 2 functions that each take a SAT solver able to encode XOR
-clauses, a target Pauli string, a list of potential Pauli strings and a set of
-qubits to consider and encode in the provided SAT solver a cover problem, either
-exact (i.e., find Pauli strings in the list of potential Pauli strings such that
-their product is exactly the target Pauli string on the considered qubits) or
-commuting (i.e., the product only needs to commute with the target Pauli string
-on the considered qubits).
-"""
-
-from __future__ import annotations
-
-from pysat.solvers import CryptoMinisat
-
-from tqecd.pauli import PauliString, pauli_literal_to_bools
-
-
-def encode_pauli_string_exact_cover_sat_problem_in_solver(
-    solver: CryptoMinisat,
-    expected_pauli_string: PauliString,
-    available_pauli_strings: list[PauliString],
-    qubits_to_consider: frozenset[int],
-) -> None:
-    """Build the SAT problem that should be solved to find an exact cover.
-
-    This function encodes the SAT problem of interest into the provided ``solver``.
-    As such, the provided ``solver`` is expected to be a newly created instance,
-    and will be mutated by this function.
-
-    The encoded problem is the following:
-
-    Find the indices of Pauli strings in ``available_pauli_strings`` such that
-    the product of all the corresponding Pauli strings is exactly equal to
-    ``expected_pauli_string`` on all the qubits in ``qubits_to_consider``.
-
-    The following post-condition should be checked after solving:
-
-    .. code-block:: python
-
-        expected_pauli_string = PauliString({})      # Fill in!
-        available_pauli_strings = [PauliString({})]  # Fill in!
-        qubits_to_consider = frozenset()             # Fill in!
-        solver = None                                # Fill in!
-
-        encode_pauli_string_exact_cover_sat_problem_in_solver(
-            solver, expected_pauli_string, available_pauli_strings, qubits_to_consider
-        )
-        indices = None # Solve the problem encoded in solver and recover the indices.
-
-        final_pauli_string = pauli_product([available_pauli_strings[i] for i in indices])
-        for q in qubits_to_consider:
-            assert final_pauli_string[q] == expected_pauli_string[q]
-
-
-    The SAT problem contains one boolean variable per Pauli string in
-    ``available_pauli_strings``, deciding if the associated Pauli string should
-    be part of the cover.
-
-    The formula is quite simple: for each qubit in ``qubits_to_consider``, the
-    following two boolean formulas should be ``True``:
-
-    1. The X part of the resulting Pauli string, obtained from XORing the X part
-       of each Pauli string weighted by the boolean variable that decides if the
-       Pauli string should be included or not, should be equal to the X part of
-       the expected final Pauli string.
-    2. Same reasoning, but for the Z part.
-
-    Note that the "should be equal to the [P] part of the expected final Pauli
-    string" is natively implemented in the SAT solver used, so we do not need to
-    adapt the formula with some binary logic tricks. The solver really solves a
-    conjunction of ``a XOR b XOR ... XOR z == [True/False]``.
-
-    This leads to a formula with ``2*len(qubits_to_consider)`` XOR clauses that
-    should all be verified for a cover to exist.
-    Each of the ``2*len(qubits_to_consider)`` XOR clauses can contain up to
-    ``len(available_pauli_strings)`` XORed items, but a simplification is made
-    that should reduce that number: if the Pauli string acts trivially w.r.t the
-    considered Pauli effect (X or Z), the boolean variable is removed from the
-    formula as we know for sure that this particular Pauli string cannot impact
-    the result on that specific qubit and Pauli effect.
-
-    Args:
-        solver: solver that will be modified in-place to encode the SAT problem.
-        expected_pauli_string: target Pauli string that should be covered by
-            strings from ``available_pauli_strings``.
-        available_pauli_strings: Pauli strings that can be used to try to cover
-            ``expected_pauli_string``.
-        qubits_to_consider: qubits on which the cover should be exactly equal to
-            ``expected_pauli_string``. Qubits not listed in this input will
-            simply be ignored and no restriction on the value of the resulting
-            cover on those qubits is added.
-    """
-
-    for qubit in qubits_to_consider:
-        expected_X, expected_Z = pauli_literal_to_bools(expected_pauli_string[qubit])
-        available_paulis_bools: list[tuple[bool, bool]] = [
-            pauli_literal_to_bools(pauli[qubit]) for pauli in available_pauli_strings
-        ]
-        # The two following lists includes the 1-based indices of Pauli strings from
-        # the provided `available_pauli_strings` input that finish with a X/Z stabilizers
-        # on the current `qubit`.
-        # For each index, we want to know if it should be included in the final cover or
-        # not, so XORing the boolean variables representing whether or not this Pauli
-        # string should be included results in the stabilizer propagated on that qubit.
-        x_clause_literals = [
-            pi + 1 for pi, (px, _) in enumerate(available_paulis_bools) if px
-        ]
-        z_clause_literals = [
-            pi + 1 for pi, (_, pz) in enumerate(available_paulis_bools) if pz
-        ]
-        solver.add_xor_clause(x_clause_literals, value=expected_X)
-        solver.add_xor_clause(z_clause_literals, value=expected_Z)
-
-
-def encode_pauli_string_commuting_cover_sat_problem_in_solver(
-    solver: CryptoMinisat,
-    expected_pauli_string: PauliString,
-    available_pauli_strings: list[PauliString],
-    qubits_to_consider: frozenset[int],
-) -> None:
-    """Build the SAT problem that should be solved to find a commuting cover.
-
-    This function encodes the SAT problem of interest into the provided ``solver``.
-    As such, the provided ``solver`` is expected to be a newly created instance,
-    and will be mutated by this function.
-
-    The encoded problem is the following:
-
-    Find the indices of Pauli strings in ``available_pauli_strings`` such that
-    the product of all the corresponding Pauli strings commutes with
-    ``expected_pauli_string`` on all the qubits in ``qubits_to_consider``.
-
-    The SAT problem contains one boolean variable per Pauli string in
-    ``available_pauli_strings``, deciding if the associated Pauli string should
-    be part of the cover.
-
-    The formula is quite simple: for each qubit in ``qubits_to_consider``, the
-    following two boolean formulas should be ``True``:
-
-    1. The X part of the resulting Pauli string, obtained from XORing the X part
-       of each Pauli string weighted by the boolean variable that decides if the
-       Pauli string should be included or not, should be equal to the X part of
-       the expected final Pauli string or be equal to I.
-    2. Same reasoning, but for the Z part.
-
-    Note that the "should be equal to the [P] part of the expected final Pauli
-    string" is natively implemented in the SAT solver used, so we do not need to
-    adapt the formula with some binary logic tricks. The solver really solves a
-    conjunction of ``a XOR b XOR ... XOR z == [True/False]``.
-
-    Also note that the "or be equal to I" part is a little bit more tricky and
-    requires some care. The logic retained for each of the possible Pauli terms
-    is explained inline within the code.
-
-    The number of XOR clauses in the formula depends on the number of each Pauli
-    terms in the input ``expected_pauli_string`` that is considered (i.e. on
-    qubits in ``qubits_to_consider``). Let [P] be the number of Pauli term equal
-    to P in ``expected_pauli_string``, the number of XOR clauses in the
-    resulting formula is ``(2 * [I] + [X] + [Y] + [Z]) * len(qubits_to_consider)``.
-    That means that there are less clauses in the SAT problem defined in this
-    function than in the one defined in
-    :func:`encode_pauli_string_exact_cover_sat_problem_in_solver`, except if
-    ``expected_pauli_string`` only contains identity gates, in which case the
-    two functions return the same number of XOR clauses.
-
-    Just like :func:`encode_pauli_string_exact_cover_sat_problem_in_solver`, the
-    number of terms in each XOR clause can vary depending on the provided Pauli
-    strings in ``available_pauli_strings``.
-
-    Also note that the problem defined above includes a trivial solution: do not
-    include any Pauli string. This lead to an identity Pauli string, that will
-    necessarily commute with the provided target. A specific clause is added to
-    the SAT problem to avoid that particular trivial solution.
-
-    Args:
-        solver: solver that will be modified in-place to encode the SAT problem.
-        expected_pauli_string: target Pauli string that should be covered by
-            strings from ``available_pauli_strings``.
-        available_pauli_strings: Pauli strings that can be used to try to cover
-            ``expected_pauli_string``.
-        qubits_to_consider: qubits on which the cover should be exactly equal to
-            ``expected_pauli_string``. Qubits not listed in this input will
-            simply be ignored and no restriction on the value of the resulting
-            cover on those qubits is added.
-    """
-    # Clause to exclude the trivial solution "include no Pauli string at all"
-    solver.add_clause([lit + 1 for lit in range(len(available_pauli_strings))])
-
-    for qubit in qubits_to_consider:
-        expected_effect = expected_pauli_string[qubit]
-        available_paulis_bools: list[tuple[bool, bool]] = [
-            pauli_literal_to_bools(pauli[qubit]) for pauli in available_pauli_strings
-        ]
-        # The two following lists includes the 1-based indices of Pauli strings from
-        # the provided `available_pauli_strings` input that finish with a X/Z stabilizers
-        # on the current `qubit`.
-        # For each index, we want to know if it should be included in the final cover or
-        # not, so XORing the boolean variables representing whether or not this Pauli
-        # string should be included results in the stabilizer propagated on that qubit.
-        x_clause_literals = [
-            pi + 1 for pi, (px, _) in enumerate(available_paulis_bools) if px
-        ]
-        z_clause_literals = [
-            pi + 1 for pi, (_, pz) in enumerate(available_paulis_bools) if pz
-        ]
-        if expected_effect == "I":
-            # Both the X and Z effect on that qubit should be OFF (i.e., identity).
-            solver.add_xor_clause(x_clause_literals, value=False)
-            solver.add_xor_clause(z_clause_literals, value=False)
-        elif expected_effect == "X":
-            # The X effect should be ON, the Z effect should be OFF.
-            # Because the identity (0, 0) is also valid, X can be either ON or OFF,
-            # so we do not need to restrict X effect.
-            solver.add_xor_clause(z_clause_literals, value=False)
-        elif expected_effect == "Y":
-            # This is an expected Y effect. This means that X and Z should be both
-            # ON (i.e., the Y effect) or both OFF (i.e., the identity effect).
-            # Rephrasing using XOR, this is equivalent to `X_effect XOR Z_effect == 0`.
-            # Note that measurements appearing twice here can be removed, as:
-            # - the XOR operation is commutative (so we can re-organise the measurements
-            #   to be sorted by index).
-            # - b XOR b = 0.
-            # - b XOR 0 = b.
-            # Also, a given measurement cannot appear more than twice, so we are fine
-            # just removing duplicated entries.
-            solver.add_xor_clause(
-                list(set(x_clause_literals) ^ set(z_clause_literals)), value=False
-            )
-        elif expected_effect == "Z":
-            # The X effect should be OFF, the Z effect should be ON.
-            # Because the identity (0, 0) is also valid, Z can be either ON or OFF,
-            # so we do not need to restrict Z effect.
-            solver.add_xor_clause(x_clause_literals, value=False)
diff --git a/src/tqecd/pauli.py b/src/tqecd/pauli.py
index 2a19e75..b1123da 100644
--- a/src/tqecd/pauli.py
+++ b/src/tqecd/pauli.py
@@ -8,8 +8,9 @@
 
 from __future__ import annotations
 
-import functools
 import operator
+from functools import reduce
+from itertools import chain
 from typing import Iterable, Literal
 
 import stim
@@ -197,6 +198,14 @@ def __hash__(self) -> int:
     def __getitem__(self, index: int) -> PAULI_STRING_TYPE:
         return self._pauli_by_qubit.get(index, "I")
 
+    def to_int(self, qubits: Iterable[int]) -> int:
+        """Convert the Pauli string to an integer representation on the provided qubits."""
+        return reduce(
+            lambda acc, bit: acc << 1 | bit,
+            chain.from_iterable(pauli_literal_to_bools(self[q]) for q in qubits),
+            0,
+        )
+
 
 def pauli_literal_to_bools(
     literal: PAULI_STRING_TYPE,
@@ -212,4 +221,4 @@ def pauli_literal_to_bools(
 
 
 def pauli_product(paulis: Iterable[PauliString]) -> PauliString:
-    return functools.reduce(operator.mul, paulis, PauliString({}))
+    return reduce(operator.mul, paulis, PauliString({}))