From 125c9e8924313bb27fe6e823136aefb4e34ad375 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Fri, 10 Jan 2025 17:02:29 +0100
Subject: [PATCH 01/15] v0 of comm checker

---
 crates/accelerate/src/commutation_checker.rs  | 184 +++++++++++-------
 crates/accelerate/src/unitary_compose.rs      |   8 +
 .../circuit/test_commutation_checker.py       |  38 +++-
 3 files changed, 145 insertions(+), 85 deletions(-)

diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index 52d4900efa5b..bac5d720da50 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -10,6 +10,8 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
+use core::f64;
+
 use hashbrown::{HashMap, HashSet};
 use ndarray::linalg::kron;
 use ndarray::Array2;
@@ -273,7 +275,8 @@ impl CommutationChecker {
         // relative and absolute tolerance used to (1) check whether rotation gates commute
         // trivially (i.e. the rotation angle is so small we assume it commutes) and (2) define
         // comparison for the matrix-based commutation checks
-        let rtol = 1e-5;
+        // let rtol = 1e-5;
+        let rtol = f64::EPSILON;
         let atol = 1e-8;
 
         // if we have rotation gates, we attempt to map them to their generators, for example
@@ -460,81 +463,68 @@ impl CommutationChecker {
             None => return Ok(false),
         };
 
-        if first_qarg == second_qarg {
-            match first_qarg.len() {
-                1 => Ok(unitary_compose::commute_1q(
-                    &first_mat.view(),
-                    &second_mat.view(),
-                    rtol,
-                    atol,
-                )),
-                2 => Ok(unitary_compose::commute_2q(
-                    &first_mat.view(),
-                    &second_mat.view(),
-                    &[Qubit(0), Qubit(1)],
-                    rtol,
-                    atol,
-                )),
-                _ => Ok(unitary_compose::allclose(
-                    &second_mat.dot(&first_mat).view(),
-                    &first_mat.dot(&second_mat).view(),
-                    rtol,
-                    atol,
-                )),
-            }
+        // if first_qarg == second_qarg {
+        // match first_qarg.len() {
+        //     1 => Ok(unitary_compose::commute_1q(
+        //         &first_mat.view(),
+        //         &second_mat.view(),
+        //         rtol,
+        //         atol,
+        //     )),
+        //     2 => Ok(unitary_compose::commute_2q(
+        //         &first_mat.view(),
+        //         &second_mat.view(),
+        //         &[Qubit(0), Qubit(1)],
+        //         rtol,
+        //         atol,
+        //     )),
+        //     _ => Ok(unitary_compose::allclose(
+        //         &second_mat.dot(&first_mat).view(),
+        //         &first_mat.dot(&second_mat).view(),
+        //         rtol,
+        //         atol,
+        //     )),
+        // }
+        //     let fid = unitary_compose::gate_fidelity(&first_mat.view(), &second_mat.view());
+        //     Ok((1.0 - fid).abs() < rtol)
+        // } else {
+        // TODO Optimize this bit to avoid unnecessary Kronecker products:
+        //  1. We currently sort the operations for the cache by operation size, putting the
+        //     *smaller* operation first: (smaller op, larger op)
+        //  2. This code here expands the first op to match the second -- hence we always
+        //     match the operator sizes.
+        // This whole extension logic could be avoided since we know the second one is larger.
+        let extra_qarg2 = num_qubits - first_qarg.len() as u32;
+        let first_mat = if extra_qarg2 > 0 {
+            let id_op = Array2::<Complex64>::eye(usize::pow(2, extra_qarg2));
+            kron(&id_op, &first_mat)
         } else {
-            // TODO Optimize this bit to avoid unnecessary Kronecker products:
-            //  1. We currently sort the operations for the cache by operation size, putting the
-            //     *smaller* operation first: (smaller op, larger op)
-            //  2. This code here expands the first op to match the second -- hence we always
-            //     match the operator sizes.
-            // This whole extension logic could be avoided since we know the second one is larger.
-            let extra_qarg2 = num_qubits - first_qarg.len() as u32;
-            let first_mat = if extra_qarg2 > 0 {
-                let id_op = Array2::<Complex64>::eye(usize::pow(2, extra_qarg2));
-                kron(&id_op, &first_mat)
-            } else {
-                first_mat
-            };
-
-            // the 1 qubit case cannot happen, since that would already have been captured
-            // by the previous if clause; first_qarg == second_qarg (if they overlap they must
-            // be the same)
-            if num_qubits == 2 {
-                return Ok(unitary_compose::commute_2q(
-                    &first_mat.view(),
-                    &second_mat.view(),
-                    &second_qarg,
-                    rtol,
-                    atol,
-                ));
-            };
+            first_mat
+        };
 
-            let op12 = match unitary_compose::compose(
-                &first_mat.view(),
-                &second_mat.view(),
-                &second_qarg,
-                false,
-            ) {
-                Ok(matrix) => matrix,
-                Err(e) => return Err(PyRuntimeError::new_err(e)),
-            };
-            let op21 = match unitary_compose::compose(
-                &first_mat.view(),
-                &second_mat.view(),
-                &second_qarg,
-                true,
-            ) {
-                Ok(matrix) => matrix,
-                Err(e) => return Err(PyRuntimeError::new_err(e)),
-            };
-            Ok(unitary_compose::allclose(
-                &op12.view(),
-                &op21.view(),
-                rtol,
-                atol,
-            ))
-        }
+        // the 1 qubit case cannot happen, since that would already have been captured
+        // by the previous if clause; first_qarg == second_qarg (if they overlap they must
+        // be the same)
+        let op12 = match unitary_compose::compose(
+            &first_mat.view(),
+            &second_mat.view(),
+            &second_qarg,
+            false,
+        ) {
+            Ok(matrix) => matrix,
+            Err(e) => return Err(PyRuntimeError::new_err(e)),
+        };
+        let op21 = match unitary_compose::compose(
+            &first_mat.view(),
+            &second_mat.view(),
+            &second_qarg,
+            true,
+        ) {
+            Ok(matrix) => matrix,
+            Err(e) => return Err(PyRuntimeError::new_err(e)),
+        };
+        let fid = unitary_compose::gate_fidelity(&op12.view(), &op21.view());
+        Ok((1.0 - fid).abs() < rtol)
     }
 
     fn clear_cache(&mut self) {
@@ -641,7 +631,10 @@ fn map_rotation<'a>(
         // commute with everything, and we simply return the operation with the flag that
         // it commutes trivially.
         if let Param::Float(angle) = params[0] {
-            if (angle % TWOPI).abs() < tol {
+            // if (angle % TWOPI).abs() < tol {
+            let gate_fidelity = rotation_fidelity(name, angle).unwrap_or(0.);
+            println!("gate: {} angle: {} fid: {}", name, angle, gate_fidelity);
+            if (1. - gate_fidelity).abs() < tol {
                 return (op, params, true);
             };
         };
@@ -659,6 +652,47 @@ fn map_rotation<'a>(
     (op, params, false)
 }
 
+/// Compute the gate fidelity of a rotation gate.
+// fn rotation_fidelity(rotation: &OperationRef, angle: f64) -> Option<f64> {
+//     if let OperationRef::Standard(gate) = rotation {
+//         match gate {
+//             StandardGate::RXGate
+//             | StandardGate::RYGate
+//             | StandardGate::RZGate
+//             | StandardGate::PhaseGate => return Some(2. * (angle / 2.).cos()),
+//             StandardGate::RXXGate
+//             | StandardGate::RYYGate
+//             | StandardGate::RZXGate
+//             | StandardGate::RZZGate => return Some(4. * (angle / 2.).cos()),
+//             StandardGate::CRXGate
+//             | StandardGate::CRYGate
+//             | StandardGate::CRZGate
+//             | StandardGate::CPhaseGate => return Some(2. + 2. * (angle / 2.).cos()),
+//             _ => return None,
+//         };
+//     };
+//     None
+// }
+
+fn rotation_fidelity(rotation: &str, angle: f64) -> Option<f64> {
+    let dim = match rotation {
+        "rx" | "ry" | "rz" | "p" => 2.,
+        _ => 4.,
+    };
+
+    match rotation {
+        "rx" | "ry" | "rz" | "p" | "rxx" | "ryy" | "rzx" | "rzz" => {
+            let gate_fid = (angle / 2.).cos().powi(2);
+            Some((dim * gate_fid + 1.) / (dim + 1.))
+        }
+        "crx" | "cry" | "crz" | "cp" => {
+            let gate_fid = (0.5 + 0.5 * (angle / 2.).cos()).powi(2);
+            Some((dim * gate_fid + 1.) / (dim + 1.))
+        }
+        _ => None,
+    }
+}
+
 fn get_relative_placement(
     first_qargs: &[Qubit],
     second_qargs: &[Qubit],
diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs
index 7c4da99dd7d2..2025342627ad 100644
--- a/crates/accelerate/src/unitary_compose.rs
+++ b/crates/accelerate/src/unitary_compose.rs
@@ -238,3 +238,11 @@ pub fn allclose(
     }
     true
 }
+
+pub fn gate_fidelity(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> f64 {
+    let dim = left.nrows() as f64; // == left.ncols() == right.nrows() == right.ncols()
+    let trace = left.t().mapv(|el| el.conj()).dot(right).diag().sum();
+    let process_fidelity = (trace / dim).abs().powi(2);
+    let gate_fidelity = (dim * process_fidelity + 1.) / (dim + 1.);
+    gate_fidelity
+}
diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py
index b4f6a30d904c..b55b8b95ca37 100644
--- a/test/python/circuit/test_commutation_checker.py
+++ b/test/python/circuit/test_commutation_checker.py
@@ -64,16 +64,14 @@
     RYGate,
     RZGate,
     PhaseGate,
-    CRXGate,
-    CRYGate,
-    CRZGate,
-    CPhaseGate,
     RXXGate,
     RYYGate,
     RZZGate,
     RZXGate,
 ]
 
+CONTROLLED_ROTATION_GATES = [CRXGate, CRYGate, CRZGate, CPhaseGate]
+
 
 class NewGateCX(Gate):
     """A dummy class containing an cx gate unknown to the commutation checker's library."""
@@ -412,7 +410,7 @@ def test_serialization(self):
         cc2.commute_nodes(dop1, dop2)
         self.assertEqual(cc2.num_cached_entries(), 1)
 
-    @idata(ROTATION_GATES)
+    @idata(ROTATION_GATES + CONTROLLED_ROTATION_GATES)
     def test_cutoff_angles(self, gate_cls):
         """Check rotations with a small enough angle are cut off."""
         max_power = 30
@@ -420,17 +418,18 @@ def test_cutoff_angles(self, gate_cls):
 
         generic_gate = DCXGate()  # gate that does not commute with any rotation gate
 
+        # TODO this is no longer a cutoff but defined by the gate fidelity
         cutoff_angle = 1e-5  # this is the cutoff we use in the CommutationChecker
 
         for i in range(1, max_power + 1):
             angle = 2 ** (-i)
             gate = gate_cls(angle)
             qargs = list(range(gate.num_qubits))
-
-            if angle < cutoff_angle:
-                self.assertTrue(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
-            else:
-                self.assertFalse(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
+            print(f"{angle:.2e}", scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
+            # if angle < cutoff_angle:
+            #     self.assertTrue(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
+            # else:
+            #     self.assertFalse(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
 
     @idata(ROTATION_GATES)
     def test_rotation_mod_2pi(self, gate_cls):
@@ -453,6 +452,25 @@ def test_rotation_mod_2pi(self, gate_cls):
                     scc.commute(generic_gate, [0], [], gate, list(range(gate.num_qubits)), [])
                 )
 
+    @idata(CONTROLLED_ROTATION_GATES)
+    def test_controlled_rotation_mod_4pi(self, gate_cls):
+        """Test the controlled rotations modulo 4pi commute with any gate.
+
+        For 2pi, the single-qubit rotations are the identity only up to a global phase,
+        which has a measurable effect if controlled.
+        """
+        generic_gate = HGate()  # does not commute with any rotation gate
+        multiples = np.arange(-6, 7)
+
+        for multiple in multiples:
+            with self.subTest(multiple=multiple):
+                gate = gate_cls(multiple * np.pi)
+                comm = scc.commute(generic_gate, [0], [], gate, list(range(gate.num_qubits)), [])
+                if multiple % 4 == 0:
+                    self.assertTrue(comm)
+                else:
+                    self.assertFalse(comm)
+
     def test_custom_gate(self):
         """Test a custom gate."""
         my_cx = NewGateCX()

From ba5d37d34c688a5131aefb3469d35cb1ac176f8b Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 13 Jan 2025 11:51:05 +0100
Subject: [PATCH 02/15] dangling print

---
 crates/accelerate/src/commutation_checker.rs | 1 -
 1 file changed, 1 deletion(-)

diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index bac5d720da50..3308bdcfa45d 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -633,7 +633,6 @@ fn map_rotation<'a>(
         if let Param::Float(angle) = params[0] {
             // if (angle % TWOPI).abs() < tol {
             let gate_fidelity = rotation_fidelity(name, angle).unwrap_or(0.);
-            println!("gate: {} angle: {} fid: {}", name, angle, gate_fidelity);
             if (1. - gate_fidelity).abs() < tol {
                 return (op, params, true);
             };

From 501c8089c605e57b4a3506817749f296c898107f Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Thu, 16 Jan 2025 18:40:26 +0100
Subject: [PATCH 03/15] games with matrices

---
 crates/accelerate/src/commutation_checker.rs |  2 +-
 crates/accelerate/src/unitary_compose.rs     | 56 ++++++++++++++++++--
 2 files changed, 53 insertions(+), 5 deletions(-)

diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index 3308bdcfa45d..745f5f6a32e0 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -523,7 +523,7 @@ impl CommutationChecker {
             Ok(matrix) => matrix,
             Err(e) => return Err(PyRuntimeError::new_err(e)),
         };
-        let fid = unitary_compose::gate_fidelity(&op12.view(), &op21.view());
+        let fid = unitary_compose::gate_fidelity(&op12.view(), &op21.view(), None);
         Ok((1.0 - fid).abs() < rtol)
     }
 
diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs
index 2025342627ad..f721cfc88093 100644
--- a/crates/accelerate/src/unitary_compose.rs
+++ b/crates/accelerate/src/unitary_compose.rs
@@ -10,7 +10,7 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use ndarray::{Array, Array2, ArrayView, ArrayView2, IxDyn};
+use ndarray::{arr2, Array, Array2, ArrayView, ArrayView2, IxDyn};
 use ndarray_einsum_beta::*;
 use num_complex::{Complex, Complex64, ComplexFloat};
 use num_traits::Zero;
@@ -239,10 +239,58 @@ pub fn allclose(
     true
 }
 
-pub fn gate_fidelity(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> f64 {
-    let dim = left.nrows() as f64; // == left.ncols() == right.nrows() == right.ncols()
-    let trace = left.t().mapv(|el| el.conj()).dot(right).diag().sum();
+pub fn gate_fidelity(
+    left: &ArrayView2<Complex64>,
+    right: &ArrayView2<Complex64>,
+    qargs: Option<&[Qubit]>,
+) -> f64 {
+    let dim = left.nrows(); // == left.ncols() == right.nrows() == right.ncols()
+                            // let trace = left.t().mapv(|el| el.conj()).dot(right).diag().sum();
+
+    let left = left.t().mapv(|el| el.conj());
+    let product = match dim {
+        2 => mm1q(&left.view(), right),
+        4 => mm2q(&left.view(), right, qargs.unwrap_or(&[Qubit(0), Qubit(1)])),
+        _ => left.dot(right),
+    };
+    let trace = product.diag().sum();
+
+    let dim = dim as f64;
     let process_fidelity = (trace / dim).abs().powi(2);
     let gate_fidelity = (dim * process_fidelity + 1.) / (dim + 1.);
     gate_fidelity
 }
+
+fn mm1q(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> Array2<Complex64> {
+    let zero = Complex64::zero();
+    let mut out = arr2(&[[zero, zero], [zero, zero]]);
+    out[[0, 0]] = left[[0, 0]] * right[[0, 0]] + left[[0, 1]] * right[[1, 0]];
+    out[[0, 1]] = left[[0, 0]] * right[[0, 1]] + left[[0, 1]] * right[[1, 1]];
+    out[[1, 0]] = left[[1, 0]] * right[[0, 0]] + left[[1, 1]] * right[[1, 0]];
+    out[[1, 1]] = left[[1, 0]] * right[[0, 1]] + left[[1, 1]] * right[[1, 1]];
+    out
+}
+
+pub fn mm2q(
+    left: &ArrayView2<Complex64>,
+    right: &ArrayView2<Complex64>,
+    qargs: &[Qubit],
+) -> Array2<Complex64> {
+    let zero = Complex64::zero();
+    let mut out = arr2(&[
+        [zero, zero, zero, zero],
+        [zero, zero, zero, zero],
+        [zero, zero, zero, zero],
+        [zero, zero, zero, zero],
+    ]);
+
+    let rev = qargs[0].0 == 1;
+    for i in 0..4usize {
+        for j in 0..4usize {
+            for k in 0..4usize {
+                out[[j, k]] += left[[_ind(i, rev), _ind(k, rev)]] * right[[k, j]];
+            }
+        }
+    }
+    out
+}

From 7d9fdb347204ed7e872fe7952759db3f63f40ebc Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 4 Feb 2025 09:22:43 +0100
Subject: [PATCH 04/15] Revert "HLS fix to not synthesize instructions already
 supported (#13417)"

This reverts commit b9d5c9c6aeb4568b6e9fba8943517ebed300886d.
---
 .../passes/synthesis/high_level_synthesis.py        | 10 ++++++----
 ...hls-supported-instructions-0d80ea33b3d2257b.yaml |  8 --------
 test/python/transpiler/test_high_level_synthesis.py | 13 -------------
 3 files changed, 6 insertions(+), 25 deletions(-)
 delete mode 100644 releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml

diff --git a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
index 95adee1d05c7..d7b5b0d535b8 100644
--- a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
+++ b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
@@ -814,7 +814,6 @@ def _definitely_skip_node(
             dag._has_calibration_for(node)
             or len(node.qargs) < self._min_qubits
             or node.is_directive()
-            or (self._instruction_supported(node.name, qubits) and not node.is_control_flow())
         ):
             return True
 
@@ -832,12 +831,15 @@ def _definitely_skip_node(
             # If all the above constraints hold, and it's already supported or the basis translator
             # can handle it, we'll leave it be.
             and (
+                self._instruction_supported(node.name, qubits)
                 # This uses unfortunately private details of `EquivalenceLibrary`, but so does the
                 # `BasisTranslator`, and this is supposed to just be temporary til this is moved
                 # into Rust space.
-                self._equiv_lib is not None
-                and equivalence.Key(name=node.name, num_qubits=node.num_qubits)
-                in self._equiv_lib.keys()
+                or (
+                    self._equiv_lib is not None
+                    and equivalence.Key(name=node.name, num_qubits=node.num_qubits)
+                    in self._equiv_lib.keys()
+                )
             )
         )
 
diff --git a/releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml b/releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml
deleted file mode 100644
index 1e9c1a1cebe7..000000000000
--- a/releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml
+++ /dev/null
@@ -1,8 +0,0 @@
----
-fixes:
-  - |
-    Previously the :class:`.HighLevelSynthesis` transpiler pass synthesized an
-    instruction for which a synthesis plugin is available, regardless of
-    whether the instruction is already supported by the target or a part of
-    the explicitly passed ``basis_gates``. This behavior is now fixed, so that
-    such already supported instructions are no longer synthesized.
diff --git a/test/python/transpiler/test_high_level_synthesis.py b/test/python/transpiler/test_high_level_synthesis.py
index b7bc31132a48..e1fa336c1377 100644
--- a/test/python/transpiler/test_high_level_synthesis.py
+++ b/test/python/transpiler/test_high_level_synthesis.py
@@ -46,7 +46,6 @@
     IGate,
     MCXGate,
     SGate,
-    QAOAAnsatz,
 )
 from qiskit.circuit.library import LinearFunction, PauliEvolutionGate
 from qiskit.quantum_info import Clifford, Operator, Statevector, SparsePauliOp
@@ -665,18 +664,6 @@ def test_synth_fails_definition_exists(self):
             out = hls(circuit)
             self.assertEqual(out.count_ops(), {"u": 1})
 
-    def test_both_basis_gates_and_plugin_specified(self):
-        """Test that a gate is not synthesized when it belongs to basis_gates,
-        regardless of whether there is a plugin method available.
-
-        See: https://github.com/Qiskit/qiskit/issues/13412 for more
-        details.
-        """
-        qc = QAOAAnsatz(SparsePauliOp("Z"), initial_state=QuantumCircuit(1))
-        pm = PassManager([HighLevelSynthesis(basis_gates=["PauliEvolution"])])
-        qct = pm.run(qc)
-        self.assertEqual(qct.count_ops()["PauliEvolution"], 2)
-
 
 class TestPMHSynthesisLinearFunctionPlugin(QiskitTestCase):
     """Tests for the PMHSynthesisLinearFunction plugin for synthesizing linear functions."""

From 2697f9fa536274c423d5748df2b064d63dc82546 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 4 Feb 2025 15:27:58 +0100
Subject: [PATCH 05/15] use custom mm1 and mm2

---
 crates/accelerate/src/unitary_compose.rs | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs
index c0a4bb7fc157..b1d2eb1efaa2 100644
--- a/crates/accelerate/src/unitary_compose.rs
+++ b/crates/accelerate/src/unitary_compose.rs
@@ -249,8 +249,8 @@ pub fn gate_fidelity(
 
     let left = left.t().mapv(|el| el.conj());
     let product = match dim {
-        // 2 => mm1q(&left.view(), right),
-        // 4 => mm2q(&left.view(), right, qargs.unwrap_or(&[Qubit(0), Qubit(1)])),
+        2 => mm1q(&left.view(), right),
+        4 => mm2q(&left.view(), right, qargs.unwrap_or(&[Qubit(0), Qubit(1)])),
         _ => left.dot(right),
     };
     let trace = product.diag().sum();

From 38e04b36f1f15073386f5a4e79dfa4e1aa36e217 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Wed, 5 Feb 2025 09:22:27 +0100
Subject: [PATCH 06/15] expose tol

---
 crates/accelerate/src/commutation_analysis.rs |   8 +-
 .../src/commutation_cancellation.rs           |   6 +-
 crates/accelerate/src/commutation_checker.rs  |  57 +++------
 crates/accelerate/src/unitary_compose.rs      | 110 +++---------------
 4 files changed, 41 insertions(+), 140 deletions(-)

diff --git a/crates/accelerate/src/commutation_analysis.rs b/crates/accelerate/src/commutation_analysis.rs
index c4fe033bfa39..57c63d2a12b6 100644
--- a/crates/accelerate/src/commutation_analysis.rs
+++ b/crates/accelerate/src/commutation_analysis.rs
@@ -54,6 +54,7 @@ pub(crate) fn analyze_commutations_inner(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
+    tol: Option<f64>,
 ) -> PyResult<(CommutationSet, NodeIndices)> {
     let mut commutation_set: CommutationSet = HashMap::new();
     let mut node_indices: NodeIndices = HashMap::new();
@@ -104,6 +105,7 @@ pub(crate) fn analyze_commutations_inner(
                             qargs2,
                             cargs2,
                             MAX_NUM_QUBITS,
+                            tol,
                         )?;
                         if !all_commute {
                             break;
@@ -134,17 +136,19 @@ pub(crate) fn analyze_commutations_inner(
 }
 
 #[pyfunction]
-#[pyo3(signature = (dag, commutation_checker))]
+#[pyo3(signature = (dag, commutation_checker, tol=None))]
 pub(crate) fn analyze_commutations(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
+    tol: Option<f64>,
 ) -> PyResult<Py<PyDict>> {
     // This returns two HashMaps:
     //   * The commuting nodes per wire: {wire: [commuting_nodes_1, commuting_nodes_2, ...]}
     //   * The index in which commutation set a given node is located on a wire: {(node, wire): index}
     // The Python dict will store both of these dictionaries in one.
-    let (commutation_set, node_indices) = analyze_commutations_inner(py, dag, commutation_checker)?;
+    let (commutation_set, node_indices) =
+        analyze_commutations_inner(py, dag, commutation_checker, tol)?;
 
     let out_dict = PyDict::new(py);
 
diff --git a/crates/accelerate/src/commutation_cancellation.rs b/crates/accelerate/src/commutation_cancellation.rs
index 7bd13c662a2e..50766c09b651 100644
--- a/crates/accelerate/src/commutation_cancellation.rs
+++ b/crates/accelerate/src/commutation_cancellation.rs
@@ -57,12 +57,13 @@ struct CancellationSetKey {
 }
 
 #[pyfunction]
-#[pyo3(signature = (dag, commutation_checker, basis_gates=None))]
+#[pyo3(signature = (dag, commutation_checker, basis_gates=None, tol=None))]
 pub(crate) fn cancel_commutations(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
     basis_gates: Option<HashSet<String>>,
+    tol: Option<f64>,
 ) -> PyResult<()> {
     let basis: HashSet<String> = if let Some(basis) = basis_gates {
         basis
@@ -97,7 +98,8 @@ pub(crate) fn cancel_commutations(
         sec_commutation_set_id), the value is the list gates that share the same gate type,
         qubits and commutation sets.
     */
-    let (commutation_set, node_indices) = analyze_commutations_inner(py, dag, commutation_checker)?;
+    let (commutation_set, node_indices) =
+        analyze_commutations_inner(py, dag, commutation_checker, tol)?;
     let mut cancellation_sets: HashMap<CancellationSetKey, Vec<NodeIndex>> = HashMap::new();
 
     (0..dag.num_qubits() as u32).for_each(|qubit| {
diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index 453c92f7c53f..c71c1e0350b8 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -130,13 +130,14 @@ impl CommutationChecker {
         }
     }
 
-    #[pyo3(signature=(op1, op2, max_num_qubits=3))]
+    #[pyo3(signature=(op1, op2, max_num_qubits=3, tol=None))]
     fn commute_nodes(
         &mut self,
         py: Python,
         op1: &DAGOpNode,
         op2: &DAGOpNode,
         max_num_qubits: u32,
+        tol: Option<f64>,
     ) -> PyResult<bool> {
         let (qargs1, qargs2) = get_bits::<Qubit>(
             py,
@@ -162,10 +163,11 @@ impl CommutationChecker {
             &qargs2,
             &cargs2,
             max_num_qubits,
+            tol,
         )
     }
 
-    #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3))]
+    #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3, tol=None))]
     #[allow(clippy::too_many_arguments)]
     fn commute(
         &mut self,
@@ -177,6 +179,7 @@ impl CommutationChecker {
         qargs2: Option<&Bound<PySequence>>,
         cargs2: Option<&Bound<PySequence>>,
         max_num_qubits: u32,
+        tol: Option<f64>,
     ) -> PyResult<bool> {
         let qargs1 = qargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?;
         let cargs1 = cargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?;
@@ -199,6 +202,7 @@ impl CommutationChecker {
             &qargs2,
             &cargs2,
             max_num_qubits,
+            tol,
         )
     }
 
@@ -269,21 +273,18 @@ impl CommutationChecker {
         qargs2: &[Qubit],
         cargs2: &[Clbit],
         max_num_qubits: u32,
+        tol: Option<f64>,
     ) -> PyResult<bool> {
-        // relative and absolute tolerance used to (1) check whether rotation gates commute
-        // trivially (i.e. the rotation angle is so small we assume it commutes) and (2) define
-        // comparison for the matrix-based commutation checks
-        // let rtol = 1e-5;
-        let rtol = f64::EPSILON;
-        let atol = 1e-8;
+        // if the average gate infidelity is below this tolerance, they commute
+        let tol = tol.unwrap_or(f64::EPSILON);
 
         // if we have rotation gates, we attempt to map them to their generators, for example
         // RX -> X or CPhase -> CZ
-        let (op1, params1, trivial1) = map_rotation(op1, params1, rtol);
+        let (op1, params1, trivial1) = map_rotation(op1, params1, tol);
         if trivial1 {
             return Ok(true);
         }
-        let (op2, params2, trivial2) = map_rotation(op2, params2, rtol);
+        let (op2, params2, trivial2) = map_rotation(op2, params2, tol);
         if trivial2 {
             return Ok(true);
         }
@@ -347,8 +348,7 @@ impl CommutationChecker {
                 second_op,
                 second_params,
                 second_qargs,
-                rtol,
-                atol,
+                tol,
             );
         }
 
@@ -383,8 +383,7 @@ impl CommutationChecker {
             second_op,
             second_params,
             second_qargs,
-            rtol,
-            atol,
+            tol,
         )?;
 
         // TODO: implement a LRU cache for this
@@ -419,8 +418,7 @@ impl CommutationChecker {
         second_op: &OperationRef,
         second_params: &[Param],
         second_qargs: &[Qubit],
-        rtol: f64,
-        atol: f64,
+        tol: f64,
     ) -> PyResult<bool> {
         // Compute relative positioning of qargs of the second gate to the first gate.
         // Since the qargs come out the same BitData, we already know there are no accidential
@@ -461,31 +459,6 @@ impl CommutationChecker {
             None => return Ok(false),
         };
 
-        // if first_qarg == second_qarg {
-        // match first_qarg.len() {
-        //     1 => Ok(unitary_compose::commute_1q(
-        //         &first_mat.view(),
-        //         &second_mat.view(),
-        //         rtol,
-        //         atol,
-        //     )),
-        //     2 => Ok(unitary_compose::commute_2q(
-        //         &first_mat.view(),
-        //         &second_mat.view(),
-        //         &[Qubit(0), Qubit(1)],
-        //         rtol,
-        //         atol,
-        //     )),
-        //     _ => Ok(unitary_compose::allclose(
-        //         &second_mat.dot(&first_mat).view(),
-        //         &first_mat.dot(&second_mat).view(),
-        //         rtol,
-        //         atol,
-        //     )),
-        // }
-        //     let fid = unitary_compose::gate_fidelity(&first_mat.view(), &second_mat.view());
-        //     Ok((1.0 - fid).abs() < rtol)
-        // } else {
         // TODO Optimize this bit to avoid unnecessary Kronecker products:
         //  1. We currently sort the operations for the cache by operation size, putting the
         //     *smaller* operation first: (smaller op, larger op)
@@ -522,7 +495,7 @@ impl CommutationChecker {
             Err(e) => return Err(PyRuntimeError::new_err(e)),
         };
         let fid = unitary_compose::gate_fidelity(&op12.view(), &op21.view(), None);
-        Ok((1.0 - fid).abs() <= rtol)
+        Ok((1.0 - fid).abs() <= tol)
     }
 
     fn clear_cache(&mut self) {
diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs
index b1d2eb1efaa2..20a681c4f55c 100644
--- a/crates/accelerate/src/unitary_compose.rs
+++ b/crates/accelerate/src/unitary_compose.rs
@@ -13,7 +13,7 @@
 use ndarray::{arr2, Array, Array2, ArrayView, ArrayView2, IxDyn};
 use ndarray_einsum_beta::*;
 use num_complex::{Complex, Complex64, ComplexFloat};
-use num_traits::Zero;
+use qiskit_circuit::util::C_ZERO;
 use qiskit_circuit::Qubit;
 
 static LOWERCASE: [u8; 26] = [
@@ -153,92 +153,6 @@ fn _einsum_matmul_index(qubits: &[u32], num_qubits: usize) -> String {
     )
 }
 
-pub fn commute_1q(
-    left: &ArrayView2<Complex64>,
-    right: &ArrayView2<Complex64>,
-    rtol: f64,
-    atol: f64,
-) -> bool {
-    // This could allow for explicit hardcoded formulas, using less FLOPS, if we only
-    // consider an absolute tolerance. But for backward compatibility we now implement the full
-    // formula including relative tolerance handling.
-    for i in 0..2usize {
-        for j in 0..2usize {
-            let mut ab = Complex64::zero();
-            let mut ba = Complex64::zero();
-            for k in 0..2usize {
-                ab += left[[i, k]] * right[[k, j]];
-                ba += right[[i, k]] * left[[k, j]];
-            }
-            let sum = ab - ba;
-            if sum.abs() > atol + ba.abs() * rtol {
-                return false;
-            }
-        }
-    }
-    true
-}
-
-pub fn commute_2q(
-    left: &ArrayView2<Complex64>,
-    right: &ArrayView2<Complex64>,
-    qargs: &[Qubit],
-    rtol: f64,
-    atol: f64,
-) -> bool {
-    let rev = qargs[0].0 == 1;
-    for i in 0..4usize {
-        for j in 0..4usize {
-            // We compute AB and BA separately, to enable checking the relative difference
-            // (AB - BA)_ij > atol + rtol * BA_ij. This is due to backward compatibility and could
-            // maybe be changed in the future to save one complex number allocation.
-            let mut ab = Complex64::zero();
-            let mut ba = Complex64::zero();
-            for k in 0..4usize {
-                ab += left[[_ind(i, rev), _ind(k, rev)]] * right[[k, j]];
-                ba += right[[i, k]] * left[[_ind(k, rev), _ind(j, rev)]];
-            }
-            let sum = ab - ba;
-            if sum.abs() > atol + ba.abs() * rtol {
-                return false;
-            }
-        }
-    }
-    true
-}
-
-#[inline]
-fn _ind(i: usize, reversed: bool) -> usize {
-    if reversed {
-        // reverse the first two bits
-        ((i & 1) << 1) + ((i & 2) >> 1)
-    } else {
-        i
-    }
-}
-
-/// For equally sized matrices, ``left`` and ``right``, check whether all entries are close
-/// by the criterion
-///
-///     |left_ij - right_ij| <= atol + rtol * right_ij
-///
-/// This is analogous to NumPy's ``allclose`` function.
-pub fn allclose(
-    left: &ArrayView2<Complex64>,
-    right: &ArrayView2<Complex64>,
-    rtol: f64,
-    atol: f64,
-) -> bool {
-    for i in 0..left.nrows() {
-        for j in 0..left.ncols() {
-            if (left[(i, j)] - right[(i, j)]).abs() > atol + rtol * right[(i, j)].abs() {
-                return false;
-            }
-        }
-    }
-    true
-}
-
 pub fn gate_fidelity(
     left: &ArrayView2<Complex64>,
     right: &ArrayView2<Complex64>,
@@ -262,8 +176,7 @@ pub fn gate_fidelity(
 }
 
 fn mm1q(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> Array2<Complex64> {
-    let zero = Complex64::zero();
-    let mut out = arr2(&[[zero, zero], [zero, zero]]);
+    let mut out = arr2(&[[C_ZERO, C_ZERO], [C_ZERO, C_ZERO]]);
     out[[0, 0]] = left[[0, 0]] * right[[0, 0]] + left[[0, 1]] * right[[1, 0]];
     out[[0, 1]] = left[[0, 0]] * right[[0, 1]] + left[[0, 1]] * right[[1, 1]];
     out[[1, 0]] = left[[1, 0]] * right[[0, 0]] + left[[1, 1]] * right[[1, 0]];
@@ -271,17 +184,26 @@ fn mm1q(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> Array2<C
     out
 }
 
+#[inline]
+fn _ind(i: usize, reversed: bool) -> usize {
+    if reversed {
+        // reverse the first two bits
+        ((i & 1) << 1) + ((i & 2) >> 1)
+    } else {
+        i
+    }
+}
+
 pub fn mm2q(
     left: &ArrayView2<Complex64>,
     right: &ArrayView2<Complex64>,
     qargs: &[Qubit],
 ) -> Array2<Complex64> {
-    let zero = Complex64::zero();
     let mut out = arr2(&[
-        [zero, zero, zero, zero],
-        [zero, zero, zero, zero],
-        [zero, zero, zero, zero],
-        [zero, zero, zero, zero],
+        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
+        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
+        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
+        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
     ]);
 
     let rev = qargs[0].0 == 1;

From c1efb19b63b5633f33b4924eeecc01e1d748d444 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Wed, 5 Feb 2025 14:18:15 +0100
Subject: [PATCH 07/15] rename tol to approximation_degree

---
 crates/accelerate/src/commutation_analysis.rs | 10 +--
 crates/accelerate/src/commutation_checker.rs  | 62 +++++--------------
 2 files changed, 22 insertions(+), 50 deletions(-)

diff --git a/crates/accelerate/src/commutation_analysis.rs b/crates/accelerate/src/commutation_analysis.rs
index 57c63d2a12b6..4aebcb50b383 100644
--- a/crates/accelerate/src/commutation_analysis.rs
+++ b/crates/accelerate/src/commutation_analysis.rs
@@ -54,7 +54,7 @@ pub(crate) fn analyze_commutations_inner(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
-    tol: Option<f64>,
+    approximation_degree: Option<f64>,
 ) -> PyResult<(CommutationSet, NodeIndices)> {
     let mut commutation_set: CommutationSet = HashMap::new();
     let mut node_indices: NodeIndices = HashMap::new();
@@ -105,7 +105,7 @@ pub(crate) fn analyze_commutations_inner(
                             qargs2,
                             cargs2,
                             MAX_NUM_QUBITS,
-                            tol,
+                            approximation_degree,
                         )?;
                         if !all_commute {
                             break;
@@ -136,19 +136,19 @@ pub(crate) fn analyze_commutations_inner(
 }
 
 #[pyfunction]
-#[pyo3(signature = (dag, commutation_checker, tol=None))]
+#[pyo3(signature = (dag, commutation_checker, approximation_degree=None))]
 pub(crate) fn analyze_commutations(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
-    tol: Option<f64>,
+    approximation_degree: Option<f64>,
 ) -> PyResult<Py<PyDict>> {
     // This returns two HashMaps:
     //   * The commuting nodes per wire: {wire: [commuting_nodes_1, commuting_nodes_2, ...]}
     //   * The index in which commutation set a given node is located on a wire: {(node, wire): index}
     // The Python dict will store both of these dictionaries in one.
     let (commutation_set, node_indices) =
-        analyze_commutations_inner(py, dag, commutation_checker, tol)?;
+        analyze_commutations_inner(py, dag, commutation_checker, approximation_degree)?;
 
     let out_dict = PyDict::new(py);
 
diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index c71c1e0350b8..d8170d487757 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -10,8 +10,6 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use core::f64;
-
 use hashbrown::{HashMap, HashSet};
 use ndarray::linalg::kron;
 use ndarray::Array2;
@@ -130,14 +128,14 @@ impl CommutationChecker {
         }
     }
 
-    #[pyo3(signature=(op1, op2, max_num_qubits=3, tol=None))]
+    #[pyo3(signature=(op1, op2, max_num_qubits=3, approximation_degree=None))]
     fn commute_nodes(
         &mut self,
         py: Python,
         op1: &DAGOpNode,
         op2: &DAGOpNode,
         max_num_qubits: u32,
-        tol: Option<f64>,
+        approximation_degree: Option<f64>,
     ) -> PyResult<bool> {
         let (qargs1, qargs2) = get_bits::<Qubit>(
             py,
@@ -163,11 +161,11 @@ impl CommutationChecker {
             &qargs2,
             &cargs2,
             max_num_qubits,
-            tol,
+            approximation_degree,
         )
     }
 
-    #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3, tol=None))]
+    #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3, approximation_degree=None))]
     #[allow(clippy::too_many_arguments)]
     fn commute(
         &mut self,
@@ -179,7 +177,7 @@ impl CommutationChecker {
         qargs2: Option<&Bound<PySequence>>,
         cargs2: Option<&Bound<PySequence>>,
         max_num_qubits: u32,
-        tol: Option<f64>,
+        approximation_degree: Option<f64>,
     ) -> PyResult<bool> {
         let qargs1 = qargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?;
         let cargs1 = cargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?;
@@ -202,7 +200,7 @@ impl CommutationChecker {
             &qargs2,
             &cargs2,
             max_num_qubits,
-            tol,
+            approximation_degree,
         )
     }
 
@@ -273,10 +271,11 @@ impl CommutationChecker {
         qargs2: &[Qubit],
         cargs2: &[Clbit],
         max_num_qubits: u32,
-        tol: Option<f64>,
+        approximation_degree: Option<f64>,
     ) -> PyResult<bool> {
-        // if the average gate infidelity is below this tolerance, they commute
-        let tol = tol.unwrap_or(f64::EPSILON);
+        // If the average gate infidelity is below this tolerance, they commute. The tolerance
+        // is set to max(machine_eps, 1 - approximation_degree).
+        let tol = f64::EPSILON.max(1. - approximation_degree.unwrap_or(1.));
 
         // if we have rotation gates, we attempt to map them to their generators, for example
         // RX -> X or CPhase -> CZ
@@ -604,7 +603,7 @@ fn map_rotation<'a>(
         // it commutes trivially.
         if let Param::Float(angle) = params[0] {
             let gate_fidelity = rotation_fidelity(name, angle).unwrap_or(0.);
-            if (1. - gate_fidelity).abs() < tol {
+            if (1. - gate_fidelity).abs() <= tol {
                 return (op, params, true);
             };
         };
@@ -622,45 +621,18 @@ fn map_rotation<'a>(
     (op, params, false)
 }
 
-/// Compute the gate fidelity of a rotation gate.
-// fn rotation_fidelity(rotation: &OperationRef, angle: f64) -> Option<f64> {
-//     if let OperationRef::Standard(gate) = rotation {
-//         match gate {
-//             StandardGate::RXGate
-//             | StandardGate::RYGate
-//             | StandardGate::RZGate
-//             | StandardGate::PhaseGate => return Some(2. * (angle / 2.).cos()),
-//             StandardGate::RXXGate
-//             | StandardGate::RYYGate
-//             | StandardGate::RZXGate
-//             | StandardGate::RZZGate => return Some(4. * (angle / 2.).cos()),
-//             StandardGate::CRXGate
-//             | StandardGate::CRYGate
-//             | StandardGate::CRZGate
-//             | StandardGate::CPhaseGate => return Some(2. + 2. * (angle / 2.).cos()),
-//             _ => return None,
-//         };
-//     };
-//     None
-// }
-
 fn rotation_fidelity(rotation: &str, angle: f64) -> Option<f64> {
     let dim = match rotation {
         "rx" | "ry" | "rz" | "p" => 2.,
         _ => 4.,
     };
 
-    match rotation {
-        "rx" | "ry" | "rz" | "p" | "rxx" | "ryy" | "rzx" | "rzz" => {
-            let gate_fid = (angle / 2.).cos().powi(2);
-            Some((dim * gate_fid + 1.) / (dim + 1.))
-        }
-        "crx" | "cry" | "crz" | "cp" => {
-            let gate_fid = (0.5 + 0.5 * (angle / 2.).cos()).powi(2);
-            Some((dim * gate_fid + 1.) / (dim + 1.))
-        }
-        _ => None,
-    }
+    let gate_fid = match rotation {
+        "rx" | "ry" | "rz" | "p" | "rxx" | "ryy" | "rzx" | "rzz" => (angle / 2.).cos().powi(2),
+        "crx" | "cry" | "crz" | "cp" => (0.5 + 0.5 * (angle / 2.).cos()).powi(2),
+        _ => return None,
+    };
+    Some((dim * gate_fid + 1.) / (dim + 1.))
 }
 
 fn get_relative_placement(

From 31164835c56326d4788e236a99fd22bff306b909 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Thu, 6 Feb 2025 11:25:43 +0100
Subject: [PATCH 08/15] note: this is not global phase correct!

---
 crates/accelerate/src/commutation_checker.rs  |  3 +-
 qiskit/circuit/commutation_checker.py         | 11 +++-
 .../passes/synthesis/high_level_synthesis.py  | 10 ++--
 ...pported-instructions-0d80ea33b3d2257b.yaml |  8 +++
 .../circuit/test_commutation_checker.py       | 52 +++++++------------
 .../transpiler/test_high_level_synthesis.py   | 13 +++++
 6 files changed, 56 insertions(+), 41 deletions(-)
 create mode 100644 releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml

diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index d8170d487757..550e3c1823e5 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -629,7 +629,8 @@ fn rotation_fidelity(rotation: &str, angle: f64) -> Option<f64> {
 
     let gate_fid = match rotation {
         "rx" | "ry" | "rz" | "p" | "rxx" | "ryy" | "rzx" | "rzz" => (angle / 2.).cos().powi(2),
-        "crx" | "cry" | "crz" | "cp" => (0.5 + 0.5 * (angle / 2.).cos()).powi(2),
+        "crx" | "cry" | "crz" => (0.5 + 0.5 * (angle / 2.).cos()).powi(2),
+        "cp" => (10. + 6. * angle.cos()) / 16.,
         _ => return None,
     };
     Some((dim * gate_fid + 1.) / (dim + 1.))
diff --git a/qiskit/circuit/commutation_checker.py b/qiskit/circuit/commutation_checker.py
index 2d474c7bac4e..994867bcf1b4 100644
--- a/qiskit/circuit/commutation_checker.py
+++ b/qiskit/circuit/commutation_checker.py
@@ -12,6 +12,7 @@
 
 """Code from commutative_analysis pass that checks commutation relations between DAG nodes."""
 
+from __future__ import annotations
 from typing import List, Union, Set, Optional
 
 from qiskit.circuit.operation import Operation
@@ -40,9 +41,10 @@ def commute_nodes(
         op1,
         op2,
         max_num_qubits: int = 3,
+        approximation_degree: float | None = None,
     ) -> bool:
         """Checks if two DAGOpNodes commute."""
-        return self.cc.commute_nodes(op1, op2, max_num_qubits)
+        return self.cc.commute_nodes(op1, op2, max_num_qubits, approximation_degree)
 
     def commute(
         self,
@@ -53,6 +55,7 @@ def commute(
         qargs2: List,
         cargs2: List,
         max_num_qubits: int = 3,
+        approximation_degree: float | None = None,
     ) -> bool:
         """
         Checks if two Operations commute. The return value of `True` means that the operations
@@ -69,11 +72,15 @@ def commute(
             cargs2: second operation's clbits.
             max_num_qubits: the maximum number of qubits to consider, the check may be skipped if
                 the number of qubits for either operation exceeds this amount.
+            approximation_degree: If the average gate fidelity in between the two operations
+                is above this number (up to machine epsilon) they are assumed to commute.
 
         Returns:
             bool: whether two operations commute.
         """
-        return self.cc.commute(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits)
+        return self.cc.commute(
+            op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits, approximation_degree
+        )
 
     def num_cached_entries(self):
         """Returns number of cached entries"""
diff --git a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
index d7b5b0d535b8..95adee1d05c7 100644
--- a/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
+++ b/qiskit/transpiler/passes/synthesis/high_level_synthesis.py
@@ -814,6 +814,7 @@ def _definitely_skip_node(
             dag._has_calibration_for(node)
             or len(node.qargs) < self._min_qubits
             or node.is_directive()
+            or (self._instruction_supported(node.name, qubits) and not node.is_control_flow())
         ):
             return True
 
@@ -831,15 +832,12 @@ def _definitely_skip_node(
             # If all the above constraints hold, and it's already supported or the basis translator
             # can handle it, we'll leave it be.
             and (
-                self._instruction_supported(node.name, qubits)
                 # This uses unfortunately private details of `EquivalenceLibrary`, but so does the
                 # `BasisTranslator`, and this is supposed to just be temporary til this is moved
                 # into Rust space.
-                or (
-                    self._equiv_lib is not None
-                    and equivalence.Key(name=node.name, num_qubits=node.num_qubits)
-                    in self._equiv_lib.keys()
-                )
+                self._equiv_lib is not None
+                and equivalence.Key(name=node.name, num_qubits=node.num_qubits)
+                in self._equiv_lib.keys()
             )
         )
 
diff --git a/releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml b/releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml
new file mode 100644
index 000000000000..1e9c1a1cebe7
--- /dev/null
+++ b/releasenotes/notes/fix-hls-supported-instructions-0d80ea33b3d2257b.yaml
@@ -0,0 +1,8 @@
+---
+fixes:
+  - |
+    Previously the :class:`.HighLevelSynthesis` transpiler pass synthesized an
+    instruction for which a synthesis plugin is available, regardless of
+    whether the instruction is already supported by the target or a part of
+    the explicitly passed ``basis_gates``. This behavior is now fixed, so that
+    such already supported instructions are no longer synthesized.
diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py
index ff6dd0af73f1..b3e58dffade9 100644
--- a/test/python/circuit/test_commutation_checker.py
+++ b/test/python/circuit/test_commutation_checker.py
@@ -75,8 +75,6 @@
     CPhaseGate,
 ]
 
-CONTROLLED_ROTATION_GATES = [CRXGate, CRYGate, CRZGate, CPhaseGate]
-
 
 class NewGateCX(Gate):
     """A dummy class containing an cx gate unknown to the commutation checker's library."""
@@ -415,7 +413,7 @@ def test_serialization(self):
         cc2.commute_nodes(dop1, dop2)
         self.assertEqual(cc2.num_cached_entries(), 1)
 
-    @idata(ROTATION_GATES + CONTROLLED_ROTATION_GATES)
+    @idata(ROTATION_GATES)
     def test_cutoff_angles(self, gate_cls):
         """Check rotations with a small enough angle are cut off."""
         max_power = 30
@@ -423,22 +421,25 @@ def test_cutoff_angles(self, gate_cls):
 
         generic_gate = DCXGate()  # gate that does not commute with any rotation gate
 
-        # TODO this is no longer a cutoff but defined by the gate fidelity
-        cutoff_angle = 1e-5  # this is the cutoff we use in the CommutationChecker
+        # the cutoff angle depends on the average gate fidelity; i.e. it is the angle
+        # for which the average gate fidelity is smaller than machine epsilon
+        if gate_cls in [CPhaseGate, CRXGate, CRYGate, CRZGate]:
+            cutoff_angle = 4.71e-8
+        else:
+            cutoff_angle = 3.65e-8
 
         for i in range(1, max_power + 1):
             angle = 2 ** (-i)
             gate = gate_cls(angle)
             qargs = list(range(gate.num_qubits))
-            print(f"{angle:.2e}", scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
-            # if angle < cutoff_angle:
-            #     self.assertTrue(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
-            # else:
-            #     self.assertFalse(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
+            if angle < cutoff_angle:
+                self.assertTrue(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
+            else:
+                self.assertFalse(scc.commute(generic_gate, [0, 1], [], gate, qargs, []))
 
     @idata(ROTATION_GATES)
-    def test_controlled_rotation_mod_4pi(self, gate_cls):
-        """Test the rotations modulo 2pi (4pi for controlled-rx/y/z) commute with any gate."""
+    def test_rotations_pi_multiples(self, gate_cls):
+        """Test the rotations modulo 2pi (crx/cry/crz modulo 4pi) commute with any gate."""
         generic_gate = HGate()  # does not commute with any rotation gate
         multiples = np.arange(-6, 7)
 
@@ -450,31 +451,18 @@ def test_controlled_rotation_mod_4pi(self, gate_cls):
                 # compute a numeric reference, that doesn't go through any special cases and
                 # uses a matrix-based commutation check
                 expected = scc.commute(
-                    generic_gate, [0], [], numeric, list(range(gate.num_qubits)), []
+                    generic_gate,
+                    [0],
+                    [],
+                    numeric,
+                    list(range(gate.num_qubits)),
+                    [],
+                    approximation_degree=1 - 1e-5,
                 )
 
                 result = scc.commute(generic_gate, [0], [], gate, list(range(gate.num_qubits)), [])
                 self.assertEqual(expected, result)
 
-    @idata(CONTROLLED_ROTATION_GATES)
-    def test_controlled_rotation_mod_4pi(self, gate_cls):
-        """Test the controlled rotations modulo 4pi commute with any gate.
-
-        For 2pi, the single-qubit rotations are the identity only up to a global phase,
-        which has a measurable effect if controlled.
-        """
-        generic_gate = HGate()  # does not commute with any rotation gate
-        multiples = np.arange(-6, 7)
-
-        for multiple in multiples:
-            with self.subTest(multiple=multiple):
-                gate = gate_cls(multiple * np.pi)
-                comm = scc.commute(generic_gate, [0], [], gate, list(range(gate.num_qubits)), [])
-                if multiple % 4 == 0:
-                    self.assertTrue(comm)
-                else:
-                    self.assertFalse(comm)
-
     def test_custom_gate(self):
         """Test a custom gate."""
         my_cx = NewGateCX()
diff --git a/test/python/transpiler/test_high_level_synthesis.py b/test/python/transpiler/test_high_level_synthesis.py
index e1fa336c1377..b7bc31132a48 100644
--- a/test/python/transpiler/test_high_level_synthesis.py
+++ b/test/python/transpiler/test_high_level_synthesis.py
@@ -46,6 +46,7 @@
     IGate,
     MCXGate,
     SGate,
+    QAOAAnsatz,
 )
 from qiskit.circuit.library import LinearFunction, PauliEvolutionGate
 from qiskit.quantum_info import Clifford, Operator, Statevector, SparsePauliOp
@@ -664,6 +665,18 @@ def test_synth_fails_definition_exists(self):
             out = hls(circuit)
             self.assertEqual(out.count_ops(), {"u": 1})
 
+    def test_both_basis_gates_and_plugin_specified(self):
+        """Test that a gate is not synthesized when it belongs to basis_gates,
+        regardless of whether there is a plugin method available.
+
+        See: https://github.com/Qiskit/qiskit/issues/13412 for more
+        details.
+        """
+        qc = QAOAAnsatz(SparsePauliOp("Z"), initial_state=QuantumCircuit(1))
+        pm = PassManager([HighLevelSynthesis(basis_gates=["PauliEvolution"])])
+        qct = pm.run(qc)
+        self.assertEqual(qct.count_ops()["PauliEvolution"], 2)
+
 
 class TestPMHSynthesisLinearFunctionPlugin(QiskitTestCase):
     """Tests for the PMHSynthesisLinearFunction plugin for synthesizing linear functions."""

From c2d77a2eee3523134df0d0c13ac8b12a9460a714 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Tue, 18 Feb 2025 17:51:12 +0100
Subject: [PATCH 09/15] don't be phase agnostic

---
 crates/accelerate/src/commutation_checker.rs | 7 +++++--
 crates/accelerate/src/unitary_compose.rs     | 9 ++++++---
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index 550e3c1823e5..bd2421619f43 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -493,8 +493,11 @@ impl CommutationChecker {
             Ok(matrix) => matrix,
             Err(e) => return Err(PyRuntimeError::new_err(e)),
         };
-        let fid = unitary_compose::gate_fidelity(&op12.view(), &op21.view(), None);
-        Ok((1.0 - fid).abs() <= tol)
+        let (fid, phase) = unitary_compose::gate_fidelity(&op12.view(), &op21.view(), None);
+
+        // we consider the gates as commuting if the process fidelity of
+        // AB (BA)^\dagger is approximately the identity and there is no global phase difference
+        Ok(phase.abs() <= tol && (1.0 - fid).abs() <= tol)
     }
 
     fn clear_cache(&mut self) {
diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs
index 20a681c4f55c..af756803a57f 100644
--- a/crates/accelerate/src/unitary_compose.rs
+++ b/crates/accelerate/src/unitary_compose.rs
@@ -157,7 +157,7 @@ pub fn gate_fidelity(
     left: &ArrayView2<Complex64>,
     right: &ArrayView2<Complex64>,
     qargs: Option<&[Qubit]>,
-) -> f64 {
+) -> (f64, f64) {
     let dim = left.nrows(); // == left.ncols() == right.nrows() == right.ncols()
                             // let trace = left.t().mapv(|el| el.conj()).dot(right).diag().sum();
 
@@ -170,9 +170,12 @@ pub fn gate_fidelity(
     let trace = product.diag().sum();
 
     let dim = dim as f64;
-    let process_fidelity = (trace / dim).abs().powi(2);
+    let normalized_trace = trace / dim;
+    let phase = normalized_trace.arg(); // compute phase difference
+
+    let process_fidelity = normalized_trace.abs().powi(2);
     let gate_fidelity = (dim * process_fidelity + 1.) / (dim + 1.);
-    gate_fidelity
+    (gate_fidelity, phase)
 }
 
 fn mm1q(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> Array2<Complex64> {

From afc4ed45b93ddd09b9321c3675658124edfd3a35 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 24 Feb 2025 14:08:44 +0100
Subject: [PATCH 10/15] additional tests

these cover cases where we accumulate roundoff errors
---
 test/python/circuit/test_commutation_checker.py      | 12 ++++++++++++
 .../transpiler/test_remove_identity_equivalent.py    |  3 +++
 2 files changed, 15 insertions(+)

diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py
index 3f32f7d65a4a..963413dc1bc6 100644
--- a/test/python/circuit/test_commutation_checker.py
+++ b/test/python/circuit/test_commutation_checker.py
@@ -142,6 +142,18 @@ def test_simple_gates(self):
         self.assertTrue(scc.commute(XGate(), [2], [], CCXGate(), [0, 1, 2], []))
         self.assertFalse(scc.commute(CCXGate(), [0, 1, 2], [], CCXGate(), [0, 2, 1], []))
 
+        self.assertTrue(scc.commute(HGate(), [0], [], HGate(), [0], []))
+
+    def test_simple_matrices(self):
+        """Test simple gates but matrix-based."""
+        had = UnitaryGate(HGate())
+        cx = UnitaryGate(CXGate())
+        x = UnitaryGate(XGate())
+        self.assertTrue(scc.commute(x, [0], [], x, [0], []))
+        self.assertTrue(scc.commute(had, [0], [], had, [0], []))
+        self.assertTrue(scc.commute(cx, [0, 1], [], cx, [0, 1], []))
+        self.assertTrue(scc.commute(x, [1], [], cx, [0, 1], []))
+
     def test_passing_quantum_registers(self):
         """Check that passing QuantumRegisters works correctly."""
         qr = QuantumRegister(4)
diff --git a/test/python/transpiler/test_remove_identity_equivalent.py b/test/python/transpiler/test_remove_identity_equivalent.py
index 736c02b765bd..495978ef585a 100644
--- a/test/python/transpiler/test_remove_identity_equivalent.py
+++ b/test/python/transpiler/test_remove_identity_equivalent.py
@@ -188,6 +188,9 @@ def to_matrix(self):
         UnitaryGate(np.array([[np.exp(1j * np.pi / 4), 0], [0, np.exp(1j * np.pi / 4)]])),
         GlobalPhaseGate(0),
         GlobalPhaseGate(np.pi / 4),
+        UnitaryGate(np.exp(-0.123j) * np.eye(2)),
+        UnitaryGate(np.exp(-0.123j) * np.eye(4)),
+        UnitaryGate(np.exp(-0.123j) * np.eye(8)),
     )
     def test_remove_identity_up_to_global_phase(self, gate):
         """Test that gates equivalent to identity up to a global phase are removed from the circuit,

From adbd9be24e1d9e7b82bbdc0d70501e20f128260c Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 3 Mar 2025 08:36:36 +0100
Subject: [PATCH 11/15] Increase threshold, better consistency

- increased the threshold for value comparisons to 16 EPS
- use same functions in RemoveIdentityEquiv and CC checker
- add reno
- add detailed docs
- add some tests for above thresholds
---
 crates/accelerate/src/commutation_analysis.rs |  6 +-
 .../src/commutation_cancellation.rs           |  6 +-
 crates/accelerate/src/commutation_checker.rs  | 47 +++++++++------
 .../accelerate/src/remove_identity_equiv.rs   | 59 +++++++++++--------
 qiskit/circuit/__init__.py                    | 10 +---
 .../optimization/commutation_analysis.py      | 26 +++++++-
 .../optimization/remove_identity_equiv.py     | 10 ++--
 .../cc-gate-fidelity-bde7a01dc8c56e29.yaml    | 24 ++++++++
 .../circuit/test_commutation_checker.py       | 53 +++++++++++------
 .../test_remove_identity_equivalent.py        |  8 +--
 10 files changed, 162 insertions(+), 87 deletions(-)
 create mode 100644 releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml

diff --git a/crates/accelerate/src/commutation_analysis.rs b/crates/accelerate/src/commutation_analysis.rs
index 4aebcb50b383..857409dc7c46 100644
--- a/crates/accelerate/src/commutation_analysis.rs
+++ b/crates/accelerate/src/commutation_analysis.rs
@@ -54,7 +54,7 @@ pub(crate) fn analyze_commutations_inner(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
-    approximation_degree: Option<f64>,
+    approximation_degree: f64,
 ) -> PyResult<(CommutationSet, NodeIndices)> {
     let mut commutation_set: CommutationSet = HashMap::new();
     let mut node_indices: NodeIndices = HashMap::new();
@@ -136,12 +136,12 @@ pub(crate) fn analyze_commutations_inner(
 }
 
 #[pyfunction]
-#[pyo3(signature = (dag, commutation_checker, approximation_degree=None))]
+#[pyo3(signature = (dag, commutation_checker, approximation_degree=1.))]
 pub(crate) fn analyze_commutations(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
-    approximation_degree: Option<f64>,
+    approximation_degree: f64,
 ) -> PyResult<Py<PyDict>> {
     // This returns two HashMaps:
     //   * The commuting nodes per wire: {wire: [commuting_nodes_1, commuting_nodes_2, ...]}
diff --git a/crates/accelerate/src/commutation_cancellation.rs b/crates/accelerate/src/commutation_cancellation.rs
index 50766c09b651..45856f132632 100644
--- a/crates/accelerate/src/commutation_cancellation.rs
+++ b/crates/accelerate/src/commutation_cancellation.rs
@@ -57,13 +57,13 @@ struct CancellationSetKey {
 }
 
 #[pyfunction]
-#[pyo3(signature = (dag, commutation_checker, basis_gates=None, tol=None))]
+#[pyo3(signature = (dag, commutation_checker, basis_gates=None, approximation_degree=1.))]
 pub(crate) fn cancel_commutations(
     py: Python,
     dag: &mut DAGCircuit,
     commutation_checker: &mut CommutationChecker,
     basis_gates: Option<HashSet<String>>,
-    tol: Option<f64>,
+    approximation_degree: f64,
 ) -> PyResult<()> {
     let basis: HashSet<String> = if let Some(basis) = basis_gates {
         basis
@@ -99,7 +99,7 @@ pub(crate) fn cancel_commutations(
         qubits and commutation sets.
     */
     let (commutation_set, node_indices) =
-        analyze_commutations_inner(py, dag, commutation_checker, tol)?;
+        analyze_commutations_inner(py, dag, commutation_checker, approximation_degree)?;
     let mut cancellation_sets: HashMap<CancellationSetKey, Vec<NodeIndex>> = HashMap::new();
 
     (0..dag.num_qubits() as u32).for_each(|qubit| {
diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index 27f324d656ae..86f7d1d6be6d 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -14,6 +14,7 @@ use hashbrown::{HashMap, HashSet};
 use ndarray::linalg::kron;
 use ndarray::Array2;
 use num_complex::Complex64;
+use num_complex::ComplexFloat;
 use once_cell::sync::Lazy;
 use smallvec::SmallVec;
 
@@ -137,14 +138,14 @@ impl CommutationChecker {
         }
     }
 
-    #[pyo3(signature=(op1, op2, max_num_qubits=3, approximation_degree=None))]
+    #[pyo3(signature=(op1, op2, max_num_qubits=3, approximation_degree=1.))]
     fn commute_nodes(
         &mut self,
         py: Python,
         op1: &DAGOpNode,
         op2: &DAGOpNode,
         max_num_qubits: u32,
-        approximation_degree: Option<f64>,
+        approximation_degree: f64,
     ) -> PyResult<bool> {
         let (qargs1, qargs2) = get_bits::<Qubit>(
             py,
@@ -174,7 +175,7 @@ impl CommutationChecker {
         )
     }
 
-    #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3, approximation_degree=None))]
+    #[pyo3(signature=(op1, qargs1, cargs1, op2, qargs2, cargs2, max_num_qubits=3, approximation_degree=1.))]
     #[allow(clippy::too_many_arguments)]
     fn commute(
         &mut self,
@@ -186,7 +187,7 @@ impl CommutationChecker {
         qargs2: Option<&Bound<PySequence>>,
         cargs2: Option<&Bound<PySequence>>,
         max_num_qubits: u32,
-        approximation_degree: Option<f64>,
+        approximation_degree: f64,
     ) -> PyResult<bool> {
         let qargs1 = qargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?;
         let cargs1 = cargs1.map_or_else(|| Ok(PyTuple::empty(py)), PySequenceMethods::to_tuple)?;
@@ -280,11 +281,12 @@ impl CommutationChecker {
         qargs2: &[Qubit],
         cargs2: &[Clbit],
         max_num_qubits: u32,
-        approximation_degree: Option<f64>,
+        approximation_degree: f64,
     ) -> PyResult<bool> {
         // If the average gate infidelity is below this tolerance, they commute. The tolerance
-        // is set to max(machine_eps, 1 - approximation_degree).
-        let tol = f64::EPSILON.max(1. - approximation_degree.unwrap_or(1.));
+        // is set to max(16 * machine_eps, 1 - approximation_degree), where the heuristic factor
+        // of 16 is set to account for round-off errors in matrix multiplication.
+        let tol = (16. * f64::EPSILON).max(1. - approximation_degree);
 
         // if we have rotation gates, we attempt to map them to their generators, for example
         // RX -> X or CPhase -> CZ
@@ -510,7 +512,10 @@ impl CommutationChecker {
 
         // we consider the gates as commuting if the process fidelity of
         // AB (BA)^\dagger is approximately the identity and there is no global phase difference
-        Ok(phase.abs() <= tol && (1.0 - fid).abs() <= tol)
+        // let dim = op12.ncols() as f64;
+        // let matrix_tol = tol * dim.powi(2);
+        let matrix_tol = tol;
+        Ok(phase.abs() <= tol && (1.0 - fid).abs() <= matrix_tol)
     }
 
     fn clear_cache(&mut self) {
@@ -618,8 +623,11 @@ fn map_rotation<'a>(
         // commute with everything, and we simply return the operation with the flag that
         // it commutes trivially.
         if let Param::Float(angle) = params[0] {
-            let gate_fidelity = rotation_fidelity(name, angle).unwrap_or(0.);
-            if (1. - gate_fidelity).abs() <= tol {
+            let (tr_over_dim, dim) = rotation_trace_and_dim(name, angle)
+                .expect("All rotation should be covered at this point");
+            let gate_fidelity = tr_over_dim.abs().powi(2);
+            let process_fidelity = (dim * gate_fidelity + 1.) / (dim + 1.);
+            if (1. - process_fidelity).abs() <= tol {
                 return (op, params, true);
             };
         };
@@ -636,19 +644,24 @@ fn map_rotation<'a>(
     (op, params, false)
 }
 
-fn rotation_fidelity(rotation: &str, angle: f64) -> Option<f64> {
+/// For a (controlled) rotation or phase gate, return a tuple ``(Tr(gate) / dim, dim)``.
+/// Returns ``None`` if the rotation gate (specified by name) is not supported.
+pub fn rotation_trace_and_dim(rotation: &str, angle: f64) -> Option<(Complex64, f64)> {
     let dim = match rotation {
-        "rx" | "ry" | "rz" | "p" => 2.,
+        "rx" | "ry" | "rz" | "p" | "u1" => 2.,
         _ => 4.,
     };
 
-    let gate_fid = match rotation {
-        "rx" | "ry" | "rz" | "p" | "rxx" | "ryy" | "rzx" | "rzz" => (angle / 2.).cos().powi(2),
-        "crx" | "cry" | "crz" => (0.5 + 0.5 * (angle / 2.).cos()).powi(2),
-        "cp" => (10. + 6. * angle.cos()) / 16.,
+    let trace_over_dim = match rotation {
+        "rx" | "ry" | "rz" | "rxx" | "ryy" | "rzx" | "rzz" => {
+            Complex64::new((angle / 2.).cos(), 0.)
+        }
+        "crx" | "cry" | "crz" => Complex64::new(0.5 + 0.5 * (angle / 2.).cos(), 0.),
+        "p" | "u1" => (1. + Complex64::new(0., angle).exp()) / 2.,
+        "cp" => (3. + Complex64::new(0., angle).exp()) / 4.,
         _ => return None,
     };
-    Some((dim * gate_fid + 1.) / (dim + 1.))
+    Some((trace_over_dim, dim))
 }
 
 fn get_relative_placement(
diff --git a/crates/accelerate/src/remove_identity_equiv.rs b/crates/accelerate/src/remove_identity_equiv.rs
index 7a65a9972019..5b8ab2610b59 100644
--- a/crates/accelerate/src/remove_identity_equiv.rs
+++ b/crates/accelerate/src/remove_identity_equiv.rs
@@ -14,6 +14,7 @@ use num_complex::ComplexFloat;
 use pyo3::prelude::*;
 use rustworkx_core::petgraph::stable_graph::NodeIndex;
 
+use crate::commutation_checker::rotation_trace_and_dim;
 use crate::nlayout::PhysicalQubit;
 use crate::target_transpiler::Target;
 use qiskit_circuit::dag_circuit::DAGCircuit;
@@ -33,12 +34,16 @@ fn remove_identity_equiv(
 ) {
     let mut remove_list: Vec<NodeIndex> = Vec::new();
     let mut global_phase_update: f64 = 0.;
+    // Minimum threshold to compare average gate fidelity to 1. Includes a heuristic factor
+    // of 16 to account for round-off errors in the trace calculations, and for consistency
+    // with the commutation checker.
+    let minimum_tol = 16. * f64::EPSILON;
 
     let get_error_cutoff = |inst: &PackedInstruction| -> f64 {
         match approx_degree {
             Some(degree) => {
                 if degree == 1.0 {
-                    f64::EPSILON
+                    minimum_tol
                 } else {
                     match target {
                         Some(target) => {
@@ -50,10 +55,10 @@ fn remove_identity_equiv(
                             let error_rate = target.get_error(inst.op.name(), qargs.as_slice());
                             match error_rate {
                                 Some(err) => err * degree,
-                                None => f64::EPSILON.max(1. - degree),
+                                None => minimum_tol.max(1. - degree),
                             }
                         }
-                        None => f64::EPSILON.max(1. - degree),
+                        None => minimum_tol.max(1. - degree),
                     }
                 }
             }
@@ -67,10 +72,10 @@ fn remove_identity_equiv(
                     let error_rate = target.get_error(inst.op.name(), qargs.as_slice());
                     match error_rate {
                         Some(err) => err,
-                        None => f64::EPSILON,
+                        None => minimum_tol,
                     }
                 }
-                None => f64::EPSILON,
+                None => minimum_tol,
             },
         }
     };
@@ -83,22 +88,24 @@ fn remove_identity_equiv(
         let view = inst.op.view();
         match view {
             OperationRef::StandardGate(gate) => {
-                let (dim, trace) = match gate {
-                    StandardGate::RXGate | StandardGate::RYGate | StandardGate::RZGate => {
-                        if let Param::Float(theta) = inst.params_view()[0] {
-                            let trace = Complex64::new((theta / 2.).cos() * 2., 0.);
-                            (2., trace)
-                        } else {
-                            continue;
-                        }
-                    }
-                    StandardGate::RXXGate
+                let (tr_over_dim, dim) = match gate {
+                    StandardGate::RXGate
+                    | StandardGate::RYGate
+                    | StandardGate::RZGate
+                    | StandardGate::PhaseGate
+                    | StandardGate::RXXGate
                     | StandardGate::RYYGate
+                    | StandardGate::RZXGate
                     | StandardGate::RZZGate
-                    | StandardGate::RZXGate => {
-                        if let Param::Float(theta) = inst.params_view()[0] {
-                            let trace = Complex64::new((theta / 2.).cos() * 4., 0.);
-                            (4., trace)
+                    | StandardGate::CRXGate
+                    | StandardGate::CRYGate
+                    | StandardGate::CRZGate
+                    | StandardGate::CPhaseGate => {
+                        let name = gate.name();
+                        if let Param::Float(angle) = inst.params_view()[0] {
+                            let (tr_over_dim, dim) =
+                                rotation_trace_and_dim(name, angle).expect("All rotation covered");
+                            (tr_over_dim, dim)
                         } else {
                             continue;
                         }
@@ -106,19 +113,19 @@ fn remove_identity_equiv(
                     _ => {
                         if let Some(matrix) = gate.matrix(inst.params_view()) {
                             let dim = matrix.shape()[0] as f64;
-                            let trace = matrix.diag().iter().sum::<Complex64>();
-                            (dim, trace)
+                            let tr_over_dim = matrix.diag().iter().sum::<Complex64>() / dim;
+                            (tr_over_dim, dim)
                         } else {
                             continue;
                         }
                     }
                 };
                 let error = get_error_cutoff(inst);
-                let f_pro = (trace / dim).abs().powi(2);
+                let f_pro = tr_over_dim.abs().powi(2);
                 let gate_fidelity = (dim * f_pro + 1.) / (dim + 1.);
                 if (1. - gate_fidelity).abs() < error {
                     remove_list.push(op_node);
-                    global_phase_update += (trace / dim).arg();
+                    global_phase_update += tr_over_dim.arg();
                 }
             }
             _ => {
@@ -127,12 +134,12 @@ fn remove_identity_equiv(
                 if let Some(matrix) = matrix {
                     let error = get_error_cutoff(inst);
                     let dim = matrix.shape()[0] as f64;
-                    let trace: Complex64 = matrix.diag().iter().sum();
-                    let f_pro = (trace / dim).abs().powi(2);
+                    let tr_over_dim = matrix.diag().iter().sum::<Complex64>() / dim;
+                    let f_pro = tr_over_dim.abs().powi(2);
                     let gate_fidelity = (dim * f_pro + 1.) / (dim + 1.);
                     if (1. - gate_fidelity).abs() < error {
                         remove_list.push(op_node);
-                        global_phase_update += (trace / dim).arg();
+                        global_phase_update += tr_over_dim.arg();
                     }
                 }
             }
diff --git a/qiskit/circuit/__init__.py b/qiskit/circuit/__init__.py
index 42cb8ddb83d1..fe30b9735815 100644
--- a/qiskit/circuit/__init__.py
+++ b/qiskit/circuit/__init__.py
@@ -806,14 +806,8 @@
     q_1: ───────────┤ X ├───┤ X ├───┤ X ├     q_1: ───┤ X ├───
                     └───┘   └───┘   └───┘             └───┘
 
-Performing these optimizations are part of the transpiler, but the tools to investigate commutations
-are available in the :class:`CommutationChecker`.
-
-.. autosummary::
-   :toctree: ../stubs/
-
-   CommutationChecker
-
+Performing these optimizations are part of the transpiler, using passes such as 
+:class:`.CommutationAnalysis`.
 
 .. _circuit-custom-gates:
 
diff --git a/qiskit/transpiler/passes/optimization/commutation_analysis.py b/qiskit/transpiler/passes/optimization/commutation_analysis.py
index d801e4775937..f2251a63f5ad 100644
--- a/qiskit/transpiler/passes/optimization/commutation_analysis.py
+++ b/qiskit/transpiler/passes/optimization/commutation_analysis.py
@@ -18,11 +18,31 @@
 
 
 class CommutationAnalysis(AnalysisPass):
-    """Analysis pass to find commutation relations between DAG nodes.
+    r"""Analysis pass to find commutation relations between DAG nodes.
 
-    ``property_set['commutation_set']`` is a dictionary that describes
-    the commutation relations on a given wire, all the gates on a wire
+    This sets ``property_set['commutation_set']`` to a dictionary that describes
+    the commutation relations on a given wire: all the gates on a wire
     are grouped into a set of gates that commute.
+
+    When possible, commutation relations are queried from a lookup table. This is the case
+    for standard gates without parameters (such as :class:`.XGate` or :class:`.HGate`) or
+    gates with free parameters (such as :class:`.RXGate` with a :class:`.ParameterExpression` as
+    angle). Otherwise, a matrix-based check is performed, where two operations are said to
+    commute, if the average gate fidelity of performing the commutation is above a certain threshold.
+    Concretely, two unitaries :math:`A` and :math:`B` on :math:`n` qubits commute if
+
+    .. math::
+
+        \frac{2^n F_{\text{process}}(AB, BA) + 1}{2^n + 1} > 1 - \varepsilon,
+
+    where
+
+    .. math::
+
+        F_{\text{process}}(U_1, U_2) = \left|\frac{\mathrm{Tr}(U_1 U_2^\dagger)}{2^n} \right|^2,
+
+    and we set :math:`\varepsilon` to ``16 * machine_eps`` to account for round-off errors on
+    few-qubit systems.
     """
 
     def __init__(self, *, _commutation_checker=None):
diff --git a/qiskit/transpiler/passes/optimization/remove_identity_equiv.py b/qiskit/transpiler/passes/optimization/remove_identity_equiv.py
index 17445eb5ad2a..53680475e5af 100644
--- a/qiskit/transpiler/passes/optimization/remove_identity_equiv.py
+++ b/qiskit/transpiler/passes/optimization/remove_identity_equiv.py
@@ -32,7 +32,7 @@ class RemoveIdentityEquivalent(TransformationPass):
 
     .. math::
 
-        \bar{F} = \frac{1 + F_{\text{process}}}{1 + d},\ 
+        \bar{F} = \frac{1 + d F_{\text{process}}}{1 + d},\ 
 
         F_{\text{process}} = \frac{|\mathrm{Tr}(G)|^2}{d^2}
 
@@ -47,9 +47,11 @@ def __init__(
         Args:
             approximation_degree: The degree to approximate for the equivalence check. This can be a
                 floating point value between 0 and 1, or ``None``. If the value is 1 this does not
-                approximate above floating point precision. For a value < 1 this is used as a scaling
-                factor for the cutoff fidelity. If the value is ``None`` this approximates up to the
-                fidelity for the gate specified in ``target``.
+                approximate above the floating point precision. For a value < 1 this is used as a
+                scaling factor for the cutoff fidelity. If the value is ``None`` this approximates up
+                to the fidelity for the gate specified in ``target``. In case no ``target`` is set
+                we approximate up to ``16 * machine_eps`` as default to account for accumulations
+                on few-qubit systems.
 
             target: If ``approximation_degree`` is set to ``None`` and a :class:`.Target` is provided
                 for this field the tolerance for determining whether an operation is equivalent to
diff --git a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
new file mode 100644
index 000000000000..d697a2044b5a
--- /dev/null
+++ b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
@@ -0,0 +1,24 @@
+---
+upgrade:
+  - |
+    Increased the minimum threshold for when gates are assumed to be the identity in
+    :class:`.RemoveIdentityEquivalent` from machine epsilon to 16 times machine epsilon to
+    account for round-off errors in the fidelity calculation and for consistency with the
+    :class:`.CommutationAnalysis`. 
+  - |
+    Updated the metric used in :class:`.CommutationAnalysis` for matrix-based commutation checks.
+    In a matrix-based check (used e.g. for non-standard gates or rotation gates with known rotation
+    angles) two gates are assumed to commute if the average gate fidelity of the commutation is
+    above ``1 - 16 * machine_eps``. This value is chosen to account for round-off errors in the
+    fidelity calculation and for consistency with :class:`.RemoveIdentityEquivalent`.
+    See the class docstring for more information.
+fixes:
+  - |
+    Fixed an inconsistency in dealing with close-to-identity gates in the transpilation process,
+    where gates that are close to the identity were assumed to commute with everything, but
+    not removed. The underlying issue were different metrics used in :class:`.RemoveIdentityEquivalent`
+    and :class:`.CommutationAnalysis` (and, by extension, :class:`.CommutativeInverseCancellation`).
+    Both now use the average gate fidelity and the same threshold to assess whether a gate
+    should be treated as identity (such as a rotation gate with very small angle). See either
+    of the aforementioned classes' docstrings for more information.
+    Fixed `#13547 <https://github.com/Qiskit/qiskit/issues/13547>`__.
diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py
index 963413dc1bc6..724bd3a64db4 100644
--- a/test/python/circuit/test_commutation_checker.py
+++ b/test/python/circuit/test_commutation_checker.py
@@ -55,6 +55,7 @@
     RZZGate,
     SGate,
     XGate,
+    YGate,
     ZGate,
     HGate,
     UnitaryGate,
@@ -113,53 +114,51 @@ def test_simple_gates(self):
         different orders of gates, different orders of qubits, different sets of
         qubits over which gates are defined, and so on."""
 
-        # should commute
+        self.assertTrue(scc.commute(HGate(), [0], [], HGate(), [0], []))
+
         self.assertTrue(scc.commute(ZGate(), [0], [], CXGate(), [0, 1], []))
-        # should not commute
         self.assertFalse(scc.commute(ZGate(), [1], [], CXGate(), [0, 1], []))
-        # should not commute
+
         self.assertFalse(scc.commute(XGate(), [0], [], CXGate(), [0, 1], []))
-        # should commute
         self.assertTrue(scc.commute(XGate(), [1], [], CXGate(), [0, 1], []))
-        # should not commute
         self.assertFalse(scc.commute(XGate(), [1], [], CXGate(), [1, 0], []))
-        # should commute
         self.assertTrue(scc.commute(XGate(), [0], [], CXGate(), [1, 0], []))
-        # should commute
         self.assertTrue(scc.commute(CXGate(), [1, 0], [], XGate(), [0], []))
-        # should not commute
         self.assertFalse(scc.commute(CXGate(), [1, 0], [], XGate(), [1], []))
-        # should commute
+
         self.assertTrue(scc.commute(CXGate(), [1, 0], [], CXGate(), [1, 0], []))
-        # should not commute
         self.assertFalse(scc.commute(CXGate(), [1, 0], [], CXGate(), [0, 1], []))
-        # should commute
         self.assertTrue(scc.commute(CXGate(), [1, 0], [], CXGate(), [1, 2], []))
-        # should not commute
         self.assertFalse(scc.commute(CXGate(), [1, 0], [], CXGate(), [2, 1], []))
-        # should commute
         self.assertTrue(scc.commute(CXGate(), [1, 0], [], CXGate(), [2, 3], []))
+
         self.assertTrue(scc.commute(XGate(), [2], [], CCXGate(), [0, 1, 2], []))
         self.assertFalse(scc.commute(CCXGate(), [0, 1, 2], [], CCXGate(), [0, 2, 1], []))
 
-        self.assertTrue(scc.commute(HGate(), [0], [], HGate(), [0], []))
+        # these would commute up to a global phase
+        self.assertFalse(scc.commute(HGate(), [0], [], YGate(), [0], []))
 
     def test_simple_matrices(self):
         """Test simple gates but matrix-based."""
+        x = UnitaryGate(XGate())
         had = UnitaryGate(HGate())
+        had2 = UnitaryGate(np.kron(HGate(), HGate()))
         cx = UnitaryGate(CXGate())
-        x = UnitaryGate(XGate())
+
         self.assertTrue(scc.commute(x, [0], [], x, [0], []))
         self.assertTrue(scc.commute(had, [0], [], had, [0], []))
+
+        self.assertTrue(scc.commute(had2, [0, 1], [], had2, [1, 0], []))
+        self.assertFalse(scc.commute(had2, [0, 1], [], cx, [1, 0], []))
         self.assertTrue(scc.commute(cx, [0, 1], [], cx, [0, 1], []))
+
+        self.assertFalse(scc.commute(x, [0], [], cx, [0, 1], []))
         self.assertTrue(scc.commute(x, [1], [], cx, [0, 1], []))
 
     def test_passing_quantum_registers(self):
         """Check that passing QuantumRegisters works correctly."""
         qr = QuantumRegister(4)
-        # should commute
         self.assertTrue(scc.commute(CXGate(), [qr[1], qr[0]], [], CXGate(), [qr[1], qr[2]], []))
-        # should not commute
         self.assertFalse(scc.commute(CXGate(), [qr[0], qr[1]], [], CXGate(), [qr[1], qr[2]], []))
 
     def test_standard_gates_commutations(self):
@@ -438,9 +437,11 @@ def test_cutoff_angles(self, gate_cls):
         # the cutoff angle depends on the average gate fidelity; i.e. it is the angle
         # for which the average gate fidelity is smaller than machine epsilon
         if gate_cls in [CPhaseGate, CRXGate, CRYGate, CRZGate]:
-            cutoff_angle = 4.71e-8
+            # cutoff_angle = 4.71e-8
+            cutoff_angle = 1.91e-7
         else:
-            cutoff_angle = 3.65e-8
+            # cutoff_angle = 3.65e-8
+            cutoff_angle = 1.34e-7
 
         for i in range(1, max_power + 1):
             angle = 2 ** (-i)
@@ -521,6 +522,20 @@ def test_2q_pauli_rot_with_non_cached(self):
         self.assertFalse(scc.commute(something_else, [0], [], RYYGate(np.pi), [1, 0], []))
         self.assertFalse(scc.commute(something_else, [1], [], RYYGate(np.pi), [1, 0], []))
 
+    def test_approximation_degree(self):
+        """Test setting the approximation degree."""
+
+        almost_identity = RZGate(1e-5)
+        other = HGate()
+
+        self.assertFalse(scc.commute(almost_identity, [0], [], other, [0], []))
+        self.assertFalse(
+            scc.commute(almost_identity, [0], [], other, [0], [], approximation_degree=1)
+        )
+        self.assertTrue(
+            scc.commute(almost_identity, [0], [], other, [0], [], approximation_degree=1 - 1e-4)
+        )
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/test/python/transpiler/test_remove_identity_equivalent.py b/test/python/transpiler/test_remove_identity_equivalent.py
index 495978ef585a..b902c7d8bff5 100644
--- a/test/python/transpiler/test_remove_identity_equivalent.py
+++ b/test/python/transpiler/test_remove_identity_equivalent.py
@@ -10,7 +10,7 @@
 # copyright notice, and modified files need to carry a notice indicating
 # that they have been altered from the originals.
 
-"""Tests for the DropNegligible transpiler pass."""
+"""Tests for the RemoveIdentityEquivalent transpiler pass."""
 
 import ddt
 import numpy as np
@@ -37,10 +37,10 @@
 
 
 @ddt.ddt
-class TestDropNegligible(QiskitTestCase):
-    """Test the DropNegligible pass."""
+class TestRemoveIdentityEquivalent(QiskitTestCase):
+    """Test the RemoveIdentityEquivalent pass."""
 
-    def test_drops_negligible_gates(self):
+    def test_remove_identity_equiv_pass(self):
         """Test that negligible gates are dropped."""
         qubits = QuantumRegister(2)
         circuit = QuantumCircuit(qubits)

From cda06f8b8aca076d3ae4e4ec361c013783764ae9 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Mon, 3 Mar 2025 09:57:48 +0100
Subject: [PATCH 12/15] Python-space approx degree oversight

---
 qiskit/circuit/commutation_checker.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/qiskit/circuit/commutation_checker.py b/qiskit/circuit/commutation_checker.py
index 994867bcf1b4..3bfaac2571f3 100644
--- a/qiskit/circuit/commutation_checker.py
+++ b/qiskit/circuit/commutation_checker.py
@@ -41,7 +41,7 @@ def commute_nodes(
         op1,
         op2,
         max_num_qubits: int = 3,
-        approximation_degree: float | None = None,
+        approximation_degree: float = 1.0,
     ) -> bool:
         """Checks if two DAGOpNodes commute."""
         return self.cc.commute_nodes(op1, op2, max_num_qubits, approximation_degree)
@@ -55,7 +55,7 @@ def commute(
         qargs2: List,
         cargs2: List,
         max_num_qubits: int = 3,
-        approximation_degree: float | None = None,
+        approximation_degree: float = 1.0,
     ) -> bool:
         """
         Checks if two Operations commute. The return value of `True` means that the operations
@@ -73,7 +73,7 @@ def commute(
             max_num_qubits: the maximum number of qubits to consider, the check may be skipped if
                 the number of qubits for either operation exceeds this amount.
             approximation_degree: If the average gate fidelity in between the two operations
-                is above this number (up to machine epsilon) they are assumed to commute.
+                is above this number (up to 16 times machine epsilon) they are assumed to commute.
 
         Returns:
             bool: whether two operations commute.

From 2c48ac4050037333d99b52133d14730d850d1a56 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Wed, 5 Mar 2025 16:59:43 +0100
Subject: [PATCH 13/15] update threshold to 1e-12

which is consistent with the 2qWeylDecomp and still stricter than Qiskit 1.x
---
 crates/accelerate/src/commutation_checker.rs           |  6 +++---
 crates/accelerate/src/remove_identity_equiv.rs         |  7 +++----
 .../notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml       | 10 +++++-----
 test/python/circuit/test_commutation_checker.py        |  6 ++----
 4 files changed, 13 insertions(+), 16 deletions(-)

diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index 86f7d1d6be6d..c5d0eeb1349e 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -284,9 +284,9 @@ impl CommutationChecker {
         approximation_degree: f64,
     ) -> PyResult<bool> {
         // If the average gate infidelity is below this tolerance, they commute. The tolerance
-        // is set to max(16 * machine_eps, 1 - approximation_degree), where the heuristic factor
-        // of 16 is set to account for round-off errors in matrix multiplication.
-        let tol = (16. * f64::EPSILON).max(1. - approximation_degree);
+        // is set to max(1e-12, 1 - approximation_degree), to account for roundoffs and for
+        // consistency with other places in Qiskit.
+        let tol = 1e-12_f64.max(1. - approximation_degree);
 
         // if we have rotation gates, we attempt to map them to their generators, for example
         // RX -> X or CPhase -> CZ
diff --git a/crates/accelerate/src/remove_identity_equiv.rs b/crates/accelerate/src/remove_identity_equiv.rs
index 5b8ab2610b59..0146faf5b3a6 100644
--- a/crates/accelerate/src/remove_identity_equiv.rs
+++ b/crates/accelerate/src/remove_identity_equiv.rs
@@ -34,10 +34,9 @@ fn remove_identity_equiv(
 ) {
     let mut remove_list: Vec<NodeIndex> = Vec::new();
     let mut global_phase_update: f64 = 0.;
-    // Minimum threshold to compare average gate fidelity to 1. Includes a heuristic factor
-    // of 16 to account for round-off errors in the trace calculations, and for consistency
-    // with the commutation checker.
-    let minimum_tol = 16. * f64::EPSILON;
+    // Minimum threshold to compare average gate fidelity to 1. This is chosen to account
+    // for roundoff errors and to be consistent with other places.
+    let minimum_tol = 1e-12_f64;
 
     let get_error_cutoff = |inst: &PackedInstruction| -> f64 {
         match approx_degree {
diff --git a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
index d697a2044b5a..115362423ee1 100644
--- a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
+++ b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
@@ -2,16 +2,16 @@
 upgrade:
   - |
     Increased the minimum threshold for when gates are assumed to be the identity in
-    :class:`.RemoveIdentityEquivalent` from machine epsilon to 16 times machine epsilon to
+    :class:`.RemoveIdentityEquivalent` from machine epsilon to ``1e-12`` to
     account for round-off errors in the fidelity calculation and for consistency with the
-    :class:`.CommutationAnalysis`. 
+    other classes, such as :class:`.CommutationAnalysis` and :class:`.TwoQubitWeylDecomposition`.
   - |
     Updated the metric used in :class:`.CommutationAnalysis` for matrix-based commutation checks.
     In a matrix-based check (used e.g. for non-standard gates or rotation gates with known rotation
     angles) two gates are assumed to commute if the average gate fidelity of the commutation is
-    above ``1 - 16 * machine_eps``. This value is chosen to account for round-off errors in the
-    fidelity calculation and for consistency with :class:`.RemoveIdentityEquivalent`.
-    See the class docstring for more information.
+    above ``1 - 1e-12``. This value is chosen to account for round-off errors in the
+    fidelity calculation and for consistency with :class:`.RemoveIdentityEquivalent` and
+    :class:`.TwoQubitWeylDecomposition`. See the class docstring for more information.
 fixes:
   - |
     Fixed an inconsistency in dealing with close-to-identity gates in the transpilation process,
diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py
index 724bd3a64db4..f46ebdbaca24 100644
--- a/test/python/circuit/test_commutation_checker.py
+++ b/test/python/circuit/test_commutation_checker.py
@@ -437,11 +437,9 @@ def test_cutoff_angles(self, gate_cls):
         # the cutoff angle depends on the average gate fidelity; i.e. it is the angle
         # for which the average gate fidelity is smaller than machine epsilon
         if gate_cls in [CPhaseGate, CRXGate, CRYGate, CRZGate]:
-            # cutoff_angle = 4.71e-8
-            cutoff_angle = 1.91e-7
+            cutoff_angle = 3.16e-6
         else:
-            # cutoff_angle = 3.65e-8
-            cutoff_angle = 1.34e-7
+            cutoff_angle = 2.2e-6
 
         for i in range(1, max_power + 1):
             angle = 2 ** (-i)

From 6bc565686ecde853d52b76573fce709558ccc76b Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Wed, 5 Mar 2025 23:11:53 +0100
Subject: [PATCH 14/15] review comments

---
 crates/accelerate/src/commutation_checker.rs  | 28 ++-----
 crates/accelerate/src/gate_metrics.rs         | 76 +++++++++++++++++++
 crates/accelerate/src/lib.rs                  |  1 +
 .../accelerate/src/remove_identity_equiv.rs   | 18 ++---
 crates/accelerate/src/unitary_compose.rs      | 43 ++---------
 qiskit/circuit/__init__.py                    | 10 ++-
 qiskit/circuit/commutation_checker.py         | 32 ++++++--
 .../optimization/commutation_analysis.py      | 20 -----
 .../cc-gate-fidelity-bde7a01dc8c56e29.yaml    | 17 ++++-
 .../circuit/test_commutation_checker.py       |  2 +-
 10 files changed, 147 insertions(+), 100 deletions(-)
 create mode 100644 crates/accelerate/src/gate_metrics.rs

diff --git a/crates/accelerate/src/commutation_checker.rs b/crates/accelerate/src/commutation_checker.rs
index 7c3109213ba3..c99b35a47391 100644
--- a/crates/accelerate/src/commutation_checker.rs
+++ b/crates/accelerate/src/commutation_checker.rs
@@ -35,6 +35,7 @@ use qiskit_circuit::operations::{
 };
 use qiskit_circuit::{BitType, Clbit, Qubit};
 
+use crate::gate_metrics;
 use crate::unitary_compose;
 use crate::QiskitError;
 
@@ -500,7 +501,7 @@ impl CommutationChecker {
             Ok(matrix) => matrix,
             Err(e) => return Err(PyRuntimeError::new_err(e)),
         };
-        let (fid, phase) = unitary_compose::gate_fidelity(&op12.view(), &op21.view(), None);
+        let (fid, phase) = gate_metrics::gate_fidelity(&op12.view(), &op21.view(), None);
 
         // we consider the gates as commuting if the process fidelity of
         // AB (BA)^\dagger is approximately the identity and there is no global phase difference
@@ -609,7 +610,10 @@ fn map_rotation<'a>(
         // commute with everything, and we simply return the operation with the flag that
         // it commutes trivially.
         if let Param::Float(angle) = params[0] {
-            let (tr_over_dim, dim) = rotation_trace_and_dim(name, angle)
+            let gate = op
+                .standard_gate()
+                .expect("Supported gates are standard gates");
+            let (tr_over_dim, dim) = gate_metrics::rotation_trace_and_dim(gate, angle)
                 .expect("All rotation should be covered at this point");
             let gate_fidelity = tr_over_dim.abs().powi(2);
             let process_fidelity = (dim * gate_fidelity + 1.) / (dim + 1.);
@@ -630,26 +634,6 @@ fn map_rotation<'a>(
     (op, params, false)
 }
 
-/// For a (controlled) rotation or phase gate, return a tuple ``(Tr(gate) / dim, dim)``.
-/// Returns ``None`` if the rotation gate (specified by name) is not supported.
-pub fn rotation_trace_and_dim(rotation: &str, angle: f64) -> Option<(Complex64, f64)> {
-    let dim = match rotation {
-        "rx" | "ry" | "rz" | "p" | "u1" => 2.,
-        _ => 4.,
-    };
-
-    let trace_over_dim = match rotation {
-        "rx" | "ry" | "rz" | "rxx" | "ryy" | "rzx" | "rzz" => {
-            Complex64::new((angle / 2.).cos(), 0.)
-        }
-        "crx" | "cry" | "crz" => Complex64::new(0.5 + 0.5 * (angle / 2.).cos(), 0.),
-        "p" | "u1" => (1. + Complex64::new(0., angle).exp()) / 2.,
-        "cp" => (3. + Complex64::new(0., angle).exp()) / 4.,
-        _ => return None,
-    };
-    Some((trace_over_dim, dim))
-}
-
 fn get_relative_placement(
     first_qargs: &[Qubit],
     second_qargs: &[Qubit],
diff --git a/crates/accelerate/src/gate_metrics.rs b/crates/accelerate/src/gate_metrics.rs
new file mode 100644
index 000000000000..d361af60b3e4
--- /dev/null
+++ b/crates/accelerate/src/gate_metrics.rs
@@ -0,0 +1,76 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2025
+//
+// 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.
+
+use ndarray::ArrayView2;
+use num_complex::{Complex64, ComplexFloat};
+use qiskit_circuit::{operations::StandardGate, Qubit};
+
+use crate::unitary_compose;
+
+/// For a (controlled) rotation or phase gate, return a tuple ``(Tr(gate) / dim, dim)``.
+/// Returns ``None`` if the rotation gate (specified by name) is not supported.
+pub fn rotation_trace_and_dim(rotation: StandardGate, angle: f64) -> Option<(Complex64, f64)> {
+    let dim = match rotation {
+        StandardGate::RXGate
+        | StandardGate::RYGate
+        | StandardGate::RZGate
+        | StandardGate::PhaseGate
+        | StandardGate::U1Gate => 2.,
+        _ => 4.,
+    };
+
+    let trace_over_dim = match rotation {
+        StandardGate::RXGate
+        | StandardGate::RYGate
+        | StandardGate::RZGate
+        | StandardGate::RXXGate
+        | StandardGate::RYYGate
+        | StandardGate::RZZGate
+        | StandardGate::RZXGate => Complex64::new((angle / 2.).cos(), 0.),
+        StandardGate::CRXGate | StandardGate::CRYGate | StandardGate::CRZGate => {
+            Complex64::new(0.5 + 0.5 * (angle / 2.).cos(), 0.)
+        }
+        StandardGate::PhaseGate | StandardGate::U1Gate => {
+            (1. + Complex64::new(0., angle).exp()) / 2.
+        }
+        StandardGate::CPhaseGate => (3. + Complex64::new(0., angle).exp()) / 4.,
+        _ => return None,
+    };
+    Some((trace_over_dim, dim))
+}
+
+pub fn gate_fidelity(
+    left: &ArrayView2<Complex64>,
+    right: &ArrayView2<Complex64>,
+    qargs: Option<&[Qubit]>,
+) -> (f64, f64) {
+    let dim = left.nrows(); // == left.ncols() == right.nrows() == right.ncols()
+                            // let trace = left.t().mapv(|el| el.conj()).dot(right).diag().sum();
+
+    let left = left.t().mapv(|el| el.conj());
+    let product = match dim {
+        2 => unitary_compose::matmul_1q(&left.view(), right),
+        4 => {
+            unitary_compose::matmul_2q(&left.view(), right, qargs.unwrap_or(&[Qubit(0), Qubit(1)]))
+        }
+        _ => left.dot(right),
+    };
+    let trace = product.diag().sum();
+
+    let dim = dim as f64;
+    let normalized_trace = trace / dim;
+    let phase = normalized_trace.arg(); // compute phase difference
+
+    let process_fidelity = normalized_trace.abs().powi(2);
+    let gate_fidelity = (dim * process_fidelity + 1.) / (dim + 1.);
+    (gate_fidelity, phase)
+}
diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs
index bcb21ee3774e..ef90ac27041e 100644
--- a/crates/accelerate/src/lib.rs
+++ b/crates/accelerate/src/lib.rs
@@ -31,6 +31,7 @@ pub mod error_map;
 pub mod euler_one_qubit_decomposer;
 pub mod filter_op_nodes;
 pub mod gate_direction;
+pub mod gate_metrics;
 pub mod gates_in_basis;
 pub mod high_level_synthesis;
 pub mod inverse_cancellation;
diff --git a/crates/accelerate/src/remove_identity_equiv.rs b/crates/accelerate/src/remove_identity_equiv.rs
index 0146faf5b3a6..33e9d630a0b6 100644
--- a/crates/accelerate/src/remove_identity_equiv.rs
+++ b/crates/accelerate/src/remove_identity_equiv.rs
@@ -14,7 +14,7 @@ use num_complex::ComplexFloat;
 use pyo3::prelude::*;
 use rustworkx_core::petgraph::stable_graph::NodeIndex;
 
-use crate::commutation_checker::rotation_trace_and_dim;
+use crate::gate_metrics::rotation_trace_and_dim;
 use crate::nlayout::PhysicalQubit;
 use crate::target_transpiler::Target;
 use qiskit_circuit::dag_circuit::DAGCircuit;
@@ -24,6 +24,8 @@ use qiskit_circuit::operations::Param;
 use qiskit_circuit::operations::StandardGate;
 use qiskit_circuit::packed_instruction::PackedInstruction;
 
+const MINIMUM_TOL: f64 = 1e-12;
+
 #[pyfunction]
 #[pyo3(signature=(dag, approx_degree=Some(1.0), target=None))]
 fn remove_identity_equiv(
@@ -36,13 +38,12 @@ fn remove_identity_equiv(
     let mut global_phase_update: f64 = 0.;
     // Minimum threshold to compare average gate fidelity to 1. This is chosen to account
     // for roundoff errors and to be consistent with other places.
-    let minimum_tol = 1e-12_f64;
 
     let get_error_cutoff = |inst: &PackedInstruction| -> f64 {
         match approx_degree {
             Some(degree) => {
                 if degree == 1.0 {
-                    minimum_tol
+                    MINIMUM_TOL
                 } else {
                     match target {
                         Some(target) => {
@@ -54,10 +55,10 @@ fn remove_identity_equiv(
                             let error_rate = target.get_error(inst.op.name(), qargs.as_slice());
                             match error_rate {
                                 Some(err) => err * degree,
-                                None => minimum_tol.max(1. - degree),
+                                None => MINIMUM_TOL.max(1. - degree),
                             }
                         }
-                        None => minimum_tol.max(1. - degree),
+                        None => MINIMUM_TOL.max(1. - degree),
                     }
                 }
             }
@@ -71,10 +72,10 @@ fn remove_identity_equiv(
                     let error_rate = target.get_error(inst.op.name(), qargs.as_slice());
                     match error_rate {
                         Some(err) => err,
-                        None => minimum_tol,
+                        None => MINIMUM_TOL,
                     }
                 }
-                None => minimum_tol,
+                None => MINIMUM_TOL,
             },
         }
     };
@@ -100,10 +101,9 @@ fn remove_identity_equiv(
                     | StandardGate::CRYGate
                     | StandardGate::CRZGate
                     | StandardGate::CPhaseGate => {
-                        let name = gate.name();
                         if let Param::Float(angle) = inst.params_view()[0] {
                             let (tr_over_dim, dim) =
-                                rotation_trace_and_dim(name, angle).expect("All rotation covered");
+                                rotation_trace_and_dim(gate, angle).expect("Since only supported rotation gates are given, the result is not None");
                             (tr_over_dim, dim)
                         } else {
                             continue;
diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs
index 242f611ab65b..d9791e392523 100644
--- a/crates/accelerate/src/unitary_compose.rs
+++ b/crates/accelerate/src/unitary_compose.rs
@@ -10,10 +10,9 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use ndarray::{arr2, Array, Array2, ArrayView, ArrayView2, IxDyn};
+use ndarray::{Array, Array2, ArrayView, ArrayView2, IxDyn};
 use ndarray_einsum::*;
-use num_complex::{Complex, Complex64, ComplexFloat};
-use qiskit_circuit::util::C_ZERO;
+use num_complex::{Complex, Complex64};
 use qiskit_circuit::Qubit;
 
 static LOWERCASE: [u8; 26] = [
@@ -153,33 +152,8 @@ fn _einsum_matmul_index(qubits: &[u32], num_qubits: usize) -> String {
     )
 }
 
-pub fn gate_fidelity(
-    left: &ArrayView2<Complex64>,
-    right: &ArrayView2<Complex64>,
-    qargs: Option<&[Qubit]>,
-) -> (f64, f64) {
-    let dim = left.nrows(); // == left.ncols() == right.nrows() == right.ncols()
-                            // let trace = left.t().mapv(|el| el.conj()).dot(right).diag().sum();
-
-    let left = left.t().mapv(|el| el.conj());
-    let product = match dim {
-        2 => mm1q(&left.view(), right),
-        4 => mm2q(&left.view(), right, qargs.unwrap_or(&[Qubit(0), Qubit(1)])),
-        _ => left.dot(right),
-    };
-    let trace = product.diag().sum();
-
-    let dim = dim as f64;
-    let normalized_trace = trace / dim;
-    let phase = normalized_trace.arg(); // compute phase difference
-
-    let process_fidelity = normalized_trace.abs().powi(2);
-    let gate_fidelity = (dim * process_fidelity + 1.) / (dim + 1.);
-    (gate_fidelity, phase)
-}
-
-fn mm1q(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> Array2<Complex64> {
-    let mut out = arr2(&[[C_ZERO, C_ZERO], [C_ZERO, C_ZERO]]);
+pub fn matmul_1q(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> Array2<Complex64> {
+    let mut out = Array2::zeros((2, 2));
     out[[0, 0]] = left[[0, 0]] * right[[0, 0]] + left[[0, 1]] * right[[1, 0]];
     out[[0, 1]] = left[[0, 0]] * right[[0, 1]] + left[[0, 1]] * right[[1, 1]];
     out[[1, 0]] = left[[1, 0]] * right[[0, 0]] + left[[1, 1]] * right[[1, 0]];
@@ -197,17 +171,12 @@ fn _ind(i: usize, reversed: bool) -> usize {
     }
 }
 
-pub fn mm2q(
+pub fn matmul_2q(
     left: &ArrayView2<Complex64>,
     right: &ArrayView2<Complex64>,
     qargs: &[Qubit],
 ) -> Array2<Complex64> {
-    let mut out = arr2(&[
-        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
-        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
-        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
-        [C_ZERO, C_ZERO, C_ZERO, C_ZERO],
-    ]);
+    let mut out = Array2::zeros((4, 4));
 
     let rev = qargs[0].0 == 1;
     for i in 0..4usize {
diff --git a/qiskit/circuit/__init__.py b/qiskit/circuit/__init__.py
index 1f37c4e82676..7040278a06ba 100644
--- a/qiskit/circuit/__init__.py
+++ b/qiskit/circuit/__init__.py
@@ -806,8 +806,14 @@
     q_1: ───────────┤ X ├───┤ X ├───┤ X ├     q_1: ───┤ X ├───
                     └───┘   └───┘   └───┘             └───┘
 
-Performing these optimizations are part of the transpiler, using passes such as 
-:class:`.CommutationAnalysis`.
+Performing these optimizations are part of the transpiler, but the tools to investigate commutations
+are available in the :class:`CommutationChecker`.
+
+.. autosummary::
+   :toctree: ../stubs/
+
+   CommutationChecker
+
 
 .. _circuit-custom-gates:
 
diff --git a/qiskit/circuit/commutation_checker.py b/qiskit/circuit/commutation_checker.py
index 3bfaac2571f3..3e4f42fe39ff 100644
--- a/qiskit/circuit/commutation_checker.py
+++ b/qiskit/circuit/commutation_checker.py
@@ -20,11 +20,31 @@
 
 
 class CommutationChecker:
-    """This code is essentially copy-pasted from commutative_analysis.py.
-    This code cleverly hashes commutativity and non-commutativity results between DAG nodes and seems
-    quite efficient for large Clifford circuits.
-    They may be other possible efficiency improvements: using rule-based commutativity analysis,
-    evicting from the cache less useful entries, etc.
+    r"""Check commutations of two operations.
+
+    Two unitaries :math:`A` and :math:`B` on :math:`n` qubits commute if
+
+    .. math::
+
+        \frac{2^n F_{\text{process}}(AB, BA) + 1}{2^n + 1} > 1 - \varepsilon,
+
+    where
+
+    .. math::
+
+        F_{\text{process}}(U_1, U_2) = \left|\frac{\mathrm{Tr}(U_1 U_2^\dagger)}{2^n} \right|^2,
+
+    and we set :math:`\varepsilon` to :math:`10^{-12}` to account for round-off errors on
+    few-qubit systems. This metric is chosen for consistency with other closeness checks in
+    Qiskit.
+
+    When possible, commutation relations are queried from a lookup table. This is the case
+    for standard gates without parameters (such as :class:`.XGate` or :class:`.HGate`) or
+    gates with free parameters (such as :class:`.RXGate` with a :class:`.ParameterExpression` as
+    angle). Otherwise, a matrix-based check is performed, where two operations are said to
+    commute, if the average gate fidelity of performing the commutation is above a certain threshold
+    (see ``approximation_degree``). The result of this commutation is then added to the
+    cached lookup table.
     """
 
     def __init__(
@@ -73,7 +93,7 @@ def commute(
             max_num_qubits: the maximum number of qubits to consider, the check may be skipped if
                 the number of qubits for either operation exceeds this amount.
             approximation_degree: If the average gate fidelity in between the two operations
-                is above this number (up to 16 times machine epsilon) they are assumed to commute.
+                is above this number (up to ``1e-12``) they are assumed to commute.
 
         Returns:
             bool: whether two operations commute.
diff --git a/qiskit/transpiler/passes/optimization/commutation_analysis.py b/qiskit/transpiler/passes/optimization/commutation_analysis.py
index f2251a63f5ad..d977f3b97514 100644
--- a/qiskit/transpiler/passes/optimization/commutation_analysis.py
+++ b/qiskit/transpiler/passes/optimization/commutation_analysis.py
@@ -23,26 +23,6 @@ class CommutationAnalysis(AnalysisPass):
     This sets ``property_set['commutation_set']`` to a dictionary that describes
     the commutation relations on a given wire: all the gates on a wire
     are grouped into a set of gates that commute.
-
-    When possible, commutation relations are queried from a lookup table. This is the case
-    for standard gates without parameters (such as :class:`.XGate` or :class:`.HGate`) or
-    gates with free parameters (such as :class:`.RXGate` with a :class:`.ParameterExpression` as
-    angle). Otherwise, a matrix-based check is performed, where two operations are said to
-    commute, if the average gate fidelity of performing the commutation is above a certain threshold.
-    Concretely, two unitaries :math:`A` and :math:`B` on :math:`n` qubits commute if
-
-    .. math::
-
-        \frac{2^n F_{\text{process}}(AB, BA) + 1}{2^n + 1} > 1 - \varepsilon,
-
-    where
-
-    .. math::
-
-        F_{\text{process}}(U_1, U_2) = \left|\frac{\mathrm{Tr}(U_1 U_2^\dagger)}{2^n} \right|^2,
-
-    and we set :math:`\varepsilon` to ``16 * machine_eps`` to account for round-off errors on
-    few-qubit systems.
     """
 
     def __init__(self, *, _commutation_checker=None):
diff --git a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
index 115362423ee1..b7534e259626 100644
--- a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
+++ b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
@@ -6,9 +6,8 @@ upgrade:
     account for round-off errors in the fidelity calculation and for consistency with the
     other classes, such as :class:`.CommutationAnalysis` and :class:`.TwoQubitWeylDecomposition`.
   - |
-    Updated the metric used in :class:`.CommutationAnalysis` for matrix-based commutation checks.
-    In a matrix-based check (used e.g. for non-standard gates or rotation gates with known rotation
-    angles) two gates are assumed to commute if the average gate fidelity of the commutation is
+    Updated the metric used to check commutations in :class:`.CommutationChecker`.
+    Two gates are assumed to commute if the average gate fidelity of the commutation is
     above ``1 - 1e-12``. This value is chosen to account for round-off errors in the
     fidelity calculation and for consistency with :class:`.RemoveIdentityEquivalent` and
     :class:`.TwoQubitWeylDecomposition`. See the class docstring for more information.
@@ -22,3 +21,15 @@ fixes:
     should be treated as identity (such as a rotation gate with very small angle). See either
     of the aforementioned classes' docstrings for more information.
     Fixed `#13547 <https://github.com/Qiskit/qiskit/issues/13547>`__.
+features_circuit:
+  - |
+    Added a new argument ``approximation_degree`` to :meth:`.CommutationChecker.commute` and
+    :meth:`.CommutationChecker.commute_nodes`, which allows to set the approximation threshold
+    for when gates are said to commute. See the docstring of :class:`.CommutationChecker` for more 
+    detail.
+features_transpiler:
+  - |
+    Added a new argument ``approximation_degree`` to :class:`.CommutationAnalysis`, which allows 
+    setting the approximation threshold for when gates are said to commute. See the class docstring 
+    for more information.
+
diff --git a/test/python/circuit/test_commutation_checker.py b/test/python/circuit/test_commutation_checker.py
index e478da6ff8e2..a2eed5d75af8 100644
--- a/test/python/circuit/test_commutation_checker.py
+++ b/test/python/circuit/test_commutation_checker.py
@@ -399,7 +399,7 @@ def test_cutoff_angles(self, gate_cls):
         generic_gate = DCXGate()  # gate that does not commute with any rotation gate
 
         # the cutoff angle depends on the average gate fidelity; i.e. it is the angle
-        # for which the average gate fidelity is smaller than machine epsilon
+        # for which the average gate fidelity is smaller than 1e-12
         if gate_cls in [CPhaseGate, CRXGate, CRYGate, CRZGate]:
             cutoff_angle = 3.16e-6
         else:

From 34e9fcca708bfe0587f69240a01288c4e1be2203 Mon Sep 17 00:00:00 2001
From: Julien Gacon <jules.gacon@googlemail.com>
Date: Wed, 5 Mar 2025 23:54:20 +0100
Subject: [PATCH 15/15] inline and rm comments

---
 crates/accelerate/src/gate_metrics.rs                     | 3 +--
 crates/accelerate/src/unitary_compose.rs                  | 2 ++
 releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml | 2 +-
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/crates/accelerate/src/gate_metrics.rs b/crates/accelerate/src/gate_metrics.rs
index d361af60b3e4..5e39d5148159 100644
--- a/crates/accelerate/src/gate_metrics.rs
+++ b/crates/accelerate/src/gate_metrics.rs
@@ -53,8 +53,7 @@ pub fn gate_fidelity(
     right: &ArrayView2<Complex64>,
     qargs: Option<&[Qubit]>,
 ) -> (f64, f64) {
-    let dim = left.nrows(); // == left.ncols() == right.nrows() == right.ncols()
-                            // let trace = left.t().mapv(|el| el.conj()).dot(right).diag().sum();
+    let dim = left.nrows();
 
     let left = left.t().mapv(|el| el.conj());
     let product = match dim {
diff --git a/crates/accelerate/src/unitary_compose.rs b/crates/accelerate/src/unitary_compose.rs
index d9791e392523..4a8805a03f95 100644
--- a/crates/accelerate/src/unitary_compose.rs
+++ b/crates/accelerate/src/unitary_compose.rs
@@ -152,6 +152,7 @@ fn _einsum_matmul_index(qubits: &[u32], num_qubits: usize) -> String {
     )
 }
 
+#[inline]
 pub fn matmul_1q(left: &ArrayView2<Complex64>, right: &ArrayView2<Complex64>) -> Array2<Complex64> {
     let mut out = Array2::zeros((2, 2));
     out[[0, 0]] = left[[0, 0]] * right[[0, 0]] + left[[0, 1]] * right[[1, 0]];
@@ -171,6 +172,7 @@ fn _ind(i: usize, reversed: bool) -> usize {
     }
 }
 
+#[inline]
 pub fn matmul_2q(
     left: &ArrayView2<Complex64>,
     right: &ArrayView2<Complex64>,
diff --git a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
index b7534e259626..e36c0c84b42b 100644
--- a/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
+++ b/releasenotes/notes/cc-gate-fidelity-bde7a01dc8c56e29.yaml
@@ -21,7 +21,7 @@ fixes:
     should be treated as identity (such as a rotation gate with very small angle). See either
     of the aforementioned classes' docstrings for more information.
     Fixed `#13547 <https://github.com/Qiskit/qiskit/issues/13547>`__.
-features_circuit:
+features_circuits:
   - |
     Added a new argument ``approximation_degree`` to :meth:`.CommutationChecker.commute` and
     :meth:`.CommutationChecker.commute_nodes`, which allows to set the approximation threshold