Support commutation check between Pauli-based gates and standard gates (After reverting #15488)

crates/pyext/src/lib.rs

@@ -72,6 +72,7 @@ fn _accelerate(m: &Bound<PyModule>) -> PyResult<()> {
     add_submodule(m, ::qiskit_transpiler::passes::sabre::sabre, "sabre")?;
     add_submodule(m, ::qiskit_accelerate::sampled_exp_val::sampled_exp_val, "sampled_exp_val")?;
     add_submodule(m, ::qiskit_quantum_info::sparse_observable::sparse_observable, "sparse_observable")?;
+    add_submodule(m, ::qiskit_quantum_info::standard_generators::standard_generators, "standard_generators")?;
     add_submodule(m, ::qiskit_quantum_info::sparse_pauli_op::sparse_pauli_op, "sparse_pauli_op")?;
     add_submodule(m, ::qiskit_transpiler::passes::scheduling_mod, "scheduling")?;
     add_submodule(m, ::qiskit_synthesis::matrix::sim::unitary_sim, "unitary_sim")?;

crates/quantum_info/src/lib.rs

@@ -14,6 +14,7 @@ pub mod clifford;
 pub mod pauli_lindblad_map;
 pub mod sparse_observable;
 pub mod sparse_pauli_op;
+pub mod standard_generators;
 pub mod unitary_compose;
 pub mod versor_u2;
 

crates/quantum_info/src/standard_generators.rs

@@ -0,0 +1,350 @@
+// This code is part of Qiskit.
+//
+// (C) Copyright IBM 2026
+//
+// 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 https://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.
+
+// AI Attribution:
+// This module was developed with assistance from GitHub Copilot integrated in VS Code.
+// The underlying model : Claude Haiku 4.5.
+// Portions of the Pauli generator mapping logic and SoA layout were generated
+// and then manually verified for mathematical correctness against the commutation logic.
+//
+// This module maps standard quantum gates to their Hamiltonian generators H_gate such that:
+//
+//   gate = exp(-i * H_gate)   (up to global phase)
+//
+// The generator is returned as a SparseObservable (a sum of Pauli tensor products).
+//
+// For single-qubit gates: X = exp(-i*(pi/2)*X), H = exp(-i*(pi/2)*(X+Z)/sqrt(2)), etc.
+// For multi-qubit controlled gates: CX = exp(-i*(pi/4)*(ZX - ZI - IX)), etc.
+// For rotation gates: RX(t) = exp(-i*(t/2)*X), RY(t) = exp(-i*(t/2)*Y), etc.
+//
+// The SoA (Struct-of-Arrays) layout is used: `bit_terms`, `indices`, and `boundaries` are
+// flattened arrays. Term `i` uses bit_terms[boundaries[i]..boundaries[i+1]] and
+// indices[boundaries[i]..boundaries[i+1]] for its Pauli operators and qubit targets.
+
+use crate::sparse_observable::SparseObservable;
+use pyo3::prelude::*;
+use pyo3::wrap_pyfunction;
+
+use num_complex::Complex64;
+use qiskit_circuit::operations::{Operation, Param, StandardGate};
+use qiskit_util::complex::{C_ZERO, c64};
+use std::f64::consts::{FRAC_PI_2, FRAC_PI_4, FRAC_PI_8, SQRT_2};
+
+const C_FRAC_PI_2: Complex64 = c64(FRAC_PI_2, 0.0);
+const C_FRAC_PI_4: Complex64 = c64(FRAC_PI_4, 0.0);
+const C_FRAC_PI_8: Complex64 = c64(FRAC_PI_8, 0.0);
+const C_FRAC_PI_2_SQRT_2: Complex64 = c64(FRAC_PI_2 / SQRT_2, 0.0);
+const C_M_FRAC_PI_4: Complex64 = c64(-FRAC_PI_4, 0.0);
+const C_M_FRAC_PI_8: Complex64 = c64(-FRAC_PI_8, 0.0);
+const C_M_FRAC_PI_2_SQRT_2: Complex64 = c64(-FRAC_PI_2 / SQRT_2, 0.0);
+
+const BETA_TOLERANCE: f64 = 1e-10;
+
+/// For parametric gates (e.g., `RX(theta)`), the generator $H$ depends on the
+/// gate parameters (e.g., $H = (\theta/2)X$). This function extracts parameter values
+/// from the `params` slice to compute the concrete coefficients for the returned
+/// `SparseObservable`.
+///
+/// If parameters are missing or symbolic (non-`Float`), it defaults to a coefficient
+/// of 1.0. This allows the commutation checker to still prove commutation in cases
+/// where the generator's Pauli structure alone is sufficient (e.g., `[theta*X, X] = 0`
+/// for any `theta`), but it avoids attempting to store parametric expressions, which
+/// `SparseObservable` does not currently support.
+pub fn standard_gate_exponent(gate: StandardGate, params: &[Param]) -> Option<SparseObservable> {
+    // let _params = params;
+    let fixed_params = params
+        .iter()
+        .map(|p| match p {
+            Param::Float(f) => *f,
+            Param::ParameterExpression(_) => 1.0,
+            Param::Obj(_) => panic!("StandardGate's don't have Param::Obj"),
+        })
+        .collect::<Vec<f64>>();
+
+    let num_qubits = gate.num_qubits();
+
+    use crate::sparse_observable::BitTerm::*;
+
+    let (coeffs, bit_terms, indices, boundaries) = match gate {
+        StandardGate::GlobalPhase => (vec![c64(-fixed_params[0], 0.)], vec![], vec![], vec![0, 0]),
+        // H = exp(-i pi/sqrt(8) (X + Z))
+        StandardGate::H => (
+            vec![C_FRAC_PI_2_SQRT_2, C_FRAC_PI_2_SQRT_2],
+            vec![X, Z],
+            vec![0, 0],
+            vec![0, 1, 2],
+        ),
+        // X = exp(-i pi/2 X), Y = exp(-i pi/2 Y), Z = exp(-i pi/2 Z)
+        StandardGate::X => (vec![C_FRAC_PI_2], vec![X], vec![0], vec![0, 1]),
+        StandardGate::Y => (vec![C_FRAC_PI_2], vec![Y], vec![0], vec![0, 1]),
+        StandardGate::Z => (vec![C_FRAC_PI_2], vec![Z], vec![0], vec![0, 1]),
+        // Identity: exp(-i 0) = I.
+        StandardGate::I => (vec![C_ZERO], vec![], vec![], vec![0, 0]),
+        // S = exp(-i pi/4 Z), Sdg = exp(-i (-pi/4) Z)
+        StandardGate::S => (vec![C_FRAC_PI_4], vec![Z], vec![0], vec![0, 1]),
+        StandardGate::Sdg => (vec![C_M_FRAC_PI_4], vec![Z], vec![0], vec![0, 1]),
+        // T = exp(-i pi/8 Z), Tdg = exp(-i (-pi/8) Z)
+        StandardGate::T => (vec![C_FRAC_PI_8], vec![Z], vec![0], vec![0, 1]),
+        StandardGate::Tdg => (vec![C_M_FRAC_PI_8], vec![Z], vec![0], vec![0, 1]),
+        // SX = exp(-i pi/4 X), SXdg = exp(-i (-pi/4) X)
+        StandardGate::SX => (vec![C_FRAC_PI_4], vec![X], vec![0], vec![0, 1]),
+        StandardGate::SXdg => (vec![C_M_FRAC_PI_4], vec![X], vec![0], vec![0, 1]),
+        // RX(t) = exp(-i t/2 X), RY, RZ=Phase, equivalently
+        StandardGate::RX => (
+            vec![c64(fixed_params[0] / 2.0, 0.0)],
+            vec![X],
+            vec![0],
+            vec![0, 1],
+        ),
+        StandardGate::RY => (
+            vec![c64(fixed_params[0] / 2.0, 0.0)],
+            vec![Y],
+            vec![0],
+            vec![0, 1],
+        ),
+        StandardGate::RZ | StandardGate::Phase | StandardGate::U1 => (
+            vec![c64(fixed_params[0] / 2.0, 0.0)],
+            vec![Z],
+            vec![0],
+            vec![0, 1],
+        ),
+        // CX = exp(-i pi/4 (ZX - ZI - IX))
+        StandardGate::CX => (
+            vec![C_M_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4],
+            vec![Z, X, Z, X],
+            vec![0, 1, 0, 1],
+            vec![0, 2, 3, 4],
+        ),
+        // CY = exp(-i pi/4 (ZY - ZI - IY))
+        StandardGate::CY => (
+            vec![C_M_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4],
+            vec![Z, Y, Z, Y],
+            vec![0, 1, 0, 1],
+            vec![0, 2, 3, 4],
+        ),
+        // CZ = exp(-i pi/4 (ZZ - ZI - IZ))
+        StandardGate::CZ => (
+            vec![C_M_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4],
+            vec![Z, Z, Z, Z],
+            vec![0, 1, 0, 1],
+            vec![0, 2, 3, 4],
+        ),
+        // CS = exp(-i (-pi/8) (ZZ - ZI - IZ))
+        StandardGate::CS => (
+            vec![C_FRAC_PI_8, C_FRAC_PI_8, C_M_FRAC_PI_8],
+            vec![Z, Z, Z, Z],
+            vec![0, 1, 0, 1],
+            vec![0, 1, 2, 4],
+        ),
+        // CSdg (inv sqrt(CZ)): factor +pi/8
+        StandardGate::CSdg => (
+            vec![C_M_FRAC_PI_8, C_M_FRAC_PI_8, C_FRAC_PI_8],
+            vec![Z, Z, Z, Z],
+            vec![0, 1, 0, 1],
+            vec![0, 1, 2, 4],
+        ),
+        // CSX (sqrt(CX)/controlled-SX): factor -pi/8
+        StandardGate::CSX => (
+            vec![C_FRAC_PI_8, C_FRAC_PI_8, C_M_FRAC_PI_8],
+            vec![Z, X, Z, X],
+            vec![0, 1, 0, 1],
+            vec![0, 1, 2, 4],
+        ),
+        // CRX(t) = exp(-i t/4 (-ZX + IX))
+        StandardGate::CRX => (
+            vec![
+                c64(-fixed_params[0] / 4.0, 0.0),
+                c64(fixed_params[0] / 4.0, 0.0),
+            ],
+            vec![Z, X, X],
+            vec![0, 1, 1],
+            vec![0, 2, 3],
+        ),
+        // CRY(t) = exp(-i t/4 (-ZY + IY))
+        StandardGate::CRY => (
+            vec![
+                c64(-fixed_params[0] / 4.0, 0.0),
+                c64(fixed_params[0] / 4.0, 0.0),
+            ],
+            vec![Z, Y, Y],
+            vec![0, 1, 1],
+            vec![0, 2, 3],
+        ),
+        // CRZ(t) = exp(-i t/4 (-ZZ + IZ))
+        StandardGate::CRZ => (
+            vec![
+                c64(-fixed_params[0] / 4.0, 0.0),
+                c64(fixed_params[0] / 4.0, 0.0),
+            ],
+            vec![Z, Z, Z],
+            vec![0, 1, 1],
+            vec![0, 2, 3],
+        ),
+        // CPhase(t) = exp(-i  t/4 (-ZZ + ZI + IZ))
+        // (same as CRZ but shifts both Z0 and Z1, not just Z1)
+        StandardGate::CPhase => (
+            vec![
+                c64(-fixed_params[0] / 4.0, 0.0),
+                c64(fixed_params[0] / 4.0, 0.0),
+                c64(fixed_params[0] / 4.0, 0.0),
+            ],
+            vec![Z, Z, Z, Z],
+            vec![0, 1, 0, 1],
+            vec![0, 2, 3, 4],
+        ),
+        // Swap = exp(-i pi/4 (XX + YY + ZZ))
+        StandardGate::Swap => (
+            vec![C_FRAC_PI_4, C_FRAC_PI_4, C_FRAC_PI_4],
+            vec![X, X, Y, Y, Z, Z],
+            vec![0, 1, 0, 1, 0, 1],
+            vec![0, 2, 4, 6],
+        ),
+        // ISwap = exp(-i (-pi/4) (XX + YY))
+        StandardGate::ISwap => (
+            vec![C_M_FRAC_PI_4, C_M_FRAC_PI_4],
+            vec![X, X, Y, Y],
+            vec![0, 1, 0, 1],
+            vec![0, 2, 4],
+        ),
+        // RXX(t) = exp(-i  t/2 XX)
+        StandardGate::RXX => (
+            vec![c64(fixed_params[0] / 2.0, 0.0)],
+            vec![X, X],
+            vec![0, 1],
+            vec![0, 2],
+        ),
+        // RYY(t) = exp(-i t/2 YY)
+        StandardGate::RYY => (
+            vec![c64(fixed_params[0] / 2.0, 0.0)],
+            vec![Y, Y],
+            vec![0, 1],
+            vec![0, 2],
+        ),
+        // RZZ(t) = exp(-i t/2 ZZ)
+        StandardGate::RZZ => (
+            vec![c64(fixed_params[0] / 2.0, 0.0)],
+            vec![Z, Z],
+            vec![0, 1],
+            vec![0, 2],
+        ),
+        // RZX(t) = exp(-i t/2 ZX)
+        StandardGate::RZX => (
+            vec![c64(fixed_params[0] / 2.0, 0.0)],
+            vec![Z, X],
+            vec![0, 1],
+            vec![0, 2],
+        ),
+        // XX+YY and XX-YY are just handled if the beta parameter is 0, in which case
+        // the generator is XX +- YY.
+        StandardGate::XXPlusYY | StandardGate::XXMinusYY => {
+            // This covers both parametric ``beta`` and if they are above tolerance
+            if fixed_params[1].abs() > BETA_TOLERANCE {
+                return None;
+            }
+
+            match gate {
+                StandardGate::XXPlusYY => (
+                    vec![
+                        c64(fixed_params[0] / 4.0, 0.0),
+                        c64(fixed_params[0] / 4.0, 0.0),
+                    ],
+                    vec![X, X, Y, Y],
+                    vec![0, 1, 0, 1],
+                    vec![0, 2, 4],
+                ),
+                StandardGate::XXMinusYY => (
+                    vec![
+                        c64(fixed_params[0] / 4.0, 0.0),
+                        c64(-fixed_params[0] / 4.0, 0.0),
+                    ],
+                    vec![X, X, Y, Y],
+                    vec![0, 1, 0, 1],
+                    vec![0, 2, 4],
+                ),
+                _ => unreachable!(),
+            }
+        }
+        // CCX = exp(-i pi/8 (ZZX - ZIX - IZX - ZZI + ZII + IZI + IIX))
+        StandardGate::CCX => (
+            vec![
+                C_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_FRAC_PI_8,
+                C_FRAC_PI_8,
+                C_FRAC_PI_8,
+            ],
+            vec![Z, Z, X, Z, X, Z, X, Z, Z, Z, Z, X],
+            vec![0, 1, 2, 0, 2, 1, 2, 0, 1, 0, 1, 2],
+            vec![0, 3, 5, 7, 9, 10, 11, 12],
+        ),
+        // Same as CCX but with Z on the target
+        StandardGate::CCZ => (
+            vec![
+                C_M_FRAC_PI_8,
+                C_FRAC_PI_8,
+                C_FRAC_PI_8,
+                C_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+            ],
+            vec![Z, Z, Z, Z, Z, Z, Z, Z, Z, Z, Z, Z],
+            vec![0, 1, 2, 0, 2, 1, 2, 0, 1, 0, 1, 2],
+            vec![0, 3, 5, 7, 9, 10, 11, 12],
+        ),
+        // CSwap = exp(-i pi/8 (ZII - ZXX - ZYY - ZZZ + IXX + IYY + IZZ))
+        StandardGate::CSwap => (
+            vec![
+                C_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_M_FRAC_PI_8,
+                C_FRAC_PI_8,
+                C_FRAC_PI_8,
+                C_FRAC_PI_8,
+            ],
+            vec![Z, Z, X, X, Z, Y, Y, Z, Z, Z, X, X, Y, Y, Z, Z],
+            vec![0, 0, 1, 2, 0, 1, 2, 0, 1, 2, 1, 2, 1, 2, 1, 2],
+            vec![0, 1, 4, 7, 10, 12, 14, 16],
+        ),
+        // ECR = exp(-i pi/sqrt(8) (IX - XY))
+        StandardGate::ECR => (
+            vec![C_FRAC_PI_2_SQRT_2, C_M_FRAC_PI_2_SQRT_2],
+            vec![X, Y, X],
+            vec![0, 0, 1],
+            vec![0, 1, 3],
+        ),
+        _ => return None,
+    };
+
+    // SAFETY: The internal data was constructed manually and is consistent.
+    Some(unsafe {
+        SparseObservable::new_unchecked(num_qubits, coeffs, bit_terms, indices, boundaries)
+    })
+}
+
+#[pyfunction(name = "_standard_gate_exponent")]
+#[pyo3(signature = (gate, params=None))]
+pub fn py_standard_gate_exponent(
+    gate: StandardGate,
+    params: Option<Vec<Param>>,
+) -> Option<SparseObservable> {
+    let params = params.unwrap_or_default();
+    standard_gate_exponent(gate, &params)
+}
+
+pub fn standard_generators(m: &Bound<PyModule>) -> PyResult<()> {
+    m.add_function(wrap_pyfunction!(py_standard_gate_exponent, m)?)?;
+    Ok(())
+}