prism-q 0.35.0

//! Gate definitions and matrix representations.
//!
//! `Gate` is 16 bytes with large payloads boxed; matrices are returned as stack arrays.

use crate::sim::unified_pauli::PauliAxis;
use num_complex::Complex64;
use smallvec::SmallVec;
use std::f64::consts::{FRAC_1_SQRT_2, PI};
use std::fmt;

/// Zero threshold on an entry's `norm_sqr` in `preserves_sparsity`.
const NEAR_ZERO_NORM_SQ: f64 = 1e-24;

/// Per-entry norm tolerance for identity, diagonal, and controlled-phase checks.
/// `RECOGNIZE_EPS` is held at this value so the matrix recognizers cannot accept a
/// gate the fusion pass's identity drop would reject.
pub(crate) const IDENTITY_EPS: f64 = 1e-12;

/// Quantum gate identifier.
///
/// Payloads larger than an `f64` are boxed to hold the enum at 16 bytes.
#[derive(Debug, Clone, PartialEq)]
#[non_exhaustive]
pub enum Gate {
    /// Identity.
    Id,
    /// Pauli-X (bit flip).
    X,
    /// Pauli-Y.
    Y,
    /// Pauli-Z (phase flip).
    Z,
    /// Hadamard.
    H,
    /// S gate (√Z).
    S,
    /// S† gate.
    Sdg,
    /// T gate (π/8).
    T,
    /// T† gate.
    Tdg,
    /// √X gate.
    SX,
    /// √X† gate.
    SXdg,

    /// Rotation about X-axis by angle (radians).
    Rx(f64),
    /// Rotation about Y-axis by angle (radians).
    Ry(f64),
    /// Rotation about Z-axis by angle (radians).
    Rz(f64),
    /// Phase gate `[[1,0],[0,e^{iθ}]]`.
    P(f64),

    /// ZZ rotation: diag(e^{-iθ/2}, e^{iθ/2}, e^{iθ/2}, e^{-iθ/2}).
    /// Qubit order: [q0, q1] (symmetric).
    Rzz(f64),

    /// Controlled-X (CNOT). Qubit order: [control, target].
    Cx,
    /// Controlled-Z. Qubit order: [q0, q1] (symmetric).
    Cz,
    /// SWAP. Qubit order: [q0, q1] (symmetric).
    Swap,

    /// Controlled-unitary on the target when the control is |1⟩. Qubit order:
    /// [control, target].
    Cu(Box<[[Complex64; 2]; 2]>),

    /// Multi-controlled unitary on the target when every control is |1⟩. Qubit
    /// order: `[ctrl_0, ctrl_1, ..., ctrl_{k-1}, target]`.
    Mcu(Box<McuData>),

    /// Pre-fused single-qubit unitary (product of consecutive gates on the same target).
    Fused(Box<[[Complex64; 2]; 2]>),

    /// Controlled-phase gates sharing a control, applied in one pass. Targets:
    /// `[control]`.
    BatchPhase(Box<BatchPhaseData>),

    /// Rzz gates applied in one pass.
    BatchRzz(Box<BatchRzzData>),

    /// A run of diagonal 1q and 2q gates applied in one sweep through a phase LUT.
    DiagonalBatch(Box<DiagonalBatchData>),

    /// Single-qubit gates on distinct qubits, applied in one tiled pass.
    MultiFused(Box<MultiFusedData>),

    /// Pre-fused two-qubit unitary, with adjacent single-qubit gates absorbed.
    Fused2q(Box<[[Complex64; 4]; 4]>),

    /// Two-qubit gates applied in one tiled pass.
    Multi2q(Box<Multi2qData>),

    /// Quantum Fourier Transform on `start..start+num`, `q[start]` the least
    /// significant bit unless `big_endian`.
    ///
    /// The flags select the gate sequence [`crate::circuit::qft_textbook_steps`]
    /// expands to: `swaps` appends the bit-reversal swaps, `big_endian` mirrors
    /// the qubit labels so the first H lands on `q[start]`, and `inverse` is the
    /// exact adjoint of that sequence. The CPU statevector runs blocks starting at
    /// `q[0]` as an FFT over each `2^num` amplitude chunk. Other blocks and other
    /// backends expand to the textbook gates before execution.
    QftBlock {
        start: u8,
        num: u8,
        inverse: bool,
        swaps: bool,
        big_endian: bool,
    },

    /// Multi-qubit Pauli rotation `exp(-i θ P / 2)`. The `PauliRotData` letter
    /// at index `i` acts on the instruction's `targets[i]`.
    ///
    /// Built by [`crate::Circuit::add_pauli_rotation`], which lowers weight-1
    /// strings to `Rx`/`Ry`/`Rz` and two-qubit `ZZ` to `Rzz`, so only the
    /// residual multi-qubit strings construct this variant. The CPU
    /// statevector applies it in one pass; backends without the kernel receive
    /// the CNOT-ladder lowering from `circuit::expand_pauli_rotations`.
    PauliRot(Box<PauliRotData>),

    /// Pauli rotations applied in order in one cache-tiled pass. Targets list every
    /// qubit the strings touch, ascending.
    MultiPauliRot(Box<MultiPauliRotData>),

    /// Dense unitary on three or more qubits, a row-major `2^k x 2^k` matrix
    /// with `targets[0]` the most significant bit of both indices, the packing
    /// [`Gate::matrix_4x4`] uses.
    ///
    /// Built only by [`Gate::unitary`], which lowers one and two qubit matrices
    /// and the structured wider ones into the existing variants. The CPU
    /// statevector applies it in one gather-scatter pass; the factored, density
    /// matrix, MPS, and tensor network backends apply it through their dense
    /// paths, and every other backend declines it by name.
    Unitary(Box<UnitaryData>),
}

/// Widest matrix [`Gate::unitary`] accepts, in qubits.
pub const MAX_UNITARY_QUBITS: usize = 4;

/// Tolerance on `max |(U^dagger U - I)_ij|` for a matrix [`Gate::unitary`]
/// accepts, matching the norm tolerance on a caller-supplied start state.
const UNITARY_EPS: f64 = 1e-9;

/// Data for a dense multi-qubit unitary.
///
/// Fields are private: the kernels read `num_qubits` and the matrix length as
/// one fact, and a disagreement between them is an out-of-bounds read rather
/// than a wrong answer. [`Gate::unitary`] is the only constructor.
#[derive(Debug, Clone, PartialEq)]
pub struct UnitaryData {
    mat: Vec<Complex64>,
    num_qubits: u8,
}

impl UnitaryData {
    /// Row-major `2^k x 2^k` entries, `targets[0]` the most significant index
    /// bit.
    pub fn matrix(&self) -> &[Complex64] {
        &self.mat
    }

    pub fn num_qubits(&self) -> usize {
        self.num_qubits as usize
    }

    fn adjoint(&self) -> UnitaryData {
        let dim = 1usize << self.num_qubits;
        let mut mat = vec![Complex64::new(0.0, 0.0); dim * dim];
        for row in 0..dim {
            for column in 0..dim {
                mat[row * dim + column] = self.mat[column * dim + row].conj();
            }
        }
        UnitaryData {
            mat,
            num_qubits: self.num_qubits,
        }
    }

    pub(crate) fn conjugated(&self) -> UnitaryData {
        UnitaryData {
            mat: self.mat.iter().map(Complex64::conj).collect(),
            num_qubits: self.num_qubits,
        }
    }
}

/// Reject a matrix [`Gate::unitary`] cannot carry: a width outside
/// `1..=MAX_UNITARY_QUBITS`, a length other than `4^k`, a non-finite entry,
/// or a deviation of `U^dagger U` from the identity past `UNITARY_EPS`.
fn check_unitary_matrix(mat: &[Complex64], num_qubits: usize) -> crate::error::Result<()> {
    let invalid = |message: String| crate::error::PrismError::InvalidParameter { message };
    if num_qubits == 0 || num_qubits > MAX_UNITARY_QUBITS {
        return Err(invalid(format!(
            "a unitary needs between 1 and {MAX_UNITARY_QUBITS} qubits, got {num_qubits}"
        )));
    }
    let dim = 1usize << num_qubits;
    if mat.len() != dim * dim {
        return Err(invalid(format!(
            "a {num_qubits}-qubit unitary needs {} entries ({dim}x{dim}), got {}",
            dim * dim,
            mat.len()
        )));
    }
    if let Some(index) = mat.iter().position(|entry| !entry.is_finite()) {
        return Err(invalid(format!(
            "unitary entry ({}, {}) is not finite",
            index / dim,
            index % dim
        )));
    }
    let mut worst = (0.0f64, 0usize, 0usize);
    for i in 0..dim {
        for j in 0..dim {
            let mut dot = Complex64::new(0.0, 0.0);
            for row in 0..dim {
                dot += mat[row * dim + i].conj() * mat[row * dim + j];
            }
            if i == j {
                dot -= Complex64::new(1.0, 0.0);
            }
            let deviation = dot.norm();
            if deviation > worst.0 {
                worst = (deviation, i, j);
            }
        }
    }
    if worst.0 > UNITARY_EPS {
        return Err(invalid(format!(
            "matrix is not unitary: U^dagger U differs from the identity by {:.3e} at ({}, {}), \
             above the {UNITARY_EPS:e} tolerance",
            worst.0, worst.1, worst.2
        )));
    }
    Ok(())
}

/// The trailing 2x2 block of `mat` when everything outside it is the identity,
/// which is the `Mcu` shape with every other qubit a control.
fn controlled_block(mat: &[Complex64], num_qubits: usize) -> Option<[[Complex64; 2]; 2]> {
    let dim = 1usize << num_qubits;
    let split = dim - 2;
    for row in 0..dim {
        for column in 0..dim {
            if row >= split && column >= split {
                continue;
            }
            let expected = if row == column { 1.0 } else { 0.0 };
            if (mat[row * dim + column] - Complex64::new(expected, 0.0)).norm() >= IDENTITY_EPS {
                return None;
            }
        }
    }
    Some([
        [mat[split * dim + split], mat[split * dim + split + 1]],
        [
            mat[(split + 1) * dim + split],
            mat[(split + 1) * dim + split + 1],
        ],
    ])
}

/// A diagonal `2^k x 2^k` matrix over `targets` as a [`Gate::DiagonalBatch`],
/// or `None` when it is not diagonal or its phases do not factor into the one
/// and two body terms [`DiagEntry`] spells (a `CCZ` does not).
///
/// The factorization is the multiplicative Moebius inversion of the diagonal
/// over the subsets of size at most two, and the product is checked against
/// every entry within `UNITARY_EPS` before the batch is trusted. The global
/// phase rides on the first target's `d0`, so the batch reproduces the matrix
/// rather than the matrix up to phase.
pub(crate) fn diagonal_batch(mat: &[Complex64], targets: &[usize]) -> Option<Gate> {
    let k = targets.len();
    let dim = 1usize << k;
    for row in 0..dim {
        for column in 0..dim {
            if row != column && mat[row * dim + column].norm() >= IDENTITY_EPS {
                return None;
            }
        }
    }
    let diag = |index: usize| mat[index * dim + index];
    let unit = |i: usize| 1usize << (k - 1 - i);
    let one = Complex64::new(1.0, 0.0);
    let g0 = diag(0);
    let g1: Vec<Complex64> = (0..k).map(|i| diag(unit(i)) / g0).collect();
    let mut g2 = vec![one; k * k];
    for i in 0..k {
        for j in i + 1..k {
            g2[i * k + j] = diag(unit(i) | unit(j)) * g0 / (diag(unit(i)) * diag(unit(j)));
        }
    }
    for index in 0..dim {
        let mut product = g0;
        for i in 0..k {
            if index & unit(i) == 0 {
                continue;
            }
            product *= g1[i];
            for j in i + 1..k {
                if index & unit(j) != 0 {
                    product *= g2[i * k + j];
                }
            }
        }
        if (product - diag(index)).norm() > UNITARY_EPS {
            return None;
        }
    }
    let mut entries = Vec::with_capacity(k + k * (k - 1) / 2);
    entries.push(DiagEntry::Phase1q {
        qubit: targets[0],
        d0: g0,
        d1: g0 * g1[0],
    });
    for i in 1..k {
        entries.push(DiagEntry::Phase1q {
            qubit: targets[i],
            d0: one,
            d1: g1[i],
        });
    }
    for i in 0..k {
        for j in i + 1..k {
            let phase = g2[i * k + j];
            if (phase - one).norm() >= IDENTITY_EPS {
                entries.push(DiagEntry::Phase2q {
                    q0: targets[i],
                    q1: targets[j],
                    phase,
                });
            }
        }
    }
    Some(Gate::DiagonalBatch(Box::new(DiagonalBatchData { entries })))
}

/// Analytic differentiation generator for a parametric gate, returned by
/// [`Gate::pauli_generator`].
///
/// Each rotation variant is `exp(-i θ/2 G)` with the named Pauli
/// generator `G` acting on the gate's `targets` (in order); `Phase` is the
/// non-Pauli phase gate `diag(1, e^{iθ})` whose generator is the projector
/// `|1⟩⟨1|` on `targets[0]`. The borrow is of the gate the generator was read
/// from, so a caller holding one cannot pair it with another gate's letters.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum GeneratorKind<'a> {
    /// `Rx(θ)`: generator `X` on `targets[0]`.
    RotX,
    /// `Ry(θ)`: generator `Y` on `targets[0]`.
    RotY,
    /// `Rz(θ)`: generator `Z` on `targets[0]`.
    RotZ,
    /// `Rzz(θ)`: generator `Z⊗Z` on `targets[0]` and `targets[1]`.
    RotZz,
    /// `P(θ)`: generator is the projector `|1⟩⟨1|` on `targets[0]`.
    Phase,
    /// `PauliRot(θ)`: generator is the Pauli string whose letter at index `i`
    /// acts on `targets[i]`.
    RotPauli(&'a [PauliAxis]),
}

/// Data for a multi-qubit Pauli rotation `exp(-i θ P / 2)`.
///
/// `axes[i]` is the letter of `P` on the owning instruction's `targets[i]`.
/// Fields are private so every construction path runs the recognizing
/// lowering in [`crate::Circuit::add_pauli_rotation`].
#[derive(Debug, Clone, PartialEq)]
pub struct PauliRotData {
    pub(crate) theta: f64,
    pub(crate) axes: Vec<PauliAxis>,
}

impl PauliRotData {
    /// Rotation angle in radians.
    pub fn theta(&self) -> f64 {
        self.theta
    }

    /// Overwrite the rotation angle, leaving the Pauli string alone. The
    /// recognizing lowering runs on the string, so rebinding an angle cannot
    /// invalidate it.
    pub fn set_theta(&mut self, theta: f64) {
        self.theta = theta;
    }

    /// Pauli letters aligned with the instruction targets.
    pub fn axes(&self) -> &[PauliAxis] {
        &self.axes
    }
}

/// Pack the Pauli string of a rotation into `(xmask, zmask, num_y)`, where
/// `xmask` covers X and Y letters and `zmask` covers Z and Y, and `axes[i]`
/// acts on `targets[i]`.
///
/// The native kernel and the adjoint gradient both derive their masks here, so
/// a change to the letter convention cannot move one without the other.
#[inline]
pub(crate) fn pauli_rot_masks(targets: &[usize], axes: &[PauliAxis]) -> (usize, usize, u32) {
    let mut xmask = 0usize;
    let mut zmask = 0usize;
    let mut num_y = 0u32;
    for (&q, axis) in targets.iter().zip(axes) {
        let bit = 1usize << q;
        match axis {
            PauliAxis::X => xmask |= bit,
            PauliAxis::Z => zmask |= bit,
            PauliAxis::Y => {
                xmask |= bit;
                zmask |= bit;
                num_y += 1;
            }
        }
    }
    (xmask, zmask, num_y)
}

/// Data for a `MultiPauliRot` batch: `(xmask, zmask, theta)` per rotation in circuit
/// order. Bit `q` of `xmask` marks an X or Y letter on qubit `q` and bit `q` of `zmask`
/// a Z or Y, so the Y count is `(xmask & zmask).count_ones()`.
#[derive(Debug, Clone, PartialEq)]
pub struct MultiPauliRotData {
    pub(crate) rotations: Vec<(usize, usize, f64)>,
}

impl MultiPauliRotData {
    pub fn rotations(&self) -> &[(usize, usize, f64)] {
        &self.rotations
    }

    /// Every qubit any rotation touches, as a bit mask.
    pub(crate) fn support(&self) -> usize {
        self.rotations.iter().fold(0, |acc, &(x, z, _)| acc | x | z)
    }

    /// Call `f` with each rotation as an angle, ascending targets, and their letters.
    #[cfg(feature = "gpu")]
    pub(crate) fn for_each_rotation(&self, mut f: impl FnMut(f64, &[usize], &[PauliAxis])) {
        let mut targets: SmallVec<[usize; 8]> = SmallVec::new();
        let mut axes: SmallVec<[PauliAxis; 8]> = SmallVec::new();
        for &(xmask, zmask, theta) in &self.rotations {
            targets.clear();
            axes.clear();
            let mut rest = xmask | zmask;
            while rest != 0 {
                let q = rest.trailing_zeros() as usize;
                rest &= rest - 1;
                targets.push(q);
                axes.push(match (xmask >> q & 1, zmask >> q & 1) {
                    (1, 1) => PauliAxis::Y,
                    (1, _) => PauliAxis::X,
                    _ => PauliAxis::Z,
                });
            }
            f(theta, &targets, &axes);
        }
    }
}

/// Data for a multi-controlled unitary gate.
#[derive(Debug, Clone, PartialEq)]
pub struct McuData {
    pub(crate) mat: [[Complex64; 2]; 2],
    pub(crate) num_controls: u8,
}

impl McuData {
    /// `mat` applied to the target when all `num_controls` controls are set.
    ///
    /// # Panics
    ///
    /// Panics if `num_controls` is below 2; [`Gate::mcu`] takes one control and
    /// returns the controlled 2x2 form instead.
    pub fn new(mat: [[Complex64; 2]; 2], num_controls: u8) -> Self {
        assert!(num_controls >= 2, "McuData needs at least 2 controls");
        Self { mat, num_controls }
    }

    /// 2x2 unitary applied to the target qubit.
    pub fn mat(&self) -> &[[Complex64; 2]; 2] {
        &self.mat
    }

    pub fn num_controls(&self) -> u8 {
        self.num_controls
    }
}

/// Data for a batched controlled-phase gate: `(target_qubit, phase)` entries, with
/// the shared control in the instruction's `targets[0]`.
#[derive(Debug, Clone, PartialEq)]
pub struct BatchPhaseData {
    pub(crate) phases: SmallVec<[(usize, Complex64); 8]>,
}

impl BatchPhaseData {
    /// # Panics
    ///
    /// Panics if there are more than [`Self::MAX_PHASES`] entries.
    pub fn new(phases: impl IntoIterator<Item = (usize, Complex64)>) -> Self {
        let phases: SmallVec<[(usize, Complex64); 8]> = phases.into_iter().collect();
        assert!(
            phases.len() <= Self::MAX_PHASES,
            "BatchPhaseData holds at most {} entries",
            Self::MAX_PHASES
        );
        Self { phases }
    }

    pub fn phases(&self) -> &[(usize, Complex64)] {
        &self.phases
    }

    /// Entry cap the fusion pass splits on, matching the kernel group tables.
    /// Targets are distinct within one payload, so a repeated `(control, target)`
    /// pair folds into a single entry rather than adding one.
    pub const MAX_PHASES: usize = 40;
}

/// Data for batched ZZ rotations: `(qubit_0, qubit_1, theta)` entries, with every
/// qubit also in the instruction's `targets`.
#[derive(Debug, Clone, PartialEq)]
pub struct BatchRzzData {
    pub(crate) edges: Vec<(usize, usize, f64)>,
}

impl BatchRzzData {
    /// # Panics
    ///
    /// Panics if there are more than [`Self::MAX_EDGES`] edges.
    pub fn new(edges: Vec<(usize, usize, f64)>) -> Self {
        assert!(
            edges.len() <= Self::MAX_EDGES,
            "BatchRzzData holds at most {} edges",
            Self::MAX_EDGES
        );
        Self { edges }
    }

    pub fn edges(&self) -> &[(usize, usize, f64)] {
        &self.edges
    }

    /// Edge cap the fusion pass splits on, matching the kernel group tables.
    pub const MAX_EDGES: usize = 32;
}

/// An individual diagonal phase contribution in a [`DiagonalBatchData`].
#[derive(Debug, Clone, PartialEq)]
pub enum DiagEntry {
    /// Diagonal on a single qubit: `state[i] *= d0` when bit 0, `*= d1` when bit 1.
    Phase1q {
        qubit: usize,
        d0: Complex64,
        d1: Complex64,
    },
    /// Phase on a qubit pair: `state[i] *= phase` when both bits are set (CZ/CPhase).
    Phase2q {
        q0: usize,
        q1: usize,
        phase: Complex64,
    },
    /// Parity-dependent phase (Rzz): `state[i] *= same` when parity is even,
    /// `state[i] *= diff` when parity is odd.
    Parity2q {
        q0: usize,
        q1: usize,
        same: Complex64,
        diff: Complex64,
    },
}

impl DiagEntry {
    /// Return the qubit and dense 2×2 matrix for a [`DiagEntry::Phase1q`]
    /// entry, or `None` for two-qubit entries.
    pub(crate) fn as_1q_matrix(&self) -> Option<(usize, [[Complex64; 2]; 2])> {
        match *self {
            DiagEntry::Phase1q { qubit, d0, d1 } => {
                let z = Complex64::new(0.0, 0.0);
                Some((qubit, [[d0, z], [z, d1]]))
            }
            _ => None,
        }
    }

    /// Return the qubit pair and dense 4×4 matrix for a [`DiagEntry::Phase2q`]
    /// or [`DiagEntry::Parity2q`] entry, or `None` for single-qubit entries.
    pub(crate) fn as_2q_matrix(&self) -> Option<(usize, usize, [[Complex64; 4]; 4])> {
        let z = Complex64::new(0.0, 0.0);
        let one = Complex64::new(1.0, 0.0);
        match *self {
            DiagEntry::Phase2q { q0, q1, phase } => Some((
                q0,
                q1,
                [
                    [one, z, z, z],
                    [z, one, z, z],
                    [z, z, one, z],
                    [z, z, z, phase],
                ],
            )),
            DiagEntry::Parity2q {
                q0, q1, same, diff, ..
            } => Some((
                q0,
                q1,
                [
                    [same, z, z, z],
                    [z, diff, z, z],
                    [z, z, diff, z],
                    [z, z, z, same],
                ],
            )),
            _ => None,
        }
    }
}

/// Data for a batched diagonal gate; the kernel builds a phase LUT indexed by the
/// bits of the distinct qubits in `entries`.
#[derive(Debug, Clone, PartialEq)]
pub struct DiagonalBatchData {
    pub(crate) entries: Vec<DiagEntry>,
}

impl DiagonalBatchData {
    pub fn new(entries: Vec<DiagEntry>) -> Self {
        Self { entries }
    }

    pub fn entries(&self) -> &[DiagEntry] {
        &self.entries
    }

    /// Widest qubit group one kernel lookup table covers. Qubits joined by a two-qubit
    /// entry share a group, so the fusion pass closes a batch before a component outgrows it.
    pub(crate) const MAX_QUBITS_PER_GROUP: usize = 10;

    /// Most lookup tables one kernel sweep combines.
    pub(crate) const MAX_GROUPS: usize = 4;
}

/// Data for a `MultiFused` batch of `(target_qubit, 2×2 matrix)` entries.
#[derive(Debug, Clone, PartialEq)]
pub struct MultiFusedData {
    pub(crate) gates: Vec<(usize, [[Complex64; 2]; 2])>,
    pub(crate) all_diagonal: bool,
}

impl MultiFusedData {
    /// Batch `(target_qubit, 2×2 matrix)` entries, computing the diagonal flag the
    /// kernels branch on.
    pub fn new(gates: Vec<(usize, [[Complex64; 2]; 2])>) -> Self {
        let all_diagonal = gates.iter().all(|(_, m)| is_diagonal_2x2(m));
        Self {
            gates,
            all_diagonal,
        }
    }

    pub fn gates(&self) -> &[(usize, [[Complex64; 2]; 2])] {
        &self.gates
    }

    /// Whether every matrix is diagonal, so the batch applies as a phase pass.
    pub fn all_diagonal(&self) -> bool {
        self.all_diagonal
    }
}

/// Data for a `Multi2q` batch of `(q0, q1, 4×4 matrix)` entries in circuit order.
#[derive(Debug, Clone, PartialEq)]
pub struct Multi2qData {
    pub(crate) gates: Vec<(usize, usize, [[Complex64; 4]; 4])>,
}

impl Multi2qData {
    pub fn new(gates: Vec<(usize, usize, [[Complex64; 4]; 4])>) -> Self {
        Self { gates }
    }

    pub fn gates(&self) -> &[(usize, usize, [[Complex64; 4]; 4])] {
        &self.gates
    }
}

/// A `Multi2q` batch runs inside one tile of `2^tile_bits` amplitudes: the
/// lowest bits of the index plus up to `tile_bits - MULTI_2Q_LOW_BITS` gathered
/// high qubits, so a batch is tileable when its gates touch at most that many
/// distinct qubits at or above [`multi_2q_low_bits`]. The tile size is per
/// platform, see [`multi_2q_tile_bits`]; [`MULTI_2Q_HIGH_BUDGET`] is the
/// widest budget any platform reaches and sizes the inline qubit lists.
pub(crate) const MULTI_2Q_MAX_TILE_BITS: usize = 18;
/// Smallest tile the knob accepts, 128 KB: half the derived floor, so a sweep can
/// try it.
pub(crate) const MULTI_2Q_MIN_TILE_BITS: usize = 13;
/// Fewest low bits a tile keeps contiguous: 64-amplitude runs, 1 KB each.
pub(crate) const MULTI_2Q_MIN_LOW_BITS: usize = 6;
/// Most distinct high qubits one `Multi2q` batch may span on any platform.
pub(crate) const MULTI_2Q_HIGH_BUDGET: usize = MULTI_2Q_MAX_TILE_BITS - MULTI_2Q_MIN_LOW_BITS;

/// Qubits from this index up count against a batch's high-qubit budget: the
/// tile keeps `2^low_bits` amplitudes contiguous per run and gathers the rest.
/// Default 6, the shortest run, since every bit taken from the run is one
/// more high qubit a pass over the state can serve: at 24 qubits a quantum
/// volume circuit went from 88 passes to 54 and ran 10% faster against 10.
/// `PRISM_MULTI_2Q_LOW_BITS` overrides it, clamped to
/// `MULTI_2Q_MIN_LOW_BITS..tile_bits`.
pub(crate) fn multi_2q_low_bits() -> usize {
    static CACHED: std::sync::OnceLock<usize> = std::sync::OnceLock::new();
    *CACHED.get_or_init(|| {
        crate::env_knobs::usize_override("PRISM_MULTI_2Q_LOW_BITS", 1)
            .unwrap_or(6)
            .clamp(MULTI_2Q_MIN_LOW_BITS, multi_2q_tile_bits() - 1)
    })
}

/// Log2 of the `Multi2q` tile in amplitudes: the tile budget from the cache topology
/// (2^14, 256 KB, on the x86 cores measured), or `PRISM_MULTI_2Q_TILE_BITS` clamped
/// to `MULTI_2Q_MIN_TILE_BITS..=MULTI_2Q_MAX_TILE_BITS`.
///
/// The tile has to sit in the core's private cache: on an i7-6700K (256 KB
/// of L2 per core, 8 MB of L3 shared by four) a 1 MB tile ran quantum volume
/// at 22 and 24 qubits 2x slower and a 2 MB tile 3x, while 512 KB read flat.
/// A wider tile takes more of a random-pair layer per pass over the state,
/// so on a core with a larger private cache the knob is worth a sweep.
pub(crate) fn multi_2q_tile_bits() -> usize {
    static CACHED: std::sync::OnceLock<usize> = std::sync::OnceLock::new();
    *CACHED.get_or_init(|| {
        crate::env_knobs::usize_override("PRISM_MULTI_2Q_TILE_BITS", 1)
            .unwrap_or_else(|| {
                (crate::backend::cache::tile_budget_bytes() / size_of::<Complex64>()).ilog2()
                    as usize
            })
            .clamp(MULTI_2Q_MIN_TILE_BITS, MULTI_2Q_MAX_TILE_BITS)
    })
}

/// [`multi_2q_tile_bits`] for a `num_qubits` statevector: cut, down to
/// [`MULTI_2Q_MIN_TILE_BITS`], so a parallel pass gathers at least
/// `cache::min_parallel_tiles` subcubes instead of leaving workers idle at the widths
/// just past the parallel threshold.
#[cfg_attr(not(feature = "parallel"), allow(unused_variables))]
pub(crate) fn multi_2q_tile_bits_for(num_qubits: usize) -> usize {
    let bits = multi_2q_tile_bits();
    #[cfg(feature = "parallel")]
    if num_qubits >= crate::backend::PARALLEL_THRESHOLD_QUBITS {
        let pool_bits = num_qubits - crate::backend::cache::min_parallel_tiles().ilog2() as usize;
        return bits.min(pool_bits).max(MULTI_2Q_MIN_TILE_BITS);
    }
    bits
}

/// Most distinct high qubits one `Multi2q` batch may span on a `num_qubits` state.
pub(crate) fn multi_2q_high_budget_for(num_qubits: usize) -> usize {
    multi_2q_tile_bits_for(num_qubits) - multi_2q_low_bits()
}

/// The high qubits `high` grows to if the pair joins a batch on a `num_qubits`
/// state, or `None` when the pair would take it past [`multi_2q_high_budget_for`].
pub(crate) fn multi_2q_join(
    high: &[usize],
    q0: usize,
    q1: usize,
    num_qubits: usize,
) -> Option<smallvec::SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]>> {
    let budget = multi_2q_high_budget_for(num_qubits);
    let low_bits = multi_2q_low_bits();
    let mut joined: smallvec::SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]> =
        high.iter().copied().collect();
    for q in [q0, q1] {
        if q >= low_bits && !joined.contains(&q) {
            if joined.len() == budget {
                return None;
            }
            joined.push(q);
        }
    }
    Some(joined)
}

/// [`multi_2q_join`] for the X and Y qubits of a Pauli string, which a
/// `MultiPauliRot` batch needs inside its tile.
pub(crate) fn pauli_rot_join(
    high: &[usize],
    qubits: impl IntoIterator<Item = usize>,
    num_qubits: usize,
) -> Option<smallvec::SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]>> {
    let budget = multi_2q_high_budget_for(num_qubits);
    let low_bits = multi_2q_low_bits();
    let mut joined: smallvec::SmallVec<[usize; MULTI_2Q_HIGH_BUDGET]> =
        high.iter().copied().collect();
    for q in qubits {
        if q >= low_bits && !joined.contains(&q) {
            if joined.len() == budget {
                return None;
            }
            joined.push(q);
        }
    }
    Some(joined)
}

/// Controlled `mat` as a 4x4 with the control on the high bit of the basis index.
#[inline]
pub(crate) fn cu_matrix_4x4(mat: &[[Complex64; 2]; 2]) -> [[Complex64; 4]; 4] {
    let z = Complex64::new(0.0, 0.0);
    let o = Complex64::new(1.0, 0.0);
    [
        [o, z, z, z],
        [z, o, z, z],
        [z, z, mat[0][0], mat[0][1]],
        [z, z, mat[1][0], mat[1][1]],
    ]
}

/// Kronecker product of two 2×2 matrices: A ⊗ B → 4×4.
///
/// Result indices: `(i*2+j, k*2+l) = A[i][k] * B[j][l]`
/// where i,k index A (targets\[0\]) and j,l index B (targets\[1\]).
#[inline]
pub(crate) fn kron_2x2(a: &[[Complex64; 2]; 2], b: &[[Complex64; 2]; 2]) -> [[Complex64; 4]; 4] {
    let mut result = [[Complex64::new(0.0, 0.0); 4]; 4];
    for i in 0..2 {
        for k in 0..2 {
            let aik = a[i][k];
            for j in 0..2 {
                for l in 0..2 {
                    result[i * 2 + j][k * 2 + l] = aik * b[j][l];
                }
            }
        }
    }
    result
}

#[inline]
pub(crate) fn mat_mul_4x4(a: &[[Complex64; 4]; 4], b: &[[Complex64; 4]; 4]) -> [[Complex64; 4]; 4] {
    let zero = Complex64::new(0.0, 0.0);
    let mut result = [[zero; 4]; 4];
    for i in 0..4 {
        for j in 0..4 {
            let mut sum = zero;
            for k in 0..4 {
                sum += a[i][k] * b[k][j];
            }
            result[i][j] = sum;
        }
    }
    result
}

fn adjoint_4x4(m: &[[Complex64; 4]; 4]) -> [[Complex64; 4]; 4] {
    let mut result = [[Complex64::new(0.0, 0.0); 4]; 4];
    for i in 0..4 {
        for j in 0..4 {
            result[i][j] = m[j][i].conj();
        }
    }
    result
}

fn adjoint_2x2(m: &[[Complex64; 2]; 2]) -> [[Complex64; 2]; 2] {
    [
        [m[0][0].conj(), m[1][0].conj()],
        [m[0][1].conj(), m[1][1].conj()],
    ]
}

#[inline]
fn count_unique_qubits<I: IntoIterator<Item = usize>>(iter: I) -> usize {
    let mut seen: SmallVec<[usize; 8]> = SmallVec::new();
    for q in iter {
        if !seen.contains(&q) {
            seen.push(q);
        }
    }
    seen.len()
}

#[inline]
fn push_unique_qubit(seen: &mut SmallVec<[usize; 8]>, qubit: usize) {
    if !seen.contains(&qubit) {
        seen.push(qubit);
    }
}

/// Combined phase factor `entries` apply to the basis state `index`.
#[inline]
pub(crate) fn diag_entries_phase(index: usize, entries: &[DiagEntry]) -> Complex64 {
    let mut combined = Complex64::new(1.0, 0.0);
    for entry in entries {
        match entry {
            DiagEntry::Phase1q { qubit, d0, d1 } => {
                combined *= if (index >> qubit) & 1 == 1 { *d1 } else { *d0 };
            }
            DiagEntry::Phase2q { q0, q1, phase } => {
                if (index >> q0) & 1 == 1 && (index >> q1) & 1 == 1 {
                    combined *= phase;
                }
            }
            DiagEntry::Parity2q { q0, q1, same, diff } => {
                let parity = ((index >> q0) ^ (index >> q1)) & 1;
                combined *= if parity == 0 { *same } else { *diff };
            }
        }
    }
    combined
}

#[inline]
fn count_unique_diag_qubits(entries: &[DiagEntry]) -> usize {
    let mut seen: SmallVec<[usize; 8]> = SmallVec::new();
    for entry in entries {
        match entry {
            DiagEntry::Phase1q { qubit, .. } => push_unique_qubit(&mut seen, *qubit),
            DiagEntry::Phase2q { q0, q1, .. } | DiagEntry::Parity2q { q0, q1, .. } => {
                push_unique_qubit(&mut seen, *q0);
                push_unique_qubit(&mut seen, *q1);
            }
        }
    }
    seen.len()
}

/// An eigenvalue angle within this of `-pi` is read at `+pi` instead.
///
/// Both spell the same eigenvalue, and the principal branch is closed at `+pi`,
/// so `pow(0.5) @ p(-pi)` gives `S` rather than `S-dagger`: a half turn written
/// the other way round is still a half turn.
const PRINCIPAL_BRANCH_SNAP: f64 = 1e-9;

/// `angle` on the principal branch `(-pi, pi]`.
pub(crate) fn principal_angle(angle: f64) -> f64 {
    let wrapped = angle.rem_euclid(std::f64::consts::TAU);
    let folded = if wrapped > std::f64::consts::PI {
        wrapped - std::f64::consts::TAU
    } else {
        wrapped
    };
    if folded <= -std::f64::consts::PI + PRINCIPAL_BRANCH_SNAP {
        std::f64::consts::PI
    } else {
        folded
    }
}

/// `e^(-i phase) U`, which has determinant 1 when `phase` is half the argument
/// of `det U`.
fn special_unitary(m: &[[Complex64; 2]; 2], phase: f64) -> [[Complex64; 2]; 2] {
    let factor = Complex64::new(0.0, -phase).exp();
    [
        [m[0][0] * factor, m[0][1] * factor],
        [m[1][0] * factor, m[1][1] * factor],
    ]
}

#[inline]
pub(crate) fn mat_mul_2x2(a: &[[Complex64; 2]; 2], b: &[[Complex64; 2]; 2]) -> [[Complex64; 2]; 2] {
    [
        [
            a[0][0] * b[0][0] + a[0][1] * b[1][0],
            a[0][0] * b[0][1] + a[0][1] * b[1][1],
        ],
        [
            a[1][0] * b[0][0] + a[1][1] * b[1][0],
            a[1][0] * b[0][1] + a[1][1] * b[1][1],
        ],
    ]
}

impl Gate {
    /// Number of qubits this gate acts on.
    #[inline]
    pub fn num_qubits(&self) -> usize {
        match self {
            Gate::Rzz(_) | Gate::Cx | Gate::Cz | Gate::Swap | Gate::Cu(_) | Gate::Fused2q(_) => 2,
            Gate::Mcu(data) => data.num_controls as usize + 1,
            Gate::BatchPhase(data) => 1 + data.phases.len(),
            Gate::QftBlock { num, .. } => *num as usize,
            Gate::BatchRzz(data) => {
                count_unique_qubits(data.edges.iter().flat_map(|&(q0, q1, _)| [q0, q1]))
            }
            Gate::DiagonalBatch(data) => count_unique_diag_qubits(&data.entries),
            Gate::PauliRot(data) => data.axes.len(),
            Gate::MultiPauliRot(data) => data.support().count_ones() as usize,
            Gate::Unitary(data) => data.num_qubits(),
            Gate::MultiFused(data) => data.gates.len(),
            Gate::Multi2q(data) => {
                count_unique_qubits(data.gates.iter().flat_map(|&(q0, q1, _)| [q0, q1]))
            }
            _ => 1,
        }
    }

    /// Returns the 2×2 unitary matrix for single-qubit gates.
    ///
    /// # Panics
    /// Panics on any variant other than the named single-qubit gates, `Fused`,
    /// and a weight-1 `PauliRot`.
    #[inline]
    pub fn matrix_2x2(&self) -> [[Complex64; 2]; 2] {
        let zero = Complex64::new(0.0, 0.0);
        let one = Complex64::new(1.0, 0.0);
        let i = Complex64::new(0.0, 1.0);
        let neg_i = Complex64::new(0.0, -1.0);
        let h = Complex64::new(FRAC_1_SQRT_2, 0.0);

        match self {
            Gate::Id => [[one, zero], [zero, one]],
            Gate::X => [[zero, one], [one, zero]],
            Gate::Y => [[zero, neg_i], [i, zero]],
            Gate::Z => [[one, zero], [zero, -one]],
            Gate::H => [[h, h], [h, -h]],
            Gate::S => [[one, zero], [zero, i]],
            Gate::Sdg => [[one, zero], [zero, neg_i]],
            Gate::T => {
                let phase = Complex64::from_polar(1.0, PI / 4.0);
                [[one, zero], [zero, phase]]
            }
            Gate::Tdg => {
                let phase = Complex64::from_polar(1.0, -PI / 4.0);
                [[one, zero], [zero, phase]]
            }
            Gate::SX => {
                let half = Complex64::new(0.5, 0.0);
                let half_i = Complex64::new(0.0, 0.5);
                [
                    [half + half_i, half - half_i],
                    [half - half_i, half + half_i],
                ]
            }
            Gate::SXdg => {
                let half = Complex64::new(0.5, 0.0);
                let half_i = Complex64::new(0.0, 0.5);
                [
                    [half - half_i, half + half_i],
                    [half + half_i, half - half_i],
                ]
            }
            Gate::Rx(theta) => {
                let c = Complex64::new((theta / 2.0).cos(), 0.0);
                let s = Complex64::new(0.0, -(theta / 2.0).sin());
                [[c, s], [s, c]]
            }
            Gate::Ry(theta) => {
                let c = Complex64::new((theta / 2.0).cos(), 0.0);
                let s = Complex64::new((theta / 2.0).sin(), 0.0);
                [[c, -s], [s, c]]
            }
            Gate::Rz(theta) => {
                let e_neg = Complex64::from_polar(1.0, -theta / 2.0);
                let e_pos = Complex64::from_polar(1.0, theta / 2.0);
                [[e_neg, zero], [zero, e_pos]]
            }
            Gate::P(theta) => {
                let phase = Complex64::from_polar(1.0, *theta);
                [[one, zero], [zero, phase]]
            }
            Gate::Fused(mat) => **mat,
            Gate::PauliRot(data) if data.axes.len() == 1 => match data.axes[0] {
                PauliAxis::X => Gate::Rx(data.theta).matrix_2x2(),
                PauliAxis::Y => Gate::Ry(data.theta).matrix_2x2(),
                PauliAxis::Z => Gate::Rz(data.theta).matrix_2x2(),
            },
            Gate::Rzz(_)
            | Gate::Cx
            | Gate::Cz
            | Gate::Swap
            | Gate::Cu(_)
            | Gate::Mcu(_)
            | Gate::BatchPhase(_)
            | Gate::QftBlock { .. }
            | Gate::BatchRzz(_)
            | Gate::DiagonalBatch(_)
            | Gate::PauliRot(_)
            | Gate::MultiPauliRot(_)
            | Gate::Unitary(_)
            | Gate::MultiFused(_)
            | Gate::Fused2q(_)
            | Gate::Multi2q(_) => {
                panic!(
                    "matrix_2x2 called on {}-qubit gate `{}`; use dedicated backend routine",
                    self.num_qubits(),
                    self.name()
                )
            }
        }
    }

    /// Returns the 4×4 unitary matrix for two-qubit gates.
    ///
    /// Matrix indices follow the convention: row/col `i*2+j` where `i` indexes
    /// `targets[0]` and `j` indexes `targets[1]`.
    ///
    /// A [`Gate::PauliRot`] is admitted at weight 2 only, as
    /// `exp(-i theta/2 * P0 x P1)` with `P0` on `targets[0]`.
    ///
    /// # Panics
    /// Panics on gates other than `Cx`, `Cz`, `Swap`, `Cu`, `Fused2q`, `Rzz`,
    /// and a weight-2 `PauliRot`.
    pub fn matrix_4x4(&self) -> [[Complex64; 4]; 4] {
        let z = Complex64::new(0.0, 0.0);
        let o = Complex64::new(1.0, 0.0);
        let m = Complex64::new(-1.0, 0.0);
        match self {
            Gate::Rzz(theta) => {
                let ps = Complex64::from_polar(1.0, -theta / 2.0);
                let pd = Complex64::from_polar(1.0, theta / 2.0);
                [[ps, z, z, z], [z, pd, z, z], [z, z, pd, z], [z, z, z, ps]]
            }
            Gate::Cx => [[o, z, z, z], [z, o, z, z], [z, z, z, o], [z, z, o, z]],
            Gate::Cz => [[o, z, z, z], [z, o, z, z], [z, z, o, z], [z, z, z, m]],
            Gate::Swap => [[o, z, z, z], [z, z, o, z], [z, o, z, z], [z, z, z, o]],
            Gate::Cu(mat) => cu_matrix_4x4(mat),
            Gate::Fused2q(mat) => **mat,
            Gate::PauliRot(data) if data.axes.len() == 2 => {
                let half = data.theta / 2.0;
                let axis = |a: &PauliAxis| match a {
                    PauliAxis::X => Gate::X.matrix_2x2(),
                    PauliAxis::Y => Gate::Y.matrix_2x2(),
                    PauliAxis::Z => Gate::Z.matrix_2x2(),
                };
                let p = kron_2x2(&axis(&data.axes[0]), &axis(&data.axes[1]));
                let cos = Complex64::new(half.cos(), 0.0);
                let minus_i_sin = Complex64::new(0.0, -half.sin());
                let mut out = [[z; 4]; 4];
                for (r, row) in out.iter_mut().enumerate() {
                    for (c, e) in row.iter_mut().enumerate() {
                        *e = minus_i_sin * p[r][c];
                    }
                    row[r] += cos;
                }
                out
            }
            _ => panic!(
                "matrix_4x4 called on non-standard-2q gate `{}`",
                self.name()
            ),
        }
    }

    /// Dense matrix over the gate's own targets, row major with `targets[0]`
    /// the most significant bit of both indices.
    ///
    /// `None` for the batched and tiled variants, whose targets the fusion
    /// passes assign rather than a caller, and which no lowering here reaches.
    pub(crate) fn dense_matrix(&self) -> Option<Vec<Complex64>> {
        let zero = Complex64::new(0.0, 0.0);
        match self {
            Gate::Mcu(data) => {
                let dim = 1usize << self.num_qubits();
                let mut out = vec![zero; dim * dim];
                for index in 0..dim - 2 {
                    out[index * dim + index] = Complex64::new(1.0, 0.0);
                }
                for row in 0..2 {
                    for column in 0..2 {
                        out[(dim - 2 + row) * dim + (dim - 2 + column)] = data.mat[row][column];
                    }
                }
                Some(out)
            }
            Gate::PauliRot(data) => {
                let width = data.axes.len();
                let dim = 1usize << width;
                let half = data.theta / 2.0;
                let mut out = vec![zero; dim * dim];
                // `P` holds one non-zero per column, so the string is walked
                // column by column rather than materialized as a product.
                for column in 0..dim {
                    let mut row = column;
                    let mut value = Complex64::new(0.0, -half.sin());
                    for (index, axis) in data.axes.iter().enumerate() {
                        let bit = width - 1 - index;
                        let set = column >> bit & 1;
                        match axis {
                            PauliAxis::X => row ^= 1 << bit,
                            PauliAxis::Y => {
                                row ^= 1 << bit;
                                value *= Complex64::new(0.0, if set == 0 { 1.0 } else { -1.0 });
                            }
                            PauliAxis::Z => {
                                if set == 1 {
                                    value = -value;
                                }
                            }
                        }
                    }
                    out[row * dim + column] += value;
                    out[column * dim + column] += Complex64::new(half.cos(), 0.0);
                }
                Some(out)
            }
            Gate::Unitary(data) => Some(data.mat.clone()),
            _ if self.num_qubits() == 1 => {
                let m = self.matrix_2x2();
                Some(vec![m[0][0], m[0][1], m[1][0], m[1][1]])
            }
            _ if self.num_qubits() == 2 => {
                let m = self.matrix_4x4();
                Some(m.iter().flat_map(|row| row.iter().copied()).collect())
            }
            _ => None,
        }
    }

    /// Lowercase gate name, the OpenQASM spelling where one exists.
    #[inline]
    pub fn name(&self) -> &'static str {
        match self {
            Gate::Id => "id",
            Gate::X => "x",
            Gate::Y => "y",
            Gate::Z => "z",
            Gate::H => "h",
            Gate::S => "s",
            Gate::Sdg => "sdg",
            Gate::T => "t",
            Gate::Tdg => "tdg",
            Gate::SX => "sx",
            Gate::SXdg => "sxdg",
            Gate::Rx(_) => "rx",
            Gate::Ry(_) => "ry",
            Gate::Rz(_) => "rz",
            Gate::P(_) => "p",
            Gate::Rzz(_) => "rzz",
            Gate::Cx => "cx",
            Gate::Cz => "cz",
            Gate::Swap => "swap",
            Gate::Cu(_) => "cu",
            Gate::Mcu(_) => "mcu",
            Gate::Fused(_) => "fused",
            Gate::BatchPhase(_) => "batch_phase",
            Gate::QftBlock { .. } => "qft_block",
            Gate::PauliRot(_) => "pauli_rot",
            Gate::MultiPauliRot(_) => "multi_pauli_rot",
            Gate::Unitary(_) => "unitary",
            Gate::BatchRzz(_) => "batch_rzz",
            Gate::DiagonalBatch(_) => "diagonal_batch",
            Gate::MultiFused(_) => "multi_fused",
            Gate::Fused2q(_) => "fused_2q",
            Gate::Multi2q(_) => "multi_2q",
        }
    }

    /// Adjoint of this gate.
    pub fn inverse(&self) -> Gate {
        match self {
            Gate::Id | Gate::X | Gate::Y | Gate::Z | Gate::H => self.clone(),
            Gate::S => Gate::Sdg,
            Gate::Sdg => Gate::S,
            Gate::T => Gate::Tdg,
            Gate::Tdg => Gate::T,
            Gate::SX => Gate::SXdg,
            Gate::SXdg => Gate::SX,
            Gate::Rx(theta) => Gate::Rx(-theta),
            Gate::Ry(theta) => Gate::Ry(-theta),
            Gate::Rz(theta) => Gate::Rz(-theta),
            Gate::P(theta) => Gate::P(-theta),
            Gate::Rzz(theta) => Gate::Rzz(-theta),
            Gate::Cx | Gate::Cz | Gate::Swap => self.clone(),
            Gate::Cu(mat) => Gate::cu(adjoint_2x2(mat)),
            Gate::Mcu(data) => Gate::mcu(adjoint_2x2(&data.mat), data.num_controls),
            Gate::Fused(mat) => Gate::Fused(Box::new(adjoint_2x2(mat))),
            Gate::BatchPhase(data) => Gate::BatchPhase(Box::new(BatchPhaseData {
                phases: data.phases.iter().map(|&(q, p)| (q, p.conj())).collect(),
            })),
            Gate::QftBlock {
                start,
                num,
                inverse,
                swaps,
                big_endian,
            } => Gate::QftBlock {
                start: *start,
                num: *num,
                inverse: !inverse,
                swaps: *swaps,
                big_endian: *big_endian,
            },
            Gate::PauliRot(data) => Gate::PauliRot(Box::new(PauliRotData {
                theta: -data.theta,
                axes: data.axes.clone(),
            })),
            Gate::MultiPauliRot(data) => Gate::MultiPauliRot(Box::new(MultiPauliRotData {
                rotations: data
                    .rotations
                    .iter()
                    .rev()
                    .map(|&(xmask, zmask, theta)| (xmask, zmask, -theta))
                    .collect(),
            })),
            Gate::Unitary(data) => Gate::Unitary(Box::new(data.adjoint())),
            Gate::BatchRzz(data) => Gate::BatchRzz(Box::new(BatchRzzData {
                edges: data
                    .edges
                    .iter()
                    .map(|&(q0, q1, theta)| (q0, q1, -theta))
                    .collect(),
            })),
            Gate::DiagonalBatch(data) => Gate::DiagonalBatch(Box::new(DiagonalBatchData {
                entries: data
                    .entries
                    .iter()
                    .map(|e| match e {
                        DiagEntry::Phase1q { qubit, d0, d1 } => DiagEntry::Phase1q {
                            qubit: *qubit,
                            d0: d0.conj(),
                            d1: d1.conj(),
                        },
                        DiagEntry::Phase2q { q0, q1, phase } => DiagEntry::Phase2q {
                            q0: *q0,
                            q1: *q1,
                            phase: phase.conj(),
                        },
                        DiagEntry::Parity2q { q0, q1, same, diff } => DiagEntry::Parity2q {
                            q0: *q0,
                            q1: *q1,
                            same: same.conj(),
                            diff: diff.conj(),
                        },
                    })
                    .collect(),
            })),
            Gate::MultiFused(data) => Gate::MultiFused(Box::new(MultiFusedData {
                gates: data
                    .gates
                    .iter()
                    .map(|&(target, mat)| (target, adjoint_2x2(&mat)))
                    .collect(),
                all_diagonal: data.all_diagonal,
            })),
            Gate::Fused2q(mat) => Gate::Fused2q(Box::new(adjoint_4x4(mat))),
            Gate::Multi2q(data) => Gate::Multi2q(Box::new(Multi2qData {
                gates: data
                    .gates
                    .iter()
                    .rev()
                    .map(|&(q0, q1, ref mat)| (q0, q1, adjoint_4x4(mat)))
                    .collect(),
            })),
        }
    }

    /// Analytic differentiation generator, or `None` for a gate with no angle in
    /// its variant (non-parametric gates, controlled unitaries built from a matrix).
    #[inline]
    pub fn pauli_generator(&self) -> Option<GeneratorKind<'_>> {
        match self {
            Gate::Rx(_) => Some(GeneratorKind::RotX),
            Gate::Ry(_) => Some(GeneratorKind::RotY),
            Gate::Rz(_) => Some(GeneratorKind::RotZ),
            Gate::Rzz(_) => Some(GeneratorKind::RotZz),
            Gate::P(_) => Some(GeneratorKind::Phase),
            Gate::PauliRot(data) => Some(GeneratorKind::RotPauli(&data.axes)),
            _ => None,
        }
    }

    /// `U^k` for a single-qubit gate; negative `k` inverts first.
    pub(crate) fn matrix_power(&self, k: i64) -> Gate {
        debug_assert_eq!(
            self.num_qubits(),
            1,
            "matrix_power only for single-qubit gates"
        );
        if k == 0 {
            return Gate::Id;
        }
        if k == 1 {
            return self.clone();
        }
        let base = if k < 0 { self.inverse() } else { self.clone() };
        let n = k.unsigned_abs() as usize;
        if n == 1 {
            return base;
        }
        let base_mat = base.matrix_2x2();
        let mut acc = base_mat;
        for _ in 1..n {
            acc = mat_mul_2x2(&base_mat, &acc);
        }
        Gate::Fused(Box::new(acc))
    }

    /// `U^t` for a real exponent, on a single-qubit gate.
    ///
    /// The principal power: `U` is written as
    /// `e^(i(phi + a)) P+ + e^(i(phi - a)) P-` through the axis-angle form
    /// every 2x2 unitary has, each eigenvalue angle is brought onto the
    /// principal branch `(-pi, pi]`, and both are scaled by `t`. Reading it
    /// from the spectrum rather than from the axis-angle pair keeps
    /// `pow(0.5) @ rx(theta)` at `rx(theta/2)` for a rotation past `pi`, where
    /// scaling the half-angle picks the other root.
    ///
    /// The two decompositions of `U` (`phi` and `phi + pi`) give the same pair
    /// of eigenvalue angles with the projectors exchanged, so the branch `arg`
    /// takes for the determinant does not reach the answer. That is what makes
    /// the sign of a zero harmless here.
    ///
    /// A whole-number exponent goes through [`Gate::matrix_power`] instead, so
    /// that `pow(2) @ t` stays the named `S` rather than becoming a matrix.
    ///
    /// # Panics
    /// Panics in debug builds if the gate is not single-qubit.
    pub(crate) fn matrix_power_real(&self, t: f64) -> Gate {
        debug_assert_eq!(
            self.num_qubits(),
            1,
            "matrix_power_real only for single-qubit gates"
        );
        if t.fract() == 0.0 && t.abs() <= i64::MAX as f64 {
            return self.matrix_power(t as i64);
        }
        let m = self.matrix_2x2();
        let det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
        let phase = det.arg() / 2.0;
        let special = special_unitary(&m, phase);
        let cosine = (special[0][0].re + special[1][1].re) / 2.0;
        let axis = [
            (special[0][1].im + special[1][0].im) / 2.0,
            (special[0][1].re - special[1][0].re) / 2.0,
            (special[0][0].im - special[1][1].im) / 2.0,
        ];
        let norm = axis.iter().map(|c| c * c).sum::<f64>().sqrt();
        let angle = norm.atan2(cosine);
        let unit = if norm > 0.0 {
            [axis[0] / norm, axis[1] / norm, axis[2] / norm]
        } else {
            [0.0, 0.0, 1.0]
        };
        let upper = Complex64::from_polar(1.0, t * principal_angle(phase + angle));
        let lower = Complex64::from_polar(1.0, t * principal_angle(phase - angle));
        let projected = [
            [
                Complex64::new(unit[2], 0.0),
                Complex64::new(unit[0], -unit[1]),
            ],
            [
                Complex64::new(unit[0], unit[1]),
                Complex64::new(-unit[2], 0.0),
            ],
        ];
        let mut powered = [[Complex64::new(0.0, 0.0); 2]; 2];
        for (row, (target, source)) in powered.iter_mut().zip(projected).enumerate() {
            for (column, (cell, entry)) in target.iter_mut().zip(source).enumerate() {
                let identity = Complex64::new(f64::from(u8::from(row == column)), 0.0);
                *cell = ((identity + entry) * upper + (identity - entry) * lower) * 0.5;
            }
        }
        Gate::Fused(Box::new(powered))
    }

    pub fn cu(mat: [[Complex64; 2]; 2]) -> Gate {
        Gate::Cu(Box::new(mat))
    }

    /// Multi-controlled unitary with `num_controls` controls. One control returns
    /// [`Gate::Cu`], so an `Mcu` always carries at least two.
    ///
    /// # Panics
    /// Panics if `num_controls` is zero; use `Gate::Fused(Box::new(mat))` for an
    /// uncontrolled matrix.
    pub fn mcu(mat: [[Complex64; 2]; 2], num_controls: u8) -> Gate {
        assert!(
            num_controls > 0,
            "Gate::mcu needs at least one control qubit"
        );
        if num_controls == 1 {
            return Gate::cu(mat);
        }
        Gate::Mcu(Box::new(McuData { mat, num_controls }))
    }

    /// Build the gate applying the `2^k x 2^k` matrix `mat` to `k` qubits.
    ///
    /// `mat` is row major with `targets[0]` the most significant bit of both
    /// indices, the packing [`Gate::matrix_4x4`] uses. The result is whichever
    /// existing variant carries the matrix: a named gate or `Fused` at one
    /// qubit; `Cx`, `Cz`, `Swap`, `Rzz`, `Cu`, or `Fused2q` at two; `Mcu` when
    /// the matrix is the identity outside its trailing 2x2 block; and
    /// [`Gate::Unitary`] only for a dense matrix on three or more qubits, so
    /// Clifford recognition, controlled-phase dispatch, and fusion keep firing
    /// on the forms they know. A diagonal on three or more qubits needs the
    /// target qubits to lower; [`crate::Circuit::add_unitary`] does that.
    ///
    /// # Errors
    /// `InvalidParameter` when `num_qubits` is zero or above
    /// [`MAX_UNITARY_QUBITS`], when `mat.len()` is not `4^num_qubits`, when an
    /// entry is not finite, or when `U^dagger U` differs from the identity by
    /// more than `1e-9` in any entry.
    pub fn unitary(mat: Vec<Complex64>, num_qubits: usize) -> crate::error::Result<Gate> {
        check_unitary_matrix(&mat, num_qubits)?;
        Ok(match num_qubits {
            1 => {
                let m = [[mat[0], mat[1]], [mat[2], mat[3]]];
                Gate::recognize_matrix(&m).unwrap_or_else(|| Gate::Fused(Box::new(m)))
            }
            2 => {
                let mut m = [[Complex64::new(0.0, 0.0); 4]; 4];
                for (row, entries) in m.iter_mut().enumerate() {
                    entries.copy_from_slice(&mat[row * 4..row * 4 + 4]);
                }
                recognize_matrix_4x4(&m)
            }
            _ => match controlled_block(&mat, num_qubits) {
                Some(block) => Gate::mcu(block, (num_qubits - 1) as u8),
                None => Gate::Unitary(Box::new(UnitaryData {
                    mat,
                    num_qubits: num_qubits as u8,
                })),
            },
        })
    }

    /// Controlled-phase gate CPhase(θ) = Cu(\[\[1,0\],\[0,e^{iθ}\]\]).
    pub fn cphase(theta: f64) -> Gate {
        let one = Complex64::new(1.0, 0.0);
        let zero = Complex64::new(0.0, 0.0);
        let phase = Complex64::from_polar(1.0, theta);
        Gate::cu([[one, zero], [zero, phase]])
    }

    /// The phase `e^{iθ}` when this is a `Cu` or `Mcu` of `[[1,0],[0,e^{iθ}]]`,
    /// which backends route to phase-only kernels.
    #[inline]
    pub fn controlled_phase(&self) -> Option<Complex64> {
        let mat = match self {
            Gate::Cu(mat) => &**mat,
            Gate::Mcu(data) => &data.mat,
            _ => return None,
        };
        if (mat[0][0].re - 1.0).abs() < IDENTITY_EPS
            && mat[0][0].im.abs() < IDENTITY_EPS
            && mat[0][1].norm() < IDENTITY_EPS
            && mat[1][0].norm() < IDENTITY_EPS
            && (mat[1][1].norm() - 1.0).abs() < IDENTITY_EPS
        {
            Some(mat[1][1])
        } else {
            None
        }
    }

    /// True if this is a diagonal single-qubit gate (matrix is `[[a,0],[0,b]]`).
    #[inline]
    pub fn is_diagonal_1q(&self) -> bool {
        match self {
            Gate::Id
            | Gate::Z
            | Gate::S
            | Gate::Sdg
            | Gate::T
            | Gate::Tdg
            | Gate::Rz(_)
            | Gate::P(_) => true,
            Gate::Fused(m) => is_diagonal_2x2(m),
            _ => false,
        }
    }

    /// True when this gate is diagonal and [`Self::push_diag_entries`] can express it.
    #[inline]
    pub(crate) fn is_diag_batchable(&self) -> bool {
        match self {
            Gate::Cz | Gate::Rzz(_) | Gate::BatchPhase(_) | Gate::BatchRzz(_) => true,
            _ if self.is_diagonal_1q() => true,
            Gate::Cu(mat) => is_diagonal_2x2(mat),
            Gate::Fused2q(mat) => is_diagonal_4x4(mat),
            Gate::MultiFused(data) => data.all_diagonal,
            Gate::Multi2q(data) => data.gates.iter().all(|(_, _, m)| is_diagonal_4x4(m)),
            _ => false,
        }
    }

    /// Append the `DiagEntry` values equivalent to this gate on `targets`. Only valid
    /// when `is_diag_batchable()` holds; off-diagonal entries are ignored.
    pub(crate) fn push_diag_entries(&self, targets: &[usize], out: &mut Vec<DiagEntry>) {
        match self {
            Gate::Cz => out.push(DiagEntry::Phase2q {
                q0: targets[0],
                q1: targets[1],
                phase: Complex64::new(-1.0, 0.0),
            }),
            Gate::Rzz(theta) => out.push(rzz_entry(targets[0], targets[1], *theta)),
            Gate::BatchRzz(data) => out.extend(
                data.edges
                    .iter()
                    .map(|&(q0, q1, theta)| rzz_entry(q0, q1, theta)),
            ),
            Gate::BatchPhase(data) => {
                out.extend(data.phases.iter().map(|&(q1, phase)| DiagEntry::Phase2q {
                    q0: targets[0],
                    q1,
                    phase,
                }))
            }
            Gate::Cu(mat) => {
                let phase = match self.controlled_phase() {
                    Some(phase) => phase,
                    None => {
                        out.push(DiagEntry::Phase1q {
                            qubit: targets[0],
                            d0: Complex64::new(1.0, 0.0),
                            d1: mat[0][0],
                        });
                        mat[1][1] / mat[0][0]
                    }
                };
                out.push(DiagEntry::Phase2q {
                    q0: targets[0],
                    q1: targets[1],
                    phase,
                });
            }
            Gate::Fused2q(mat) => push_diag_4x4_entries(targets[0], targets[1], mat, out),
            Gate::Multi2q(data) => {
                for (q0, q1, mat) in &data.gates {
                    push_diag_4x4_entries(*q0, *q1, mat, out);
                }
            }
            Gate::MultiFused(data) => {
                out.extend(data.gates.iter().map(|(qubit, m)| DiagEntry::Phase1q {
                    qubit: *qubit,
                    d0: m[0][0],
                    d1: m[1][1],
                }))
            }
            _ => {
                let mat = self.matrix_2x2();
                out.push(DiagEntry::Phase1q {
                    qubit: targets[0],
                    d0: mat[0][0],
                    d1: mat[1][1],
                });
            }
        }
    }

    #[inline]
    pub(crate) fn is_self_inverse_2q(&self) -> bool {
        matches!(self, Gate::Cx | Gate::Cz | Gate::Swap)
    }

    /// True if this gate maps basis states to basis states up to a phase (diagonal
    /// and permutation gates), so it keeps the count of non-zero amplitudes.
    #[inline]
    pub(crate) fn preserves_sparsity(&self) -> bool {
        match self {
            Gate::Id | Gate::X | Gate::Y | Gate::Z => true,
            Gate::S | Gate::Sdg | Gate::T | Gate::Tdg => true,
            Gate::Rz(_) | Gate::P(_) => true,
            Gate::Rzz(_) | Gate::Cx | Gate::Cz | Gate::Swap => true,
            Gate::Cu(mat) | Gate::Fused(mat) => {
                let is_diag = mat[0][1].norm_sqr() < NEAR_ZERO_NORM_SQ
                    && mat[1][0].norm_sqr() < NEAR_ZERO_NORM_SQ;
                let is_antidiag = mat[0][0].norm_sqr() < NEAR_ZERO_NORM_SQ
                    && mat[1][1].norm_sqr() < NEAR_ZERO_NORM_SQ;
                is_diag || is_antidiag
            }
            Gate::Mcu(data) => {
                let m = &data.mat;
                let is_diag = m[0][1].norm_sqr() < NEAR_ZERO_NORM_SQ
                    && m[1][0].norm_sqr() < NEAR_ZERO_NORM_SQ;
                let is_antidiag = m[0][0].norm_sqr() < NEAR_ZERO_NORM_SQ
                    && m[1][1].norm_sqr() < NEAR_ZERO_NORM_SQ;
                is_diag || is_antidiag
            }
            Gate::BatchPhase(_) | Gate::BatchRzz(_) | Gate::DiagonalBatch(_) => true,
            Gate::PauliRot(data) => data.axes.iter().all(|axis| *axis == PauliAxis::Z),
            _ => false,
        }
    }

    /// Recognize a 2x2 unitary as a named gate, so passes downstream of fusion see
    /// the Cliffords it produced (T·T → S).
    ///
    /// The match is exact, not up to a global phase: a named gate carries no
    /// scalar, so recognizing `e^{iπ/4}·X` as `X` would change the amplitudes.
    /// Such a product stays a [`Gate::Fused`].
    pub(crate) fn recognize_matrix(mat: &[[Complex64; 2]; 2]) -> Option<Gate> {
        NAMED_1Q_CANDIDATES
            .iter()
            .find(|candidate| matrices_equal(mat, &candidate.matrix_2x2(), RECOGNIZE_EPS))
            .cloned()
    }

    /// Recognize a 2x2 unitary as a named gate or an `Rz`, up to a global phase.
    ///
    /// Returns the gate and the scalar `c` with `mat = c * gate.matrix_2x2()`.
    /// Diagonal matrices outside the named set come back as `Rz(theta)` with
    /// `theta` the argument of `mat[1][1] / mat[0][0]`.
    ///
    /// The phase is dropped on purpose here and nowhere else: fusion keeps the
    /// exact rule above because dense backends export amplitudes, while the
    /// Pauli engines and the stabilizer-rank expansion only ever return
    /// expectation values, probabilities, and samples, all of which are
    /// invariant under a global phase on an uncontrolled gate. A controlled use
    /// restores `c` as a phase rotation on the control.
    pub(crate) fn recognize_matrix_up_to_phase(
        mat: &[[Complex64; 2]; 2],
    ) -> Option<(Gate, Complex64)> {
        for candidate in NAMED_1Q_CANDIDATES {
            let named = candidate.matrix_2x2();
            let (r, c) = if named[0][0].norm_sqr() > 0.25 {
                (0, 0)
            } else {
                (1, 0)
            };
            let phase = mat[r][c] / named[r][c];
            if (phase.norm() - 1.0).abs() > RECOGNIZE_EPS {
                continue;
            }
            let scaled = [
                [named[0][0] * phase, named[0][1] * phase],
                [named[1][0] * phase, named[1][1] * phase],
            ];
            if matrices_equal(mat, &scaled, RECOGNIZE_EPS) {
                return Some((candidate.clone(), phase));
            }
        }
        if is_diagonal_2x2(mat) {
            let theta = (mat[1][1] * mat[0][0].conj()).arg();
            let phase = mat[0][0] * Complex64::from_polar(1.0, theta / 2.0);
            return Some((Gate::Rz(theta), phase));
        }
        None
    }

    /// True for the named Clifford variants; a Clifford carried as a matrix or an
    /// angle (`Fused`, `Rz(π/2)`) reports false.
    #[inline]
    pub fn is_clifford(&self) -> bool {
        matches!(
            self,
            Gate::Id
                | Gate::X
                | Gate::Y
                | Gate::Z
                | Gate::H
                | Gate::S
                | Gate::Sdg
                | Gate::SX
                | Gate::SXdg
                | Gate::Cx
                | Gate::Cz
                | Gate::Swap
        )
    }
}

/// Whether a 2x2 matrix is diagonal (both off-diagonal norms below `IDENTITY_EPS`).
#[inline]
pub(crate) fn is_diagonal_2x2(mat: &[[Complex64; 2]; 2]) -> bool {
    mat[0][1].norm() < IDENTITY_EPS && mat[1][0].norm() < IDENTITY_EPS
}

/// Whether a 2x2 matrix is antidiagonal (both diagonal norms below `IDENTITY_EPS`).
#[inline]
pub(crate) fn is_antidiagonal_2x2(mat: &[[Complex64; 2]; 2]) -> bool {
    mat[0][0].norm() < IDENTITY_EPS && mat[1][1].norm() < IDENTITY_EPS
}

/// Whether a 4x4 matrix is diagonal (all off-diagonal norms below `IDENTITY_EPS`).
#[inline]
pub(crate) fn is_diagonal_4x4(mat: &[[Complex64; 4]; 4]) -> bool {
    for (r, row) in mat.iter().enumerate() {
        for (c, value) in row.iter().enumerate() {
            if r != c && value.norm() >= IDENTITY_EPS {
                return false;
            }
        }
    }
    true
}

fn rzz_entry(q0: usize, q1: usize, theta: f64) -> DiagEntry {
    DiagEntry::Parity2q {
        q0,
        q1,
        same: Complex64::from_polar(1.0, -theta / 2.0),
        diff: Complex64::from_polar(1.0, theta / 2.0),
    }
}

/// Split a diagonal two-qubit unitary on `[q0, q1]` (row index `2 * b0 + b1`) into a
/// phase on each qubit and one on `|11>`. Unitarity keeps every divisor at unit modulus.
fn push_diag_4x4_entries(
    q0: usize,
    q1: usize,
    mat: &[[Complex64; 4]; 4],
    out: &mut Vec<DiagEntry>,
) {
    let (d00, d01, d10, d11) = (mat[0][0], mat[1][1], mat[2][2], mat[3][3]);
    out.push(DiagEntry::Phase1q {
        qubit: q0,
        d0: d00,
        d1: d10,
    });
    out.push(DiagEntry::Phase1q {
        qubit: q1,
        d0: Complex64::new(1.0, 0.0),
        d1: d01 / d00,
    });
    out.push(DiagEntry::Phase2q {
        q0,
        q1,
        phase: d11 * d00 / (d10 * d01),
    });
}

/// Per-entry magnitude tolerance for `recognize_matrix` and
/// `recognize_matrix_up_to_phase`.
///
/// Taken from `IDENTITY_EPS` so neither can accept a matrix the fusion pass
/// would refuse to drop as an identity: a disc of this radius on each entry
/// sits inside the box the identity check puts on the real and imaginary
/// parts. A tolerance on `norm_sqr` would admit the square root of this, which
/// is large enough to swallow a whole rotation: `Rx(1e-5)` differs from the
/// identity by 5e-6 per entry but by only 2.5e-11 squared.
const RECOGNIZE_EPS: f64 = IDENTITY_EPS;

/// Named single-qubit gates the matrix recognizers try, in order.
const NAMED_1Q_CANDIDATES: &[Gate] = &[
    Gate::Id,
    Gate::H,
    Gate::X,
    Gate::Y,
    Gate::Z,
    Gate::S,
    Gate::Sdg,
    Gate::T,
    Gate::Tdg,
    Gate::SX,
    Gate::SXdg,
];

fn matrices_equal(a: &[[Complex64; 2]; 2], b: &[[Complex64; 2]; 2], eps: f64) -> bool {
    (0..2).all(|i| (0..2).all(|j| (a[i][j] - b[i][j]).norm() <= eps))
}

fn matrices_equal_4x4(a: &[[Complex64; 4]; 4], b: &[[Complex64; 4]; 4], eps: f64) -> bool {
    (0..4).all(|i| (0..4).all(|j| (a[i][j] - b[i][j]).norm() <= eps))
}

/// The named two-qubit gate `mat` equals within `RECOGNIZE_EPS`, then `Rzz`,
/// then `Cu` when the matrix is the identity outside its trailing block, else
/// `Fused2q`. Exact, not up to a global phase, for the reason
/// [`Gate::recognize_matrix`] gives.
fn recognize_matrix_4x4(mat: &[[Complex64; 4]; 4]) -> Gate {
    for candidate in [Gate::Cx, Gate::Cz, Gate::Swap] {
        if matrices_equal_4x4(mat, &candidate.matrix_4x4(), RECOGNIZE_EPS) {
            return candidate;
        }
    }
    if is_diagonal_4x4(mat) {
        let theta = 2.0 * mat[1][1].arg();
        if matrices_equal_4x4(mat, &Gate::Rzz(theta).matrix_4x4(), RECOGNIZE_EPS) {
            return Gate::Rzz(theta);
        }
    }
    let flat: Vec<Complex64> = mat.iter().flat_map(|row| row.iter().copied()).collect();
    match controlled_block(&flat, 2) {
        Some(block) => Gate::cu(block),
        None => Gate::Fused2q(Box::new(*mat)),
    }
}

fn format_angle(theta: f64) -> String {
    const FRACTIONS: &[(f64, &str)] = &[
        (1.0, "π"),
        (-1.0, "-π"),
        (0.5, "π/2"),
        (-0.5, "-π/2"),
        (0.25, "π/4"),
        (-0.25, "-π/4"),
        (1.0 / 3.0, "π/3"),
        (-1.0 / 3.0, "-π/3"),
        (2.0 / 3.0, "2π/3"),
        (-2.0 / 3.0, "-2π/3"),
        (1.0 / 6.0, "π/6"),
        (-1.0 / 6.0, "-π/6"),
        (5.0 / 6.0, "5π/6"),
        (-5.0 / 6.0, "-5π/6"),
        (1.0 / 8.0, "π/8"),
        (-1.0 / 8.0, "-π/8"),
        (3.0 / 8.0, "3π/8"),
        (-3.0 / 8.0, "-3π/8"),
        (1.5, "3π/2"),
        (-1.5, "-3π/2"),
        (2.0, "2π"),
        (-2.0, "-2π"),
    ];
    let ratio = theta / std::f64::consts::PI;
    for &(frac, label) in FRACTIONS {
        if (ratio - frac).abs() < 1e-10 {
            return label.to_string();
        }
    }
    format!("{:.4}", theta)
}

impl fmt::Display for Gate {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Gate::Id => f.write_str("I"),
            Gate::X => f.write_str("X"),
            Gate::Y => f.write_str("Y"),
            Gate::Z => f.write_str("Z"),
            Gate::H => f.write_str("H"),
            Gate::S => f.write_str("S"),
            Gate::Sdg => f.write_str("Sdg"),
            Gate::T => f.write_str("T"),
            Gate::Tdg => f.write_str("Tdg"),
            Gate::SX => f.write_str("SX"),
            Gate::SXdg => f.write_str("SXdg"),
            Gate::Rx(t) => write!(f, "Rx({})", format_angle(*t)),
            Gate::Ry(t) => write!(f, "Ry({})", format_angle(*t)),
            Gate::Rz(t) => write!(f, "Rz({})", format_angle(*t)),
            Gate::P(t) => write!(f, "P({})", format_angle(*t)),
            Gate::Rzz(t) => write!(f, "Rzz({})", format_angle(*t)),
            Gate::Cx => f.write_str("CX"),
            Gate::Cz => f.write_str("CZ"),
            Gate::Swap => f.write_str("SWAP"),
            Gate::Cu(_) => f.write_str("CU"),
            Gate::Mcu(data) => write!(f, "MCU({}ctrl)", data.num_controls),
            Gate::Fused(_) => f.write_str("U"),
            Gate::Fused2q(_) => f.write_str("U2"),
            Gate::Unitary(data) => write!(f, "U{}", data.num_qubits()),
            Gate::MultiFused(data) => write!(f, "MF[{}]", data.gates.len()),
            Gate::BatchPhase(data) => write!(f, "BP[{}]", data.phases.len()),
            Gate::QftBlock {
                start,
                num,
                inverse,
                swaps,
                big_endian,
            } => {
                let end = *start as usize + *num as usize;
                let prefix = if *inverse { "I" } else { "" };
                let no_swaps = if *swaps { "" } else { ",ns" };
                let endian = if *big_endian { ",be" } else { "" };
                write!(f, "{prefix}QFT[{start}..{end}{no_swaps}{endian}]")
            }
            Gate::PauliRot(data) => {
                f.write_str("R[")?;
                for axis in &data.axes {
                    f.write_str(match axis {
                        PauliAxis::X => "X",
                        PauliAxis::Y => "Y",
                        PauliAxis::Z => "Z",
                    })?;
                }
                write!(f, "]({})", format_angle(data.theta))
            }
            Gate::MultiPauliRot(data) => write!(f, "MR[{}]", data.rotations.len()),
            Gate::BatchRzz(data) => write!(f, "BZZ[{}]", data.edges.len()),
            Gate::DiagonalBatch(data) => write!(f, "BD[{}]", data.entries.len()),
            Gate::Multi2q(data) => write!(f, "M2[{}]", data.gates.len()),
        }
    }
}

pub(crate) mod spectral;

#[cfg(test)]
mod tests;