Add Rust matrix support for PauliProductRotation

crates/circuit/src/gate_matrix.rs

@@ -13,6 +13,10 @@
 use num_complex::Complex64;
 use std::f64::consts::FRAC_1_SQRT_2;
 
+// Ensure necessary imports are present
+use ndarray::Array2;
+
+
 use crate::util::{
     C_M_ONE, C_ONE, C_ZERO, GateArray0Q, GateArray1Q, GateArray2Q, GateArray3Q, GateArray4Q, IM,
     M_IM, c64,
@@ -516,3 +520,91 @@ pub fn xx_plus_yy_gate(theta: f64, beta: f64) -> GateArray2Q {
         [C_ZERO, C_ZERO, C_ZERO, C_ONE],
     ]
 }
+/// Compute the unitary matrix for a PauliProductRotation gate.
+///
+/// The gate implements: exp(-i * theta/2 * P) = cos(theta/2)*I - i*sin(theta/2)*P
+/// where P is the tensor product of Pauli operators given by `pauli_str`.
+///
+/// # Arguments
+/// * `theta`     - Rotation angle in radians (f64)
+/// * `pauli_str` - Pauli string label e.g. "XYZ", "ZZ", "X" (rightmost = qubit 0)
+///
+/// # Returns
+/// Dense 2^n × 2^n complex unitary matrix as `Array2<Complex64>`
+pub fn pauli_product_rotation_matrix(
+    angle: f64,
+    z: &[bool],
+    x: &[bool],
+) -> Array2<Complex64> {
+
+    let n = z.len();
+    let dim = 1usize << n;
+
+    let pauli_i = Array2::from_shape_vec(
+        (2, 2),
+        vec![c64(1.,0.), c64(0.,0.), c64(0.,0.), c64(1.,0.)],
+    ).unwrap();
+    let pauli_x = Array2::from_shape_vec(
+        (2, 2),
+        vec![c64(0.,0.), c64(1.,0.), c64(1.,0.), c64(0.,0.)],
+    ).unwrap();
+    let pauli_y = Array2::from_shape_vec(
+        (2, 2),
+        vec![c64(0.,0.), c64(0.,-1.), c64(0.,1.), c64(0.,0.)],
+    ).unwrap();
+    let pauli_z = Array2::from_shape_vec(
+        (2, 2),
+        vec![c64(1.,0.), c64(0.,0.), c64(0.,0.), c64(-1.,0.)],
+    ).unwrap();
+
+    // Build full tensor product using z/x symplectic representation:
+    // z=false x=false → I
+    // z=false x=true  → X
+    // z=true  x=true  → Y
+    // z=true  x=false → Z
+    let first = match (z[0], x[0]) {
+        (false, false) => pauli_i.clone(),
+        (false, true)  => pauli_x.clone(),
+        (true,  true)  => pauli_y.clone(),
+        (true,  false) => pauli_z.clone(),
+    };
+
+    let full_pauli = (1..n).fold(first, |acc, i| {
+        let single = match (z[i], x[i]) {
+            (false, false) => pauli_i.clone(),
+            (false, true)  => pauli_x.clone(),
+            (true,  true)  => pauli_y.clone(),
+            (true,  false) => pauli_z.clone(),
+        };
+        kron(&single, &acc)
+    });
+
+    let cos_val = (angle / 2.0).cos();
+    let sin_val = (angle / 2.0).sin();
+
+    let identity: Array2<Complex64> = Array2::eye(dim)
+        .mapv(|x: f64| c64(x, 0.));
+
+    identity.mapv(|v| v * c64(cos_val, 0.))
+        + full_pauli.mapv(|v| v * c64(0., -sin_val))
+
+}
+
+/// Kronecker (tensor) product of two complex matrices.
+fn kron(a: &Array2<Complex64>, b: &Array2<Complex64>) -> Array2<Complex64> {
+    let (ra, ca) = a.dim();
+    let (rb, cb) = b.dim();
+    let mut out = Array2::zeros((ra * rb, ca * cb));
+    for i in 0..ra {
+        for j in 0..ca {
+            let aij = a[[i, j]];
+            for p in 0..rb {
+                for q in 0..cb {
+                    out[[i * rb + p, j * cb + q]] = aij * b[[p, q]];
+                }
+            }
+        }
+    }
+    out
+}
+

crates/circuit/src/operations.rs

@@ -3335,22 +3335,18 @@ pub struct PauliProductRotation {
     pub angle: Param,
 }
 
-impl Operation for PauliProductRotation {
-    fn name(&self) -> &str {
-        "pauli_product_rotation"
-    }
-    fn num_qubits(&self) -> u32 {
-        self.z.len() as u32
-    }
-    fn num_clbits(&self) -> u32 {
-        0
-    }
-    fn num_params(&self) -> u32 {
-        1
-    }
-    fn directive(&self) -> bool {
-        false
-    }
+
+// Helper function to convert x/z bool vectors into a Pauli string in "IXYZ" notation.
+fn pauli_product_to_string(zs: &[bool], xs: &[bool]) -> String {
+    zs.iter()
+        .zip(xs.iter())
+        .map(|(&z, &x)| match (z, x) {
+            (false, false) => 'I',
+            (false, true) => 'X',
+            (true, false) => 'Z',
+            (true, true) => 'Y',
+        })
+        .collect()
 }
 
 impl PauliProductRotation {
@@ -3372,8 +3368,43 @@ impl PauliProductRotation {
             )?;
         Ok(gate.unbind())
     }
+
+    pub fn matrix(&self) -> Option<Array2<Complex64>> {
+        let angle = match self.angle {
+            Param::Float(f) => f,
+            _ => return None,
+        };
+        Some(gate_matrix::pauli_product_rotation_matrix(
+            angle,
+            &self.z,
+            &self.x,
+        ))
+    }
 }
 
+impl Operation for PauliProductRotation {
+    fn name(&self) -> &str {
+        "pauli_product_rotation"
+    }
+
+    fn num_qubits(&self) -> u32 {
+        self.z.len() as u32
+    }
+
+    fn num_clbits(&self) -> u32 {
+        0
+    }
+
+    fn num_params(&self) -> u32 {
+        1
+    }
+
+    fn directive(&self) -> bool {
+        false
+    }
+}
+
+
 impl PartialEq for PauliProductRotation {
     fn eq(&self, other: &Self) -> bool {
         self.x == other.x

crates/circuit/src/packed_instruction.rs

@@ -23,6 +23,7 @@ use crate::operations::{
     PyOperationTypes, PythonOperation, StandardGate, StandardInstruction, UnitaryGate,
 };
 use crate::{Block, Clbit, Qubit};
+
 use hashbrown::HashMap;
 use nalgebra::Matrix2;
 use ndarray::{Array2, CowArray, Ix2};