fix(compile): canonicalize KAK Weyl boundary representatives

crates/binding-c/tests/compile_pass2_test.rs

@@ -22,6 +22,124 @@ fn zeroed<T>() -> T {
     unsafe { std::mem::zeroed() }
 }
 
+#[test]
+fn test_weyl_boundary_kak_synthesis_and_unitary_decomposition() {
+    use binding_c::circuit::{circuit_to_matrix, circuit_to_matrix_len};
+    use cqlib_core::circuit::StandardGate;
+    use ndarray::{Array2, linalg::kron};
+    // Reported boundary unitary in the row-major order accepted by the C ABI.
+    let flat = [
+        Complex64::new(-0.14755155196621783, 0.08535345227001921),
+        Complex64::new(-0.015112015359785388, -0.026124324914469213),
+        Complex64::new(-4.3021142204224816e-16, 0.06601980439593577),
+        Complex64::new(-0.9826870609532827, 2.7755575615628914e-16),
+        Complex64::new(-0.006212290419590961, -0.010739262071974568),
+        Complex64::new(-0.3589338503302527, 0.2076307761864266),
+        Complex64::new(0.9070757873156866, -5.412337245047638e-16),
+        Complex64::new(-6.245004513516506e-16, -0.07152302867497204),
+        Complex64::new(-0.9706933646378622, 2.7755575615628914e-17),
+        Complex64::new(2.6020852139652106e-18, -0.15594957493615416),
+        Complex64::new(0.03569705171577718, -0.061709927849413246),
+        Complex64::new(0.1457506851638859, 0.08431171332109447),
+        Complex64::new(7.112366251504909e-17, -0.16894909068671718),
+        Complex64::new(-0.8960049266516804, -3.608224830031759e-16),
+        Complex64::new(-0.35455306241796475, -0.20509664240745978),
+        Complex64::new(-0.014674445936026376, 0.02536789332498595),
+    ];
+    let matrix = Array2::from_shape_vec((4, 4), flat.to_vec()).unwrap();
+    let kak = kak_decompose(flat.as_ptr());
+    assert!(!kak.is_null());
+    let (mut a, mut b, mut c) = (0.0, 0.0, 0.0);
+    assert_eq!(
+        kak_decomposition_coordinates(kak, &mut a, &mut b, &mut c),
+        0
+    );
+    assert_eq!(a, std::f64::consts::FRAC_PI_4);
+    assert!(c > 0.0);
+    let mut local = Vec::new();
+    for tag in [
+        KAK_FACTOR_K1L,
+        KAK_FACTOR_K1R,
+        KAK_FACTOR_K2L,
+        KAK_FACTOR_K2R,
+    ] {
+        let mut values = vec![Complex64::new(0.0, 0.0); 4];
+        assert_eq!(
+            kak_decomposition_local_factor(kak, tag, values.as_mut_ptr()),
+            0
+        );
+        local.push(Array2::from_shape_vec((2, 2), values).unwrap());
+    }
+    let rxx = StandardGate::RXX.matrix(&[-2.0 * a]).unwrap().into_owned();
+    let ryy = StandardGate::RYY.matrix(&[-2.0 * b]).unwrap().into_owned();
+    let rzz = StandardGate::RZZ.matrix(&[-2.0 * c]).unwrap().into_owned();
+    let core = rxx.dot(&ryy).dot(&rzz);
+    let phase = Complex64::from_polar(1.0, kak_decomposition_global_phase(kak));
+    let reconstructed = kron(&local[0], &local[1])
+        .dot(&core)
+        .dot(&kron(&local[2], &local[3]));
+    assert!(
+        reconstructed
+            .iter()
+            .zip(matrix.iter())
+            .all(|(actual, expected)| (phase * actual - expected).norm() < 1e-12)
+    );
+    kak_decomposition_free(kak);
+
+    for basis in [
+        TWO_QUBIT_BASIS_CX,
+        TWO_QUBIT_BASIS_CY,
+        TWO_QUBIT_BASIS_CZ,
+        TWO_QUBIT_BASIS_RZZ,
+        TWO_QUBIT_BASIS_PAULI_ROTATIONS,
+    ] {
+        let synthesis = synthesize_numeric_2q_unitary(flat.as_ptr(), 1, 0, basis);
+        assert!(!synthesis.is_null(), "basis={basis}");
+        assert!(two_qubit_synthesis_num_operations(synthesis) > 0);
+        assert!(two_qubit_synthesis_global_phase(synthesis).is_finite());
+        two_qubit_synthesis_free(synthesis);
+    }
+
+    let source = circuit_new(2);
+    let targets = [1, 0];
+    assert_eq!(
+        circuit_unitary(
+            source,
+            name("boundary").as_ptr(),
+            2,
+            flat.as_ptr(),
+            targets.as_ptr(),
+            2
+        ),
+        0
+    );
+    let output = decompose_unitaries(source, decompose_unitary_config_default());
+    assert!(!output.is_null());
+    let matrix_of = |circuit| {
+        let len = circuit_to_matrix_len(circuit, std::ptr::null(), 0);
+        let mut values = vec![0.0; 2 * len];
+        assert_eq!(
+            circuit_to_matrix(
+                circuit,
+                std::ptr::null(),
+                0,
+                values.as_mut_ptr(),
+                values.len()
+            ),
+            len
+        );
+        values
+    };
+    assert!(
+        matrix_of(output)
+            .iter()
+            .zip(matrix_of(source))
+            .all(|(a, b)| (a - b).abs() < 1e-9)
+    );
+    circuit_free(output);
+    circuit_free(source);
+}
+
 // =====  Decomposition and resynthesis  =====
 
 #[test]

crates/binding-python/tests/compile/transform/test_decompose.py

@@ -51,6 +51,141 @@ def operation_names(circuit: Circuit) -> list[str]:
     return [operation.instruction.instruction.name for operation in circuit.operations]
 
 
+@pytest.fixture
+def weyl_boundary_matrix() -> np.ndarray:
+    # Preserve the reported matrix, including its floating-point roundoff.
+    return np.array(
+        [
+            [
+                complex(-0.14755155196621783, 0.08535345227001921),
+                complex(-0.015112015359785388, -0.026124324914469213),
+                complex(-4.3021142204224816e-16, 0.06601980439593577),
+                complex(-0.9826870609532827, 2.7755575615628914e-16),
+            ],
+            [
+                complex(-0.006212290419590961, -0.010739262071974568),
+                complex(-0.3589338503302527, 0.2076307761864266),
+                complex(0.9070757873156866, -5.412337245047638e-16),
+                complex(-6.245004513516506e-16, -0.07152302867497204),
+            ],
+            [
+                complex(-0.9706933646378622, 2.7755575615628914e-17),
+                complex(2.6020852139652106e-18, -0.15594957493615416),
+                complex(0.03569705171577718, -0.061709927849413246),
+                complex(0.1457506851638859, 0.08431171332109447),
+            ],
+            [
+                complex(7.112366251504909e-17, -0.16894909068671718),
+                complex(-0.8960049266516804, -3.608224830031759e-16),
+                complex(-0.35455306241796475, -0.20509664240745978),
+                complex(-0.014674445936026376, 0.02536789332498595),
+            ],
+        ],
+        dtype=np.complex128,
+    )
+
+
+def reconstructed_kak(decomposition: KakDecomposition) -> np.ndarray:
+    core = np.eye(4, dtype=np.complex128)
+    paulis = (
+        np.array([[0, 1], [1, 0]], dtype=np.complex128),
+        np.array([[0, -1j], [1j, 0]], dtype=np.complex128),
+        np.diag([1, -1]),
+    )
+    for angle, pauli in zip(
+        (decomposition.a, decomposition.b, decomposition.c), paulis, strict=True
+    ):
+        core = core @ (
+            np.cos(angle) * np.eye(4) + 1j * np.sin(angle) * np.kron(pauli, pauli)
+        )
+    return (
+        np.exp(1j * decomposition.global_phase)
+        * np.kron(decomposition.k1l, decomposition.k1r)
+        @ core
+        @ np.kron(decomposition.k2l, decomposition.k2r)
+    )
+
+
+@pytest.mark.parametrize("phase", [0.0, -1.1, np.pi / 4, np.pi / 2, np.pi, 2.7])
+def test_weyl_boundary_kak_preserves_phase_and_su2_factors(
+    weyl_boundary_matrix: np.ndarray, phase: float
+) -> None:
+    source = np.exp(1j * phase) * weyl_boundary_matrix
+    decomposition = kak_decompose(source)
+    assert decomposition.a == np.pi / 4
+    assert decomposition.c > 0
+    np.testing.assert_allclose(
+        reconstructed_kak(decomposition), source, atol=1e-12, rtol=0
+    )
+    for factor in (
+        decomposition.k1l,
+        decomposition.k1r,
+        decomposition.k2l,
+        decomposition.k2r,
+    ):
+        np.testing.assert_allclose(
+            factor.conj().T @ factor, np.eye(2), atol=1e-12, rtol=0
+        )
+        np.testing.assert_allclose(np.linalg.det(factor), 1, atol=1e-12, rtol=0)
+
+
+@pytest.mark.parametrize(
+    "basis",
+    [
+        TwoQubitUnitaryDecomposeBasis.cx(),
+        TwoQubitUnitaryDecomposeBasis.cy(),
+        TwoQubitUnitaryDecomposeBasis.cz(),
+        TwoQubitUnitaryDecomposeBasis.rzz(),
+        TwoQubitUnitaryDecomposeBasis.pauli_rotations(),
+    ],
+)
+@pytest.mark.parametrize("qubits", [[0, 1], [1, 0]])
+def test_weyl_boundary_synthesis_and_cached_unitary_lowering(
+    weyl_boundary_matrix: np.ndarray,
+    basis: TwoQubitUnitaryDecomposeBasis,
+    qubits: list[int],
+) -> None:
+    matrix = np.exp(0.37j) * weyl_boundary_matrix
+    gate = UnitaryGate("boundary", 2).with_matrix(matrix)
+    source = Circuit(2)
+    source.append_unitary_gate(gate, qubits)
+    synthesis = synthesize_numeric_2q_unitary(matrix, *qubits, basis)
+    output = Circuit.from_operations([Qubit(0), Qubit(1)], synthesis.operations)
+    output.set_global_phase(synthesis.global_phase)
+    np.testing.assert_allclose(
+        output.to_matrix(), source.to_matrix(), atol=1e-9, rtol=0
+    )
+
+    source.append_unitary_gate(gate, qubits)
+    config = UnitaryDecomposeConfig(two_qubit_basis=basis)
+    result, stats = decompose_unitaries_with_rule_stats(source, config)
+    np.testing.assert_allclose(
+        result.circuit.to_matrix(), source.to_matrix(), atol=1e-9, rtol=0
+    )
+    assert (stats.hits, stats.misses, stats.inserts) == (1, 1, 1)
+
+
+def test_weyl_boundary_interior_negative_c_is_not_rejected() -> None:
+    identity = np.eye(4, dtype=np.complex128)
+    x = np.array([[0, 1], [1, 0]], dtype=np.complex128)
+    y = np.array([[0, -1j], [1j, 0]], dtype=np.complex128)
+    z = np.diag([1, -1])
+    for distance in (1e-12, 1e-10, 5e-9, 2e-8):
+        source = identity.copy()
+        for angle, pauli in zip(
+            (np.pi / 4 - distance, 0.2, -0.03), (x, y, z), strict=True
+        ):
+            source = source @ (
+                np.cos(angle) * identity + 1j * np.sin(angle) * np.kron(pauli, pauli)
+            )
+        decomposition = kak_decompose(source)
+        assert decomposition.a < np.pi / 4
+        assert decomposition.c < 0
+        np.testing.assert_allclose(
+            reconstructed_kak(decomposition), source, atol=1e-12, rtol=0
+        )
+
+
 def test_decompose_module_and_types_are_registered() -> None:
     assert "cqlib._native.compile.transform.decompose" in sys.modules
     assert decompose.expand_definitions is expand_definitions
@@ -259,23 +394,7 @@ def test_kak_decomposition_reconstructs_matrix_and_returns_owned_arrays() -> Non
     )
 
     decomposition = kak_decompose(source)
-    identity = np.eye(4, dtype=np.complex128)
-    x = np.array([[0, 1], [1, 0]], dtype=np.complex128)
-    y = np.array([[0, -1j], [1j, 0]], dtype=np.complex128)
-    z = np.array([[1, 0], [0, -1]], dtype=np.complex128)
-    xx, yy, zz = np.kron(x, x), np.kron(y, y), np.kron(z, z)
-    cartan = (
-        np.cos(decomposition.a) * identity + 1j * np.sin(decomposition.a) * xx
-    ) @ (np.cos(decomposition.b) * identity + 1j * np.sin(decomposition.b) * yy)
-    cartan = cartan @ (
-        np.cos(decomposition.c) * identity + 1j * np.sin(decomposition.c) * zz
-    )
-    reconstructed = (
-        np.exp(1j * decomposition.global_phase)
-        * np.kron(decomposition.k1l, decomposition.k1r)
-        @ cartan
-        @ np.kron(decomposition.k2l, decomposition.k2r)
-    )
+    reconstructed = reconstructed_kak(decomposition)
 
     np.testing.assert_allclose(reconstructed, source, atol=1e-7)
     local = decomposition.k1l

crates/cqlib-core/src/compile/test_utils.rs

@@ -25,6 +25,38 @@ use std::sync::Arc;
 
 const EPSILON: f64 = 1e-9;
 
+/// The reported unitary whose Cartan coordinates lie on the a=pi/4 wall.
+pub(crate) fn kak_boundary_matrix() -> ndarray::Array2<num_complex::Complex64> {
+    use num_complex::Complex64;
+
+    ndarray::array![
+        [
+            Complex64::new(-0.14755155196621783, 0.08535345227001921),
+            Complex64::new(-0.015112015359785388, -0.026124324914469213),
+            Complex64::new(-4.3021142204224816e-16, 0.06601980439593577),
+            Complex64::new(-0.9826870609532827, 2.7755575615628914e-16),
+        ],
+        [
+            Complex64::new(-0.006212290419590961, -0.010739262071974568),
+            Complex64::new(-0.3589338503302527, 0.2076307761864266),
+            Complex64::new(0.9070757873156866, -5.412337245047638e-16),
+            Complex64::new(-6.245004513516506e-16, -0.07152302867497204),
+        ],
+        [
+            Complex64::new(-0.9706933646378622, 2.7755575615628914e-17),
+            Complex64::new(2.6020852139652106e-18, -0.15594957493615416),
+            Complex64::new(0.03569705171577718, -0.061709927849413246),
+            Complex64::new(0.1457506851638859, 0.08431171332109447),
+        ],
+        [
+            Complex64::new(7.112366251504909e-17, -0.16894909068671718),
+            Complex64::new(-0.8960049266516804, -3.608224830031759e-16),
+            Complex64::new(-0.35455306241796475, -0.20509664240745978),
+            Complex64::new(-0.014674445936026376, 0.02536789332498595),
+        ],
+    ]
+}
+
 pub(super) fn build_device_synthesis_context(
     device: &Device,
     circuit: &Circuit,