Implement some standard gates in rust

.github/workflows/tests.yml

@@ -36,6 +36,15 @@ jobs:
           python -m pip install -U -r requirements.txt -c constraints.txt
           python -m pip install -U -r requirements-dev.txt -c constraints.txt
           python -m pip install -c constraints.txt -e .
+        if: matrix.python-version == '3.10'
+        env:
+          QISKIT_NO_CACHE_GATES: 1
+      - name: 'Install dependencies'
+        run: |
+          python -m pip install -U -r requirements.txt -c constraints.txt
+          python -m pip install -U -r requirements-dev.txt -c constraints.txt
+          python -m pip install -c constraints.txt -e .
+        if: matrix.python-version == '3.12'
       - name: 'Install optionals'
         run: |
           python -m pip install -r requirements-optional.txt -c constraints.txt

CONTRIBUTING.md

@@ -135,6 +135,18 @@ Note that in order to run `python setup.py ...` commands you need have build
 dependency packages installed in your environment, which are listed in the
 `pyproject.toml` file under the `[build-system]` section.
 
+### Compile time options
+
+When building qiskit from source there are options available to control how
+Qiskit is build. Right now the only option is if you set the environment
+variable `QISKIT_NO_CACHE_GATES=1` this will disable runtime caching of
+Python gate objects when accessing them from a `QuantumCircuit` or `DAGCircuit`.
+This makes a tradeoff between runtime performance for Python access and memory
+overhead. Caching gates will result in better runtime for users of Python at
+the cost of increased memory consumption. If you're working with any custom
+transpiler passes written in python or are otherwise using a workflow that
+repeatedly accesses the `operation` attribute of a `CircuitInstruction` or `op`
+attribute of `DAGOpNode` enabling caching is recommended.
 
 ## Issues and pull requests
 

Cargo.toml

@@ -16,6 +16,11 @@ license = "Apache-2.0"
 [workspace.dependencies]
 indexmap.version = "2.2.6"
 hashbrown.version = "0.14.0"
+num-complex = "0.4"
+ndarray = "^0.15.6"
+numpy = "0.21.0"
+smallvec = "1.13"
+
 # Most of the crates don't need the feature `extension-module`, since only `qiskit-pyext` builds an
 # actual C extension (the feature disables linking in `libpython`, which is forbidden in Python
 # distributions).  We only activate that feature when building the C extension module; we still need

crates/accelerate/Cargo.toml

@@ -11,29 +11,29 @@ doctest = false
 
 [dependencies]
 rayon = "1.10"
-numpy = "0.21.0"
+numpy.workspace = true
 rand = "0.8"
 rand_pcg = "0.3"
 rand_distr = "0.4.3"
 ahash = "0.8.11"
 num-traits = "0.2"
-num-complex = "0.4"
+num-complex.workspace = true
 num-bigint = "0.4"
 rustworkx-core = "0.14"
 faer = "0.19.0"
 itertools = "0.13.0"
 qiskit-circuit.workspace = true
 
 [dependencies.smallvec]
-version = "1.13"
+workspace = true
 features = ["union"]
 
 [dependencies.pyo3]
 workspace = true
 features = ["hashbrown", "indexmap", "num-complex", "num-bigint", "smallvec"]
 
 [dependencies.ndarray]
-version = "^0.15.6"
+workspace = true
 features = ["rayon", "approx-0_5"]
 
 [dependencies.approx]

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]

crates/accelerate/src/euler_one_qubit_decomposer.rs

@@ -31,8 +31,12 @@ use ndarray::prelude::*;
 use numpy::PyReadonlyArray2;
 use pyo3::pybacked::PyBackedStr;
 
+use qiskit_circuit::util::M_IM;
 use qiskit_circuit::SliceOrInt;
 
+// Macro c64!
+use qiskit_circuit::c64;
+
 pub const ANGLE_ZERO_EPSILON: f64 = 1e-12;
 
 #[pyclass(module = "qiskit._accelerate.euler_one_qubit_decomposer")]
@@ -855,16 +859,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 = (M_IM * det.ln()).re / 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());

crates/accelerate/src/isometry.rs

@@ -23,7 +23,8 @@ use itertools::Itertools;
 use ndarray::prelude::*;
 use numpy::{IntoPyArray, PyReadonlyArray1, PyReadonlyArray2};
 
-use crate::two_qubit_decompose::ONE_QUBIT_IDENTITY;
+use qiskit_circuit::gate_matrix::ONE_QUBIT_IDENTITY;
+use qiskit_circuit::util::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], ZERO], [ZERO, diag[2 * i + 1]]]).dot(&gate);
             res.into_pyarray_bound(py).into()
         })
         .collect()

crates/accelerate/src/pauli_exp_val.rs

@@ -18,6 +18,9 @@ use pyo3::prelude::*;
 use pyo3::wrap_pyfunction;
 use rayon::prelude::*;
 
+// Macro c64!
+use qiskit_circuit::c64;
+
 use crate::getenv_use_multiple_threads;
 
 const PARALLEL_THRESHOLD: usize = 19;
@@ -88,15 +91,15 @@ 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
                     - data_arr[index_1].re * data_arr[index_0].im,
             ))
         .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

crates/accelerate/src/sampled_exp_val.rs

@@ -17,6 +17,9 @@ use numpy::PyReadonlyArray1;
 use pyo3::prelude::*;
 use pyo3::wrap_pyfunction;
 
+// Macro c64!
+use qiskit_circuit::c64;
+
 use crate::pauli_exp_val::fast_sum;
 
 const OPER_TABLE_SIZE: usize = (b'Z' as usize) + 1;
@@ -81,7 +84,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)
 }

crates/accelerate/src/sparse_pauli_op.rs

@@ -22,9 +22,14 @@ use numpy::{PyArray1, PyArray2, PyReadonlyArray1, PyReadonlyArray2, PyUntypedArr
 use hashbrown::HashMap;
 use ndarray::{s, Array1, Array2, ArrayView1, ArrayView2, Axis};
 use num_complex::Complex64;
+// For method is_zero
 use num_traits::Zero;
 use rayon::prelude::*;
 
+use qiskit_circuit::util::{ONE, ZERO};
+// Macro c64!
+use qiskit_circuit::c64;
+
 use crate::rayon_ext::*;
 
 /// Find the unique elements of an array.
@@ -257,9 +262,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 +316,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(ZERO) += coeff;
         }
         for ((x, z), coeff) in hash_table {
-            if coeff == Complex64::new(0.0, 0.0) {
+            if coeff == ZERO {
                 continue;
             }
             self.x_like.push(x);
@@ -347,7 +352,7 @@ pub fn decompose_dense(
     let mut coeffs = vec![];
     if num_qubits > 0 {
         decompose_dense_inner(
-            Complex64::new(1.0, 0.0),
+            ONE,
             num_qubits,
             &[],
             operator.as_array(),
@@ -485,7 +490,7 @@ fn decompose_dense_inner(
     );
     recurse_if_nonzero(
         Pauli::Y,
-        0.5 * Complex64::i() * factor,
+        c64!(0, 0.5) * factor,
         &block.slice(s![..mid, mid..]) - &block.slice(s![mid.., ..mid]),
     );
     recurse_if_nonzero(
@@ -532,7 +537,7 @@ fn to_matrix_dense_inner(paulis: &MatrixCompressedPaulis, parallel: bool) -> Vec
         // Doing the initialisation 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(ZERO);
         for ((&x_like, &z_like), &coeff) in paulis
             .x_like
             .iter()
@@ -667,7 +672,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 = 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 +753,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 = 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 +849,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),
             ],
         }
     }