diff --git a/crates/accelerate/src/convert_2q_block_matrix.rs b/crates/accelerate/src/convert_2q_block_matrix.rs
index e311c129b11b..9c179397d641 100644
--- a/crates/accelerate/src/convert_2q_block_matrix.rs
+++ b/crates/accelerate/src/convert_2q_block_matrix.rs
@@ -20,10 +20,7 @@ use numpy::ndarray::{aview2, Array2, ArrayView2};
 use numpy::{IntoPyArray, PyArray2, PyReadonlyArray2};
 use smallvec::SmallVec;
 
-static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] = [
-    [Complex64::new(1., 0.), Complex64::new(0., 0.)],
-    [Complex64::new(0., 0.), Complex64::new(1., 0.)],
-];
+use qiskit_circuit::gate_matrix::ONE_QUBIT_IDENTITY;
 
 /// Return the matrix Operator resulting from a block of Instructions.
 #[pyfunction]
diff --git a/crates/accelerate/src/euler_one_qubit_decomposer.rs b/crates/accelerate/src/euler_one_qubit_decomposer.rs
index 1fd5fd7834ff..9f10f76de467 100644
--- a/crates/accelerate/src/euler_one_qubit_decomposer.rs
+++ b/crates/accelerate/src/euler_one_qubit_decomposer.rs
@@ -31,6 +31,7 @@ use ndarray::prelude::*;
 use numpy::PyReadonlyArray2;
 use pyo3::pybacked::PyBackedStr;
 
+use qiskit_circuit::util::c64;
 use qiskit_circuit::SliceOrInt;
 
 pub const ANGLE_ZERO_EPSILON: f64 = 1e-12;
@@ -855,16 +856,16 @@ pub fn params_xyx(unitary: PyReadonlyArray2<Complex64>) -> [f64; 4] {
 
 fn params_xzx_inner(umat: ArrayView2<Complex64>) -> [f64; 4] {
     let det = det_one_qubit(umat);
-    let phase = (Complex64::new(0., -1.) * det.ln()).re / 2.;
+    let phase = det.ln().im / 2.;
     let sqrt_det = det.sqrt();
     let mat_zyz = arr2(&[
         [
-            Complex64::new((umat[[0, 0]] / sqrt_det).re, (umat[[1, 0]] / sqrt_det).im),
-            Complex64::new((umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im),
+            c64((umat[[0, 0]] / sqrt_det).re, (umat[[1, 0]] / sqrt_det).im),
+            c64((umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im),
         ],
         [
-            Complex64::new(-(umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im),
-            Complex64::new((umat[[0, 0]] / sqrt_det).re, -(umat[[1, 0]] / sqrt_det).im),
+            c64(-(umat[[1, 0]] / sqrt_det).re, (umat[[0, 0]] / sqrt_det).im),
+            c64((umat[[0, 0]] / sqrt_det).re, -(umat[[1, 0]] / sqrt_det).im),
         ],
     ]);
     let [theta, phi, lam, phase_zxz] = params_zxz_inner(mat_zyz.view());
diff --git a/crates/accelerate/src/isometry.rs b/crates/accelerate/src/isometry.rs
index a3a8be38dae2..ceaba2946b3a 100644
--- a/crates/accelerate/src/isometry.rs
+++ b/crates/accelerate/src/isometry.rs
@@ -24,6 +24,7 @@ use ndarray::prelude::*;
 use numpy::{IntoPyArray, PyReadonlyArray1, PyReadonlyArray2};
 
 use qiskit_circuit::gate_matrix::ONE_QUBIT_IDENTITY;
+use qiskit_circuit::util::C_ZERO;
 
 /// Find special unitary matrix that maps [c0,c1] to [r,0] or [0,r] if basis_state=0 or
 /// basis_state=1 respectively
@@ -315,11 +316,7 @@ pub fn merge_ucgate_and_diag(
         .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);
+            let res = aview2(&[[diag[2 * i], C_ZERO], [C_ZERO, diag[2 * i + 1]]]).dot(&gate);
             res.into_pyarray_bound(py).into()
         })
         .collect()
diff --git a/crates/accelerate/src/pauli_exp_val.rs b/crates/accelerate/src/pauli_exp_val.rs
index 52a2fc07f81d..8ee4b019b3e0 100644
--- a/crates/accelerate/src/pauli_exp_val.rs
+++ b/crates/accelerate/src/pauli_exp_val.rs
@@ -19,6 +19,7 @@ use pyo3::wrap_pyfunction;
 use rayon::prelude::*;
 
 use crate::getenv_use_multiple_threads;
+use qiskit_circuit::util::c64;
 
 const PARALLEL_THRESHOLD: usize = 19;
 
@@ -88,7 +89,7 @@ pub fn expval_pauli_with_x(
         let index_0 = ((i << 1) & mask_u) | (i & mask_l);
         let index_1 = index_0 ^ x_mask;
         let val_0 = (phase
-            * Complex64::new(
+            * c64(
                 data_arr[index_1].re * data_arr[index_0].re
                     + data_arr[index_1].im * data_arr[index_0].im,
                 data_arr[index_1].im * data_arr[index_0].re
@@ -96,7 +97,7 @@ pub fn expval_pauli_with_x(
             ))
         .re;
         let val_1 = (phase
-            * Complex64::new(
+            * c64(
                 data_arr[index_0].re * data_arr[index_1].re
                     + data_arr[index_0].im * data_arr[index_1].im,
                 data_arr[index_0].im * data_arr[index_1].re
diff --git a/crates/accelerate/src/sampled_exp_val.rs b/crates/accelerate/src/sampled_exp_val.rs
index b51ca3c98f0e..0b8836a94165 100644
--- a/crates/accelerate/src/sampled_exp_val.rs
+++ b/crates/accelerate/src/sampled_exp_val.rs
@@ -18,6 +18,7 @@ use pyo3::prelude::*;
 use pyo3::wrap_pyfunction;
 
 use crate::pauli_exp_val::fast_sum;
+use qiskit_circuit::util::c64;
 
 const OPER_TABLE_SIZE: usize = (b'Z' as usize) + 1;
 const fn generate_oper_table() -> [[f64; 2]; OPER_TABLE_SIZE] {
@@ -81,7 +82,7 @@ pub fn sampled_expval_complex(
     let out: Complex64 = oper_strs
         .into_iter()
         .enumerate()
-        .map(|(idx, string)| coeff_arr[idx] * Complex64::new(bitstring_expval(&dist, string), 0.))
+        .map(|(idx, string)| coeff_arr[idx] * c64(bitstring_expval(&dist, string), 0.))
         .sum();
     Ok(out.re)
 }
diff --git a/crates/accelerate/src/sparse_pauli_op.rs b/crates/accelerate/src/sparse_pauli_op.rs
index e0c80f716161..8a51d8ee781c 100644
--- a/crates/accelerate/src/sparse_pauli_op.rs
+++ b/crates/accelerate/src/sparse_pauli_op.rs
@@ -23,6 +23,7 @@ use hashbrown::HashMap;
 use ndarray::{s, Array1, Array2, ArrayView1, ArrayView2, Axis};
 use num_complex::Complex64;
 use num_traits::Zero;
+use qiskit_circuit::util::{c64, C_ONE, C_ZERO};
 use rayon::prelude::*;
 
 use crate::rayon_ext::*;
@@ -257,9 +258,9 @@ impl<'py> ZXPaulisView<'py> {
                 let ys = (xs & zs).count_ones();
                 match (phase as u32 + ys) % 4 {
                     0 => coeff,
-                    1 => Complex64::new(coeff.im, -coeff.re),
-                    2 => Complex64::new(-coeff.re, -coeff.im),
-                    3 => Complex64::new(-coeff.im, coeff.re),
+                    1 => c64(coeff.im, -coeff.re),
+                    2 => c64(-coeff.re, -coeff.im),
+                    3 => c64(-coeff.im, coeff.re),
                     _ => unreachable!(),
                 }
             })
@@ -311,10 +312,10 @@ impl MatrixCompressedPaulis {
             .zip(self.z_like.drain(..))
             .zip(self.coeffs.drain(..))
         {
-            *hash_table.entry(key).or_insert(Complex64::new(0.0, 0.0)) += coeff;
+            *hash_table.entry(key).or_insert(C_ZERO) += coeff;
         }
         for ((x, z), coeff) in hash_table {
-            if coeff == Complex64::new(0.0, 0.0) {
+            if coeff.is_zero() {
                 continue;
             }
             self.x_like.push(x);
@@ -347,7 +348,7 @@ pub fn decompose_dense(
     let mut coeffs = vec![];
     if num_qubits > 0 {
         decompose_dense_inner(
-            Complex64::new(1.0, 0.0),
+            C_ONE,
             num_qubits,
             &[],
             operator.as_array(),
@@ -532,7 +533,7 @@ fn to_matrix_dense_inner(paulis: &MatrixCompressedPaulis, parallel: bool) -> Vec
         // Doing the initialization here means that when we're in parallel contexts, we do the
         // zeroing across the whole threadpool.  This also seems to give a speed-up in serial
         // contexts, but I don't understand that. ---Jake
-        row.fill(Complex64::new(0.0, 0.0));
+        row.fill(C_ZERO);
         for ((&x_like, &z_like), &coeff) in paulis
             .x_like
             .iter()
@@ -667,7 +668,7 @@ macro_rules! impl_to_matrix_sparse {
                     ((i_row as $uint_ty) ^ (paulis.x_like[a] as $uint_ty))
                         .cmp(&((i_row as $uint_ty) ^ (paulis.x_like[b] as $uint_ty)))
                 });
-                let mut running = Complex64::new(0.0, 0.0);
+                let mut running = C_ZERO;
                 let mut prev_index = i_row ^ (paulis.x_like[order[0]] as usize);
                 for (x_like, z_like, coeff) in order
                     .iter()
@@ -748,7 +749,7 @@ macro_rules! impl_to_matrix_sparse {
                             (i_row as $uint_ty ^ paulis.x_like[a] as $uint_ty)
                                 .cmp(&(i_row as $uint_ty ^ paulis.x_like[b] as $uint_ty))
                         });
-                        let mut running = Complex64::new(0.0, 0.0);
+                        let mut running = C_ZERO;
                         let mut prev_index = i_row ^ (paulis.x_like[order[0]] as usize);
                         for (x_like, z_like, coeff) in order
                             .iter()
@@ -844,11 +845,11 @@ mod tests {
             // Deliberately using multiples of small powers of two so the floating-point addition
             // of them is associative.
             coeffs: vec![
-                Complex64::new(0.25, 0.5),
-                Complex64::new(0.125, 0.25),
-                Complex64::new(0.375, 0.125),
-                Complex64::new(-0.375, 0.0625),
-                Complex64::new(-0.5, -0.25),
+                c64(0.25, 0.5),
+                c64(0.125, 0.25),
+                c64(0.375, 0.125),
+                c64(-0.375, 0.0625),
+                c64(-0.5, -0.25),
             ],
         }
     }
diff --git a/crates/accelerate/src/two_qubit_decompose.rs b/crates/accelerate/src/two_qubit_decompose.rs
index e8c572b04039..8637cb03c735 100644
--- a/crates/accelerate/src/two_qubit_decompose.rs
+++ b/crates/accelerate/src/two_qubit_decompose.rs
@@ -52,67 +52,28 @@ use rand_distr::StandardNormal;
 use rand_pcg::Pcg64Mcg;
 
 use qiskit_circuit::gate_matrix::{CX_GATE, H_GATE, ONE_QUBIT_IDENTITY, SX_GATE, X_GATE};
+use qiskit_circuit::util::{c64, GateArray1Q, GateArray2Q, C_M_ONE, C_ONE, C_ZERO, IM, M_IM};
 use qiskit_circuit::SliceOrInt;
 
-const PI2: f64 = PI / 2.0;
-const PI4: f64 = PI / 4.0;
+const PI2: f64 = PI / 2.;
+const PI4: f64 = PI / 4.;
 const PI32: f64 = 3.0 * PI2;
 const TWO_PI: f64 = 2.0 * PI;
 
 const C1: c64 = c64 { re: 1.0, im: 0.0 };
 
-static B_NON_NORMALIZED: [[Complex64; 4]; 4] = [
-    [
-        Complex64::new(1.0, 0.),
-        Complex64::new(0., 1.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-    ],
-    [
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 1.),
-        Complex64::new(1.0, 0.0),
-    ],
-    [
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 1.),
-        Complex64::new(-1., 0.),
-    ],
-    [
-        Complex64::new(1., 0.),
-        Complex64::new(0., -1.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-    ],
+static B_NON_NORMALIZED: GateArray2Q = [
+    [C_ONE, IM, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, IM, C_ONE],
+    [C_ZERO, C_ZERO, IM, C_M_ONE],
+    [C_ONE, M_IM, C_ZERO, C_ZERO],
 ];
 
-static B_NON_NORMALIZED_DAGGER: [[Complex64; 4]; 4] = [
-    [
-        Complex64::new(0.5, 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0.5, 0.0),
-    ],
-    [
-        Complex64::new(0., -0.5),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.),
-        Complex64::new(0., 0.5),
-    ],
-    [
-        Complex64::new(0., 0.),
-        Complex64::new(0., -0.5),
-        Complex64::new(0., -0.5),
-        Complex64::new(0., 0.),
-    ],
-    [
-        Complex64::new(0., 0.),
-        Complex64::new(0.5, 0.),
-        Complex64::new(-0.5, 0.),
-        Complex64::new(0., 0.),
-    ],
+static B_NON_NORMALIZED_DAGGER: GateArray2Q = [
+    [c64(0.5, 0.), C_ZERO, C_ZERO, c64(0.5, 0.)],
+    [c64(0., -0.5), C_ZERO, C_ZERO, c64(0., 0.5)],
+    [C_ZERO, c64(0., -0.5), c64(0., -0.5), C_ZERO],
+    [C_ZERO, c64(0.5, 0.), c64(-0.5, 0.), C_ZERO],
 ];
 
 enum MagicBasisTransform {
@@ -318,29 +279,26 @@ fn closest_partial_swap(a: f64, b: f64, c: f64) -> f64 {
 
 fn rx_matrix(theta: f64) -> Array2<Complex64> {
     let half_theta = theta / 2.;
-    let cos = Complex64::new(half_theta.cos(), 0.);
-    let isin = Complex64::new(0., -half_theta.sin());
+    let cos = c64(half_theta.cos(), 0.);
+    let isin = c64(0., -half_theta.sin());
     array![[cos, isin], [isin, cos]]
 }
 
 fn ry_matrix(theta: f64) -> Array2<Complex64> {
     let half_theta = theta / 2.;
-    let cos = Complex64::new(half_theta.cos(), 0.);
-    let sin = Complex64::new(half_theta.sin(), 0.);
+    let cos = c64(half_theta.cos(), 0.);
+    let sin = c64(half_theta.sin(), 0.);
     array![[cos, -sin], [sin, cos]]
 }
 
 fn rz_matrix(theta: f64) -> Array2<Complex64> {
-    let ilam2 = Complex64::new(0., 0.5 * theta);
-    array![
-        [(-ilam2).exp(), Complex64::new(0., 0.)],
-        [Complex64::new(0., 0.), ilam2.exp()]
-    ]
+    let ilam2 = c64(0., 0.5 * theta);
+    array![[(-ilam2).exp(), C_ZERO], [C_ZERO, ilam2.exp()]]
 }
 
 fn compute_unitary(sequence: &TwoQubitSequenceVec, global_phase: f64) -> Array2<Complex64> {
     let identity = aview2(&ONE_QUBIT_IDENTITY);
-    let phase = Complex64::new(0., global_phase).exp();
+    let phase = c64(0., global_phase).exp();
     let mut matrix = Array2::from_diag(&arr1(&[phase, phase, phase, phase]));
     sequence
         .iter()
@@ -375,7 +333,6 @@ fn compute_unitary(sequence: &TwoQubitSequenceVec, global_phase: f64) -> Array2<
 }
 
 const DEFAULT_FIDELITY: f64 = 1.0 - 1.0e-9;
-const C1_IM: Complex64 = Complex64::new(0.0, 1.0);
 
 #[derive(Clone, Debug, Copy)]
 #[pyclass(module = "qiskit._accelerate.two_qubit_decompose")]
@@ -500,18 +457,9 @@ impl TwoQubitWeylDecomposition {
     }
 }
 
-static IPZ: [[Complex64; 2]; 2] = [
-    [C1_IM, Complex64::new(0., 0.)],
-    [Complex64::new(0., 0.), Complex64::new(0., -1.)],
-];
-static IPY: [[Complex64; 2]; 2] = [
-    [Complex64::new(0., 0.), Complex64::new(1., 0.)],
-    [Complex64::new(-1., 0.), Complex64::new(0., 0.)],
-];
-static IPX: [[Complex64; 2]; 2] = [
-    [Complex64::new(0., 0.), C1_IM],
-    [C1_IM, Complex64::new(0., 0.)],
-];
+static IPZ: GateArray1Q = [[IM, C_ZERO], [C_ZERO, M_IM]];
+static IPY: GateArray1Q = [[C_ZERO, C_ONE], [C_M_ONE, C_ZERO]];
+static IPX: GateArray1Q = [[C_ZERO, IM], [IM, C_ZERO]];
 
 #[pymethods]
 impl TwoQubitWeylDecomposition {
@@ -671,7 +619,7 @@ impl TwoQubitWeylDecomposition {
         temp.diag_mut()
             .iter_mut()
             .enumerate()
-            .for_each(|(index, x)| *x = (C1_IM * d[index]).exp());
+            .for_each(|(index, x)| *x = (IM * d[index]).exp());
         let k1 = magic_basis_transform(u_p.dot(&p).dot(&temp).view(), MagicBasisTransform::Into);
         let k2 = magic_basis_transform(p.t(), MagicBasisTransform::Into);
 
@@ -737,7 +685,7 @@ impl TwoQubitWeylDecomposition {
         let is_close = |ap: f64, bp: f64, cp: f64| -> bool {
             let [da, db, dc] = [a - ap, b - bp, c - cp];
             let tr = 4.
-                * Complex64::new(
+                * c64(
                     da.cos() * db.cos() * dc.cos(),
                     da.sin() * db.sin() * dc.sin(),
                 );
@@ -1016,13 +964,13 @@ impl TwoQubitWeylDecomposition {
                 b - specialized.b,
                 -c - specialized.c,
             ];
-            4. * Complex64::new(
+            4. * c64(
                 da.cos() * db.cos() * dc.cos(),
                 da.sin() * db.sin() * dc.sin(),
             )
         } else {
             let [da, db, dc] = [a - specialized.a, b - specialized.b, c - specialized.c];
-            4. * Complex64::new(
+            4. * c64(
                 da.cos() * db.cos() * dc.cos(),
                 da.sin() * db.sin() * dc.sin(),
             )
@@ -1597,20 +1545,14 @@ impl TwoQubitBasisDecomposer {
     }
 }
 
-static K12R_ARR: [[Complex64; 2]; 2] = [
-    [
-        Complex64::new(0., FRAC_1_SQRT_2),
-        Complex64::new(FRAC_1_SQRT_2, 0.),
-    ],
-    [
-        Complex64::new(-FRAC_1_SQRT_2, 0.),
-        Complex64::new(0., -FRAC_1_SQRT_2),
-    ],
+static K12R_ARR: GateArray1Q = [
+    [c64(0., FRAC_1_SQRT_2), c64(FRAC_1_SQRT_2, 0.)],
+    [c64(-FRAC_1_SQRT_2, 0.), c64(0., -FRAC_1_SQRT_2)],
 ];
 
-static K12L_ARR: [[Complex64; 2]; 2] = [
-    [Complex64::new(0.5, 0.5), Complex64::new(0.5, 0.5)],
-    [Complex64::new(-0.5, 0.5), Complex64::new(0.5, -0.5)],
+static K12L_ARR: GateArray1Q = [
+    [c64(0.5, 0.5), c64(0.5, 0.5)],
+    [c64(-0.5, 0.5), c64(0.5, -0.5)],
 ];
 
 fn decomp0_inner(target: &TwoQubitWeylDecomposition) -> SmallVec<[Array2<Complex64>; 8]> {
@@ -1650,90 +1592,71 @@ impl TwoQubitBasisDecomposer {
         // Create some useful matrices U1, U2, U3 are equivalent to the basis,
         // expand as Ui = Ki1.Ubasis.Ki2
         let b = basis_decomposer.b;
-        let temp = Complex64::new(0.5, -0.5);
+        let temp = c64(0.5, -0.5);
         let k11l = array![
-            [
-                temp * (Complex64::new(0., -1.) * Complex64::new(0., -b).exp()),
-                temp * Complex64::new(0., -b).exp()
-            ],
-            [
-                temp * (Complex64::new(0., -1.) * Complex64::new(0., b).exp()),
-                temp * -(Complex64::new(0., b).exp())
-            ],
+            [temp * (M_IM * c64(0., -b).exp()), temp * c64(0., -b).exp()],
+            [temp * (M_IM * c64(0., b).exp()), temp * -(c64(0., b).exp())],
         ];
         let k11r = array![
             [
-                FRAC_1_SQRT_2 * (Complex64::new(0., 1.) * Complex64::new(0., -b).exp()),
-                FRAC_1_SQRT_2 * -Complex64::new(0., -b).exp()
+                FRAC_1_SQRT_2 * (IM * c64(0., -b).exp()),
+                FRAC_1_SQRT_2 * -c64(0., -b).exp()
             ],
             [
-                FRAC_1_SQRT_2 * Complex64::new(0., b).exp(),
-                FRAC_1_SQRT_2 * (Complex64::new(0., -1.) * Complex64::new(0., b).exp())
+                FRAC_1_SQRT_2 * c64(0., b).exp(),
+                FRAC_1_SQRT_2 * (M_IM * c64(0., b).exp())
             ],
         ];
         let k12l = aview2(&K12L_ARR);
         let k12r = aview2(&K12R_ARR);
         let k32l_k21l = array![
             [
-                FRAC_1_SQRT_2 * Complex64::new(1., (2. * b).cos()),
-                FRAC_1_SQRT_2 * (Complex64::new(0., 1.) * (2. * b).sin())
+                FRAC_1_SQRT_2 * c64(1., (2. * b).cos()),
+                FRAC_1_SQRT_2 * (IM * (2. * b).sin())
             ],
             [
-                FRAC_1_SQRT_2 * (Complex64::new(0., 1.) * (2. * b).sin()),
-                FRAC_1_SQRT_2 * Complex64::new(1., -(2. * b).cos())
+                FRAC_1_SQRT_2 * (IM * (2. * b).sin()),
+                FRAC_1_SQRT_2 * c64(1., -(2. * b).cos())
             ],
         ];
-        let temp = Complex64::new(0.5, 0.5);
+        let temp = c64(0.5, 0.5);
         let k21r = array![
             [
-                temp * (Complex64::new(0., -1.) * Complex64::new(0., -2. * b).exp()),
-                temp * Complex64::new(0., -2. * b).exp()
+                temp * (M_IM * c64(0., -2. * b).exp()),
+                temp * c64(0., -2. * b).exp()
             ],
             [
-                temp * (Complex64::new(0., 1.) * Complex64::new(0., 2. * b).exp()),
-                temp * Complex64::new(0., 2. * b).exp()
+                temp * (IM * c64(0., 2. * b).exp()),
+                temp * c64(0., 2. * b).exp()
             ],
         ];
-        const K22L_ARR: [[Complex64; 2]; 2] = [
-            [
-                Complex64::new(FRAC_1_SQRT_2, 0.),
-                Complex64::new(-FRAC_1_SQRT_2, 0.),
-            ],
-            [
-                Complex64::new(FRAC_1_SQRT_2, 0.),
-                Complex64::new(FRAC_1_SQRT_2, 0.),
-            ],
+        const K22L_ARR: GateArray1Q = [
+            [c64(FRAC_1_SQRT_2, 0.), c64(-FRAC_1_SQRT_2, 0.)],
+            [c64(FRAC_1_SQRT_2, 0.), c64(FRAC_1_SQRT_2, 0.)],
         ];
         let k22l = aview2(&K22L_ARR);
-        let k22r_arr: [[Complex64; 2]; 2] = [
-            [Complex64::zero(), Complex64::new(1., 0.)],
-            [Complex64::new(-1., 0.), Complex64::zero()],
-        ];
+        let k22r_arr: GateArray1Q = [[Complex64::zero(), C_ONE], [C_M_ONE, Complex64::zero()]];
         let k22r = aview2(&k22r_arr);
         let k31l = array![
             [
-                FRAC_1_SQRT_2 * Complex64::new(0., -b).exp(),
-                FRAC_1_SQRT_2 * Complex64::new(0., -b).exp()
+                FRAC_1_SQRT_2 * c64(0., -b).exp(),
+                FRAC_1_SQRT_2 * c64(0., -b).exp()
             ],
             [
-                FRAC_1_SQRT_2 * -Complex64::new(0., b).exp(),
-                FRAC_1_SQRT_2 * Complex64::new(0., b).exp()
+                FRAC_1_SQRT_2 * -c64(0., b).exp(),
+                FRAC_1_SQRT_2 * c64(0., b).exp()
             ],
         ];
-        let temp = Complex64::new(0., 1.);
         let k31r = array![
-            [temp * Complex64::new(0., b).exp(), Complex64::zero()],
-            [Complex64::zero(), temp * -Complex64::new(0., -b).exp()],
+            [IM * c64(0., b).exp(), Complex64::zero()],
+            [Complex64::zero(), M_IM * c64(0., -b).exp()],
         ];
-        let temp = Complex64::new(0.5, 0.5);
+        let temp = c64(0.5, 0.5);
         let k32r = array![
+            [temp * c64(0., b).exp(), temp * -c64(0., -b).exp()],
             [
-                temp * Complex64::new(0., b).exp(),
-                temp * -Complex64::new(0., -b).exp()
-            ],
-            [
-                temp * (Complex64::new(0., -1.) * Complex64::new(0., b).exp()),
-                temp * (Complex64::new(0., -1.) * Complex64::new(0., -b).exp())
+                temp * (M_IM * c64(0., b).exp()),
+                temp * (M_IM * c64(0., -b).exp())
             ],
         ];
         let k1ld = transpose_conjugate(basis_decomposer.K1l.view());
@@ -1793,11 +1716,11 @@ impl TwoQubitBasisDecomposer {
 
     fn traces(&self, target: &TwoQubitWeylDecomposition) -> [Complex64; 4] {
         [
-            4. * Complex64::new(
+            4. * c64(
                 target.a.cos() * target.b.cos() * target.c.cos(),
                 target.a.sin() * target.b.sin() * target.c.sin(),
             ),
-            4. * Complex64::new(
+            4. * c64(
                 (PI4 - target.a).cos()
                     * (self.basis_decomposer.b - target.b).cos()
                     * target.c.cos(),
@@ -1805,8 +1728,8 @@ impl TwoQubitBasisDecomposer {
                     * (self.basis_decomposer.b - target.b).sin()
                     * target.c.sin(),
             ),
-            Complex64::new(4. * target.c.cos(), 0.),
-            Complex64::new(4., 0.),
+            c64(4. * target.c.cos(), 0.),
+            c64(4., 0.),
         ]
     }
 
diff --git a/crates/accelerate/src/uc_gate.rs b/crates/accelerate/src/uc_gate.rs
index 3a5f74a6f0b1..21fd7fa04656 100644
--- a/crates/accelerate/src/uc_gate.rs
+++ b/crates/accelerate/src/uc_gate.rs
@@ -21,14 +21,14 @@ use ndarray::prelude::*;
 use numpy::{IntoPyArray, PyReadonlyArray2};
 
 use crate::euler_one_qubit_decomposer::det_one_qubit;
+use qiskit_circuit::util::{c64, C_ZERO, IM};
 
-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);
+const RZ_PI2_11: Complex64 = c64(FRAC_1_SQRT_2, -FRAC_1_SQRT_2);
+const RZ_PI2_00: Complex64 = c64(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.
@@ -48,10 +48,10 @@ fn demultiplex_single_uc(
     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 r1 = (IM / 2. * (PI / 2. - phi / 2. - x11.arg())).exp();
+    let r2 = (IM / 2. * (PI / 2. - phi / 2. + x11.arg() + PI)).exp();
 
-    let r = array![[r1, Complex64::new(0., 0.)], [Complex64::new(0., 0.), r2],];
+    let r = array![[r1, C_ZERO], [C_ZERO, r2],];
 
     let decomp = r
         .dot(&x)
@@ -67,7 +67,7 @@ fn demultiplex_single_uc(
     // 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 {
+    if (diag[0] + IM).abs() < EPS {
         diag = diag.slice(s![..;-1]).to_owned();
         u = u.slice(s![.., ..;-1]).to_owned();
     }
diff --git a/crates/circuit/src/gate_matrix.rs b/crates/circuit/src/gate_matrix.rs
index 2a3fcdf88289..2f085ea79c0a 100644
--- a/crates/circuit/src/gate_matrix.rs
+++ b/crates/circuit/src/gate_matrix.rs
@@ -10,22 +10,16 @@
 // copyright notice, and modified files need to carry a notice indicating
 // that they have been altered from the originals.
 
-use num_complex::Complex64;
 use std::f64::consts::FRAC_1_SQRT_2;
 
-// num-complex exposes an equivalent function but it's not a const function
-// so it's not compatible with static definitions. This is a const func and
-// just reduces the amount of typing we need.
-#[inline(always)]
-const fn c64(re: f64, im: f64) -> Complex64 {
-    Complex64::new(re, im)
-}
+use crate::util::{
+    c64, GateArray0Q, GateArray1Q, GateArray2Q, GateArray3Q, C_M_ONE, C_ONE, C_ZERO, IM, M_IM,
+};
 
-pub static ONE_QUBIT_IDENTITY: [[Complex64; 2]; 2] =
-    [[c64(1., 0.), c64(0., 0.)], [c64(0., 0.), c64(1., 0.)]];
+pub static ONE_QUBIT_IDENTITY: GateArray1Q = [[C_ONE, C_ZERO], [C_ZERO, C_ONE]];
 
 #[inline]
-pub fn r_gate(theta: f64, phi: f64) -> [[Complex64; 2]; 2] {
+pub fn r_gate(theta: f64, phi: f64) -> GateArray1Q {
     let half_theta = theta / 2.;
     let cost = c64(half_theta.cos(), 0.);
     let sint = half_theta.sin();
@@ -38,7 +32,7 @@ pub fn r_gate(theta: f64, phi: f64) -> [[Complex64; 2]; 2] {
 }
 
 #[inline]
-pub fn rx_gate(theta: f64) -> [[Complex64; 2]; 2] {
+pub fn rx_gate(theta: f64) -> GateArray1Q {
     let half_theta = theta / 2.;
     let cos = c64(half_theta.cos(), 0.);
     let isin = c64(0., -half_theta.sin());
@@ -46,7 +40,7 @@ pub fn rx_gate(theta: f64) -> [[Complex64; 2]; 2] {
 }
 
 #[inline]
-pub fn ry_gate(theta: f64) -> [[Complex64; 2]; 2] {
+pub fn ry_gate(theta: f64) -> GateArray1Q {
     let half_theta = theta / 2.;
     let cos = c64(half_theta.cos(), 0.);
     let sin = c64(half_theta.sin(), 0.);
@@ -54,213 +48,150 @@ pub fn ry_gate(theta: f64) -> [[Complex64; 2]; 2] {
 }
 
 #[inline]
-pub fn rz_gate(theta: f64) -> [[Complex64; 2]; 2] {
+pub fn rz_gate(theta: f64) -> GateArray1Q {
     let ilam2 = c64(0., 0.5 * theta);
-    [[(-ilam2).exp(), c64(0., 0.)], [c64(0., 0.), ilam2.exp()]]
+    [[(-ilam2).exp(), C_ZERO], [C_ZERO, ilam2.exp()]]
 }
 
-pub static H_GATE: [[Complex64; 2]; 2] = [
+pub static H_GATE: GateArray1Q = [
     [c64(FRAC_1_SQRT_2, 0.), c64(FRAC_1_SQRT_2, 0.)],
     [c64(FRAC_1_SQRT_2, 0.), c64(-FRAC_1_SQRT_2, 0.)],
 ];
 
-pub static CX_GATE: [[Complex64; 4]; 4] = [
-    [c64(1., 0.), c64(0., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(0., 0.), c64(1., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(1., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(1., 0.), c64(0., 0.), c64(0., 0.)],
+pub static CX_GATE: GateArray2Q = [
+    [C_ONE, C_ZERO, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ZERO, C_ONE],
+    [C_ZERO, C_ZERO, C_ONE, C_ZERO],
+    [C_ZERO, C_ONE, C_ZERO, C_ZERO],
 ];
 
-pub static SX_GATE: [[Complex64; 2]; 2] = [
+pub static SX_GATE: GateArray1Q = [
     [c64(0.5, 0.5), c64(0.5, -0.5)],
     [c64(0.5, -0.5), c64(0.5, 0.5)],
 ];
 
-pub static SXDG_GATE: [[Complex64; 2]; 2] = [
+pub static SXDG_GATE: GateArray1Q = [
     [c64(0.5, -0.5), c64(0.5, 0.5)],
     [c64(0.5, 0.5), c64(0.5, -0.5)],
 ];
 
-pub static X_GATE: [[Complex64; 2]; 2] = [[c64(0., 0.), c64(1., 0.)], [c64(1., 0.), c64(0., 0.)]];
+pub static X_GATE: GateArray1Q = [[C_ZERO, C_ONE], [C_ONE, C_ZERO]];
 
-pub static Z_GATE: [[Complex64; 2]; 2] = [[c64(1., 0.), c64(0., 0.)], [c64(0., 0.), c64(-1., 0.)]];
+pub static Z_GATE: GateArray1Q = [[C_ONE, C_ZERO], [C_ZERO, C_M_ONE]];
 
-pub static Y_GATE: [[Complex64; 2]; 2] = [[c64(0., 0.), c64(0., -1.)], [c64(0., 1.), c64(0., 0.)]];
+pub static Y_GATE: GateArray1Q = [[C_ZERO, M_IM], [IM, C_ZERO]];
 
-pub static CZ_GATE: [[Complex64; 4]; 4] = [
-    [c64(1., 0.), c64(0., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(1., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(1., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(0., 0.), c64(-1., 0.)],
+pub static CZ_GATE: GateArray2Q = [
+    [C_ONE, C_ZERO, C_ZERO, C_ZERO],
+    [C_ZERO, C_ONE, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ONE, C_ZERO],
+    [C_ZERO, C_ZERO, C_ZERO, C_M_ONE],
 ];
 
-pub static CY_GATE: [[Complex64; 4]; 4] = [
-    [c64(1., 0.), c64(0., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(0., 0.), c64(0., -1.)],
-    [c64(0., 0.), c64(0., 0.), c64(1., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 1.), c64(0., 0.), c64(0., 0.)],
+pub static CY_GATE: GateArray2Q = [
+    [C_ONE, C_ZERO, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ZERO, M_IM],
+    [C_ZERO, C_ZERO, C_ONE, C_ZERO],
+    [C_ZERO, IM, C_ZERO, C_ZERO],
 ];
 
-pub static CCX_GATE: [[Complex64; 8]; 8] = [
+pub static CCX_GATE: GateArray3Q = [
     [
-        c64(1., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
+        C_ONE, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO,
     ],
     [
-        c64(0., 0.),
-        c64(1., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
+        C_ZERO, C_ONE, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO,
     ],
     [
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(1., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
+        C_ZERO, C_ZERO, C_ONE, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO,
     ],
     [
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(1., 0.),
+        C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ONE,
     ],
     [
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(1., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
+        C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ONE, C_ZERO, C_ZERO, C_ZERO,
     ],
     [
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(1., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
+        C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ONE, C_ZERO, C_ZERO,
     ],
     [
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(1., 0.),
-        c64(0., 0.),
+        C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ZERO, C_ONE, C_ZERO,
     ],
     [
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(1., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
-        c64(0., 0.),
+        C_ZERO, C_ZERO, C_ZERO, C_ONE, C_ZERO, C_ZERO, C_ZERO, C_ZERO,
     ],
 ];
 
-pub static ECR_GATE: [[Complex64; 4]; 4] = [
+pub static ECR_GATE: GateArray2Q = [
     [
-        c64(0., 0.),
+        C_ZERO,
         c64(FRAC_1_SQRT_2, 0.),
-        c64(0., 0.),
+        C_ZERO,
         c64(0., FRAC_1_SQRT_2),
     ],
     [
         c64(FRAC_1_SQRT_2, 0.),
-        c64(0., 0.),
+        C_ZERO,
         c64(0., -FRAC_1_SQRT_2),
-        c64(0., 0.),
+        C_ZERO,
     ],
     [
-        c64(0., 0.),
+        C_ZERO,
         c64(0., FRAC_1_SQRT_2),
-        c64(0., 0.),
+        C_ZERO,
         c64(FRAC_1_SQRT_2, 0.),
     ],
     [
         c64(0., -FRAC_1_SQRT_2),
-        c64(0., 0.),
+        C_ZERO,
         c64(FRAC_1_SQRT_2, 0.),
-        c64(0., 0.),
+        C_ZERO,
     ],
 ];
 
-pub static SWAP_GATE: [[Complex64; 4]; 4] = [
-    [c64(1., 0.), c64(0., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(1., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(1., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(0., 0.), c64(1., 0.)],
+pub static SWAP_GATE: GateArray2Q = [
+    [C_ONE, C_ZERO, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ONE, C_ZERO],
+    [C_ZERO, C_ONE, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ZERO, C_ONE],
 ];
-pub static ISWAP_GATE: [[Complex64; 4]; 4] = [
-    [c64(1., 0.), c64(0., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(0., 1.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 1.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(0., 0.), c64(1., 0.)],
+pub static ISWAP_GATE: GateArray2Q = [
+    [C_ONE, C_ZERO, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, IM, C_ZERO],
+    [C_ZERO, IM, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ZERO, C_ONE],
 ];
 
-pub static S_GATE: [[Complex64; 2]; 2] = [[c64(1., 0.), c64(0., 0.)], [c64(0., 0.), c64(0., 1.)]];
+pub static S_GATE: GateArray1Q = [[C_ONE, C_ZERO], [C_ZERO, IM]];
 
-pub static SDG_GATE: [[Complex64; 2]; 2] =
-    [[c64(1., 0.), c64(0., 0.)], [c64(0., 0.), c64(0., -1.)]];
+pub static SDG_GATE: GateArray1Q = [[C_ONE, C_ZERO], [C_ZERO, M_IM]];
 
-pub static T_GATE: [[Complex64; 2]; 2] = [
-    [c64(1., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(FRAC_1_SQRT_2, FRAC_1_SQRT_2)],
-];
+pub static T_GATE: GateArray1Q = [[C_ONE, C_ZERO], [C_ZERO, c64(FRAC_1_SQRT_2, FRAC_1_SQRT_2)]];
 
-pub static TDG_GATE: [[Complex64; 2]; 2] = [
-    [c64(1., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(FRAC_1_SQRT_2, -FRAC_1_SQRT_2)],
+pub static TDG_GATE: GateArray1Q = [
+    [C_ONE, C_ZERO],
+    [C_ZERO, c64(FRAC_1_SQRT_2, -FRAC_1_SQRT_2)],
 ];
 
-pub static DCX_GATE: [[Complex64; 4]; 4] = [
-    [c64(1., 0.), c64(0., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(0., 0.), c64(1., 0.)],
-    [c64(0., 0.), c64(1., 0.), c64(0., 0.), c64(0., 0.)],
-    [c64(0., 0.), c64(0., 0.), c64(1., 0.), c64(0., 0.)],
+pub static DCX_GATE: GateArray2Q = [
+    [C_ONE, C_ZERO, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ZERO, C_ONE],
+    [C_ZERO, C_ONE, C_ZERO, C_ZERO],
+    [C_ZERO, C_ZERO, C_ONE, C_ZERO],
 ];
 
 #[inline]
-pub fn global_phase_gate(theta: f64) -> [[Complex64; 1]; 1] {
+pub fn global_phase_gate(theta: f64) -> GateArray0Q {
     [[c64(0., theta).exp()]]
 }
 
 #[inline]
-pub fn phase_gate(lam: f64) -> [[Complex64; 2]; 2] {
-    [
-        [c64(1., 0.), c64(0., 0.)],
-        [c64(0., 0.), c64(0., lam).exp()],
-    ]
+pub fn phase_gate(lam: f64) -> GateArray1Q {
+    [[C_ONE, C_ZERO], [C_ZERO, c64(0., lam).exp()]]
 }
 
 #[inline]
-pub fn u_gate(theta: f64, phi: f64, lam: f64) -> [[Complex64; 2]; 2] {
+pub fn u_gate(theta: f64, phi: f64, lam: f64) -> GateArray1Q {
     let cos = (theta / 2.).cos();
     let sin = (theta / 2.).sin();
     [
@@ -270,37 +201,34 @@ pub fn u_gate(theta: f64, phi: f64, lam: f64) -> [[Complex64; 2]; 2] {
 }
 
 #[inline]
-pub fn xx_minus_yy_gate(theta: f64, beta: f64) -> [[Complex64; 4]; 4] {
+pub fn xx_minus_yy_gate(theta: f64, beta: f64) -> GateArray2Q {
     let cos = (theta / 2.).cos();
     let sin = (theta / 2.).sin();
     [
         [
             c64(cos, 0.),
-            c64(0., 0.),
-            c64(0., 0.),
+            C_ZERO,
+            C_ZERO,
             c64(0., -sin) * c64(0., -beta).exp(),
         ],
-        [c64(0., 0.), c64(1., 0.), c64(0., 0.), c64(0., 0.)],
-        [c64(0., 0.), c64(0., 0.), c64(1., 0.), c64(0., 0.)],
+        [C_ZERO, C_ONE, C_ZERO, C_ZERO],
+        [C_ZERO, C_ZERO, C_ONE, C_ZERO],
         [
             c64(0., -sin) * c64(0., beta).exp(),
-            c64(0., 0.),
-            c64(0., 0.),
+            C_ZERO,
+            C_ZERO,
             c64(cos, 0.),
         ],
     ]
 }
 
 #[inline]
-pub fn u1_gate(lam: f64) -> [[Complex64; 2]; 2] {
-    [
-        [c64(1., 0.), c64(0., 0.)],
-        [c64(0., 0.), c64(0., lam).exp()],
-    ]
+pub fn u1_gate(lam: f64) -> GateArray1Q {
+    [[C_ONE, C_ZERO], [C_ZERO, c64(0., lam).exp()]]
 }
 
 #[inline]
-pub fn u2_gate(phi: f64, lam: f64) -> [[Complex64; 2]; 2] {
+pub fn u2_gate(phi: f64, lam: f64) -> GateArray1Q {
     [
         [
             c64(FRAC_1_SQRT_2, 0.),
@@ -314,7 +242,7 @@ pub fn u2_gate(phi: f64, lam: f64) -> [[Complex64; 2]; 2] {
 }
 
 #[inline]
-pub fn u3_gate(theta: f64, phi: f64, lam: f64) -> [[Complex64; 2]; 2] {
+pub fn u3_gate(theta: f64, phi: f64, lam: f64) -> GateArray1Q {
     let cos = (theta / 2.).cos();
     let sin = (theta / 2.).sin();
     [
@@ -324,23 +252,23 @@ pub fn u3_gate(theta: f64, phi: f64, lam: f64) -> [[Complex64; 2]; 2] {
 }
 
 #[inline]
-pub fn xx_plus_yy_gate(theta: f64, beta: f64) -> [[Complex64; 4]; 4] {
+pub fn xx_plus_yy_gate(theta: f64, beta: f64) -> GateArray2Q {
     let cos = (theta / 2.).cos();
     let sin = (theta / 2.).sin();
     [
-        [c64(1., 0.), c64(0., 0.), c64(0., 0.), c64(0., 0.)],
+        [C_ONE, C_ZERO, C_ZERO, C_ZERO],
         [
-            c64(0., 0.),
+            C_ZERO,
             c64(cos, 0.),
             c64(0., -sin) * c64(0., -beta).exp(),
-            c64(0., 0.),
+            C_ZERO,
         ],
         [
-            c64(0., 0.),
+            C_ZERO,
             c64(0., -sin) * c64(0., beta).exp(),
             c64(cos, 0.),
-            c64(0., 0.),
+            C_ZERO,
         ],
-        [c64(0., 0.), c64(0., 0.), c64(0., 0.), c64(1., 0.)],
+        [C_ZERO, C_ZERO, C_ZERO, C_ONE],
     ]
 }
diff --git a/crates/circuit/src/lib.rs b/crates/circuit/src/lib.rs
index d7f285911750..9fcaa36480cf 100644
--- a/crates/circuit/src/lib.rs
+++ b/crates/circuit/src/lib.rs
@@ -17,6 +17,7 @@ pub mod gate_matrix;
 pub mod imports;
 pub mod operations;
 pub mod parameter_table;
+pub mod util;
 
 mod bit_data;
 mod interner;
diff --git a/crates/circuit/src/operations.rs b/crates/circuit/src/operations.rs
index e0e93726735d..ff730744c80f 100644
--- a/crates/circuit/src/operations.rs
+++ b/crates/circuit/src/operations.rs
@@ -24,9 +24,6 @@ use pyo3::prelude::*;
 use pyo3::{intern, IntoPy, Python};
 use smallvec::smallvec;
 
-const PI2: f64 = PI / 2.0;
-const PI4: f64 = PI / 4.0;
-
 /// Valid types for an operation field in a CircuitInstruction
 ///
 /// These are basically the types allowed in a QuantumCircuit
@@ -563,7 +560,11 @@ impl Operation for StandardGate {
                         1,
                         [(
                             Self::UGate,
-                            smallvec![Param::Float(PI), Param::Float(PI2), Param::Float(PI2),],
+                            smallvec![
+                                Param::Float(PI),
+                                Param::Float(PI / 2.),
+                                Param::Float(PI / 2.),
+                            ],
                             smallvec![Qubit(0)],
                         )],
                         FLOAT_ZERO,
@@ -732,7 +733,7 @@ impl Operation for StandardGate {
                         1,
                         [(
                             Self::UGate,
-                            smallvec![Param::Float(PI2), Param::Float(0.), Param::Float(PI)],
+                            smallvec![Param::Float(PI / 2.), Param::Float(0.), Param::Float(PI)],
                             smallvec![Qubit(0)],
                         )],
                         FLOAT_ZERO,
@@ -763,7 +764,7 @@ impl Operation for StandardGate {
                         1,
                         [(
                             Self::PhaseGate,
-                            smallvec![Param::Float(PI2)],
+                            smallvec![Param::Float(PI / 2.)],
                             smallvec![Qubit(0)],
                         )],
                         FLOAT_ZERO,
@@ -793,7 +794,7 @@ impl Operation for StandardGate {
                         1,
                         [(
                             Self::PhaseGate,
-                            smallvec![Param::Float(-PI2)],
+                            smallvec![Param::Float(-PI / 2.)],
                             smallvec![Qubit(0)],
                         )],
                         FLOAT_ZERO,
@@ -823,7 +824,7 @@ impl Operation for StandardGate {
                         1,
                         [(
                             Self::PhaseGate,
-                            smallvec![Param::Float(PI4)],
+                            smallvec![Param::Float(PI / 4.)],
                             smallvec![Qubit(0)],
                         )],
                         FLOAT_ZERO,
@@ -853,7 +854,7 @@ impl Operation for StandardGate {
                         1,
                         [(
                             Self::PhaseGate,
-                            smallvec![Param::Float(-PI4)],
+                            smallvec![Param::Float(-PI / 4.)],
                             smallvec![Qubit(0)],
                         )],
                         FLOAT_ZERO,
@@ -895,9 +896,9 @@ impl Operation for StandardGate {
                                 smallvec![multiply_param(beta, -1.0, py)],
                                 q1.clone(),
                             ),
-                            (Self::RZGate, smallvec![Param::Float(-PI2)], q0.clone()),
+                            (Self::RZGate, smallvec![Param::Float(-PI / 2.)], q0.clone()),
                             (Self::SXGate, smallvec![], q0.clone()),
-                            (Self::RZGate, smallvec![Param::Float(PI2)], q0.clone()),
+                            (Self::RZGate, smallvec![Param::Float(PI / 2.)], q0.clone()),
                             (Self::SGate, smallvec![], q1.clone()),
                             (Self::CXGate, smallvec![], q0_1.clone()),
                             (
@@ -912,9 +913,9 @@ impl Operation for StandardGate {
                             ),
                             (Self::CXGate, smallvec![], q0_1),
                             (Self::SdgGate, smallvec![], q1.clone()),
-                            (Self::RZGate, smallvec![Param::Float(-PI2)], q0.clone()),
+                            (Self::RZGate, smallvec![Param::Float(-PI / 2.)], q0.clone()),
                             (Self::SXdgGate, smallvec![], q0.clone()),
-                            (Self::RZGate, smallvec![Param::Float(PI2)], q0),
+                            (Self::RZGate, smallvec![Param::Float(PI / 2.)], q0),
                             (Self::RZGate, smallvec![beta.clone()], q1),
                         ],
                         FLOAT_ZERO,
@@ -934,9 +935,9 @@ impl Operation for StandardGate {
                         2,
                         [
                             (Self::RZGate, smallvec![beta.clone()], q0.clone()),
-                            (Self::RZGate, smallvec![Param::Float(-PI2)], q1.clone()),
+                            (Self::RZGate, smallvec![Param::Float(-PI / 2.)], q1.clone()),
                             (Self::SXGate, smallvec![], q1.clone()),
-                            (Self::RZGate, smallvec![Param::Float(PI2)], q1.clone()),
+                            (Self::RZGate, smallvec![Param::Float(PI / 2.)], q1.clone()),
                             (Self::SGate, smallvec![], q0.clone()),
                             (Self::CXGate, smallvec![], q1_0.clone()),
                             (
@@ -951,9 +952,9 @@ impl Operation for StandardGate {
                             ),
                             (Self::CXGate, smallvec![], q1_0),
                             (Self::SdgGate, smallvec![], q0.clone()),
-                            (Self::RZGate, smallvec![Param::Float(-PI2)], q1.clone()),
+                            (Self::RZGate, smallvec![Param::Float(-PI / 2.)], q1.clone()),
                             (Self::SXdgGate, smallvec![], q1.clone()),
-                            (Self::RZGate, smallvec![Param::Float(PI2)], q1),
+                            (Self::RZGate, smallvec![Param::Float(PI / 2.)], q1),
                             (Self::RZGate, smallvec![multiply_param(beta, -1.0, py)], q0),
                         ],
                         FLOAT_ZERO,
@@ -964,7 +965,7 @@ impl Operation for StandardGate {
             Self::CRXGate | Self::CRYGate | Self::CRZGate => todo!(),
             Self::RGate => Python::with_gil(|py| -> Option<CircuitData> {
                 let theta_expr = clone_param(&params[0], py);
-                let phi_expr1 = add_param(&params[1], -PI2, py);
+                let phi_expr1 = add_param(&params[1], -PI / 2., py);
                 let phi_expr2 = multiply_param(&phi_expr1, -1.0, py);
                 let defparams = smallvec![theta_expr, phi_expr1, phi_expr2];
                 Some(
diff --git a/crates/circuit/src/util.rs b/crates/circuit/src/util.rs
new file mode 100644
index 000000000000..11562b0a48cd
--- /dev/null
+++ b/crates/circuit/src/util.rs
@@ -0,0 +1,48 @@
+// 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;
+
+// This is a very conservative version of an abbreviation for constructing new Complex64.
+// A couple of alternatives to this function are
+// `c64<T: Into<f64>, V: Into<f64>>(re: T, im: V) -> Complex64`
+// Disadvantages are:
+//  1. Some people don't like that this allows things like `c64(1, 0)`. Presumably,
+//     they prefer a more explicit construction.
+//  2. This will not work in `const` and `static` constructs.
+// Another alternative is
+//   macro_rules! c64 {
+//       ($re: expr, $im: expr $(,)*) => {
+//           Complex64::new($re as f64, $im as f64)
+//       };
+// Advantages: This allows things like `c64!(1, 2.0)`, including in
+// `static` and `const` constructs.
+//  Disadvantages:
+//  1. Three characters `c64!` rather than two `c64`.
+//  2. Some people prefer the opposite of the advantages, i.e. more explicitness.
+/// Create a new [`Complex<f64>`]
+#[inline(always)]
+pub const fn c64(re: f64, im: f64) -> Complex64 {
+    Complex64::new(re, im)
+}
+
+pub type GateArray0Q = [[Complex64; 1]; 1];
+pub type GateArray1Q = [[Complex64; 2]; 2];
+pub type GateArray2Q = [[Complex64; 4]; 4];
+pub type GateArray3Q = [[Complex64; 8]; 8];
+
+// Use prefix `C_` to distinguish from real, for example
+pub const C_ZERO: Complex64 = c64(0., 0.);
+pub const C_ONE: Complex64 = c64(1., 0.);
+pub const C_M_ONE: Complex64 = c64(-1., 0.);
+pub const IM: Complex64 = c64(0., 1.);
+pub const M_IM: Complex64 = c64(0., -1.);