diff --git a/crates/binding-c/tests/compile_pass2_test.rs b/crates/binding-c/tests/compile_pass2_test.rs index 30edee93..1f2c4ffa 100644 --- a/crates/binding-c/tests/compile_pass2_test.rs +++ b/crates/binding-c/tests/compile_pass2_test.rs @@ -22,6 +22,124 @@ fn zeroed() -> 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] diff --git a/crates/binding-python/tests/compile/transform/test_decompose.py b/crates/binding-python/tests/compile/transform/test_decompose.py index 106606f0..b3ee211b 100644 --- a/crates/binding-python/tests/compile/transform/test_decompose.py +++ b/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 diff --git a/crates/cqlib-core/src/compile/test_utils.rs b/crates/cqlib-core/src/compile/test_utils.rs index 295a0cbd..d112da1e 100644 --- a/crates/cqlib-core/src/compile/test_utils.rs +++ b/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 { + 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, diff --git a/crates/cqlib-core/src/compile/transform/decompose/unitary/decompose.rs b/crates/cqlib-core/src/compile/transform/decompose/unitary/decompose.rs index 6bf63b8e..4df35a39 100644 --- a/crates/cqlib-core/src/compile/transform/decompose/unitary/decompose.rs +++ b/crates/cqlib-core/src/compile/transform/decompose/unitary/decompose.rs @@ -731,6 +731,167 @@ mod tests { } } + #[test] + fn weyl_boundary_unitary_gate_decomposition_preserves_phase_and_device_constraints() { + use crate::compile::test_utils::{build_device_synthesis_context, kak_boundary_matrix}; + use crate::compile::transform::decompose::unitary::DeviceSynthesisPlacement; + use crate::device::Device; + + let gate = UnitaryGate::new("boundary", 2, 0) + .with_matrix(kak_boundary_matrix()) + .unwrap(); + let mut source = Circuit::new(2); + source + .unitary(gate.clone(), vec![Qubit::new(1), Qubit::new(0)]) + .unwrap(); + source + .unitary(gate, vec![Qubit::new(1), Qubit::new(0)]) + .unwrap(); + source.set_global_phase(Parameter::from(0.37)); + let expected = circuit_to_matrix(&source, None).unwrap(); + + for gate in [ + StandardGate::CX, + StandardGate::CY, + StandardGate::CZ, + StandardGate::RZZ, + ] { + let config = config_for_native_2q(gate); + let (generic, stats) = + decompose_unitaries_with_rule_stats(&source, config.clone()).unwrap(); + assert_eq!(stats.hits, 1); + assert_eq!(stats.misses, 1); + assert_abs_diff_eq!( + expected, + circuit_to_matrix(&generic.circuit, None).unwrap(), + epsilon = 1e-9 + ); + + let device = Device::bidirectional_line("boundary", 2) + .unwrap() + .with_native_gates( + config + .two_qubit_target + .native_1q() + .iter() + .copied() + .chain([gate]) + .map(Instruction::Standard) + .collect(), + ) + .unwrap(); + // Matrix-backed gates are synthesized before physical routing. + let context = build_device_synthesis_context( + &device, + &source, + DeviceSynthesisPlacement::PreLayoutEnvelope, + ) + .unwrap(); + let result = DecomposeUnitaries::new_device_aware(config.clone(), context) + .transform_resolved(&source, None) + .unwrap_or_else(|error| panic!("gate={gate:?}: {error}")); + assert!(result.changed); + assert!( + result + .circuit + .operations() + .iter() + .all(|op| match op.instruction { + Instruction::Standard(g) => g.num_qubits() == 1 || g == gate, + _ => false, + }) + ); + assert_abs_diff_eq!( + expected, + circuit_to_matrix(&result.circuit, None).unwrap(), + epsilon = 1e-9 + ); + } + } + + #[test] + fn weyl_boundary_does_not_bypass_unsupported_device_interactions() { + use crate::compile::test_utils::{build_device_synthesis_context, kak_boundary_matrix}; + use crate::compile::transform::decompose::unitary::DeviceSynthesisPlacement; + use crate::device::Device; + let mut source = Circuit::new(2); + source + .unitary( + UnitaryGate::new("boundary", 2, 0) + .with_matrix(kak_boundary_matrix()) + .unwrap(), + vec![Qubit::new(0), Qubit::new(1)], + ) + .unwrap(); + let device = Device::bidirectional_line("no-entangler", 2) + .unwrap() + .with_native_gates(vec![Instruction::Standard(StandardGate::U)]) + .unwrap(); + let context = build_device_synthesis_context( + &device, + &source, + DeviceSynthesisPlacement::PreLayoutEnvelope, + ) + .unwrap(); + let error = DecomposeUnitaries::new_device_aware(Default::default(), context) + .transform_resolved(&source, None) + .unwrap_err(); + assert!( + error + .to_string() + .contains("no exact device-aware two-qubit synthesis candidate") + ); + } + + #[test] + fn weyl_boundary_phase_stays_inside_classical_control_body() { + use crate::compile::test_utils::kak_boundary_matrix; + let gate = UnitaryGate::new("boundary", 2, 0) + .with_matrix(kak_boundary_matrix()) + .unwrap(); + let mut expected_body = Circuit::new(2); + expected_body + .unitary(gate.clone(), vec![Qubit::new(0), Qubit::new(1)]) + .unwrap(); + let mut source = Circuit::new(2); + source + .if_else( + ClassicalExpr::bool_literal(true), + |body| { + body.append( + Instruction::UnitaryGate(Box::new(gate)), + [Qubit::new(0), Qubit::new(1)], + [], + None, + ) + }, + |_| Ok(()), + ) + .unwrap(); + let output = decompose_unitaries(&source, Default::default()).unwrap(); + assert!(output.global_phase().is_zero()); + let Instruction::ClassicalControl(ClassicalControlOp::If(control)) = + &output.operations()[0].instruction + else { + panic!("expected if body") + }; + let mut body = Circuit::new(2); + for op in control.then_body().operations() { + let params = op + .params + .iter() + .map(|param| output.parameter_value(param).unwrap()) + .collect::>(); + body.append(op.instruction.clone(), op.qubits.clone(), params, None) + .unwrap(); + } + assert_abs_diff_eq!( + circuit_to_matrix(&body, None).unwrap(), + circuit_to_matrix(&expected_body, None).unwrap(), + epsilon = 1e-9 + ); + } + fn operation_parameter(circuit: &Circuit, operation: &Operation, position: usize) -> Parameter { match operation.params.get(position) { Some(CircuitParam::Fixed(value)) => Parameter::from(*value), diff --git a/crates/cqlib-core/src/compile/transform/decompose/unitary/two_qubit_kak.rs b/crates/cqlib-core/src/compile/transform/decompose/unitary/two_qubit_kak.rs index c57807f3..3c6f7ebe 100644 --- a/crates/cqlib-core/src/compile/transform/decompose/unitary/two_qubit_kak.rs +++ b/crates/cqlib-core/src/compile/transform/decompose/unitary/two_qubit_kak.rs @@ -24,7 +24,9 @@ //! //! The returned coordinates satisfy the Weyl chamber convention //! `pi/4 >= a >= b >= |c|`, with the additional boundary convention that -//! `c >= 0` when `a = pi/4`. The compiler-facing lowering lives in +//! `c >= 0` when `a = pi/4`. Coordinates within floating-point roundoff of +//! that wall are normalized to `a = pi/4` before selecting its representative. +//! The compiler-facing lowering lives in //! `unitary_2q.rs`; this file owns only the numerical decomposition contract. //! //! The implementation validates finite 4x4 input, the intermediate symmetric @@ -49,6 +51,10 @@ const UNITARY_EPS: f64 = 1e-8; const RECONSTRUCTION_EPS: f64 = 1e-7; const AUTONNE_EPS: f64 = 10.0 * UNITARY_EPS; const KRON_RANK_EPS: f64 = 1e-8; +// Angle extraction and folding accumulate a small number of f64 rounding +// errors. This allowance is deliberately independent of input unitarity and +// reconstruction tolerances: genuine interior coordinates may have c < 0. +const WEYL_BOUNDARY_EPS: f64 = 32.0 * f64::EPSILON; static B_RAW: [[Complex64; 4]; 4] = [ [c(1.0, 0.0), c(0.0, 1.0), c(0.0, 0.0), c(0.0, 0.0)], @@ -301,7 +307,7 @@ pub fn kak_decompose(matrix: &Array2) -> Result) -> Result Result<(), CompilerError> { let eps = UNITARY_EPS; if a < -eps || b < -eps { @@ -337,9 +364,11 @@ fn validate_weyl_chamber(a: f64, b: f64, c: f64) -> Result<(), CompilerError> { c.abs() ))); } - if (a - FRAC_PI_4).abs() < eps && c < -eps { + // Boundary normalization sets a exactly to this constant. An epsilon band + // here would incorrectly reject valid interior points with negative c. + if a == FRAC_PI_4 && c < 0.0 { return Err(kak_failed(format!( - "Weyl chamber violation: a=pi/4 but c={c} < 0" + "Weyl chamber violation: a={a} on pi/4 wall but c={c} < 0" ))); } Ok(()) @@ -768,6 +797,7 @@ fn decompose_two_qubit_product_gate( #[cfg(test)] mod tests { use super::*; + use crate::compile::test_utils::kak_boundary_matrix; use approx::assert_abs_diff_eq; use ndarray::linalg::kron; @@ -813,6 +843,76 @@ mod tests { decomp.c, decomp.b ); + if decomp.a == FRAC_PI_4 { + assert!( + decomp.c >= 0.0, + "negative c={} on the a=pi/4 wall", + decomp.c + ); + } + } + + #[test] + fn weyl_boundary_reported_matrix_preserves_full_reconstruction() { + for phase in [0.0, -1.1, FRAC_PI_4, FRAC_PI_2, PI, 2.7] { + let matrix = + kak_boundary_matrix().mapv(|value| value * Complex64::from_polar(1.0, phase)); + let decomp = kak_decompose(&matrix).unwrap(); + assert_eq!(decomp.a, FRAC_PI_4); + assert!(decomp.c > 0.0); + assert_weyl_chamber(&decomp); + validate_local_su2(&decomp).unwrap(); + assert!(max_entry_diff(&matrix, &reconstructed_from_kak(&decomp)) < 1e-12); + } + } + + #[test] + fn weyl_boundary_distinguishes_roundoff_from_interior_coordinates() { + let left = kron( + &StandardGate::RY.matrix(&[0.7]).unwrap(), + &StandardGate::RZ.matrix(&[-0.4]).unwrap(), + ); + let right = kron( + &StandardGate::RX.matrix(&[-1.1]).unwrap(), + &StandardGate::RY.matrix(&[0.3]).unwrap(), + ); + for delta in [0.0, 4.0 * f64::EPSILON, 1e-12, 1e-10, 5e-9, 2e-8] { + for sign in [-1.0, 1.0] { + for c in [-0.03, 0.03, -1e-12, 1e-12] { + let a = FRAC_PI_4 + sign * delta; + let rxx = StandardGate::RXX.matrix(&[-2.0 * a]).unwrap(); + let ryy = StandardGate::RYY.matrix(&[-0.4]).unwrap(); + let rzz = StandardGate::RZZ.matrix(&[-2.0 * c]).unwrap(); + let matrix = left + .dot(&rxx.dot(&ryy.dot(&rzz.dot(&right)))) + .mapv(|value| value * Complex64::from_polar(1.0, 0.37)); + let decomp = kak_decompose(&matrix).unwrap(); + assert_weyl_chamber(&decomp); + validate_local_su2(&decomp).unwrap(); + if delta <= 4.0 * f64::EPSILON { + assert_eq!(decomp.a, FRAC_PI_4); + assert!(decomp.c >= 0.0); + } else { + assert!(decomp.a < FRAC_PI_4); + let expected_c = if sign < 0.0 { c } else { -c }; + assert_eq!(decomp.c.is_sign_negative(), expected_c.is_sign_negative()); + assert!((decomp.a - (FRAC_PI_4 - delta)).abs() < 1e-14); + } + assert!(max_entry_diff(&matrix, &reconstructed_from_kak(&decomp)) < 1e-12); + } + } + } + } + + #[test] + fn weyl_boundary_validator_accepts_negative_c_in_the_interior() { + for distance in [1e-12, 1e-10, 5e-9, 2e-8] { + validate_weyl_chamber(FRAC_PI_4 - distance, 0.2, -0.03).unwrap(); + } + assert!(validate_weyl_chamber(FRAC_PI_4, 0.2, -0.03).is_err()); + assert!(validate_weyl_chamber(FRAC_PI_4, 0.2, -1e-12).is_err()); + assert!(validate_weyl_chamber(FRAC_PI_4 + 1e-6, 0.2, 0.03).is_err()); + assert!(validate_weyl_chamber(0.3, 0.2, 0.25).is_err()); } #[test] diff --git a/crates/cqlib-core/src/compile/transform/decompose/unitary/unitary_2q_test.rs b/crates/cqlib-core/src/compile/transform/decompose/unitary/unitary_2q_test.rs index bebe41fe..96056123 100644 --- a/crates/cqlib-core/src/compile/transform/decompose/unitary/unitary_2q_test.rs +++ b/crates/cqlib-core/src/compile/transform/decompose/unitary/unitary_2q_test.rs @@ -18,6 +18,66 @@ use super::unitary_2q::{ }; use crate::circuit::{Instruction, Qubit, StandardGate, ValueInstruction}; +#[test] +fn weyl_boundary_all_backends_preserve_phase_and_qubit_order() { + use super::unitary_2q::{TwoQubitSynthesisCandidate, synthesize_numeric_2q_unitary}; + use crate::compile::test_utils::kak_boundary_matrix; + use num_complex::Complex64; + for phase in [0.0, 0.37] { + let matrix = kak_boundary_matrix().mapv(|value| value * Complex64::from_polar(1.0, phase)); + for qubits in [ + [Qubit::new(0), Qubit::new(1)], + [Qubit::new(3), Qubit::new(1)], + ] { + for backend in [ + TwoQubitUnitaryDecomposeBasis::Cx, + TwoQubitUnitaryDecomposeBasis::Cy, + TwoQubitUnitaryDecomposeBasis::Cz, + TwoQubitUnitaryDecomposeBasis::Rzz, + TwoQubitUnitaryDecomposeBasis::PauliRotations, + ] { + let result = synthesize_numeric_2q_unitary(&matrix, qubits, backend).unwrap(); + let candidate = TwoQubitSynthesisCandidate { + backend, + operations: result.operations, + global_phase: result.global_phase, + cost: Default::default(), + }; + validate_numeric_2q_candidate(&matrix, qubits, &candidate).unwrap(); + assert!( + candidate + .operations + .iter() + .all(|op| op.qubits.iter().all(|q| qubits.contains(q))) + ); + assert!( + candidate + .operations + .iter() + .filter(|op| op.qubits.len() == 2) + .all(|op| { + let ValueInstruction::Instruction(Instruction::Standard(gate)) = + op.instruction + else { + return false; + }; + match backend { + TwoQubitUnitaryDecomposeBasis::Cx => gate == StandardGate::CX, + TwoQubitUnitaryDecomposeBasis::Cy => gate == StandardGate::CY, + TwoQubitUnitaryDecomposeBasis::Cz => gate == StandardGate::CZ, + TwoQubitUnitaryDecomposeBasis::Rzz => gate == StandardGate::RZZ, + TwoQubitUnitaryDecomposeBasis::PauliRotations => matches!( + gate, + StandardGate::RXX | StandardGate::RYY | StandardGate::RZZ + ), + } + }) + ); + } + } + } +} + fn count_cx(candidate: &super::unitary_2q::TwoQubitSynthesisCandidate) -> usize { candidate .operations diff --git a/crates/cqlib-core/src/compile/transform/resynthesis/resynthesizer_test.rs b/crates/cqlib-core/src/compile/transform/resynthesis/resynthesizer_test.rs index 7b7422ad..b98ea601 100644 --- a/crates/cqlib-core/src/compile/transform/resynthesis/resynthesizer_test.rs +++ b/crates/cqlib-core/src/compile/transform/resynthesis/resynthesizer_test.rs @@ -112,6 +112,71 @@ fn two_qubit_op_count(circuit: &Circuit) -> usize { .count() } +#[test] +fn weyl_boundary_resynthesis_improves_a_valid_interior_block() { + use crate::circuit::test_utils::assert_matrix_approx_eq; + use crate::circuit::{Parameter, circuit_to_matrix}; + + let q0 = Qubit::new(0); + let q1 = Qubit::new(1); + let mut circuit = Circuit::new(2); + circuit + .rxx(q0, q1, -2.0 * (std::f64::consts::FRAC_PI_4 - 5e-9)) + .unwrap(); + circuit.ryy(q0, q1, -0.4).unwrap(); + circuit.rzz(q0, q1, 0.06).unwrap(); + for _ in 0..2 { + circuit.cx(q0, q1).unwrap(); + circuit.cx(q0, q1).unwrap(); + } + circuit.set_global_phase(Parameter::from(0.37)); + let expected = circuit_to_matrix(&circuit, None).unwrap(); + let (generic, stats) = + resynthesize_two_qubit_blocks_with_cache_budget(&circuit, cx_config(), 4096).unwrap(); + assert!(generic.changed); + assert!(two_qubit_op_count(&generic.circuit) < two_qubit_op_count(&circuit)); + assert_eq!(stats.kak_failed_hits, 0); + assert_matrix_approx_eq( + &circuit_to_matrix(&generic.circuit, None).unwrap(), + &expected, + 1e-9, + ); + + let device = Device::bidirectional_line("boundary", 2) + .unwrap() + // Cost source and replacement on a device that can execute both. + .with_native_gates(vec![ + Instruction::Standard(StandardGate::U), + Instruction::Standard(StandardGate::CX), + Instruction::Standard(StandardGate::RXX), + Instruction::Standard(StandardGate::RYY), + Instruction::Standard(StandardGate::RZZ), + ]) + .unwrap(); + for placement in [ + DeviceSynthesisPlacement::PreLayoutEnvelope, + DeviceSynthesisPlacement::ExactPhysical, + ] { + let context = build_device_synthesis_context(&device, &circuit, placement).unwrap(); + let (result, stats) = resynthesize_two_qubit_blocks_with_device_cache_budget( + &circuit, + cx_config(), + context, + 4096, + ) + .unwrap(); + assert!(result.changed, "placement={placement:?}, stats={stats:?}"); + assert!(two_qubit_op_count(&result.circuit) < two_qubit_op_count(&circuit)); + assert_eq!(stats.kak_failed_hits, 0); + assert_matrix_approx_eq( + &circuit_to_matrix(&result.circuit, None).unwrap(), + &expected, + 1e-9, + ); + device.validate_circuit(&result.circuit).unwrap(); + } +} + #[test] fn cancels_adjacent_cx_pair() { let q0 = Qubit::new(0); diff --git a/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache.rs b/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache.rs index 7c6987b4..fa933682 100644 --- a/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache.rs +++ b/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache.rs @@ -52,7 +52,7 @@ use super::commutation::OperationView; use super::cost::ResynthesisCost; pub(super) const RESYNTHESIS_SYNTHESIS_CACHE_BUDGET: usize = 4096; -const SYNTHESIS_ALGORITHM_REVISION: u32 = 2; +const SYNTHESIS_ALGORITHM_REVISION: u32 = 3; #[derive(Clone, Debug, PartialEq, Eq)] struct SynthesisCacheNamespace { diff --git a/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache_test.rs b/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache_test.rs index 4712d112..bfdabcdc 100644 --- a/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache_test.rs +++ b/crates/cqlib-core/src/compile/transform/resynthesis/synthesis_cache_test.rs @@ -28,6 +28,20 @@ fn matrix_with(real: f64) -> Array2 { matrix } +#[test] +fn weyl_boundary_kak_cache_stores_a_successful_decomposition() { + let matrix = crate::compile::test_utils::kak_boundary_matrix(); + let mut cache = TwoQubitSynthesisCache::new_native_session(); + let first = cache.kak_decomposition(&matrix).unwrap().unwrap(); + let second = cache.kak_decomposition(&matrix).unwrap().unwrap(); + assert!(Arc::ptr_eq(&first, &second)); + assert_eq!(first.a, std::f64::consts::FRAC_PI_4); + assert!(first.c > 0.0); + assert_eq!(cache.stats().kak_misses, 1); + assert_eq!(cache.stats().kak_hits, 1); + assert_eq!(cache.stats().kak_failed_hits, 0); +} + #[test] fn exact_key_distinguishes_qargs_bits_and_signed_zero() { let q0 = Qubit::new(0); diff --git a/crates/cqlib-core/src/compile/workflow_test.rs b/crates/cqlib-core/src/compile/workflow_test.rs index 7a87dcd7..3b3e5d4f 100644 --- a/crates/cqlib-core/src/compile/workflow_test.rs +++ b/crates/cqlib-core/src/compile/workflow_test.rs @@ -54,6 +54,61 @@ fn compile_config(mode: CompileMode) -> CompileConfig { } } +#[test] +fn weyl_boundary_workflow_preserves_phase_for_logical_basis_and_device_targets() { + use crate::circuit::{circuit_to_matrix, test_utils::assert_matrix_approx_eq}; + use crate::compile::test_utils::kak_boundary_matrix; + let mut circuit = Circuit::new(2); + circuit + .unitary( + UnitaryGate::new("boundary", 2, 0) + .with_matrix(kak_boundary_matrix()) + .unwrap(), + vec![Qubit::new(1), Qubit::new(0)], + ) + .unwrap(); + circuit.set_global_phase(Parameter::from(-0.41)); + let expected = circuit_to_matrix(&circuit, None).unwrap(); + let native = vec![ + Instruction::Standard(StandardGate::U), + Instruction::Standard(StandardGate::CX), + ]; + let device = Device::bidirectional_line("boundary", 2) + .unwrap() + .with_native_gates(native.clone()) + .unwrap(); + for mode in [CompileMode::Normal, CompileMode::Enhanced] { + for target in [ + CompileTarget::Logical, + CompileTarget::Basis(native.clone()), + CompileTarget::Device(DeviceCompileTarget { + device: device.clone(), + initial_layout: Some(Layout::from_pairs(&[(0, 0), (1, 1)], 2).unwrap()), + seed: Some(41), + }), + ] { + let is_device = matches!(&target, CompileTarget::Device(_)); + let result = compile( + &circuit, + CompileConfig { + target, + ..compile_config(mode) + }, + ) + .unwrap(); + assert!(result.step_changed("decompose.unitary")); + assert_matrix_approx_eq( + &circuit_to_matrix(&result.circuit, None).unwrap(), + &expected, + 1e-9, + ); + if is_device { + device.validate_circuit(&result.circuit).unwrap(); + } + } + } +} + #[test] fn xy_basis_compiles_generic_u_to_three_rotations_in_both_modes() { let mut circuit = Circuit::new(1);