cqlib-core 0.1.0-beta.4

// This code is part of Cqlib.
//
// (C) Copyright China Telecom Quantum Group 2025-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 http://www.apache.org/licenses/LICENSE-2.0.
//
// Any modifications or derivative works of this code must retain this
// copyright notice, and modified files need to carry a notice indicating
// that they have been altered from the originals.

//! Structural algebraic commutation checks.
//!
//! This module contains symbolic proofs that do not require building matrices.
//! It recognizes gates that can be modeled as Pauli-axis products, rotations
//! around tensor-product axes, controlled single-axis operations, diagonal
//! operators, and a few symmetric two-qubit gate families.  The implementation
//! favors soundness over coverage: when an axis relation or parameter identity
//! cannot be proven from the available symbolic expressions, it returns `None`
//! and lets the high-level checker try later oracles.
//!
//! # Examples
//!
//! ```rust
//! use cqlib_core::circuit::{Instruction, Parameter, Qubit, StandardGate};
//! use cqlib_core::compile::commutation::{
//!     Commutation, CommutationChecker, CommutationConfig,
//! };
//!
//! let checker = CommutationChecker::with_config(CommutationConfig {
//!     enable_rule_oracle: false,
//!     enable_matrix_fallback: false,
//!     ..CommutationConfig::default()
//! });
//! let result = checker.check(
//!     &Instruction::Standard(StandardGate::RXX),
//!     &[Qubit::new(0), Qubit::new(1)],
//!     &[Parameter::symbol("a")],
//!     &Instruction::Standard(StandardGate::RZZ),
//!     &[Qubit::new(0), Qubit::new(1)],
//!     &[Parameter::symbol("b")],
//! );
//!
//! assert_eq!(result, Some(Commutation::Exact));
//! ```

use crate::circuit::{Instruction, MCGate, Parameter, Qubit, StandardGate};
use crate::compile::PARAMETER_EQ_TOLERANCE;
use crate::compile::commutation::checker::{Commutation, CommutationResult};
use smallvec::{SmallVec, smallvec};
use std::f64::consts::{FRAC_PI_2, PI};

const TWO_PI: f64 = 2.0 * PI;
const FOUR_PI: f64 = 4.0 * PI;

/// Cartesian Pauli axes used by Pauli products and rotations.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum PauliAxis {
    /// Pauli X axis.
    X,
    /// Pauli Y axis.
    Y,
    /// Pauli Z axis.
    Z,
}

/// Single-qubit axis understood by the algebraic prover.
///
/// `Planar(phi)` represents an axis in the XY plane.  Symbolically, it is
/// compared modulo `PI` because opposite axes generate the same product up to a
/// sign and the same rotation family under the commutation tests used here.
#[derive(Debug, Clone)]
enum Axis {
    /// One of the three Cartesian Pauli axes.
    Pauli(PauliAxis),
    /// Axis in the XY plane, parameterized by its angle from X.
    Planar(Parameter),
}

/// One tensor factor in a Pauli-like axis string.
#[derive(Debug, Clone)]
struct AxisTerm {
    /// Concrete circuit qubit on which this tensor factor acts.
    qubit: Qubit,
    /// Axis acting on `qubit`.
    axis: Axis,
}

/// Compact tensor-product axis representation.
///
/// Most gates handled here touch one or two qubits, so `SmallVec` avoids heap
/// allocation on the common paths.
type AxisString = SmallVec<[AxisTerm; 2]>;

/// Whether an axis string denotes a product or an arbitrary rotation.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum AxisOpKind {
    /// A Pauli-like product, including rotations known to be product-valued.
    Product,
    /// A rotation generated by the axis string.
    Rotation,
}

/// Axis-string abstraction used for parity-based commutation checks.
///
/// Two Pauli strings commute when the number of qubits on which their axes
/// anti-commute is even.  Product/product pairs with odd parity still commute
/// up to the global phase `PI`; arbitrary rotations do not receive that proof.
#[derive(Debug, Clone)]
struct AxisOp {
    /// Tensor factors participating in the axis string.
    terms: AxisString,
    /// Product versus rotation semantics for odd anti-commutation parity.
    kind: AxisOpKind,
}

/// Controlled single-target operation represented by its control set and axis.
///
/// Controls are treated as Z-basis predicates.  The target operation is modeled
/// only by its single target axis, so multi-target controlled bases are outside
/// this structural proof.
#[derive(Debug, Clone)]
struct ControlledAxis {
    /// Control qubits, in the concrete operation order.
    controls: SmallVec<[Qubit; 4]>,
    /// Target qubit of the controlled axis operation.
    target: Qubit,
    /// Axis applied on the target when controls are satisfied.
    axis: Axis,
}

/// Proven relationship between two single-qubit axes.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum AxisRelation {
    /// The axes are equal up to the equivalence used by this prover.
    Same,
    /// The axes anti-commute on the same qubit.
    Orthogonal,
    /// The prover cannot establish either relation symbolically.
    Unknown,
}

/// Symmetric two-qubit families with pair-order-independent commutation facts.
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
enum TwoQubitFamily {
    /// The SWAP permutation.
    Swap,
    /// The parameterized fSim family.
    Fsim,
    /// Symmetric Pauli interactions such as RXX and RYY.
    SymmetricPauliRotation,
    /// Symmetric diagonal interactions such as RZZ and CZ.
    SymmetricDiagonal,
}

/// Attempts symbolic algebraic commutation proofs for one operation pair.
///
/// This function is called after the high-level checker has already handled
/// arity validation, identity, same-application, and disjoint-support facts.
/// Returning `None` keeps the result conservative and allows rule or matrix
/// fallbacks to run.
pub fn algebraic_commutation(
    lhs_inst: &Instruction,
    lhs_qubits: &[Qubit],
    lhs_params: &[Parameter],
    rhs_inst: &Instruction,
    rhs_qubits: &[Qubit],
    rhs_params: &[Parameter],
) -> CommutationResult {
    if instruction_is_global_phase(lhs_inst, lhs_params)
        || instruction_is_global_phase(rhs_inst, rhs_params)
    {
        return Some(Commutation::Exact);
    }

    if instruction_is_diagonal(lhs_inst) && instruction_is_diagonal(rhs_inst) {
        return Some(Commutation::Exact);
    }

    if let (Some(lhs), Some(rhs)) = (
        axis_op_for_instruction(lhs_inst, lhs_qubits, lhs_params),
        axis_op_for_instruction(rhs_inst, rhs_qubits, rhs_params),
    ) && let Some(result) = axis_ops_commutation(&lhs, &rhs)
    {
        return Some(result);
    }

    if let Some(result) = controlled_commutation(
        lhs_inst, lhs_qubits, lhs_params, rhs_inst, rhs_qubits, rhs_params,
    ) {
        return Some(result);
    }

    special_two_qubit_commutation(lhs_inst, lhs_qubits, rhs_inst, rhs_qubits)
}

/// Returns whether an instruction application is a scalar global phase.
///
/// For multi-controlled gates, a scalar base gate generally becomes a relative
/// phase conditioned on controls, so only base identities are treated as
/// globally harmless.
fn instruction_is_global_phase(instruction: &Instruction, params: &[Parameter]) -> bool {
    match instruction {
        Instruction::Standard(gate) => standard_is_global_phase(*gate, params),
        Instruction::McGate(gate) if gate.num_qubits() == gate.base_gate().num_qubits() => {
            standard_is_global_phase(*gate.base_gate(), params)
        }
        Instruction::McGate(gate) => standard_is_identity(*gate.base_gate(), params),
        _ => false,
    }
}

/// Returns whether a standard gate is represented by a scalar matrix.
fn standard_is_global_phase(gate: StandardGate, params: &[Parameter]) -> bool {
    match gate {
        StandardGate::I | StandardGate::GPhase => true,
        StandardGate::RX
        | StandardGate::RY
        | StandardGate::RZ
        | StandardGate::RXX
        | StandardGate::RYY
        | StandardGate::RZZ
        | StandardGate::RZX
        | StandardGate::RXY => angle_eq_mod(&params[0], 0.0, TWO_PI),
        StandardGate::Phase => angle_eq_mod(&params[0], 0.0, TWO_PI),
        StandardGate::CRX | StandardGate::CRY | StandardGate::CRZ => {
            angle_eq_mod(&params[0], 0.0, FOUR_PI)
        }
        StandardGate::FSIM => {
            angle_eq_mod(&params[0], 0.0, TWO_PI) && angle_eq_mod(&params[1], 0.0, TWO_PI)
        }
        _ => false,
    }
}

/// Returns whether a standard gate is exactly the identity.
///
/// This differs from [`standard_is_global_phase`] for gates that may be scalar
/// but not identity, because controlled scalar phases are not global phases of
/// the whole operation.
fn standard_is_identity(gate: StandardGate, params: &[Parameter]) -> bool {
    match gate {
        StandardGate::I => true,
        StandardGate::GPhase => angle_eq_mod(&params[0], 0.0, TWO_PI),
        StandardGate::RX
        | StandardGate::RY
        | StandardGate::RZ
        | StandardGate::RXX
        | StandardGate::RYY
        | StandardGate::RZZ
        | StandardGate::RZX
        | StandardGate::RXY
        | StandardGate::CRX
        | StandardGate::CRY
        | StandardGate::CRZ => angle_eq_mod(&params[0], 0.0, FOUR_PI),
        StandardGate::Phase => angle_eq_mod(&params[0], 0.0, TWO_PI),
        StandardGate::FSIM => {
            angle_eq_mod(&params[0], 0.0, TWO_PI) && angle_eq_mod(&params[1], 0.0, TWO_PI)
        }
        _ => false,
    }
}

/// Returns whether the instruction is diagonal in the computational basis.
///
/// Two diagonal operators commute exactly even when their parameters are
/// symbolic.
fn instruction_is_diagonal(instruction: &Instruction) -> bool {
    match instruction {
        Instruction::Standard(gate) => gate.is_diagonal(),
        Instruction::McGate(gate) => gate.base_gate().is_diagonal(),
        _ => false,
    }
}

/// Applies Pauli-string parity commutation to two axis operations.
fn axis_ops_commutation(lhs: &AxisOp, rhs: &AxisOp) -> CommutationResult {
    let anti_commutations = axis_string_anti_commutations(&lhs.terms, &rhs.terms)?;
    if anti_commutations % 2 == 0 {
        return Some(Commutation::Exact);
    }

    if lhs.kind == AxisOpKind::Product && rhs.kind == AxisOpKind::Product {
        return Some(Commutation::UpToGlobalPhase(Parameter::from(PI)));
    }

    None
}

/// Counts same-qubit axis pairs that are proven to anti-commute.
///
/// Returns `None` as soon as any overlapping pair has an unknown relation,
/// preserving the conservative contract of the algebraic prover.
fn axis_string_anti_commutations(lhs: &AxisString, rhs: &AxisString) -> Option<usize> {
    let mut anti_commutations = 0usize;
    for lhs_term in lhs {
        for rhs_term in rhs {
            if lhs_term.qubit != rhs_term.qubit {
                continue;
            }

            match axis_relation(&lhs_term.axis, &rhs_term.axis) {
                AxisRelation::Same => {}
                AxisRelation::Orthogonal => anti_commutations += 1,
                AxisRelation::Unknown => return None,
            }
        }
    }
    Some(anti_commutations)
}

/// Converts a supported instruction application into an axis operation.
///
/// Multi-controlled gates with no extra controls are treated as their base gate;
/// genuinely controlled gates are handled by the controlled-axis path instead.
fn axis_op_for_instruction(
    instruction: &Instruction,
    qubits: &[Qubit],
    params: &[Parameter],
) -> Option<AxisOp> {
    match instruction {
        Instruction::Standard(gate) => axis_op_for_standard(*gate, qubits, params),
        Instruction::McGate(gate) if gate.num_qubits() == gate.base_gate().num_qubits() => {
            axis_op_for_standard(*gate.base_gate(), qubits, params)
        }
        _ => None,
    }
}

/// Converts a supported standard gate application into an axis operation.
///
/// The conversion records only the commutation-relevant generator axis and
/// whether the application behaves as a product or as a general rotation.
fn axis_op_for_standard(
    gate: StandardGate,
    qubits: &[Qubit],
    params: &[Parameter],
) -> Option<AxisOp> {
    match gate {
        StandardGate::X if qubits.len() == 1 => Some(AxisOp {
            terms: smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::X),
            }],
            kind: AxisOpKind::Product,
        }),
        StandardGate::Y if qubits.len() == 1 => Some(AxisOp {
            terms: smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::Y),
            }],
            kind: AxisOpKind::Product,
        }),
        StandardGate::Z if qubits.len() == 1 => Some(AxisOp {
            terms: smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::Z),
            }],
            kind: AxisOpKind::Product,
        }),
        StandardGate::X2P | StandardGate::X2M if qubits.len() == 1 => Some(AxisOp {
            terms: smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::X),
            }],
            kind: AxisOpKind::Rotation,
        }),
        StandardGate::Y2P | StandardGate::Y2M if qubits.len() == 1 => Some(AxisOp {
            terms: smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::Y),
            }],
            kind: AxisOpKind::Rotation,
        }),
        StandardGate::S | StandardGate::SDG | StandardGate::T | StandardGate::TDG
            if qubits.len() == 1 =>
        {
            Some(AxisOp {
                terms: smallvec![AxisTerm {
                    qubit: qubits[0],
                    axis: Axis::Pauli(PauliAxis::Z),
                }],
                kind: AxisOpKind::Rotation,
            })
        }
        StandardGate::RX if qubits.len() == 1 => Some(rotation_or_product(
            &params[0],
            smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::X),
            }],
        )),
        StandardGate::RY if qubits.len() == 1 => Some(rotation_or_product(
            &params[0],
            smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::Y),
            }],
        )),
        StandardGate::RZ | StandardGate::Phase if qubits.len() == 1 => Some(rotation_or_product(
            &params[0],
            smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Pauli(PauliAxis::Z),
            }],
        )),
        StandardGate::RXX if qubits.len() == 2 => Some(rotation_or_product(
            &params[0],
            smallvec![
                AxisTerm {
                    qubit: qubits[0],
                    axis: Axis::Pauli(PauliAxis::X),
                },
                AxisTerm {
                    qubit: qubits[1],
                    axis: Axis::Pauli(PauliAxis::X),
                }
            ],
        )),
        StandardGate::RYY if qubits.len() == 2 => Some(rotation_or_product(
            &params[0],
            smallvec![
                AxisTerm {
                    qubit: qubits[0],
                    axis: Axis::Pauli(PauliAxis::Y),
                },
                AxisTerm {
                    qubit: qubits[1],
                    axis: Axis::Pauli(PauliAxis::Y),
                }
            ],
        )),
        StandardGate::RZZ if qubits.len() == 2 => Some(rotation_or_product(
            &params[0],
            smallvec![
                AxisTerm {
                    qubit: qubits[0],
                    axis: Axis::Pauli(PauliAxis::Z),
                },
                AxisTerm {
                    qubit: qubits[1],
                    axis: Axis::Pauli(PauliAxis::Z),
                }
            ],
        )),
        StandardGate::RZX if qubits.len() == 2 => Some(rotation_or_product(
            &params[0],
            smallvec![
                AxisTerm {
                    qubit: qubits[0],
                    axis: Axis::Pauli(PauliAxis::Z),
                },
                AxisTerm {
                    qubit: qubits[1],
                    axis: Axis::Pauli(PauliAxis::X),
                }
            ],
        )),
        StandardGate::RXY if qubits.len() == 1 => Some(rotation_or_product(
            &params[0],
            smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Planar(params[1].clone()),
            }],
        )),
        StandardGate::XY if qubits.len() == 1 => Some(AxisOp {
            terms: smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Planar(params[0].clone()),
            }],
            kind: AxisOpKind::Product,
        }),
        StandardGate::XY2P | StandardGate::XY2M if qubits.len() == 1 => Some(AxisOp {
            terms: smallvec![AxisTerm {
                qubit: qubits[0],
                axis: Axis::Planar(params[0].clone()),
            }],
            kind: AxisOpKind::Rotation,
        }),
        _ => None,
    }
}

fn rotation_or_product(angle: &Parameter, terms: AxisString) -> AxisOp {
    // A zero rotation has empty support for commutation purposes.  A PI
    // rotation is product-valued up to global phase, so odd Pauli parity can
    // still produce a valid global-phase proof.
    if angle_eq_mod(angle, 0.0, TWO_PI) {
        return AxisOp {
            terms: SmallVec::new(),
            kind: AxisOpKind::Rotation,
        };
    }

    AxisOp {
        terms,
        kind: if angle_eq_mod(angle, PI, TWO_PI) {
            AxisOpKind::Product
        } else {
            AxisOpKind::Rotation
        },
    }
}

/// Attempts commutation proofs involving controlled single-axis operations.
///
/// Controlled gates commute with compatible operations on controls (Z axis) and
/// with compatible operations on their target axis.  Cases where a control is
/// also the other operation's target are left to later proof sources because
/// the control predicate and target action no longer separate cleanly.
fn controlled_commutation(
    lhs_inst: &Instruction,
    lhs_qubits: &[Qubit],
    lhs_params: &[Parameter],
    rhs_inst: &Instruction,
    rhs_qubits: &[Qubit],
    rhs_params: &[Parameter],
) -> CommutationResult {
    let lhs_controlled = controlled_axis_for_instruction(lhs_inst, lhs_qubits, lhs_params);
    let rhs_controlled = controlled_axis_for_instruction(rhs_inst, rhs_qubits, rhs_params);

    match (lhs_controlled, rhs_controlled) {
        (Some(lhs), Some(rhs)) => controlled_pair_commutation(&lhs, &rhs),
        (Some(controlled), None) => {
            let rhs_axis = axis_op_for_instruction(rhs_inst, rhs_qubits, rhs_params)?;
            controlled_axis_op_commutation(&controlled, &rhs_axis)
        }
        (None, Some(controlled)) => {
            let lhs_axis = axis_op_for_instruction(lhs_inst, lhs_qubits, lhs_params)?;
            controlled_axis_op_commutation(&controlled, &lhs_axis)
        }
        (None, None) => None,
    }
}

/// Checks two controlled-axis operations against each other.
///
/// If neither operation targets the other's control qubits, then distinct
/// targets commute.  Shared targets require the same target axis.
fn controlled_pair_commutation(lhs: &ControlledAxis, rhs: &ControlledAxis) -> CommutationResult {
    if lhs.controls.contains(&rhs.target) || rhs.controls.contains(&lhs.target) {
        return None;
    }

    if lhs.target == rhs.target {
        if axis_relation(&lhs.axis, &rhs.axis) == AxisRelation::Same {
            return Some(Commutation::Exact);
        }
        return None;
    }

    Some(Commutation::Exact)
}

/// Checks a controlled-axis operation against a non-controlled axis operation.
///
/// The non-controlled operation may touch controls only through Z-like axes and
/// may touch the target only through the controlled target axis.
fn controlled_axis_op_commutation(
    controlled: &ControlledAxis,
    axis_op: &AxisOp,
) -> CommutationResult {
    for term in &axis_op.terms {
        if controlled.controls.contains(&term.qubit)
            && axis_relation(&term.axis, &Axis::Pauli(PauliAxis::Z)) != AxisRelation::Same
        {
            return None;
        }

        if term.qubit == controlled.target
            && axis_relation(&term.axis, &controlled.axis) != AxisRelation::Same
        {
            return None;
        }
    }

    Some(Commutation::Exact)
}

/// Extracts a controlled-axis model from a supported instruction application.
fn controlled_axis_for_instruction(
    instruction: &Instruction,
    qubits: &[Qubit],
    params: &[Parameter],
) -> Option<ControlledAxis> {
    match instruction {
        Instruction::Standard(gate) => controlled_axis_for_standard(*gate, qubits),
        Instruction::McGate(gate) => controlled_axis_for_mc_gate(gate, qubits, params),
        _ => None,
    }
}

/// Extracts the controlled-axis model for built-in controlled standard gates.
fn controlled_axis_for_standard(gate: StandardGate, qubits: &[Qubit]) -> Option<ControlledAxis> {
    match gate {
        StandardGate::CX if qubits.len() == 2 => Some(ControlledAxis {
            controls: smallvec![qubits[0]],
            target: qubits[1],
            axis: Axis::Pauli(PauliAxis::X),
        }),
        StandardGate::CY if qubits.len() == 2 => Some(ControlledAxis {
            controls: smallvec![qubits[0]],
            target: qubits[1],
            axis: Axis::Pauli(PauliAxis::Y),
        }),
        StandardGate::CZ if qubits.len() == 2 => Some(ControlledAxis {
            controls: smallvec![qubits[0]],
            target: qubits[1],
            axis: Axis::Pauli(PauliAxis::Z),
        }),
        StandardGate::CCX if qubits.len() == 3 => Some(ControlledAxis {
            controls: smallvec![qubits[0], qubits[1]],
            target: qubits[2],
            axis: Axis::Pauli(PauliAxis::X),
        }),
        StandardGate::CRX if qubits.len() == 2 => Some(ControlledAxis {
            controls: smallvec![qubits[0]],
            target: qubits[1],
            axis: Axis::Pauli(PauliAxis::X),
        }),
        StandardGate::CRY if qubits.len() == 2 => Some(ControlledAxis {
            controls: smallvec![qubits[0]],
            target: qubits[1],
            axis: Axis::Pauli(PauliAxis::Y),
        }),
        StandardGate::CRZ if qubits.len() == 2 => Some(ControlledAxis {
            controls: smallvec![qubits[0]],
            target: qubits[1],
            axis: Axis::Pauli(PauliAxis::Z),
        }),
        _ => None,
    }
}

/// Extracts a controlled-axis model from an [`MCGate`].
///
/// Additional controls are prepended to the base gate's own controls.  When the
/// base gate is not itself controlled, it must reduce to a single-axis target
/// operation; multi-qubit base generators are outside this model.
fn controlled_axis_for_mc_gate(
    gate: &MCGate,
    qubits: &[Qubit],
    params: &[Parameter],
) -> Option<ControlledAxis> {
    let added_controls = gate
        .num_qubits()
        .checked_sub(gate.base_gate().num_qubits())?;
    if added_controls == 0 {
        return controlled_axis_for_standard(*gate.base_gate(), qubits);
    }

    let base_qubits = &qubits[added_controls..];
    let mut controlled =
        if let Some(controlled) = controlled_axis_for_standard(*gate.base_gate(), base_qubits) {
            controlled
        } else {
            let axis_op = axis_op_for_standard(*gate.base_gate(), base_qubits, params)?;
            if axis_op.terms.len() != 1 {
                return None;
            }
            let term = axis_op.terms[0].clone();
            ControlledAxis {
                controls: SmallVec::new(),
                target: term.qubit,
                axis: term.axis,
            }
        };
    for &control in qubits[..added_controls].iter().rev() {
        controlled.controls.insert(0, control);
    }
    Some(controlled)
}

/// Handles pair-order-independent commutation facts for selected two-qubit gates.
///
/// These families are symmetric on the unordered qubit pair, so operations on
/// `(a, b)` and `(b, a)` can still be recognized as acting on the same support.
fn special_two_qubit_commutation(
    lhs_inst: &Instruction,
    lhs_qubits: &[Qubit],
    rhs_inst: &Instruction,
    rhs_qubits: &[Qubit],
) -> CommutationResult {
    let lhs = two_qubit_family(lhs_inst, lhs_qubits)?;
    let rhs = two_qubit_family(rhs_inst, rhs_qubits)?;
    if lhs.0 != rhs.0 {
        return None;
    }

    match (lhs.1, rhs.1) {
        (TwoQubitFamily::Swap, TwoQubitFamily::Swap)
        | (TwoQubitFamily::Fsim, TwoQubitFamily::Fsim)
        | (TwoQubitFamily::Swap, TwoQubitFamily::Fsim)
        | (TwoQubitFamily::Fsim, TwoQubitFamily::Swap)
        | (TwoQubitFamily::Swap, TwoQubitFamily::SymmetricPauliRotation)
        | (TwoQubitFamily::SymmetricPauliRotation, TwoQubitFamily::Swap)
        | (TwoQubitFamily::Swap, TwoQubitFamily::SymmetricDiagonal)
        | (TwoQubitFamily::SymmetricDiagonal, TwoQubitFamily::Swap)
        | (TwoQubitFamily::Fsim, TwoQubitFamily::SymmetricDiagonal)
        | (TwoQubitFamily::SymmetricDiagonal, TwoQubitFamily::Fsim) => Some(Commutation::Exact),
        _ => None,
    }
}

/// Classifies a supported two-qubit instruction by unordered qubit pair and family.
fn two_qubit_family(
    instruction: &Instruction,
    qubits: &[Qubit],
) -> Option<((Qubit, Qubit), TwoQubitFamily)> {
    let gate = match instruction {
        Instruction::Standard(gate) => *gate,
        Instruction::McGate(gate) if gate.num_qubits() == gate.base_gate().num_qubits() => {
            *gate.base_gate()
        }
        _ => return None,
    };
    if qubits.len() != 2 {
        return None;
    }

    let pair = if qubits[0] <= qubits[1] {
        (qubits[0], qubits[1])
    } else {
        (qubits[1], qubits[0])
    };
    let family = match gate {
        StandardGate::SWAP => TwoQubitFamily::Swap,
        StandardGate::FSIM => TwoQubitFamily::Fsim,
        StandardGate::RXX | StandardGate::RYY => TwoQubitFamily::SymmetricPauliRotation,
        StandardGate::RZZ | StandardGate::CZ => TwoQubitFamily::SymmetricDiagonal,
        _ => return None,
    };
    Some((pair, family))
}

/// Determines whether two axes are equal, anti-commuting, or not provably related.
///
/// Planar axes are compared modulo `PI`; a shift by `PI / 2` proves
/// orthogonality.  Symbolic expressions that cannot be reduced to one of those
/// relationships are marked unknown.
fn axis_relation(lhs: &Axis, rhs: &Axis) -> AxisRelation {
    match (lhs, rhs) {
        (Axis::Pauli(lhs), Axis::Pauli(rhs)) if lhs == rhs => AxisRelation::Same,
        (Axis::Pauli(_), Axis::Pauli(_)) => AxisRelation::Orthogonal,
        (Axis::Planar(_), Axis::Pauli(PauliAxis::Z))
        | (Axis::Pauli(PauliAxis::Z), Axis::Planar(_)) => AxisRelation::Orthogonal,
        (Axis::Planar(phi), Axis::Pauli(axis)) | (Axis::Pauli(axis), Axis::Planar(phi)) => {
            match *axis {
                PauliAxis::X if angle_eq_mod(phi, 0.0, PI) => AxisRelation::Same,
                PauliAxis::X if angle_eq_mod(phi, FRAC_PI_2, PI) => AxisRelation::Orthogonal,
                PauliAxis::Y if angle_eq_mod(phi, FRAC_PI_2, PI) => AxisRelation::Same,
                PauliAxis::Y if angle_eq_mod(phi, 0.0, PI) => AxisRelation::Orthogonal,
                PauliAxis::Z => AxisRelation::Orthogonal,
                _ => AxisRelation::Unknown,
            }
        }
        (Axis::Planar(lhs), Axis::Planar(rhs)) => {
            if lhs.provably_equal_modulo(rhs, &Parameter::from(PI), PARAMETER_EQ_TOLERANCE) {
                return AxisRelation::Same;
            }

            let rhs_plus_half_pi = rhs.clone() + FRAC_PI_2;
            if lhs.provably_equal_modulo(
                &rhs_plus_half_pi,
                &Parameter::from(PI),
                PARAMETER_EQ_TOLERANCE,
            ) {
                return AxisRelation::Orthogonal;
            }

            AxisRelation::Unknown
        }
    }
}

/// Returns whether `angle == expected (mod modulus)` is provable symbolically.
fn angle_eq_mod(angle: &Parameter, expected: f64, modulus: f64) -> bool {
    angle.provably_equal_modulo(
        &Parameter::from(expected),
        &Parameter::from(modulus),
        PARAMETER_EQ_TOLERANCE,
    )
}