From 6dca86e2200f1c27bd7e8aef47901a69e10317d5 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Tue, 16 Apr 2024 13:48:41 -0400
Subject: [PATCH 01/20] Oxidize the numeric code in the Isometry gate class

This commit ports the numeric portion of the Isometry gate class to
rust. While this will likely improve the performance slightly this move
is more to make isolate this code from blas/lapack in numpy. We're
hitting some stability issues on arm64 mac in CI and moving this code to
rust should hopefully fix this issue. As this is more for functional
reasons no real performance tuning was done on this port, there are
likely several opportunities to improve the runtime performance of the
code.
---
 Cargo.lock                                    |  10 +
 crates/accelerate/Cargo.toml                  |   1 +
 crates/accelerate/src/isometry.rs             | 375 ++++++++++++++++++
 crates/accelerate/src/lib.rs                  |   2 +
 crates/accelerate/src/two_qubit_decompose.rs  |   2 +-
 qiskit/__init__.py                            |   1 +
 .../library/generalized_gates/isometry.py     | 160 ++------
 7 files changed, 425 insertions(+), 126 deletions(-)
 create mode 100644 crates/accelerate/src/isometry.rs

diff --git a/Cargo.lock b/Cargo.lock
index 02bf4b2bb934..f4cd4458a301 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -589,6 +589,15 @@ dependencies = [
  "either",
 ]
 
+[[package]]
+name = "itertools"
+version = "0.12.1"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ba291022dbbd398a455acf126c1e341954079855bc60dfdda641363bd6922569"
+dependencies = [
+ "either",
+]
+
 [[package]]
 name = "jod-thread"
 version = "0.1.2"
@@ -1097,6 +1106,7 @@ dependencies = [
  "faer-ext",
  "hashbrown 0.14.3",
  "indexmap 2.2.6",
+ "itertools 0.12.1",
  "ndarray",
  "num-bigint",
  "num-complex",
diff --git a/crates/accelerate/Cargo.toml b/crates/accelerate/Cargo.toml
index 36fc5d408435..d7266fccc754 100644
--- a/crates/accelerate/Cargo.toml
+++ b/crates/accelerate/Cargo.toml
@@ -26,6 +26,7 @@ num-complex = "0.4"
 num-bigint = "0.4"
 rustworkx-core = "0.14"
 faer = "0.18.2"
+itertools = "0.12.1"
 
 [dependencies.smallvec]
 version = "1.13"
diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
new file mode 100644
index 000000000000..88d66e0136d5
--- /dev/null
+++ b/crates/accelerate/src/isometry.rs
@@ -0,0 +1,375 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2024
+//
+// 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 approx::abs_diff_eq;
+use num_complex::{Complex64, ComplexFloat};
+use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+use pyo3::Python;
+
+use hashbrown::HashSet;
+use itertools::Itertools;
+use ndarray::prelude::*;
+use numpy::{IntoPyArray, PyReadonlyArray2};
+use rayon::prelude::*;
+
+use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY;
+
+/// Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or
+/// basis_state=1 respectively
+#[pyfunction]
+pub fn reverse_qubit_state(
+    py: Python,
+    state: Vec<Complex64>,
+    basis_state: usize,
+    epsilon: f64,
+) -> PyObject {
+    reverse_qubit_state_inner(&state, basis_state, epsilon)
+        .into_pyarray_bound(py)
+        .into()
+}
+
+#[inline(always)]
+fn l2_norm(vec: &[Complex64]) -> f64 {
+    vec.iter()
+        .fold(0., |acc, elem| acc + elem.abs().powi(2))
+        .sqrt()
+}
+
+fn reverse_qubit_state_inner(
+    state: &[Complex64],
+    basis_state: usize,
+    epsilon: f64,
+) -> Array2<Complex64> {
+    let r = l2_norm(state);
+    if r < epsilon {
+        Array2::eye(2)
+    } else if basis_state == 0 {
+        array![
+            [state[0].conj() / r, state[1].conj() / r],
+            [-state[1] / r, state[0] / r],
+        ]
+    } else {
+        array![
+            [-state[1] / r, state[0] / r],
+            [state[0].conj() / r, state[1].conj() / r],
+        ]
+    }
+}
+
+#[pyfunction]
+pub fn find_squs_for_disentangling(
+    py: Python,
+    v: PyReadonlyArray2<Complex64>,
+    k: usize,
+    s: usize,
+    epsilon: f64,
+    n: usize,
+) -> Vec<PyObject> {
+    let v = v.as_array();
+    let k_prime = 0;
+    let i_start = if b(k, s + 1) == 0 {
+        a(k, s + 1)
+    } else {
+        a(k, s + 1) + 1
+    };
+    let mut output: Vec<Array2<Complex64>> = (0..i_start).map(|_| Array2::eye(2)).collect();
+    let mut squs: Vec<Array2<Complex64>> = (i_start..2_usize.pow((n - s - 1) as u32))
+        .map(|i| {
+            reverse_qubit_state_inner(
+                &[
+                    v[[2 * i * 2_usize.pow(s as u32) + b(k, s), k_prime]],
+                    v[[(2 * i + 1) * 2_usize.pow(s as u32) + b(k, s), k_prime]],
+                ],
+                k_s(k, s),
+                epsilon,
+            )
+        })
+        .collect();
+    output.append(&mut squs);
+    output
+        .into_iter()
+        .map(|x| x.into_pyarray_bound(py).into())
+        .collect()
+}
+
+#[pyfunction]
+pub fn apply_ucg(
+    py: Python,
+    m: PyReadonlyArray2<Complex64>,
+    k: usize,
+    single_qubit_gates: Vec<PyReadonlyArray2<Complex64>>,
+) -> PyObject {
+    let mut m = m.as_array().to_owned();
+    let shape = m.shape();
+    let num_qubits = shape[0].ilog2();
+    let num_col = shape[1];
+    let spacing: usize = 2_usize.pow(num_qubits - k as u32 - 1);
+    for j in 0..2_usize.pow(num_qubits - 1) {
+        let i = (j / spacing) * spacing + j;
+        let gate_index = i / (2_usize.pow(num_qubits - k as u32));
+        for col in 0..num_col {
+            let gate = single_qubit_gates[gate_index].as_array();
+            let a = m[[i, col]];
+            let b = m[[i + spacing, col]];
+            m[[i, col]] = gate[[0, 0]] * a + gate[[0, 1]] * b;
+            m[[i + spacing, col]] = gate[[1, 0]] * a + gate[[1, 1]] * b;
+        }
+    }
+    m.into_pyarray_bound(py).into()
+}
+
+#[inline(always)]
+fn bin_to_int(bin: &[u8]) -> usize {
+    bin.iter()
+        .fold(0_usize, |acc, digit| (acc << 1) + *digit as usize)
+}
+
+#[pyfunction]
+pub fn apply_diagonal_gate(
+    py: Python,
+    m: PyReadonlyArray2<Complex64>,
+    action_qubit_labels: Vec<usize>,
+    diag: Vec<Complex64>,
+) -> PyObject {
+    let mut m = m.as_array().to_owned();
+    let shape = m.shape();
+    let num_qubits = shape[0].ilog2();
+    let num_col = shape[1];
+    for state in std::iter::repeat([0_u8, 1_u8])
+        .take(num_qubits as usize)
+        .multi_cartesian_product()
+    {
+        let state_on_action_qubits: Vec<u8> =
+            action_qubit_labels.iter().map(|i| state[*i]).collect();
+        let diag_index = bin_to_int(&state_on_action_qubits);
+        let i = bin_to_int(&state);
+        for j in 0..num_col {
+            m[[i, j]] = diag[diag_index] * m[[i, j]]
+        }
+    }
+    m.into_pyarray_bound(py).into()
+}
+
+#[pyfunction]
+pub fn apply_diagonal_gate_to_diag(
+    mut m_diagonal: Vec<Complex64>,
+    action_qubit_labels: Vec<usize>,
+    diag: Vec<Complex64>,
+    num_qubits: usize,
+) -> Vec<Complex64> {
+    if m_diagonal.is_empty() {
+        return m_diagonal;
+    }
+    for state in std::iter::repeat([0_u8, 1_u8])
+        .take(num_qubits)
+        .multi_cartesian_product()
+        .take(m_diagonal.len())
+    {
+        let state_on_action_qubits: Vec<u8> =
+            action_qubit_labels.iter().map(|i| state[*i]).collect();
+        let diag_index = bin_to_int(&state_on_action_qubits);
+        let i = bin_to_int(&state);
+        m_diagonal[i] *= diag[diag_index]
+    }
+    m_diagonal
+}
+
+/// Helper method for _apply_multi_controlled_gate. This constructs the basis states the MG gate
+/// is acting on for a specific state state_free of the qubits we neither control nor act on
+fn construct_basis_states(
+    state_free: &[u8],
+    control_labels: &[usize],
+    target_label: usize,
+) -> [usize; 2] {
+    let size = state_free.len() + control_labels.len() + 1;
+    let mut e1: Vec<u8> = Vec::with_capacity(size);
+    let mut e2: Vec<u8> = Vec::with_capacity(size);
+    let mut j = 0;
+    let control_set: HashSet<usize> = control_labels.iter().copied().collect();
+    for i in 0..size {
+        if control_set.contains(&i) {
+            e1.push(1);
+            e2.push(1);
+        } else if i == target_label {
+            e1.push(0);
+            e2.push(1);
+        } else {
+            assert!(j <= 1);
+            e1.push(state_free[j]);
+            e2.push(state_free[j]);
+            j += 1
+        }
+    }
+    [bin_to_int(&e1), bin_to_int(&e2)]
+}
+
+#[pyfunction]
+pub fn apply_multi_controlled_gate(
+    py: Python,
+    m: PyReadonlyArray2<Complex64>,
+    mut control_labels: Vec<usize>,
+    target_label: usize,
+    gate: PyReadonlyArray2<Complex64>,
+) -> PyObject {
+    let mut m = m.as_array().to_owned();
+    let gate = gate.as_array();
+    let shape = m.shape();
+    let num_qubits = shape[0].ilog2();
+    let num_col = shape[1];
+    control_labels.par_sort();
+    let free_qubits = num_qubits as usize - control_labels.len() - 1;
+    if free_qubits == 0 {
+        let [e1, e2] = construct_basis_states(&[], &control_labels, target_label);
+        for i in 0..num_col {
+            let temp: Vec<_> = gate
+                .dot(&aview2(&[[m[[e1, i]]], [m[[e2, i]]]]))
+                .into_iter()
+                .take(2)
+                .collect();
+            m[[e1, i]] = temp[0];
+            m[[e2, i]] = temp[1];
+        }
+        return m.into_pyarray_bound(py).into();
+    }
+    for state_free in std::iter::repeat([0_u8, 1_u8])
+        .take(free_qubits)
+        .multi_cartesian_product()
+    {
+        let [e1, e2] = construct_basis_states(&state_free, &control_labels, target_label);
+        for i in 0..num_col {
+            let temp: Vec<_> = gate
+                .dot(&aview2(&[[m[[e1, i]]], [m[[e2, i]]]]))
+                .into_iter()
+                .take(2)
+                .collect();
+            m[[e1, i]] = temp[0];
+            m[[e2, i]] = temp[1];
+        }
+    }
+    m.into_pyarray_bound(py).into()
+}
+
+#[pyfunction]
+pub fn ucg_is_identity_up_to_global_phase(
+    single_qubit_gates: Vec<PyReadonlyArray2<Complex64>>,
+    epsilon: f64,
+) -> bool {
+    let global_phase: Complex64 = if single_qubit_gates[0].as_array()[[0, 0]].abs() >= epsilon {
+        single_qubit_gates[0].as_array()[[0, 0]].finv()
+    } else {
+        return false;
+    };
+    for raw_gate in single_qubit_gates {
+        let gate = raw_gate.as_array();
+        if !abs_diff_eq!(
+            gate.mapv(|x| x * global_phase),
+            aview2(&ONE_QUBIT_IDENTITY),
+            epsilon = 1e-8 // Default tolerance from numpy for allclose()
+        ) {
+            return false;
+        }
+    }
+    true
+}
+
+#[pyfunction]
+fn diag_is_identity_up_to_global_phase(diag: Vec<Complex64>, epsilon: f64) -> bool {
+    let global_phase: Complex64 = if diag[0].abs() >= epsilon {
+        diag[0].finv()
+    } else {
+        return false;
+    };
+    for d in diag {
+        if (global_phase * d - 1.0).abs() >= epsilon {
+            return false;
+        }
+    }
+    true
+}
+
+#[pyfunction]
+pub fn merge_ucgate_and_diag(
+    py: Python,
+    single_qubit_gates: Vec<PyReadonlyArray2<Complex64>>,
+    diag: Vec<Complex64>,
+) -> Vec<PyObject> {
+    single_qubit_gates
+        .iter()
+        .enumerate()
+        .map(|(i, raw_gate)| {
+            let gate = raw_gate.as_array();
+            let res = aview2(&[
+                [diag[2 * i], Complex64::new(0., 0.)],
+                [Complex64::new(0., 0.), diag[2 * i + 1]],
+            ])
+            .dot(&gate);
+            res.into_pyarray_bound(py).into()
+        })
+        .collect()
+}
+
+#[inline(always)]
+fn get_binary_rep_as_list(n: usize, num_digits: usize) -> Vec<u8> {
+    // TODO: Come up with a better implementation of this which doesn't involve
+    // strings or doubly reversing
+    let binary_str = format!("{:0num_digits$b}", n, num_digits = num_digits);
+    let mut res: Vec<u8> = binary_str
+        .chars()
+        .map(|c| c.to_digit(10).unwrap() as u8)
+        .rev()
+        .take(num_digits)
+        .collect();
+    res.reverse();
+    res
+}
+
+#[inline(always)]
+#[pyfunction]
+fn k_s(k: usize, s: usize) -> usize {
+    if k == 0 {
+        0
+    } else {
+        let num_digits = s + 1;
+        let res = get_binary_rep_as_list(k, num_digits);
+        res[0] as usize
+    }
+}
+
+#[inline(always)]
+#[pyfunction]
+fn a(k: usize, s: usize) -> usize {
+    k / 2_usize.pow(s as u32)
+}
+
+#[inline(always)]
+#[pyfunction]
+fn b(k: usize, s: usize) -> usize {
+    k - (a(k, s) * 2_usize.pow(s as u32))
+}
+
+#[pymodule]
+pub fn isometry(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_wrapped(wrap_pyfunction!(diag_is_identity_up_to_global_phase))?;
+    m.add_wrapped(wrap_pyfunction!(find_squs_for_disentangling))?;
+    m.add_wrapped(wrap_pyfunction!(reverse_qubit_state))?;
+    m.add_wrapped(wrap_pyfunction!(apply_ucg))?;
+    m.add_wrapped(wrap_pyfunction!(apply_diagonal_gate))?;
+    m.add_wrapped(wrap_pyfunction!(apply_diagonal_gate_to_diag))?;
+    m.add_wrapped(wrap_pyfunction!(apply_multi_controlled_gate))?;
+    m.add_wrapped(wrap_pyfunction!(ucg_is_identity_up_to_global_phase))?;
+    m.add_wrapped(wrap_pyfunction!(merge_ucgate_and_diag))?;
+    m.add_wrapped(wrap_pyfunction!(a))?;
+    m.add_wrapped(wrap_pyfunction!(b))?;
+    m.add_wrapped(wrap_pyfunction!(k_s))?;
+    Ok(())
+}
diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs
index 106764d4c0fb..6d8ee08e1e45 100644
--- a/crates/accelerate/src/lib.rs
+++ b/crates/accelerate/src/lib.rs
@@ -20,6 +20,7 @@ mod dense_layout;
 mod edge_collections;
 mod error_map;
 mod euler_one_qubit_decomposer;
+mod isometry;
 mod nlayout;
 mod optimize_1q_gates;
 mod pauli_exp_val;
@@ -62,6 +63,7 @@ fn _accelerate(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pymodule!(vf2_layout::vf2_layout))?;
     m.add_wrapped(wrap_pymodule!(two_qubit_decompose::two_qubit_decompose))?;
     m.add_wrapped(wrap_pymodule!(utils::utils))?;
+    m.add_wrapped(wrap_pymodule!(isometry::isometry))?;
     m.add_wrapped(wrap_pymodule!(
         euler_one_qubit_decomposer::euler_one_qubit_decomposer
     ))?;
diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs
index de3b80970f3c..7cebebc42dd6 100644
--- a/crates/accelerate/src/two_qubit_decompose.rs
+++ b/crates/accelerate/src/two_qubit_decompose.rs
@@ -58,7 +58,7 @@ const TWO_PI: f64 = 2.0 * PI;
 
 const C1: c64 = c64 { re: 1.0, im: 0.0 };
 
-static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [
+pub static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [
     [Complex64::new(1., 0.), Complex64::new(0., 0.)],
     [Complex64::new(0., 0.), Complex64::new(1., 0.)],
 ];
diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 9ef3487361c9..18ed22e208fa 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -71,6 +71,7 @@
 sys.modules["qiskit._accelerate.sampled_exp_val"] = qiskit._accelerate.sampled_exp_val
 sys.modules["qiskit._accelerate.vf2_layout"] = qiskit._accelerate.vf2_layout
 sys.modules["qiskit._accelerate.error_map"] = qiskit._accelerate.error_map
+sys.modules["qiskit._accelerate.isometry"] = qiskit._accelerate.isometry
 sys.modules["qiskit._accelerate.euler_one_qubit_decomposer"] = (
     qiskit._accelerate.euler_one_qubit_decomposer
 )
diff --git a/qiskit/circuit/library/generalized_gates/isometry.py b/qiskit/circuit/library/generalized_gates/isometry.py
index 1294feb26342..2bef0dcba543 100644
--- a/qiskit/circuit/library/generalized_gates/isometry.py
+++ b/qiskit/circuit/library/generalized_gates/isometry.py
@@ -30,6 +30,7 @@
 from qiskit.circuit.quantumregister import QuantumRegister
 from qiskit.exceptions import QiskitError
 from qiskit.quantum_info.operators.predicates import is_isometry
+from qiskit._accelerate import isometry as isometry_rs
 
 from .diagonal import Diagonal
 from .uc import UCGate
@@ -157,7 +158,9 @@ def _gates_to_uncompute(self):
         # correspond to the firstfew columns of the identity matrix up to diag, and hence we only
         # have to save a list containing them.
         for column_index in range(2**m):
-            self._decompose_column(circuit, q, diag, remaining_isometry, column_index)
+            remaining_isometry, diag = self._decompose_column(
+                circuit, q, diag, remaining_isometry, column_index
+            )
             # extract phase of the state that was sent to the basis state ket(column_index)
             diag.append(remaining_isometry[column_index, 0])
             # remove first column (which is now stored in diag)
@@ -173,7 +176,10 @@ def _decompose_column(self, circuit, q, diag, remaining_isometry, column_index):
         """
         n = int(math.log2(self.iso_data.shape[0]))
         for s in range(n):
-            self._disentangle(circuit, q, diag, remaining_isometry, column_index, s)
+            remaining_isometry, diag = self._disentangle(
+                circuit, q, diag, remaining_isometry, column_index, s
+            )
+        return remaining_isometry, diag
 
     def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
         """
@@ -205,12 +211,14 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
                 circuit, q, gate, control_labels, target_label
             )
             # apply the MCG to the remaining isometry
-            _apply_multi_controlled_gate(v, control_labels, target_label, gate)
+            v = _apply_multi_controlled_gate(v, control_labels, target_label, gate)
             # correct for the implementation "up to diagonal"
             diag_mcg_inverse = np.conj(diagonal_mcg).tolist()
-            _apply_diagonal_gate(v, control_labels + [target_label], diag_mcg_inverse)
+            v = _apply_diagonal_gate(v, control_labels + [target_label], diag_mcg_inverse)
             # update the diag according to the applied diagonal gate
-            _apply_diagonal_gate_to_diag(diag, control_labels + [target_label], diag_mcg_inverse, n)
+            diag = _apply_diagonal_gate_to_diag(
+                diag, control_labels + [target_label], diag_mcg_inverse, n
+            )
 
         # UCGate to disentangle a qubit:
         # Find the UCGate, decompose it and apply it to the remaining isometry
@@ -224,38 +232,24 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
             diagonal_ucg_inverse = np.conj(diagonal_ucg).tolist()
             single_qubit_gates = _merge_UCGate_and_diag(single_qubit_gates, diagonal_ucg_inverse)
             # apply the UCGate (with the merged diagonal gate) to the remaining isometry
-            _apply_ucg(v, len(control_labels), single_qubit_gates)
+            v = _apply_ucg(v, len(control_labels), single_qubit_gates)
             # update the diag according to the applied diagonal gate
-            _apply_diagonal_gate_to_diag(
+            diag = _apply_diagonal_gate_to_diag(
                 diag, control_labels + [target_label], diagonal_ucg_inverse, n
             )
             # # correct for the implementation "up to diagonal"
             # diag_inv = np.conj(diag).tolist()
             # _apply_diagonal_gate(v, control_labels + [target_label], diag_inv)
+        return v, diag
 
     # This method finds the single-qubit gates for a UCGate to disentangle a qubit:
     # we consider the n-qubit state v[:,0] starting with k zeros (in the computational basis).
     # The qubit with label n-s-1 is disentangled into the basis state k_s(k,s).
     def _find_squs_for_disentangling(self, v, k, s):
-        k_prime = 0
-        n = int(math.log2(self.iso_data.shape[0]))
-        if _b(k, s + 1) == 0:
-            i_start = _a(k, s + 1)
-        else:
-            i_start = _a(k, s + 1) + 1
-        id_list = [np.eye(2, 2) for _ in range(i_start)]
-        squs = [
-            _reverse_qubit_state(
-                [
-                    v[2 * i * 2**s + _b(k, s), k_prime],
-                    v[(2 * i + 1) * 2**s + _b(k, s), k_prime],
-                ],
-                _k_s(k, s),
-                self._epsilon,
-            )
-            for i in range(i_start, 2 ** (n - s - 1))
-        ]
-        return id_list + squs
+        res = isometry_rs.find_squs_for_disentangling(
+            v, k, s, self._epsilon, n=int(math.log2(self.iso_data.shape[0]))
+        )
+        return res
 
     # Append a UCGate up to diagonal to the circuit circ.
     def _append_ucg_up_to_diagonal(self, circ, q, single_qubit_gates, control_labels, target_label):
@@ -341,15 +335,8 @@ def inv_gate(self):
 # Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or
 # basis_state=1 respectively
 def _reverse_qubit_state(state, basis_state, epsilon):
-    state = np.array(state)
-    r = np.linalg.norm(state)
-    if r < epsilon:
-        return np.eye(2, 2)
-    if basis_state == 0:
-        m = np.array([[np.conj(state[0]), np.conj(state[1])], [-state[1], state[0]]]) / r
-    else:
-        m = np.array([[-state[1], state[0]], [np.conj(state[0]), np.conj(state[1])]]) / r
-    return m
+    res = isometry_rs.reverse_qubit_state(state, basis_state, epsilon)
+    return res
 
 
 # Methods for applying gates to matrices (should be moved to Qiskit AER)
@@ -373,17 +360,8 @@ def _reverse_qubit_state(state, basis_state, epsilon):
 def _apply_ucg(m, k, single_qubit_gates):
     # ToDo: Improve efficiency by parallelizing the gate application. A generalized version of
     # ToDo: this method should be implemented by the state vector simulator in Qiskit AER.
-    num_qubits = int(math.log2(m.shape[0]))
-    num_col = m.shape[1]
-    spacing = 2 ** (num_qubits - k - 1)
-    for j in range(2 ** (num_qubits - 1)):
-        i = (j // spacing) * spacing + j
-        gate_index = i // (2 ** (num_qubits - k))
-        for col in range(num_col):
-            m[np.array([i, i + spacing]), np.array([col, col])] = np.ndarray.flatten(
-                single_qubit_gates[gate_index].dot(np.array([[m[i, col]], [m[i + spacing, col]]]))
-            ).tolist()
-    return m
+    res = isometry_rs.apply_ucg(m, k, single_qubit_gates)
+    return res
 
 
 # Apply a diagonal gate with diagonal entries liste in diag and acting on qubits with labels
@@ -396,16 +374,8 @@ def _apply_ucg(m, k, single_qubit_gates):
 def _apply_diagonal_gate(m, action_qubit_labels, diag):
     # ToDo: Improve efficiency by parallelizing the gate application. A generalized version of
     # ToDo: this method should be implemented by the state vector simulator in Qiskit AER.
-    num_qubits = int(math.log2(m.shape[0]))
-    num_cols = m.shape[1]
-    basis_states = list(itertools.product([0, 1], repeat=num_qubits))
-    for state in basis_states:
-        state_on_action_qubits = [state[i] for i in action_qubit_labels]
-        diag_index = _bin_to_int(state_on_action_qubits)
-        i = _bin_to_int(state)
-        for j in range(num_cols):
-            m[i, j] = diag[diag_index] * m[i, j]
-    return m
+    res = isometry_rs.apply_diagonal_gate(m, action_qubit_labels, diag)
+    return res
 
 
 # Special case of the method _apply_diagonal_gate, where the input m is a diagonal matrix on the
@@ -417,15 +387,8 @@ def _apply_diagonal_gate(m, action_qubit_labels, diag):
 
 
 def _apply_diagonal_gate_to_diag(m_diagonal, action_qubit_labels, diag, num_qubits):
-    if not m_diagonal:
-        return m_diagonal
-    basis_states = list(itertools.product([0, 1], repeat=num_qubits))
-    for state in basis_states[: len(m_diagonal)]:
-        state_on_action_qubits = [state[i] for i in action_qubit_labels]
-        diag_index = _bin_to_int(state_on_action_qubits)
-        i = _bin_to_int(state)
-        m_diagonal[i] *= diag[diag_index]
-    return m_diagonal
+    res = isometry_rs.apply_diagonal_gate_to_diag(m_diagonal, action_qubit_labels, diag, num_qubits)
+    return res
 
 
 # Apply a MC single-qubit gate (given by the 2*2 unitary input: gate) with controlling on
@@ -436,43 +399,9 @@ def _apply_diagonal_gate_to_diag(m_diagonal, action_qubit_labels, diag, num_qubi
 
 
 def _apply_multi_controlled_gate(m, control_labels, target_label, gate):
-    # ToDo: This method should be integrated into the state vector simulator in Qiskit AER.
-    num_qubits = int(math.log2(m.shape[0]))
-    num_cols = m.shape[1]
-    control_labels.sort()
-    free_qubits = num_qubits - len(control_labels) - 1
-    basis_states_free = list(itertools.product([0, 1], repeat=free_qubits))
-    for state_free in basis_states_free:
-        (e1, e2) = _construct_basis_states(state_free, control_labels, target_label)
-        for i in range(num_cols):
-            m[np.array([e1, e2]), np.array([i, i])] = np.ndarray.flatten(
-                gate.dot(np.array([[m[e1, i]], [m[e2, i]]]))
-            ).tolist()
-    return m
-
-
-# Helper method for _apply_multi_controlled_gate. This constructs the basis states the MG gate
-# is acting on for a specific state state_free of the qubits we neither control nor act on.
-
-
-def _construct_basis_states(state_free, control_labels, target_label):
-    e1 = []
-    e2 = []
-    j = 0
-    for i in range(len(state_free) + len(control_labels) + 1):
-        if i in control_labels:
-            e1.append(1)
-            e2.append(1)
-        elif i == target_label:
-            e1.append(0)
-            e2.append(1)
-        else:
-            e1.append(state_free[j])
-            e2.append(state_free[j])
-            j += 1
-    out1 = _bin_to_int(e1)
-    out2 = _bin_to_int(e2)
-    return out1, out2
+    # ToDo: This method should be integrated into the state vector simulator in Qiskit AER
+    res = isometry_rs.apply_multi_controlled_gate(m, control_labels, target_label, gate)
+    return res
 
 
 # Some helper methods:
@@ -500,10 +429,6 @@ def _bin_to_int(binary_digits_as_list):
     return int("".join(str(x) for x in binary_digits_as_list), 2)
 
 
-def _ct(m):
-    return np.transpose(np.conjugate(m))
-
-
 def _get_binary_rep_as_list(n, num_digits):
     binary_string = np.binary_repr(n).zfill(num_digits)
     binary = []
@@ -517,9 +442,8 @@ def _get_binary_rep_as_list(n, num_digits):
 
 
 def _merge_UCGate_and_diag(single_qubit_gates, diag):
-    for i, gate in enumerate(single_qubit_gates):
-        single_qubit_gates[i] = np.array([[diag[2 * i], 0.0], [0.0, diag[2 * i + 1]]]).dot(gate)
-    return single_qubit_gates
+    res = isometry_rs.merge_ucgate_and_diag(single_qubit_gates, diag)
+    return res
 
 
 # Helper variables/functions for the column-by-column decomposition
@@ -553,22 +477,8 @@ def _k_s(k, s):
 
 
 def _ucg_is_identity_up_to_global_phase(single_qubit_gates, epsilon):
-    if not np.abs(single_qubit_gates[0][0, 0]) < epsilon:
-        global_phase = 1.0 / (single_qubit_gates[0][0, 0])
-    else:
-        return False
-    for gate in single_qubit_gates:
-        if not np.allclose(global_phase * gate, np.eye(2, 2)):
-            return False
-    return True
+    return isometry_rs.ucg_is_identity_up_to_global_phase(single_qubit_gates, epsilon)
 
 
 def _diag_is_identity_up_to_global_phase(diag, epsilon):
-    if not np.abs(diag[0]) < epsilon:
-        global_phase = 1.0 / (diag[0])
-    else:
-        return False
-    for d in diag:
-        if not np.abs(global_phase * d - 1) < epsilon:
-            return False
-    return True
+    return isometry_rs.diag_is_identity_up_to_global_phase(diag, epsilon)

From 5c21069de6504d2e11a2dd89a6cbb05012b26f27 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Wed, 17 Apr 2024 15:54:22 -0400
Subject: [PATCH 02/20] Remove unused import

---
 qiskit/circuit/library/generalized_gates/isometry.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/qiskit/circuit/library/generalized_gates/isometry.py b/qiskit/circuit/library/generalized_gates/isometry.py
index 2bef0dcba543..eb6e1cf988b6 100644
--- a/qiskit/circuit/library/generalized_gates/isometry.py
+++ b/qiskit/circuit/library/generalized_gates/isometry.py
@@ -21,7 +21,6 @@
 
 from __future__ import annotations
 
-import itertools
 import math
 import numpy as np
 from qiskit.circuit.exceptions import CircuitError

From cf3c5ff11591005790adc2d9684e66356e9b8529 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Wed, 17 Apr 2024 15:56:48 -0400
Subject: [PATCH 03/20] Use rust impl for small utility functions too

---
 qiskit/circuit/library/generalized_gates/isometry.py | 10 +++-------
 1 file changed, 3 insertions(+), 7 deletions(-)

diff --git a/qiskit/circuit/library/generalized_gates/isometry.py b/qiskit/circuit/library/generalized_gates/isometry.py
index eb6e1cf988b6..92b9170e75fd 100644
--- a/qiskit/circuit/library/generalized_gates/isometry.py
+++ b/qiskit/circuit/library/generalized_gates/isometry.py
@@ -453,11 +453,11 @@ def _merge_UCGate_and_diag(single_qubit_gates, diag):
 
 
 def _a(k, s):
-    return k // 2**s
+    return isometry_rs.a(k, s)
 
 
 def _b(k, s):
-    return k - (_a(k, s) * 2**s)
+    return isometry_rs.b(k, s)
 
 
 # given a binary representation of k with binary digits [k_{n-1},..,k_1,k_0],
@@ -465,11 +465,7 @@ def _b(k, s):
 
 
 def _k_s(k, s):
-    if k == 0:
-        return 0
-    else:
-        num_digits = s + 1
-        return _get_binary_rep_as_list(k, num_digits)[0]
+    return isometry_rs.k_s(k, s)
 
 
 # Check if a gate of a special form is equal to the identity gate up to global phase

From 8107fb250625a30fcbd62c973336df821a14d20c Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Wed, 17 Apr 2024 17:36:34 -0400
Subject: [PATCH 04/20] Oxidize the linalg in UCGate too

The UCGate class is used almost exclusively by the Isometry class to
build up the definition of the isometry circuit. There were also some
linear algebra inside the function which could also be the source of the
stability issues we were seeing on arm64. This commit ports this
function as part of the larger isometry migration.
---
 crates/accelerate/src/lib.rs                  |  2 +
 crates/accelerate/src/uc_gate.rs              | 85 +++++++++++++++++++
 qiskit/__init__.py                            |  1 +
 .../circuit/library/generalized_gates/uc.py   | 25 +-----
 4 files changed, 90 insertions(+), 23 deletions(-)
 create mode 100644 crates/accelerate/src/uc_gate.rs

diff --git a/crates/accelerate/src/lib.rs b/crates/accelerate/src/lib.rs
index 6d8ee08e1e45..8e5f50868154 100644
--- a/crates/accelerate/src/lib.rs
+++ b/crates/accelerate/src/lib.rs
@@ -31,6 +31,7 @@ mod sampled_exp_val;
 mod sparse_pauli_op;
 mod stochastic_swap;
 mod two_qubit_decompose;
+mod uc_gate;
 mod utils;
 mod vf2_layout;
 
@@ -64,6 +65,7 @@ fn _accelerate(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pymodule!(two_qubit_decompose::two_qubit_decompose))?;
     m.add_wrapped(wrap_pymodule!(utils::utils))?;
     m.add_wrapped(wrap_pymodule!(isometry::isometry))?;
+    m.add_wrapped(wrap_pymodule!(uc_gate::uc_gate))?;
     m.add_wrapped(wrap_pymodule!(
         euler_one_qubit_decomposer::euler_one_qubit_decomposer
     ))?;
diff --git a/crates/accelerate/src/uc_gate.rs b/crates/accelerate/src/uc_gate.rs
new file mode 100644
index 000000000000..8a68af11c986
--- /dev/null
+++ b/crates/accelerate/src/uc_gate.rs
@@ -0,0 +1,85 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2024
+//
+// 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 num_complex::{Complex64, ComplexFloat};
+use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+use pyo3::Python;
+use std::f64::consts::PI;
+
+use ndarray::prelude::*;
+use numpy::{IntoPyArray, PyReadonlyArray2};
+use faer_ext::{IntoFaerComplex, IntoNdarrayComplex};
+
+use crate::euler_one_qubit_decomposer::det_one_qubit;
+
+const PI2: f64 = PI / 2.;
+const EPS: f64 = 1e-10;
+
+/// This method implements the decomposition given in equation (3) in
+/// https://arxiv.org/pdf/quant-ph/0410066.pdf.
+///
+/// The decomposition is used recursively to decompose uniformly controlled gates.
+///
+/// a,b = single qubit unitaries
+/// v,u,r = outcome of the decomposition given in the reference mentioned above
+///
+/// (see there for the details).
+#[pyfunction]
+pub fn demultiplex_single_uc(
+    py: Python,
+    a: PyReadonlyArray2<Complex64>,
+    b: PyReadonlyArray2<Complex64>,
+) -> (PyObject, PyObject, PyObject) {
+    let a = a.as_array();
+    let b = b.as_array();
+    let x = a.dot(&b.mapv(|x| x.conj()).t());
+    let det_x = det_one_qubit(x.view());
+    let x11 = x[[0, 0]] / det_x.sqrt();
+    let phi = det_x.arg();
+
+    let r1 = (Complex64::new(0., 1.) / 2. * (PI2 - phi / 2. - x11.arg())).exp();
+    let r2 = (Complex64::new(0., 1.) / 2. * (PI2 - phi / 2. + x11.arg() + PI)).exp();
+
+    let r = array![
+        [r1, Complex64::new(0., 0.)],
+        [Complex64::new(0., 0.), r2],
+    ];
+
+    let decomp = r.dot(&x).dot(&r).view().into_faer_complex().complex_eigendecomposition();
+    let mut u: Array2<Complex64> = decomp.u().into_ndarray_complex().to_owned();
+    let s = decomp.s().column_vector();
+    let mut diag: Array1<Complex64> = Array1::from_shape_fn(u.shape()[0], |i| {
+        s[i].to_num_complex()
+    });
+
+    // If d is not equal to diag(i,-i), then we put it into this "standard" form
+    // (see eq. (13) in https://arxiv.org/pdf/quant-ph/0410066.pdf) by interchanging
+    // the eigenvalues and eigenvectors
+    if (diag[0] + Complex64::new(0., 1.)).abs() < EPS {
+        diag = diag.slice(s![..;-1]).to_owned();
+        u = u.slice(s![.., ..;-1]).to_owned();
+//        d = np.flip(d, 0);
+//        u = np.flip(u, 1);
+    }
+    diag.mapv_inplace(|x| x.sqrt());
+    let d = Array2::from_diag(&diag);
+    let v = d.dot(&u.mapv(|x| x.conj()).t()).dot(&r.mapv(|x| x.conj()).t()).dot(&b);
+    (v.into_pyarray_bound(py).into(), u.into_pyarray_bound(py).into(), r.into_pyarray_bound(py).into())
+}
+
+
+#[pymodule]
+pub fn uc_gate(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_wrapped(wrap_pyfunction!(demultiplex_single_uc))?;
+    Ok(())
+}
diff --git a/qiskit/__init__.py b/qiskit/__init__.py
index 18ed22e208fa..1fa3e4f564e4 100644
--- a/qiskit/__init__.py
+++ b/qiskit/__init__.py
@@ -72,6 +72,7 @@
 sys.modules["qiskit._accelerate.vf2_layout"] = qiskit._accelerate.vf2_layout
 sys.modules["qiskit._accelerate.error_map"] = qiskit._accelerate.error_map
 sys.modules["qiskit._accelerate.isometry"] = qiskit._accelerate.isometry
+sys.modules["qiskit._accelerate.uc_gate"] = qiskit._accelerate.uc_gate
 sys.modules["qiskit._accelerate.euler_one_qubit_decomposer"] = (
     qiskit._accelerate.euler_one_qubit_decomposer
 )
diff --git a/qiskit/circuit/library/generalized_gates/uc.py b/qiskit/circuit/library/generalized_gates/uc.py
index 2d650e98466a..19280f835556 100644
--- a/qiskit/circuit/library/generalized_gates/uc.py
+++ b/qiskit/circuit/library/generalized_gates/uc.py
@@ -21,7 +21,6 @@
 
 from __future__ import annotations
 
-import cmath
 import math
 
 import numpy as np
@@ -33,14 +32,11 @@
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.exceptions import QiskitError
-
-# pylint: disable=cyclic-import
-from qiskit.synthesis.one_qubit.one_qubit_decompose import OneQubitEulerDecomposer
+from qiskit._accelerate import uc_gate
 
 from .diagonal import Diagonal
 
 _EPS = 1e-10  # global variable used to chop very small numbers to zero
-_DECOMPOSER1Q = OneQubitEulerDecomposer("U3")
 
 
 class UCGate(Gate):
@@ -274,24 +270,7 @@ def _demultiplex_single_uc(self, a, b):
         v,u,r = outcome of the decomposition given in the reference mentioned above
         (see there for the details).
         """
-        # The notation is chosen as in https://arxiv.org/pdf/quant-ph/0410066.pdf.
-        x = a.dot(UCGate._ct(b))
-        det_x = np.linalg.det(x)
-        x11 = x.item((0, 0)) / cmath.sqrt(det_x)
-        phi = cmath.phase(det_x)
-        r1 = cmath.exp(1j / 2 * (np.pi / 2 - phi / 2 - cmath.phase(x11)))
-        r2 = cmath.exp(1j / 2 * (np.pi / 2 - phi / 2 + cmath.phase(x11) + np.pi))
-        r = np.array([[r1, 0], [0, r2]], dtype=complex)
-        d, u = np.linalg.eig(r.dot(x).dot(r))
-        # If d is not equal to diag(i,-i), then we put it into this "standard" form
-        # (see eq. (13) in https://arxiv.org/pdf/quant-ph/0410066.pdf) by interchanging
-        # the eigenvalues and eigenvectors.
-        if abs(d[0] + 1j) < _EPS:
-            d = np.flip(d, 0)
-            u = np.flip(u, 1)
-        d = np.diag(np.sqrt(d))
-        v = d.dot(UCGate._ct(u)).dot(UCGate._ct(r)).dot(b)
-        return v, u, r
+        return uc_gate.demultiplex_single_uc(a, b)
 
     @staticmethod
     def _ct(m):

From b0b7a3382913aa203c3c9f13e73208fc1aa136bc Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Wed, 17 Apr 2024 18:57:38 -0400
Subject: [PATCH 05/20] Migrate another numeric helper method of UCGate

---
 crates/accelerate/src/uc_gate.rs              | 95 ++++++++++++++++++-
 .../circuit/library/generalized_gates/uc.py   | 62 +-----------
 2 files changed, 91 insertions(+), 66 deletions(-)

diff --git a/crates/accelerate/src/uc_gate.rs b/crates/accelerate/src/uc_gate.rs
index 47f80494a24d..8005a21d19cc 100644
--- a/crates/accelerate/src/uc_gate.rs
+++ b/crates/accelerate/src/uc_gate.rs
@@ -42,6 +42,18 @@ pub fn demultiplex_single_uc(
 ) -> (PyObject, PyObject, PyObject) {
     let a = a.as_array();
     let b = b.as_array();
+    let [v, u, r] = demultiplex_single_uc_inner(a, b);
+    (
+        v.into_pyarray_bound(py).into(),
+        u.into_pyarray_bound(py).into(),
+        r.into_pyarray_bound(py).into(),
+    )
+}
+
+fn demultiplex_single_uc_inner(
+    a: ArrayView2<Complex64>,
+    b: ArrayView2<Complex64>,
+) -> [Array2<Complex64>; 3] {
     let x = a.dot(&b.mapv(|x| x.conj()).t());
     let det_x = det_one_qubit(x.view());
     let x11 = x[[0, 0]] / det_x.sqrt();
@@ -69,8 +81,6 @@ pub fn demultiplex_single_uc(
     if (diag[0] + Complex64::new(0., 1.)).abs() < EPS {
         diag = diag.slice(s![..;-1]).to_owned();
         u = u.slice(s![.., ..;-1]).to_owned();
-        //        d = np.flip(d, 0);
-        //        u = np.flip(u, 1);
     }
     diag.mapv_inplace(|x| x.sqrt());
     let d = Array2::from_diag(&diag);
@@ -78,15 +88,90 @@ pub fn demultiplex_single_uc(
         .dot(&u.mapv(|x| x.conj()).t())
         .dot(&r.mapv(|x| x.conj()).t())
         .dot(&b);
+    [v, u, r]
+}
+
+#[pyfunction]
+pub fn dec_ucg_help(
+    py: Python,
+    sq_gates: Vec<PyReadonlyArray2<Complex64>>,
+    num_qubits: u32,
+) -> (Vec<PyObject>, PyObject) {
+    let mut single_qubit_gates: Vec<Array2<Complex64>> = sq_gates
+        .into_iter()
+        .map(|x| x.as_array().to_owned())
+        .collect();
+    let mut diag: Array1<Complex64> = Array1::ones(2_usize.pow(num_qubits));
+    let num_controls = num_qubits - 1;
+    for dec_step in 0..num_controls {
+        let num_ucgs = 2_u32.pow(dec_step);
+        // The decomposition works recursively and the followign loop goes over the different
+        // UCGates that arise in the decomposition
+        for ucg_index in 0..num_ucgs {
+            let len_ucg = 2_u32.pow(num_controls - dec_step);
+            for i in 0..len_ucg / 2 {
+                let shift = ucg_index * len_ucg;
+                let a = single_qubit_gates[(shift + i) as usize].view();
+                let b = single_qubit_gates[(shift + len_ucg / 2 + i) as usize].view();
+                // Apply the decomposition for UCGates given in equation (3) in
+                // https://arxiv.org/pdf/quant-ph/0410066.pdf
+                // to demultiplex one control of all the num_ucgs uniformly-controlled gates
+                // with log2(len_ucg) uniform controls
+                let [v, u, r] = demultiplex_single_uc_inner(a, b);
+                // replace the single-qubit gates with v,u (the already existing ones
+                // are not needed any more)
+                single_qubit_gates[(shift + i) as usize] = v;
+                single_qubit_gates[(shift + len_ucg / 2 + i) as usize] = u;
+                // Now we decompose the gates D as described in Figure 4 in
+                // https://arxiv.org/pdf/quant-ph/0410066.pdf and merge some of the gates
+                // into the UCGates and the diagonal at the end of the circuit
+
+                // Remark: The Rz(pi/2) rotation acting on the target qubit and the Hadamard
+                // gates arising in the decomposition of D are ignored for the moment (they will
+                // be added together with the C-NOT gates at the end of the decomposition
+                // (in the method dec_ucg()))
+                let rz_11 = (-Complex64::new(0., 0.5 * PI2)).exp();
+                let rz_00 = Complex64::new(0., 0.5 * PI2).exp();
+                let r_conj_t = r.mapv(|x| x.conj()).t().to_owned();
+                if ucg_index < num_ucgs - 1 {
+                    // Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and
+                    // merge the UC-Rz rotation with the following UCGate,
+                    // which hasn't been decomposed yet
+                    let k = (shift + len_ucg + i) as usize;
+
+                    single_qubit_gates[k] = single_qubit_gates[k].dot(&r_conj_t);
+                    single_qubit_gates[k].mapv_inplace(|x| x * rz_00);
+                    let k = k + len_ucg as usize / 2;
+                    single_qubit_gates[k] = single_qubit_gates[k].dot(&r);
+                    single_qubit_gates[k].mapv_inplace(|x| x * rz_11);
+                } else {
+                    // Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and merge
+                    // the trailing UC-Rz rotation into a diagonal gate at the end of the circuit
+                    for ucg_index_2 in 0..num_ucgs {
+                        let shift_2 = ucg_index_2 * len_ucg;
+                        let k = (2 * (i + shift_2)) as usize;
+                        diag[k] *= r_conj_t[[0, 0]] * rz_00;
+                        diag[k + 1] *= r_conj_t[[1, 1]] * rz_00;
+                        let k = len_ucg as usize + k;
+                        diag[k] *= r[[0, 0]] * rz_11;
+                        diag[k + 1] *= r[[1, 1]] * rz_11;
+                    }
+                }
+            }
+        }
+    }
     (
-        v.into_pyarray_bound(py).into(),
-        u.into_pyarray_bound(py).into(),
-        r.into_pyarray_bound(py).into(),
+        single_qubit_gates
+            .into_iter()
+            .map(|x| x.into_pyarray_bound(py).into())
+            .collect(),
+        diag.into_pyarray_bound(py).into(),
     )
 }
 
 #[pymodule]
 pub fn uc_gate(m: &Bound<PyModule>) -> PyResult<()> {
     m.add_wrapped(wrap_pyfunction!(demultiplex_single_uc))?;
+    m.add_wrapped(wrap_pyfunction!(dec_ucg_help))?;
     Ok(())
 }
diff --git a/qiskit/circuit/library/generalized_gates/uc.py b/qiskit/circuit/library/generalized_gates/uc.py
index 19280f835556..f38fa5f729c9 100644
--- a/qiskit/circuit/library/generalized_gates/uc.py
+++ b/qiskit/circuit/library/generalized_gates/uc.py
@@ -199,67 +199,7 @@ def _dec_ucg_help(self):
         https://arxiv.org/pdf/quant-ph/0410066.pdf.
         """
         single_qubit_gates = [gate.astype(complex) for gate in self.params]
-        diag = np.ones(2**self.num_qubits, dtype=complex)
-        num_contr = self.num_qubits - 1
-        for dec_step in range(num_contr):
-            num_ucgs = 2**dec_step
-            # The decomposition works recursively and the following loop goes over the different
-            # UCGates that arise in the decomposition
-            for ucg_index in range(num_ucgs):
-                len_ucg = 2 ** (num_contr - dec_step)
-                for i in range(int(len_ucg / 2)):
-                    shift = ucg_index * len_ucg
-                    a = single_qubit_gates[shift + i]
-                    b = single_qubit_gates[shift + len_ucg // 2 + i]
-                    # Apply the decomposition for UCGates given in equation (3) in
-                    # https://arxiv.org/pdf/quant-ph/0410066.pdf
-                    # to demultiplex one control of all the num_ucgs uniformly-controlled gates
-                    #  with log2(len_ucg) uniform controls
-                    v, u, r = self._demultiplex_single_uc(a, b)
-                    #  replace the single-qubit gates with v,u (the already existing ones
-                    #  are not needed any more)
-                    single_qubit_gates[shift + i] = v
-                    single_qubit_gates[shift + len_ucg // 2 + i] = u
-                    # Now we decompose the gates D as described in Figure 4  in
-                    # https://arxiv.org/pdf/quant-ph/0410066.pdf and merge some of the gates
-                    # into the UCGates and the diagonal at the end of the circuit
-
-                    # Remark: The Rz(pi/2) rotation acting on the target qubit and the Hadamard
-                    # gates arising in the decomposition of D are ignored for the moment (they will
-                    # be added together with the C-NOT gates at the end of the decomposition
-                    # (in the method dec_ucg()))
-                    if ucg_index < num_ucgs - 1:
-                        # Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and
-                        # merge the UC-Rz rotation with the following UCGate,
-                        # which hasn't been decomposed yet.
-                        k = shift + len_ucg + i
-                        single_qubit_gates[k] = single_qubit_gates[k].dot(
-                            UCGate._ct(r)
-                        ) * UCGate._rz(np.pi / 2).item((0, 0))
-                        k = k + len_ucg // 2
-                        single_qubit_gates[k] = single_qubit_gates[k].dot(r) * UCGate._rz(
-                            np.pi / 2
-                        ).item((1, 1))
-                    else:
-                        # Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and merge
-                        # the trailing UC-Rz rotation into a diagonal gate at the end of the circuit
-                        for ucg_index_2 in range(num_ucgs):
-                            shift_2 = ucg_index_2 * len_ucg
-                            k = 2 * (i + shift_2)
-                            diag[k] = (
-                                diag[k]
-                                * UCGate._ct(r).item((0, 0))
-                                * UCGate._rz(np.pi / 2).item((0, 0))
-                            )
-                            diag[k + 1] = (
-                                diag[k + 1]
-                                * UCGate._ct(r).item((1, 1))
-                                * UCGate._rz(np.pi / 2).item((0, 0))
-                            )
-                            k = len_ucg + k
-                            diag[k] *= r.item((0, 0)) * UCGate._rz(np.pi / 2).item((1, 1))
-                            diag[k + 1] *= r.item((1, 1)) * UCGate._rz(np.pi / 2).item((1, 1))
-        return single_qubit_gates, diag
+        return uc_gate.dec_ucg_help(single_qubit_gates, self.num_qubits)
 
     def _demultiplex_single_uc(self, a, b):
         """

From eaf43f023ba7adf1a38dc7315a9041c32e8a4abb Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Wed, 17 Apr 2024 20:05:47 -0400
Subject: [PATCH 06/20] Remove now unused code paths

---
 crates/accelerate/src/uc_gate.rs              | 21 ++-----------------
 .../library/generalized_gates/isometry.py     |  4 ----
 .../circuit/library/generalized_gates/uc.py   | 15 -------------
 3 files changed, 2 insertions(+), 38 deletions(-)

diff --git a/crates/accelerate/src/uc_gate.rs b/crates/accelerate/src/uc_gate.rs
index 8005a21d19cc..111ae97c68b4 100644
--- a/crates/accelerate/src/uc_gate.rs
+++ b/crates/accelerate/src/uc_gate.rs
@@ -34,23 +34,7 @@ const EPS: f64 = 1e-10;
 /// v,u,r = outcome of the decomposition given in the reference mentioned above
 ///
 /// (see there for the details).
-#[pyfunction]
-pub fn demultiplex_single_uc(
-    py: Python,
-    a: PyReadonlyArray2<Complex64>,
-    b: PyReadonlyArray2<Complex64>,
-) -> (PyObject, PyObject, PyObject) {
-    let a = a.as_array();
-    let b = b.as_array();
-    let [v, u, r] = demultiplex_single_uc_inner(a, b);
-    (
-        v.into_pyarray_bound(py).into(),
-        u.into_pyarray_bound(py).into(),
-        r.into_pyarray_bound(py).into(),
-    )
-}
-
-fn demultiplex_single_uc_inner(
+fn demultiplex_single_uc(
     a: ArrayView2<Complex64>,
     b: ArrayView2<Complex64>,
 ) -> [Array2<Complex64>; 3] {
@@ -117,7 +101,7 @@ pub fn dec_ucg_help(
                 // https://arxiv.org/pdf/quant-ph/0410066.pdf
                 // to demultiplex one control of all the num_ucgs uniformly-controlled gates
                 // with log2(len_ucg) uniform controls
-                let [v, u, r] = demultiplex_single_uc_inner(a, b);
+                let [v, u, r] = demultiplex_single_uc(a, b);
                 // replace the single-qubit gates with v,u (the already existing ones
                 // are not needed any more)
                 single_qubit_gates[(shift + i) as usize] = v;
@@ -171,7 +155,6 @@ pub fn dec_ucg_help(
 
 #[pymodule]
 pub fn uc_gate(m: &Bound<PyModule>) -> PyResult<()> {
-    m.add_wrapped(wrap_pyfunction!(demultiplex_single_uc))?;
     m.add_wrapped(wrap_pyfunction!(dec_ucg_help))?;
     Ok(())
 }
diff --git a/qiskit/circuit/library/generalized_gates/isometry.py b/qiskit/circuit/library/generalized_gates/isometry.py
index 92b9170e75fd..024ce4fca2f2 100644
--- a/qiskit/circuit/library/generalized_gates/isometry.py
+++ b/qiskit/circuit/library/generalized_gates/isometry.py
@@ -424,10 +424,6 @@ def _reverse_qubit_oder(qubits):
 # Convert list of binary digits to integer
 
 
-def _bin_to_int(binary_digits_as_list):
-    return int("".join(str(x) for x in binary_digits_as_list), 2)
-
-
 def _get_binary_rep_as_list(n, num_digits):
     binary_string = np.binary_repr(n).zfill(num_digits)
     binary = []
diff --git a/qiskit/circuit/library/generalized_gates/uc.py b/qiskit/circuit/library/generalized_gates/uc.py
index f38fa5f729c9..f54567123e02 100644
--- a/qiskit/circuit/library/generalized_gates/uc.py
+++ b/qiskit/circuit/library/generalized_gates/uc.py
@@ -201,21 +201,6 @@ def _dec_ucg_help(self):
         single_qubit_gates = [gate.astype(complex) for gate in self.params]
         return uc_gate.dec_ucg_help(single_qubit_gates, self.num_qubits)
 
-    def _demultiplex_single_uc(self, a, b):
-        """
-        This method implements the decomposition given in equation (3) in
-        https://arxiv.org/pdf/quant-ph/0410066.pdf.
-        The decomposition is used recursively to decompose uniformly controlled gates.
-        a,b = single qubit unitaries
-        v,u,r = outcome of the decomposition given in the reference mentioned above
-        (see there for the details).
-        """
-        return uc_gate.demultiplex_single_uc(a, b)
-
-    @staticmethod
-    def _ct(m):
-        return np.transpose(np.conjugate(m))
-
     @staticmethod
     def _rz(alpha):
         return np.array([[np.exp(1j * alpha / 2), 0], [0, np.exp(-1j * alpha / 2)]])

From 23f5e2e61b6a1faf81682c8e7dd72025f8d6298e Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 19 Apr 2024 07:48:19 -0400
Subject: [PATCH 07/20] Remove bitstring usage with bitwise ops

This commit removes the use of bit string manipulations that were
faithfully ported from the original python logic (but left a bad taste
in my mouth) into more efficient bitwise operations (which were
possible in the original python too).
---
 crates/accelerate/src/isometry.rs | 22 ++++------------------
 1 file changed, 4 insertions(+), 18 deletions(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 88d66e0136d5..579fbc6db9e5 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.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 std::ops::BitAnd;
+
 use approx::abs_diff_eq;
 use num_complex::{Complex64, ComplexFloat};
 use pyo3::prelude::*;
@@ -318,30 +320,14 @@ pub fn merge_ucgate_and_diag(
         .collect()
 }
 
-#[inline(always)]
-fn get_binary_rep_as_list(n: usize, num_digits: usize) -> Vec<u8> {
-    // TODO: Come up with a better implementation of this which doesn't involve
-    // strings or doubly reversing
-    let binary_str = format!("{:0num_digits$b}", n, num_digits = num_digits);
-    let mut res: Vec<u8> = binary_str
-        .chars()
-        .map(|c| c.to_digit(10).unwrap() as u8)
-        .rev()
-        .take(num_digits)
-        .collect();
-    res.reverse();
-    res
-}
-
 #[inline(always)]
 #[pyfunction]
 fn k_s(k: usize, s: usize) -> usize {
     if k == 0 {
         0
     } else {
-        let num_digits = s + 1;
-        let res = get_binary_rep_as_list(k, num_digits);
-        res[0] as usize
+        let filter = 1 << s;
+        k.bitand(filter) >> s
     }
 }
 

From cb8fc67d9bf888dd99a394a70d82c868b8a0a3c5 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 19 Apr 2024 08:05:06 -0400
Subject: [PATCH 08/20] Mostly replace Vec<u8> usage with bitwise operations

The use of intermediate Vec<u8> as proxy bitstrings was originally
ported nearly exactly from the python implementation. But since
everything is working now this commit switches to use bitwise operations
where it makes sense as this will be more efficient.
---
 crates/accelerate/src/isometry.rs | 35 ++++++++++++++++++-------------
 1 file changed, 20 insertions(+), 15 deletions(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 579fbc6db9e5..5f87edab97a5 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -151,9 +151,9 @@ pub fn apply_diagonal_gate(
         .take(num_qubits as usize)
         .multi_cartesian_product()
     {
-        let state_on_action_qubits: Vec<u8> =
-            action_qubit_labels.iter().map(|i| state[*i]).collect();
-        let diag_index = bin_to_int(&state_on_action_qubits);
+        let diag_index = action_qubit_labels
+            .iter()
+            .fold(0_usize, |acc, i| (acc << 1) + state[*i] as usize);
         let i = bin_to_int(&state);
         for j in 0..num_col {
             m[[i, j]] = diag[diag_index] * m[[i, j]]
@@ -177,9 +177,9 @@ pub fn apply_diagonal_gate_to_diag(
         .multi_cartesian_product()
         .take(m_diagonal.len())
     {
-        let state_on_action_qubits: Vec<u8> =
-            action_qubit_labels.iter().map(|i| state[*i]).collect();
-        let diag_index = bin_to_int(&state_on_action_qubits);
+        let diag_index = action_qubit_labels
+            .iter()
+            .fold(0_usize, |acc, i| (acc << 1) + state[*i] as usize);
         let i = bin_to_int(&state);
         m_diagonal[i] *= diag[diag_index]
     }
@@ -194,25 +194,30 @@ fn construct_basis_states(
     target_label: usize,
 ) -> [usize; 2] {
     let size = state_free.len() + control_labels.len() + 1;
-    let mut e1: Vec<u8> = Vec::with_capacity(size);
-    let mut e2: Vec<u8> = Vec::with_capacity(size);
+    let mut e1: usize = 0;
+    let mut e2: usize = 0;
     let mut j = 0;
     let control_set: HashSet<usize> = control_labels.iter().copied().collect();
     for i in 0..size {
         if control_set.contains(&i) {
-            e1.push(1);
-            e2.push(1);
+            e1 <<= 1;
+            e1 += 1;
+            e2 <<= 1;
+            e2 += 1;
         } else if i == target_label {
-            e1.push(0);
-            e2.push(1);
+            e1 <<= 1;
+            e2 <<= 1;
+            e2 += 1;
         } else {
             assert!(j <= 1);
-            e1.push(state_free[j]);
-            e2.push(state_free[j]);
+            e1 <<= 1;
+            e1 += state_free[j] as usize;
+            e2 <<= 1;
+            e2 += state_free[j] as usize;
             j += 1
         }
     }
-    [bin_to_int(&e1), bin_to_int(&e2)]
+    [e1, e2]
 }
 
 #[pyfunction]

From 6a6f3bcb99a4acbfa6bd2bace69ed913e7b43041 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 06:15:56 -0400
Subject: [PATCH 09/20] Apply suggestions from code review

Co-authored-by: Jake Lishman <jake@binhbar.com>
---
 crates/accelerate/src/isometry.rs | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 5f87edab97a5..f3821520c484 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -68,6 +68,10 @@ fn reverse_qubit_state_inner(
     }
 }
 
+/// This method finds the single-qubit gates for a UCGate to disentangle a qubit:
+/// we consider the n-qubit state v[:,0] starting with k zeros (in the computational basis).
+/// The qubit with label n-s-1 is disentangled into the basis state k_s(k,s).
+
 #[pyfunction]
 pub fn find_squs_for_disentangling(
     py: Python,

From 6e86efe470ba617e014aabd4c350c543ef5ecaf5 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 06:40:41 -0400
Subject: [PATCH 10/20] Remove python side call sites

---
 .../library/generalized_gates/isometry.py     | 153 +++---------------
 1 file changed, 24 insertions(+), 129 deletions(-)

diff --git a/qiskit/circuit/library/generalized_gates/isometry.py b/qiskit/circuit/library/generalized_gates/isometry.py
index 024ce4fca2f2..2c597b7507e0 100644
--- a/qiskit/circuit/library/generalized_gates/isometry.py
+++ b/qiskit/circuit/library/generalized_gates/isometry.py
@@ -164,7 +164,9 @@ def _gates_to_uncompute(self):
             diag.append(remaining_isometry[column_index, 0])
             # remove first column (which is now stored in diag)
             remaining_isometry = remaining_isometry[:, 1:]
-        if len(diag) > 1 and not _diag_is_identity_up_to_global_phase(diag, self._epsilon):
+        if len(diag) > 1 and not isometry_rs.diag_is_identity_up_to_global_phase(
+            diag, self._epsilon
+        ):
             diagonal = Diagonal(np.conj(diag))
             circuit.append(diagonal, q_input)
         return circuit
@@ -194,13 +196,19 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
         n = int(math.log2(self.iso_data.shape[0]))
 
         # MCG to set one entry to zero (preparation for disentangling with UCGate):
-        index1 = 2 * _a(k, s + 1) * 2**s + _b(k, s + 1)
-        index2 = (2 * _a(k, s + 1) + 1) * 2**s + _b(k, s + 1)
+        index1 = 2 * isometry_rs.a(k, s + 1) * 2**s + isometry_rs.b(k, s + 1)
+        index2 = (2 * isometry_rs.a(k, s + 1) + 1) * 2**s + isometry_rs.b(k, s + 1)
         target_label = n - s - 1
         # Check if a MCG is required
-        if _k_s(k, s) == 0 and _b(k, s + 1) != 0 and np.abs(v[index2, k_prime]) > self._epsilon:
+        if (
+            isometry_rs.k_s(k, s) == 0
+            and isometry_rs.b(k, s + 1) != 0
+            and np.abs(v[index2, k_prime]) > self._epsilon
+        ):
             # Find the MCG, decompose it and apply it to the remaining isometry
-            gate = _reverse_qubit_state([v[index1, k_prime], v[index2, k_prime]], 0, self._epsilon)
+            gate = isometry_rs.reverse_qubit_state(
+                [v[index1, k_prime], v[index2, k_prime]], 0, self._epsilon
+            )
             control_labels = [
                 i
                 for i, x in enumerate(_get_binary_rep_as_list(k, n))
@@ -210,30 +218,34 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
                 circuit, q, gate, control_labels, target_label
             )
             # apply the MCG to the remaining isometry
-            v = _apply_multi_controlled_gate(v, control_labels, target_label, gate)
+            v = isometry_rs.apply_multi_controlled_gate(v, control_labels, target_label, gate)
             # correct for the implementation "up to diagonal"
             diag_mcg_inverse = np.conj(diagonal_mcg).tolist()
-            v = _apply_diagonal_gate(v, control_labels + [target_label], diag_mcg_inverse)
+            v = isometry_rs.apply_diagonal_gate(
+                v, control_labels + [target_label], diag_mcg_inverse
+            )
             # update the diag according to the applied diagonal gate
-            diag = _apply_diagonal_gate_to_diag(
+            diag = isometry_rs.apply_diagonal_gate_to_diag(
                 diag, control_labels + [target_label], diag_mcg_inverse, n
             )
 
         # UCGate to disentangle a qubit:
         # Find the UCGate, decompose it and apply it to the remaining isometry
         single_qubit_gates = self._find_squs_for_disentangling(v, k, s)
-        if not _ucg_is_identity_up_to_global_phase(single_qubit_gates, self._epsilon):
+        if not isometry_rs.ucg_is_identity_up_to_global_phase(single_qubit_gates, self._epsilon):
             control_labels = list(range(target_label))
             diagonal_ucg = self._append_ucg_up_to_diagonal(
                 circuit, q, single_qubit_gates, control_labels, target_label
             )
             # merge the diagonal into the UCGate for efficient application of both together
             diagonal_ucg_inverse = np.conj(diagonal_ucg).tolist()
-            single_qubit_gates = _merge_UCGate_and_diag(single_qubit_gates, diagonal_ucg_inverse)
+            single_qubit_gates = isometry_rs.merge_ucgate_and_diag(
+                single_qubit_gates, diagonal_ucg_inverse
+            )
             # apply the UCGate (with the merged diagonal gate) to the remaining isometry
-            v = _apply_ucg(v, len(control_labels), single_qubit_gates)
+            v = isometry_rs.apply_ucg(v, len(control_labels), single_qubit_gates)
             # update the diag according to the applied diagonal gate
-            diag = _apply_diagonal_gate_to_diag(
+            diag = isometry_rs.apply_diagonal_gate_to_diag(
                 diag, control_labels + [target_label], diagonal_ucg_inverse, n
             )
             # # correct for the implementation "up to diagonal"
@@ -331,81 +343,6 @@ def inv_gate(self):
         return self._inverse
 
 
-# Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or
-# basis_state=1 respectively
-def _reverse_qubit_state(state, basis_state, epsilon):
-    res = isometry_rs.reverse_qubit_state(state, basis_state, epsilon)
-    return res
-
-
-# Methods for applying gates to matrices (should be moved to Qiskit AER)
-
-# Input: matrix m with 2^n rows (and arbitrary many columns). Think of the columns as states
-#  on n qubits. The method applies a uniformly controlled gate (UCGate) to all the columns, where
-#  the UCGate is specified by the inputs k and single_qubit_gates:
-
-#  k =  number of controls. We assume that the controls are on the k most significant qubits
-#       (and the target is on the (k+1)th significant qubit)
-#  single_qubit_gates =     [u_0,...,u_{2^k-1}], where the u_i's are 2*2 unitaries
-#                           (provided as numpy arrays)
-
-# The order of the single-qubit unitaries is such that the first unitary u_0 is applied to the
-# (k+1)th significant qubit if the control qubits are in the state ket(0...00), the gate u_1 is
-# applied if the control qubits are in the state ket(0...01), and so on.
-
-# The input matrix m and the single-qubit gates have to be of dtype=complex.
-
-
-def _apply_ucg(m, k, single_qubit_gates):
-    # ToDo: Improve efficiency by parallelizing the gate application. A generalized version of
-    # ToDo: this method should be implemented by the state vector simulator in Qiskit AER.
-    res = isometry_rs.apply_ucg(m, k, single_qubit_gates)
-    return res
-
-
-# Apply a diagonal gate with diagonal entries liste in diag and acting on qubits with labels
-#  action_qubit_labels to a matrix m.
-# The input matrix m has to be of dtype=complex
-# The qubit labels are such that label 0 corresponds to the most significant qubit, label 1 to
-#  the second most significant qubit, and so on ...
-
-
-def _apply_diagonal_gate(m, action_qubit_labels, diag):
-    # ToDo: Improve efficiency by parallelizing the gate application. A generalized version of
-    # ToDo: this method should be implemented by the state vector simulator in Qiskit AER.
-    res = isometry_rs.apply_diagonal_gate(m, action_qubit_labels, diag)
-    return res
-
-
-# Special case of the method _apply_diagonal_gate, where the input m is a diagonal matrix on the
-# log2(len(m_diagonal)) least significant qubits (this method is more efficient in this case
-# than _apply_diagonal_gate). The input m_diagonal is provided as a list of diagonal entries.
-# The diagonal diag is applied on the qubits with labels listed in action_qubit_labels. The input
-# num_qubits gives the total number of considered qubits (this input is required to interpret the
-# action_qubit_labels in relation to the least significant qubits).
-
-
-def _apply_diagonal_gate_to_diag(m_diagonal, action_qubit_labels, diag, num_qubits):
-    res = isometry_rs.apply_diagonal_gate_to_diag(m_diagonal, action_qubit_labels, diag, num_qubits)
-    return res
-
-
-# Apply a MC single-qubit gate (given by the 2*2 unitary input: gate) with controlling on
-# the qubits with label control_labels and acting on the qubit with label target_label
-# to a matrix m. The input matrix m and the gate have to be of dtype=complex. The qubit labels are
-# such that label 0 corresponds to the most significant qubit, label 1 to the second most
-# significant qubit, and so on ...
-
-
-def _apply_multi_controlled_gate(m, control_labels, target_label, gate):
-    # ToDo: This method should be integrated into the state vector simulator in Qiskit AER
-    res = isometry_rs.apply_multi_controlled_gate(m, control_labels, target_label, gate)
-    return res
-
-
-# Some helper methods:
-
-
 # Get the qubits in the list qubits corresponding to the labels listed in labels. The total number
 # of qubits is given by num_qubits (and determines the convention for the qubit labeling)
 
@@ -431,45 +368,3 @@ def _get_binary_rep_as_list(n, num_digits):
         for c in line:
             binary.append(int(c))
     return binary[-num_digits:]
-
-
-# absorb a diagonal gate into a UCGate
-
-
-def _merge_UCGate_and_diag(single_qubit_gates, diag):
-    res = isometry_rs.merge_ucgate_and_diag(single_qubit_gates, diag)
-    return res
-
-
-# Helper variables/functions for the column-by-column decomposition
-
-
-# a(k,s) and b(k,s) are positive integers such that k = a(k,s)2^s + b(k,s)
-# (with the maximal choice of a(k,s))
-
-
-def _a(k, s):
-    return isometry_rs.a(k, s)
-
-
-def _b(k, s):
-    return isometry_rs.b(k, s)
-
-
-# given a binary representation of k with binary digits [k_{n-1},..,k_1,k_0],
-# the method k_s(k, s) returns k_s
-
-
-def _k_s(k, s):
-    return isometry_rs.k_s(k, s)
-
-
-# Check if a gate of a special form is equal to the identity gate up to global phase
-
-
-def _ucg_is_identity_up_to_global_phase(single_qubit_gates, epsilon):
-    return isometry_rs.ucg_is_identity_up_to_global_phase(single_qubit_gates, epsilon)
-
-
-def _diag_is_identity_up_to_global_phase(diag, epsilon):
-    return isometry_rs.diag_is_identity_up_to_global_phase(diag, epsilon)

From 3676469a1000c5e980ae4013c97981f3c6bbf0bf Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 06:48:21 -0400
Subject: [PATCH 11/20] Fix integer typing in uc_gate.rs

---
 crates/accelerate/src/uc_gate.rs | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/crates/accelerate/src/uc_gate.rs b/crates/accelerate/src/uc_gate.rs
index 111ae97c68b4..2639c012276a 100644
--- a/crates/accelerate/src/uc_gate.rs
+++ b/crates/accelerate/src/uc_gate.rs
@@ -88,15 +88,15 @@ pub fn dec_ucg_help(
     let mut diag: Array1<Complex64> = Array1::ones(2_usize.pow(num_qubits));
     let num_controls = num_qubits - 1;
     for dec_step in 0..num_controls {
-        let num_ucgs = 2_u32.pow(dec_step);
+        let num_ucgs = 2_usize.pow(dec_step);
         // The decomposition works recursively and the followign loop goes over the different
         // UCGates that arise in the decomposition
         for ucg_index in 0..num_ucgs {
-            let len_ucg = 2_u32.pow(num_controls - dec_step);
+            let len_ucg = 2_usize.pow(num_controls - dec_step);
             for i in 0..len_ucg / 2 {
                 let shift = ucg_index * len_ucg;
-                let a = single_qubit_gates[(shift + i) as usize].view();
-                let b = single_qubit_gates[(shift + len_ucg / 2 + i) as usize].view();
+                let a = single_qubit_gates[shift + i].view();
+                let b = single_qubit_gates[shift + len_ucg / 2 + i].view();
                 // Apply the decomposition for UCGates given in equation (3) in
                 // https://arxiv.org/pdf/quant-ph/0410066.pdf
                 // to demultiplex one control of all the num_ucgs uniformly-controlled gates
@@ -104,8 +104,8 @@ pub fn dec_ucg_help(
                 let [v, u, r] = demultiplex_single_uc(a, b);
                 // replace the single-qubit gates with v,u (the already existing ones
                 // are not needed any more)
-                single_qubit_gates[(shift + i) as usize] = v;
-                single_qubit_gates[(shift + len_ucg / 2 + i) as usize] = u;
+                single_qubit_gates[shift + i] = v;
+                single_qubit_gates[shift + len_ucg / 2 + i] = u;
                 // Now we decompose the gates D as described in Figure 4 in
                 // https://arxiv.org/pdf/quant-ph/0410066.pdf and merge some of the gates
                 // into the UCGates and the diagonal at the end of the circuit
@@ -121,11 +121,11 @@ pub fn dec_ucg_help(
                     // Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and
                     // merge the UC-Rz rotation with the following UCGate,
                     // which hasn't been decomposed yet
-                    let k = (shift + len_ucg + i) as usize;
+                    let k = shift + len_ucg + i;
 
                     single_qubit_gates[k] = single_qubit_gates[k].dot(&r_conj_t);
                     single_qubit_gates[k].mapv_inplace(|x| x * rz_00);
-                    let k = k + len_ucg as usize / 2;
+                    let k = k + len_ucg / 2;
                     single_qubit_gates[k] = single_qubit_gates[k].dot(&r);
                     single_qubit_gates[k].mapv_inplace(|x| x * rz_11);
                 } else {
@@ -133,10 +133,10 @@ pub fn dec_ucg_help(
                     // the trailing UC-Rz rotation into a diagonal gate at the end of the circuit
                     for ucg_index_2 in 0..num_ucgs {
                         let shift_2 = ucg_index_2 * len_ucg;
-                        let k = (2 * (i + shift_2)) as usize;
+                        let k = 2 * (i + shift_2);
                         diag[k] *= r_conj_t[[0, 0]] * rz_00;
                         diag[k + 1] *= r_conj_t[[1, 1]] * rz_00;
-                        let k = len_ucg as usize + k;
+                        let k = len_ucg + k;
                         diag[k] *= r[[0, 0]] * rz_11;
                         diag[k + 1] *= r[[1, 1]] * rz_11;
                     }

From 29ffc994fd51c2f5d4fa8162bfa3a9966bbc70e5 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 06:54:21 -0400
Subject: [PATCH 12/20] Simplify basis state bitshift loop logic

---
 crates/accelerate/src/isometry.rs | 8 ++------
 1 file changed, 2 insertions(+), 6 deletions(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index f3821520c484..877ad7884828 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -203,20 +203,16 @@ fn construct_basis_states(
     let mut j = 0;
     let control_set: HashSet<usize> = control_labels.iter().copied().collect();
     for i in 0..size {
-        if control_set.contains(&i) {
             e1 <<= 1;
-            e1 += 1;
             e2 <<= 1;
+        if control_set.contains(&i) {
+            e1 += 1;
             e2 += 1;
         } else if i == target_label {
-            e1 <<= 1;
-            e2 <<= 1;
             e2 += 1;
         } else {
             assert!(j <= 1);
-            e1 <<= 1;
             e1 += state_free[j] as usize;
-            e2 <<= 1;
             e2 += state_free[j] as usize;
             j += 1
         }

From 366a934f9a33c4528bb375dedbf6749be2409b70 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 06:57:00 -0400
Subject: [PATCH 13/20] Build set of control labels outside
 construct_basis_states

---
 crates/accelerate/src/isometry.rs | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 877ad7884828..976219de8dda 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -22,7 +22,6 @@ use hashbrown::HashSet;
 use itertools::Itertools;
 use ndarray::prelude::*;
 use numpy::{IntoPyArray, PyReadonlyArray2};
-use rayon::prelude::*;
 
 use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY;
 
@@ -194,17 +193,16 @@ pub fn apply_diagonal_gate_to_diag(
 /// is acting on for a specific state state_free of the qubits we neither control nor act on
 fn construct_basis_states(
     state_free: &[u8],
-    control_labels: &[usize],
+    control_set: &HashSet<usize>,
     target_label: usize,
 ) -> [usize; 2] {
-    let size = state_free.len() + control_labels.len() + 1;
+    let size = state_free.len() + control_set.len() + 1;
     let mut e1: usize = 0;
     let mut e2: usize = 0;
     let mut j = 0;
-    let control_set: HashSet<usize> = control_labels.iter().copied().collect();
     for i in 0..size {
-            e1 <<= 1;
-            e2 <<= 1;
+        e1 <<= 1;
+        e2 <<= 1;
         if control_set.contains(&i) {
             e1 += 1;
             e2 += 1;
@@ -224,7 +222,7 @@ fn construct_basis_states(
 pub fn apply_multi_controlled_gate(
     py: Python,
     m: PyReadonlyArray2<Complex64>,
-    mut control_labels: Vec<usize>,
+    control_labels: Vec<usize>,
     target_label: usize,
     gate: PyReadonlyArray2<Complex64>,
 ) -> PyObject {
@@ -233,10 +231,10 @@ pub fn apply_multi_controlled_gate(
     let shape = m.shape();
     let num_qubits = shape[0].ilog2();
     let num_col = shape[1];
-    control_labels.par_sort();
     let free_qubits = num_qubits as usize - control_labels.len() - 1;
+    let control_set: HashSet<usize> = control_labels.into_iter().collect();
     if free_qubits == 0 {
-        let [e1, e2] = construct_basis_states(&[], &control_labels, target_label);
+        let [e1, e2] = construct_basis_states(&[], &control_set, target_label);
         for i in 0..num_col {
             let temp: Vec<_> = gate
                 .dot(&aview2(&[[m[[e1, i]]], [m[[e2, i]]]]))
@@ -252,7 +250,7 @@ pub fn apply_multi_controlled_gate(
         .take(free_qubits)
         .multi_cartesian_product()
     {
-        let [e1, e2] = construct_basis_states(&state_free, &control_labels, target_label);
+        let [e1, e2] = construct_basis_states(&state_free, &control_set, target_label);
         for i in 0..num_col {
             let temp: Vec<_> = gate
                 .dot(&aview2(&[[m[[e1, i]]], [m[[e2, i]]]]))

From 4843ef5c41652afa6488d49d60cd75cfbd5a52e1 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 07:03:59 -0400
Subject: [PATCH 14/20] Use 2 element array for reverse_qubit_state

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

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 976219de8dda..2eb45cc36542 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -30,7 +30,7 @@ use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY;
 #[pyfunction]
 pub fn reverse_qubit_state(
     py: Python,
-    state: Vec<Complex64>,
+    state: [Complex64; 2],
     basis_state: usize,
     epsilon: f64,
 ) -> PyObject {
@@ -47,7 +47,7 @@ fn l2_norm(vec: &[Complex64]) -> f64 {
 }
 
 fn reverse_qubit_state_inner(
-    state: &[Complex64],
+    state: &[Complex64; 2],
     basis_state: usize,
     epsilon: f64,
 ) -> Array2<Complex64> {

From d26bd4aadbac6001317ed8ecbb3e16c648cd2483 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 07:04:13 -0400
Subject: [PATCH 15/20] Micro optimize reverse_qubit_state

---
 crates/accelerate/src/isometry.rs | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 2eb45cc36542..ab7dbeebcc4c 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -52,17 +52,18 @@ fn reverse_qubit_state_inner(
     epsilon: f64,
 ) -> Array2<Complex64> {
     let r = l2_norm(state);
+    let r_inv = 1. / r;
     if r < epsilon {
         Array2::eye(2)
     } else if basis_state == 0 {
         array![
-            [state[0].conj() / r, state[1].conj() / r],
-            [-state[1] / r, state[0] / r],
+            [state[0].conj() * r_inv, state[1].conj() * r_inv],
+            [-state[1] * r_inv, state[0] * r_inv],
         ]
     } else {
         array![
-            [-state[1] / r, state[0] / r],
-            [state[0].conj() / r, state[1].conj() / r],
+            [-state[1] * r_inv, state[0] * r_inv],
+            [state[0].conj() * r_inv, state[1].conj() * r_inv],
         ]
     }
 }

From 091d229957271de68904f5b3141364915585f8d1 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 07:23:18 -0400
Subject: [PATCH 16/20] Use 1d numpy arrays for diagonal inputs

---
 crates/accelerate/src/isometry.rs              | 18 ++++++++++--------
 .../library/generalized_gates/isometry.py      |  4 ++--
 2 files changed, 12 insertions(+), 10 deletions(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index ab7dbeebcc4c..c4f40bc8db6a 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -21,7 +21,7 @@ use pyo3::Python;
 use hashbrown::HashSet;
 use itertools::Itertools;
 use ndarray::prelude::*;
-use numpy::{IntoPyArray, PyReadonlyArray2};
+use numpy::{IntoPyArray, PyReadonlyArray2, PyReadonlyArray1};
 
 use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY;
 
@@ -145,8 +145,9 @@ pub fn apply_diagonal_gate(
     py: Python,
     m: PyReadonlyArray2<Complex64>,
     action_qubit_labels: Vec<usize>,
-    diag: Vec<Complex64>,
-) -> PyObject {
+    diag: PyReadonlyArray1<Complex64>,
+) -> PyResult<PyObject> {
+    let diag = diag.as_slice()?;
     let mut m = m.as_array().to_owned();
     let shape = m.shape();
     let num_qubits = shape[0].ilog2();
@@ -163,18 +164,19 @@ pub fn apply_diagonal_gate(
             m[[i, j]] = diag[diag_index] * m[[i, j]]
         }
     }
-    m.into_pyarray_bound(py).into()
+    Ok(m.into_pyarray_bound(py).into())
 }
 
 #[pyfunction]
 pub fn apply_diagonal_gate_to_diag(
     mut m_diagonal: Vec<Complex64>,
     action_qubit_labels: Vec<usize>,
-    diag: Vec<Complex64>,
+    diag: PyReadonlyArray1<Complex64>,
     num_qubits: usize,
-) -> Vec<Complex64> {
+) -> PyResult<Vec<Complex64>> {
+    let diag = diag.as_slice()?;
     if m_diagonal.is_empty() {
-        return m_diagonal;
+        return Ok(m_diagonal);
     }
     for state in std::iter::repeat([0_u8, 1_u8])
         .take(num_qubits)
@@ -187,7 +189,7 @@ pub fn apply_diagonal_gate_to_diag(
         let i = bin_to_int(&state);
         m_diagonal[i] *= diag[diag_index]
     }
-    m_diagonal
+    Ok(m_diagonal)
 }
 
 /// Helper method for _apply_multi_controlled_gate. This constructs the basis states the MG gate
diff --git a/qiskit/circuit/library/generalized_gates/isometry.py b/qiskit/circuit/library/generalized_gates/isometry.py
index 2c597b7507e0..c180e7a14484 100644
--- a/qiskit/circuit/library/generalized_gates/isometry.py
+++ b/qiskit/circuit/library/generalized_gates/isometry.py
@@ -220,7 +220,7 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
             # apply the MCG to the remaining isometry
             v = isometry_rs.apply_multi_controlled_gate(v, control_labels, target_label, gate)
             # correct for the implementation "up to diagonal"
-            diag_mcg_inverse = np.conj(diagonal_mcg).tolist()
+            diag_mcg_inverse = np.conj(diagonal_mcg).astype(complex, copy=False)
             v = isometry_rs.apply_diagonal_gate(
                 v, control_labels + [target_label], diag_mcg_inverse
             )
@@ -238,7 +238,7 @@ def _disentangle(self, circuit, q, diag, remaining_isometry, column_index, s):
                 circuit, q, single_qubit_gates, control_labels, target_label
             )
             # merge the diagonal into the UCGate for efficient application of both together
-            diagonal_ucg_inverse = np.conj(diagonal_ucg).tolist()
+            diagonal_ucg_inverse = np.conj(diagonal_ucg).astype(complex, copy=False)
             single_qubit_gates = isometry_rs.merge_ucgate_and_diag(
                 single_qubit_gates, diagonal_ucg_inverse
             )

From 604e33aafe7c9fd4ea193bd90b527d2a3b94a141 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Fri, 26 Apr 2024 07:43:09 -0400
Subject: [PATCH 17/20] Fix lint

---
 crates/accelerate/src/isometry.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index c4f40bc8db6a..472874e9e9c1 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -21,7 +21,7 @@ use pyo3::Python;
 use hashbrown::HashSet;
 use itertools::Itertools;
 use ndarray::prelude::*;
-use numpy::{IntoPyArray, PyReadonlyArray2, PyReadonlyArray1};
+use numpy::{IntoPyArray, PyReadonlyArray1, PyReadonlyArray2};
 
 use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY;
 

From 3bd43586629c100f180915fb32e7217db5942c8b Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Mon, 29 Apr 2024 13:48:49 -0400
Subject: [PATCH 18/20] Update crates/accelerate/src/isometry.rs

Co-authored-by: John Lapeyre <jlapeyre@users.noreply.github.com>
---
 crates/accelerate/src/isometry.rs | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 472874e9e9c1..906c9c1c870e 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -42,8 +42,7 @@ pub fn reverse_qubit_state(
 #[inline(always)]
 fn l2_norm(vec: &[Complex64]) -> f64 {
     vec.iter()
-        .fold(0., |acc, elem| acc + elem.abs().powi(2))
-        .sqrt()
+        .fold(0., |acc, elem| acc + elem.norm_sqr())
 }
 
 fn reverse_qubit_state_inner(

From b8b82a6bf187cc1fa5d07009b2a5fdbff7811027 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Mon, 29 Apr 2024 15:15:49 -0400
Subject: [PATCH 19/20] Add back sqrt() accidentally removed by inline
 suggestion

---
 crates/accelerate/src/isometry.rs | 1 +
 1 file changed, 1 insertion(+)

diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index 906c9c1c870e..8d0761666bb6 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -43,6 +43,7 @@ pub fn reverse_qubit_state(
 fn l2_norm(vec: &[Complex64]) -> f64 {
     vec.iter()
         .fold(0., |acc, elem| acc + elem.norm_sqr())
+        .sqrt()
 }
 
 fn reverse_qubit_state_inner(

From fcf587f661fe83970a1313e712aa6ccb189aac72 Mon Sep 17 00:00:00 2001
From: Matthew Treinish <mtreinish@kortar.org>
Date: Tue, 30 Apr 2024 15:49:28 -0400
Subject: [PATCH 20/20] Use a constant for rz pi/2 elements

---
 crates/accelerate/src/uc_gate.rs | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/crates/accelerate/src/uc_gate.rs b/crates/accelerate/src/uc_gate.rs
index 2639c012276a..3a5f74a6f0b1 100644
--- a/crates/accelerate/src/uc_gate.rs
+++ b/crates/accelerate/src/uc_gate.rs
@@ -14,7 +14,7 @@ use num_complex::{Complex64, ComplexFloat};
 use pyo3::prelude::*;
 use pyo3::wrap_pyfunction;
 use pyo3::Python;
-use std::f64::consts::PI;
+use std::f64::consts::{FRAC_1_SQRT_2, PI};
 
 use faer_ext::{IntoFaerComplex, IntoNdarrayComplex};
 use ndarray::prelude::*;
@@ -25,6 +25,11 @@ use crate::euler_one_qubit_decomposer::det_one_qubit;
 const PI2: f64 = PI / 2.;
 const EPS: f64 = 1e-10;
 
+// These constants are the non-zero elements of an RZ gate's unitary with an
+// angle of pi / 2
+const RZ_PI2_11: Complex64 = Complex64::new(FRAC_1_SQRT_2, -FRAC_1_SQRT_2);
+const RZ_PI2_00: Complex64 = Complex64::new(FRAC_1_SQRT_2, FRAC_1_SQRT_2);
+
 /// This method implements the decomposition given in equation (3) in
 /// https://arxiv.org/pdf/quant-ph/0410066.pdf.
 ///
@@ -114,8 +119,6 @@ pub fn dec_ucg_help(
                 // gates arising in the decomposition of D are ignored for the moment (they will
                 // be added together with the C-NOT gates at the end of the decomposition
                 // (in the method dec_ucg()))
-                let rz_11 = (-Complex64::new(0., 0.5 * PI2)).exp();
-                let rz_00 = Complex64::new(0., 0.5 * PI2).exp();
                 let r_conj_t = r.mapv(|x| x.conj()).t().to_owned();
                 if ucg_index < num_ucgs - 1 {
                     // Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and
@@ -124,21 +127,21 @@ pub fn dec_ucg_help(
                     let k = shift + len_ucg + i;
 
                     single_qubit_gates[k] = single_qubit_gates[k].dot(&r_conj_t);
-                    single_qubit_gates[k].mapv_inplace(|x| x * rz_00);
+                    single_qubit_gates[k].mapv_inplace(|x| x * RZ_PI2_00);
                     let k = k + len_ucg / 2;
                     single_qubit_gates[k] = single_qubit_gates[k].dot(&r);
-                    single_qubit_gates[k].mapv_inplace(|x| x * rz_11);
+                    single_qubit_gates[k].mapv_inplace(|x| x * RZ_PI2_11);
                 } else {
                     // Absorb the Rz(pi/2) rotation on the control into the UC-Rz gate and merge
                     // the trailing UC-Rz rotation into a diagonal gate at the end of the circuit
                     for ucg_index_2 in 0..num_ucgs {
                         let shift_2 = ucg_index_2 * len_ucg;
                         let k = 2 * (i + shift_2);
-                        diag[k] *= r_conj_t[[0, 0]] * rz_00;
-                        diag[k + 1] *= r_conj_t[[1, 1]] * rz_00;
+                        diag[k] *= r_conj_t[[0, 0]] * RZ_PI2_00;
+                        diag[k + 1] *= r_conj_t[[1, 1]] * RZ_PI2_00;
                         let k = len_ucg + k;
-                        diag[k] *= r[[0, 0]] * rz_11;
-                        diag[k + 1] *= r[[1, 1]] * rz_11;
+                        diag[k] *= r[[0, 0]] * RZ_PI2_11;
+                        diag[k + 1] *= r[[1, 1]] * RZ_PI2_11;
                     }
                 }
             }