diff --git a/qiskit/circuit/library/generalized_gates/linear_function.py b/qiskit/circuit/library/generalized_gates/linear_function.py
index 15ab05ce426b..29c26b004f91 100644
--- a/qiskit/circuit/library/generalized_gates/linear_function.py
+++ b/qiskit/circuit/library/generalized_gates/linear_function.py
@@ -16,6 +16,7 @@
 import numpy as np
 from qiskit.circuit import QuantumCircuit, Gate
 from qiskit.circuit.exceptions import CircuitError
+from qiskit.synthesis.linear import check_invertible_binary_matrix
 
 
 class LinearFunction(Gate):
@@ -110,8 +111,7 @@ def __init__(
 
             # Optionally, check that the matrix is invertible
             if validate_input:
-                det = np.linalg.det(linear) % 2
-                if not np.allclose(det, 1):
+                if not check_invertible_binary_matrix(linear):
                     raise CircuitError(
                         "A linear function must be represented by an invertible matrix."
                     )
@@ -134,9 +134,9 @@ def synthesize(self):
         Returns:
             QuantumCircuit: A circuit implementing the evolution.
         """
-        from qiskit.transpiler.synthesis import cnot_synth
+        from qiskit.synthesis.linear import synth_cnot_count_full_pmh
 
-        return cnot_synth(self.linear)
+        return synth_cnot_count_full_pmh(self.linear)
 
     @property
     def linear(self):
diff --git a/qiskit/quantum_info/operators/dihedral/random.py b/qiskit/quantum_info/operators/dihedral/random.py
index 6ba2a354f4fd..d531eb9a43c6 100644
--- a/qiskit/quantum_info/operators/dihedral/random.py
+++ b/qiskit/quantum_info/operators/dihedral/random.py
@@ -48,10 +48,9 @@ def random_cnotdihedral(num_qubits, seed=None):
 
     # Random affine function
     # Random invertible binary matrix
-    det = 0
-    while np.allclose(det, 0) or np.allclose(det, 2):
-        linear = rng.integers(2, size=(num_qubits, num_qubits))
-        det = np.linalg.det(linear) % 2
+    from qiskit.synthesis.linear import random_invertible_binary_matrix
+
+    linear = random_invertible_binary_matrix(num_qubits, seed=rng)
     elem.linear = linear
 
     # Random shift
diff --git a/qiskit/synthesis/linear/__init__.py b/qiskit/synthesis/linear/__init__.py
new file mode 100644
index 000000000000..2229eda1f987
--- /dev/null
+++ b/qiskit/synthesis/linear/__init__.py
@@ -0,0 +1,21 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2018.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Module containing cnot circuits and cnot-phase circuit synthesize."""
+
+
+from .graysynth import graysynth, synth_cnot_count_full_pmh
+from .linear_matrix_utils import (
+    random_invertible_binary_matrix,
+    calc_inverse_matrix,
+    check_invertible_binary_matrix,
+)
diff --git a/qiskit/transpiler/synthesis/graysynth.py b/qiskit/synthesis/linear/graysynth.py
similarity index 99%
rename from qiskit/transpiler/synthesis/graysynth.py
rename to qiskit/synthesis/linear/graysynth.py
index 033e55dfe250..bb61f4bfb715 100644
--- a/qiskit/transpiler/synthesis/graysynth.py
+++ b/qiskit/synthesis/linear/graysynth.py
@@ -176,11 +176,11 @@ def graysynth(cnots, angles, section_size=2):
         else:
             sta.append([cnots1, list(set(ilist).difference([j])), qubit])
         sta.append([cnots0, list(set(ilist).difference([j])), qubit])
-    qcir &= cnot_synth(state, section_size).inverse()
+    qcir &= synth_cnot_count_full_pmh(state, section_size).inverse()
     return qcir
 
 
-def cnot_synth(state, section_size=2):
+def synth_cnot_count_full_pmh(state, section_size=2):
     """
     This function is an implementation of the Patel–Markov–Hayes algorithm
     for optimal synthesis of linear reversible circuits, as specified by an
diff --git a/qiskit/synthesis/linear/linear_circuits_utils.py b/qiskit/synthesis/linear/linear_circuits_utils.py
new file mode 100644
index 000000000000..375c91e9e35f
--- /dev/null
+++ b/qiskit/synthesis/linear/linear_circuits_utils.py
@@ -0,0 +1,102 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Utility functions for handling linear reversible circuits."""
+
+import copy
+from typing import Callable
+import numpy as np
+from qiskit import QuantumCircuit
+from qiskit.exceptions import QiskitError
+from qiskit.circuit.exceptions import CircuitError
+from . import calc_inverse_matrix, check_invertible_binary_matrix
+
+
+def transpose_cx_circ(qc: QuantumCircuit):
+    """Takes a circuit having only CX gates, and calculates its transpose.
+    This is done by recursively replacing CX(i, j) with CX(j, i) in all instructions.
+
+    Args:
+        qc: a QuantumCircuit containing only CX gates.
+
+    Returns:
+        QuantumCircuit: the transposed circuit.
+
+    Raises:
+        CircuitError: if qc has a non-CX gate.
+    """
+    transposed_circ = QuantumCircuit(qc.qubits, qc.clbits, name=qc.name + "_transpose")
+    for instruction in reversed(qc.data):
+        if instruction.operation.name != "cx":
+            raise CircuitError("The circuit contains non-CX gates.")
+        transposed_circ._append(instruction.replace(qubits=reversed(instruction.qubits)))
+    return transposed_circ
+
+
+def optimize_cx_4_options(function: Callable, mat: np.ndarray, optimize_count: bool = True):
+    """Get the best implementation of a circuit implementing a binary invertible matrix M,
+    by considering all four options: M,M^(-1),M^T,M^(-1)^T.
+    Optimizing either the CX count or the depth.
+
+    Args:
+        function: the synthesis function.
+        mat: a binary invertible matrix.
+        optimize_count: True if the number of CX gates in optimize, False if the depth is optimized.
+
+    Returns:
+        QuantumCircuit: an optimized QuantumCircuit, has the best depth or CX count of the four options.
+
+    Raises:
+        QiskitError: if mat is not an invertible matrix.
+    """
+    if not check_invertible_binary_matrix(mat):
+        raise QiskitError("The matrix is not invertible.")
+
+    qc = function(mat)
+    best_qc = qc
+    best_depth = qc.depth()
+    best_count = qc.count_ops()["cx"]
+
+    for i in range(1, 4):
+        mat_cpy = copy.deepcopy(mat)
+        # i=1 inverse, i=2 transpose, i=3 transpose and inverse
+        if i == 1:
+            mat_cpy = calc_inverse_matrix(mat_cpy)
+            qc = function(mat_cpy)
+            qc = qc.inverse()
+        elif i == 2:
+            mat_cpy = np.transpose(mat_cpy)
+            qc = function(mat_cpy)
+            qc = transpose_cx_circ(qc)
+        elif i == 3:
+            mat_cpy = calc_inverse_matrix(np.transpose(mat_cpy))
+            qc = function(mat_cpy)
+            qc = transpose_cx_circ(qc)
+            qc = qc.inverse()
+
+        new_depth = qc.depth()
+        new_count = qc.count_ops()["cx"]
+        # Prioritize count, and if it has the same count, then also consider depth
+        better_count = (optimize_count and best_count > new_count) or (
+            not optimize_count and best_depth == new_depth and best_count > new_count
+        )
+        # Prioritize depth, and if it has the same depth, then also consider count
+        better_depth = (not optimize_count and best_depth > new_depth) or (
+            optimize_count and best_count == new_count and best_depth > new_depth
+        )
+
+        if better_count or better_depth:
+            best_count = new_count
+            best_depth = new_depth
+            best_qc = qc
+
+    return best_qc
diff --git a/qiskit/synthesis/linear/linear_matrix_utils.py b/qiskit/synthesis/linear/linear_matrix_utils.py
new file mode 100644
index 000000000000..96e00d722b33
--- /dev/null
+++ b/qiskit/synthesis/linear/linear_matrix_utils.py
@@ -0,0 +1,149 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Utility functions for handling binary matrices."""
+
+from typing import Optional, Union
+import numpy as np
+from qiskit.exceptions import QiskitError
+
+
+def check_invertible_binary_matrix(mat: np.ndarray):
+    """Check that a binary matrix is invertible.
+
+    Args:
+        mat: a binary matrix.
+
+    Returns:
+        bool: True if mat in invertible and False otherwise.
+    """
+    if len(mat.shape) != 2 or mat.shape[0] != mat.shape[1]:
+        return False
+
+    mat = _gauss_elimination(mat)
+    rank = _compute_rank_after_gauss_elim(mat)
+    return rank == mat.shape[0]
+
+
+def random_invertible_binary_matrix(
+    num_qubits: int, seed: Optional[Union[np.random.Generator, int]] = None
+):
+    """Generates a random invertible n x n binary matrix.
+
+    Args:
+        num_qubits: the matrix size.
+        seed: a random seed.
+
+    Returns:
+        np.ndarray: A random invertible binary matrix of size num_qubits.
+    """
+    if isinstance(seed, np.random.Generator):
+        rng = seed
+    else:
+        rng = np.random.default_rng(seed)
+
+    rank = 0
+    while rank != num_qubits:
+        mat = rng.integers(2, size=(num_qubits, num_qubits))
+        mat_gauss = mat.copy()
+        mat_gauss = _gauss_elimination(mat_gauss)
+        rank = _compute_rank_after_gauss_elim(mat_gauss)
+    return mat
+
+
+def _gauss_elimination(mat, ncols=None, full_elim=False):
+    """Gauss elimination of a matrix mat with m rows and n columns.
+    If full_elim = True, it allows full elimination of mat[:, 0 : ncols]
+    Mutates and returns the matrix mat."""
+
+    # Treat the matrix A as containing integer values
+    mat = np.array(mat, dtype=int, copy=False)
+
+    m = mat.shape[0]  # no. of rows
+    n = mat.shape[1]  # no. of columns
+    if ncols is not None:
+        n = min(n, ncols)  # no. of active columns
+
+    r = 0  # current rank
+    k = 0  # current pivot column
+    while (r < m) and (k < n):
+        is_non_zero = False
+        new_r = r
+        for j in range(k, n):
+            for i in range(r, m):
+                if mat[i][j]:
+                    is_non_zero = True
+                    k = j
+                    new_r = i
+                    break
+            if is_non_zero:
+                break
+        if not is_non_zero:
+            return mat  # A is in the canonical form
+
+        if new_r != r:
+            mat[[r, new_r]] = mat[[new_r, r]]
+
+        if full_elim:
+            for i in range(0, r):
+                if mat[i][k]:
+                    mat[i] = mat[i] ^ mat[r]
+
+        for i in range(r + 1, m):
+            if mat[i][k]:
+                mat[i] = mat[i] ^ mat[r]
+        r += 1
+
+    return mat
+
+
+def calc_inverse_matrix(mat: np.ndarray, verify: bool = False):
+    """Given a square numpy(dtype=int) matrix mat, tries to compute its inverse.
+
+    Args:
+        mat: a boolean square matrix.
+        verify: if True asserts that the multiplication of mat and its inverse is the identity matrix.
+
+    Returns:
+        np.ndarray: the inverse matrix.
+
+    Raises:
+         QiskitError: if the matrix is not square.
+         QiskitError: if the matrix is not invertible.
+    """
+
+    if mat.shape[0] != mat.shape[1]:
+        raise QiskitError("Matrix to invert is a non-square matrix.")
+
+    n = mat.shape[0]
+    # concatenate the matrix and identity
+    mat1 = np.concatenate((mat, np.eye(n, dtype=int)), axis=1)
+    mat1 = _gauss_elimination(mat1, None, full_elim=True)
+
+    r = _compute_rank_after_gauss_elim(mat1[:, 0:n])
+
+    if r < n:
+        raise QiskitError("The matrix is not invertible.")
+
+    matinv = mat1[:, n : 2 * n]
+
+    if verify:
+        mat2 = np.dot(mat, matinv) % 2
+        assert np.array_equal(mat2, np.eye(n))
+
+    return matinv
+
+
+def _compute_rank_after_gauss_elim(mat):
+    """Given a matrix A after Gaussian elimination, computes its rank
+    (i.e. simply the number of nonzero rows)"""
+    return np.sum(mat.any(axis=1))
diff --git a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
index 0edb47fab932..168139d42333 100644
--- a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
+++ b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
@@ -19,7 +19,7 @@
 from qiskit.dagcircuit.dagcircuit import DAGCircuit
 from qiskit.transpiler.exceptions import TranspilerError
 from qiskit.quantum_info import decompose_clifford
-from qiskit.transpiler.synthesis import cnot_synth
+from qiskit.synthesis.linear import synth_cnot_count_full_pmh
 from .plugin import HighLevelSynthesisPluginManager, HighLevelSynthesisPlugin
 
 
@@ -168,5 +168,5 @@ class DefaultSynthesisLinearFunction(HighLevelSynthesisPlugin):
 
     def run(self, high_level_object, **options):
         """Run synthesis for the given LinearFunction."""
-        decomposition = cnot_synth(high_level_object.linear)
+        decomposition = synth_cnot_count_full_pmh(high_level_object.linear)
         return decomposition
diff --git a/qiskit/transpiler/synthesis/__init__.py b/qiskit/transpiler/synthesis/__init__.py
index 575daf2625b0..bf4e7e4888a5 100644
--- a/qiskit/transpiler/synthesis/__init__.py
+++ b/qiskit/transpiler/synthesis/__init__.py
@@ -11,6 +11,3 @@
 # that they have been altered from the originals.
 
 """Module containing transpiler synthesize."""
-
-
-from .graysynth import graysynth, cnot_synth
diff --git a/releasenotes/notes/linear_function_synthesis_utils-f2f96924ca45e1fb.yaml b/releasenotes/notes/linear_function_synthesis_utils-f2f96924ca45e1fb.yaml
new file mode 100644
index 000000000000..a509ab6c03ec
--- /dev/null
+++ b/releasenotes/notes/linear_function_synthesis_utils-f2f96924ca45e1fb.yaml
@@ -0,0 +1,23 @@
+---
+features:
+  - |
+    Add new utility functions to handle linear functions and circuits:
+    * A utility function that checks whether a NxN binary matrix is invertible.
+    * A utility function that computes an inverse of a NxN invertible binary matrix.
+    * A utility function that calculates the transpose of a linear reversible circuit.
+    * Implement the following "meta" idea for synthesizing linear functions:
+    given an invertible binary matrix A and a synthesis algorithm f: A -> QuantumCircuit,
+    in addition to applying f to A, we can also apply f to A^{-1} and invert the computed circuit,
+    apply f to A^{T} and transpose the computed circuit, and to apply f to A^{T}^{-1} and to invert
+    and transpose the computed circuit, thus leading to 4 generally different synthesized circuits for A,
+    allowing to pick the best one in terms of depth or the number of gates.
+upgrade:
+  - |
+    Improve code structure by moving the internal code graysynth.py from qiskit/transpiler/synthesis
+    to qiskit/synthesis/linear, and renaming cnot_synth to synth_cnot_count_full_pmh.
+fixes:
+  - |
+    Provide an internal utility function to generate a random NxN invertible binary matrix.
+    This function is completely robust, as opposed to the previously implemented method based on
+    computing determinant (via np.linalg.det) and checking that it's close to 1,
+    where we already saw cases of floating point errors leading to an incorrect result (for N>100).
diff --git a/test/python/circuit/library/test_linear_function.py b/test/python/circuit/library/test_linear_function.py
index ee96fe6ba70d..b7ca592ace80 100644
--- a/test/python/circuit/library/test_linear_function.py
+++ b/test/python/circuit/library/test_linear_function.py
@@ -22,6 +22,7 @@
 from qiskit.circuit.library.standard_gates import CXGate, SwapGate
 from qiskit.circuit.library.generalized_gates import LinearFunction
 from qiskit.circuit.exceptions import CircuitError
+from qiskit.synthesis.linear import random_invertible_binary_matrix
 
 from qiskit.quantum_info.operators import Operator
 
@@ -51,23 +52,6 @@ def random_linear_circuit(num_qubits, num_gates, seed=None):
     return circ
 
 
-def random_invertible_binary_matrix(num_qubits, seed=None):
-    """Generates a random invertible n x n binary matrix."""
-
-    # This code is adapted from random_cnotdihedral
-    if isinstance(seed, np.random.Generator):
-        rng = seed
-    else:
-        rng = np.random.default_rng(seed)
-
-    det = 0
-    while np.allclose(det, 0) or np.allclose(det, 2):
-        binary_matrix = rng.integers(2, size=(num_qubits, num_qubits))
-        det = np.linalg.det(binary_matrix) % 2
-
-    return binary_matrix
-
-
 @ddt
 class TestLinearFunctions(QiskitTestCase):
     """Tests for clifford append gate functions."""
@@ -132,7 +116,7 @@ def test_conversion_to_linear_function_and_back(self, num_qubits):
     def test_patel_markov_hayes(self):
         """Checks the explicit example from Patel-Markov-Hayes's paper."""
 
-        # This code is adapted from test_synthesis.py
+        # This code is adapted from test_gray_synthesis.py
         binary_matrix = [
             [1, 1, 0, 0, 0, 0],
             [1, 0, 0, 1, 1, 0],
diff --git a/test/python/synthesis/__init__.py b/test/python/synthesis/__init__.py
new file mode 100644
index 000000000000..f7c02dacd381
--- /dev/null
+++ b/test/python/synthesis/__init__.py
@@ -0,0 +1,13 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Qiskit synthesis unit tests."""
diff --git a/test/python/transpiler/test_synthesis.py b/test/python/synthesis/test_gray_synthesis.py
similarity index 98%
rename from test/python/transpiler/test_synthesis.py
rename to test/python/synthesis/test_gray_synthesis.py
index 1f7fb297d5a4..738567df0d86 100644
--- a/test/python/transpiler/test_synthesis.py
+++ b/test/python/synthesis/test_gray_synthesis.py
@@ -10,12 +10,14 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""Test synthesis algorithms"""
+"""Test cnot circuit and cnot-phase circuit synthesis algorithms"""
+
+import unittest
 
 from qiskit.circuit import QuantumCircuit, QuantumRegister
 from qiskit.quantum_info.operators import Operator
 from qiskit.extensions.unitary import UnitaryGate
-from qiskit.transpiler.synthesis import graysynth, cnot_synth
+from qiskit.synthesis.linear import graysynth, synth_cnot_count_full_pmh
 from qiskit.test import QiskitTestCase
 
 
@@ -237,7 +239,7 @@ def test_patel_markov_hayes(self):
             [1, 1, 0, 1, 1, 1],
             [0, 0, 1, 1, 1, 0],
         ]
-        c_patel = cnot_synth(state)
+        c_patel = synth_cnot_count_full_pmh(state)
         unitary_patel = UnitaryGate(Operator(c_patel))
 
         # Create the circuit displayed above:
@@ -262,3 +264,7 @@ def test_patel_markov_hayes(self):
 
         # Check if the two circuits are equivalent
         self.assertEqual(unitary_patel, unitary_compare)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/python/synthesis/test_linear_synthesis.py b/test/python/synthesis/test_linear_synthesis.py
new file mode 100644
index 000000000000..04c0ad50c5f5
--- /dev/null
+++ b/test/python/synthesis/test_linear_synthesis.py
@@ -0,0 +1,114 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2022.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Test linear reversible circuits synthesis functions."""
+
+import unittest
+
+import numpy as np
+from qiskit import QuantumCircuit
+from qiskit.circuit.library import LinearFunction
+from qiskit.synthesis.linear import (
+    synth_cnot_count_full_pmh,
+    random_invertible_binary_matrix,
+    check_invertible_binary_matrix,
+    calc_inverse_matrix,
+)
+from qiskit.synthesis.linear.linear_circuits_utils import transpose_cx_circ, optimize_cx_4_options
+from qiskit.test import QiskitTestCase
+
+
+class TestLinearSynth(QiskitTestCase):
+    """Test the linear reversible circuit synthesis functions."""
+
+    def test_lnn_circuit(self):
+        """Test the synthesis of a CX circuit with LNN connectivity."""
+
+        n = 5
+        qc = QuantumCircuit(n)
+        for i in range(n - 1):
+            qc.cx(i, i + 1)
+        mat = LinearFunction(qc).linear
+
+        for optimized in [True, False]:
+            optimized_qc = optimize_cx_4_options(
+                synth_cnot_count_full_pmh, mat, optimize_count=optimized
+            )
+            self.assertEqual(optimized_qc.depth(), 4)
+            self.assertEqual(optimized_qc.count_ops()["cx"], 4)
+
+    def test_full_circuit(self):
+        """Test the synthesis of a CX circuit with full connectivity."""
+
+        n = 5
+        qc = QuantumCircuit(n)
+        for i in range(n):
+            for j in range(i + 1, n):
+                qc.cx(i, j)
+        mat = LinearFunction(qc).linear
+
+        for optimized in [True, False]:
+            optimized_qc = optimize_cx_4_options(
+                synth_cnot_count_full_pmh, mat, optimize_count=optimized
+            )
+            self.assertEqual(optimized_qc.depth(), 4)
+            self.assertEqual(optimized_qc.count_ops()["cx"], 4)
+
+    def test_transpose_circ(self):
+        """Test the transpose_cx_circ() function."""
+        n = 5
+        mat = random_invertible_binary_matrix(n, seed=1234)
+        qc = synth_cnot_count_full_pmh(mat)
+        transposed_qc = transpose_cx_circ(qc)
+        transposed_mat = LinearFunction(transposed_qc).linear.astype(int)
+        self.assertTrue((mat.transpose() == transposed_mat).all())
+
+    def test_example_circuit(self):
+        """Test the synthesis of an example CX circuit which provides different CX count
+        and depth for different optimization methods."""
+
+        qc = QuantumCircuit(9)
+        qc.swap(8, 7)
+        qc.swap(7, 6)
+        qc.cx(5, 6)
+        qc.cx(6, 5)
+        qc.swap(4, 5)
+        qc.cx(3, 4)
+        qc.cx(4, 3)
+        qc.swap(2, 3)
+        qc.cx(1, 2)
+        qc.cx(2, 1)
+        qc.cx(0, 1)
+        qc.cx(1, 0)
+        mat = LinearFunction(qc).linear
+
+        optimized_qc = optimize_cx_4_options(synth_cnot_count_full_pmh, mat, optimize_count=True)
+        self.assertEqual(optimized_qc.depth(), 17)
+        self.assertEqual(optimized_qc.count_ops()["cx"], 20)
+
+        optimized_qc = optimize_cx_4_options(synth_cnot_count_full_pmh, mat, optimize_count=False)
+        self.assertEqual(optimized_qc.depth(), 15)
+        self.assertEqual(optimized_qc.count_ops()["cx"], 23)
+
+    def test_invertible_matrix(self):
+        """Test the functions for generating a random invertible matrix and inverting it."""
+        n = 5
+        mat = random_invertible_binary_matrix(n, seed=1234)
+        out = check_invertible_binary_matrix(mat)
+        mat_inv = calc_inverse_matrix(mat, verify=True)
+        mat_out = np.dot(mat, mat_inv) % 2
+        self.assertTrue(np.array_equal(mat_out, np.eye(n)))
+        self.assertTrue(out)
+
+
+if __name__ == "__main__":
+    unittest.main()