From 30f7e9dd096627ae5257c9778bf1dac141f9e9c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=81=A5=E5=81=A5?= Date: Thu, 8 Oct 2026 11:49:42 +0800 Subject: [PATCH 1/4] fix(qis): correct simulation edge cases and numerical errors Validate state dimensions and finite inputs before constructing simulators. Respect original qubit labels and identity-gate noise during noisy simulation. Stabilize Renyi entropy, validate real observable expectations and variance, retain small evolution coefficients, and preserve global phase during ZNE folding. Include regression tests and matching Rust, Python, and C API documentation. Validation: 16 core regression tests, 4 C regression tests, and 1200 Python QIS tests passed against the staged source snapshot. Rust formatting and Python formatting checks passed. --- crates/binding-c/src/qis/density_matrix.rs | 31 +- .../binding-c/src/qis/density_matrix_noise.rs | 9 +- crates/binding-c/src/qis/statevector.rs | 9 +- .../tests/qis_constructor_regressions.rs | 51 +++ crates/binding-c/tests/qis_p2_regressions.rs | 49 +++ crates/binding-python/cqlib/qis/entropy.pyi | 11 +- .../binding-python/cqlib/qis/hamiltonian.pyi | 16 +- crates/binding-python/cqlib/qis/pauli.pyi | 13 +- .../cqlib/qis/state/density_matrix.pyi | 6 +- .../cqlib/qis/state/density_matrix_noise.pyi | 6 +- .../cqlib/qis/state/statevector.pyi | 6 +- crates/binding-python/src/qis/entropy.rs | 11 +- crates/binding-python/src/qis/hamiltonian.rs | 16 + crates/binding-python/src/qis/pauli.rs | 15 +- .../src/qis/state/density_matrix.rs | 23 +- .../src/qis/state/density_matrix_noise.rs | 15 +- .../src/qis/state/statevector.rs | 14 +- .../tests/qis/test_p1_regressions.py | 124 ++++++++ .../tests/qis/test_p2_regressions.py | 240 ++++++++++++++ crates/cqlib-core/src/device/noise.rs | 2 +- .../src/error_mitigation/zne_mitigation.rs | 2 + crates/cqlib-core/src/qis/entropy.rs | 80 +++-- crates/cqlib-core/src/qis/evolution.rs | 4 +- crates/cqlib-core/src/qis/hamiltonian.rs | 13 +- crates/cqlib-core/src/qis/observable.rs | 35 +- crates/cqlib-core/src/qis/pauli.rs | 9 + .../src/qis/state/aligned_buffer.rs | 9 + .../src/qis/state/density_matrix.rs | 102 ++++-- .../src/qis/state/density_matrix_noise.rs | 253 +++++++++------ .../src/qis/state/density_matrix_test.rs | 14 + crates/cqlib-core/src/qis/state/mod.rs | 16 + .../cqlib-core/src/qis/state/statevector.rs | 36 ++- crates/cqlib-core/tests/p1_regressions.rs | 232 ++++++++++++++ crates/cqlib-core/tests/p2_regressions.rs | 299 ++++++++++++++++++ docs/api/cn/c/3_qis/1_statevector.md | 2 + docs/api/cn/c/3_qis/2_density_matrix.md | 6 + docs/api/cn/c/3_qis/3_density_matrix_noise.md | 2 + docs/api/cn/python/3_qis/1_statevector.md | 2 + docs/api/cn/python/3_qis/2_density_matrix.md | 2 + .../cn/python/3_qis/3_density_matrix_noise.md | 2 + docs/api/cn/python/3_qis/7_hamiltonian.md | 2 + docs/api/cn/python/3_qis/9_metrics_entropy.md | 6 +- docs/api/cn/rust/3_qis/1_statevector.md | 3 + docs/api/cn/rust/3_qis/2_density_matrix.md | 5 + .../cn/rust/3_qis/3_density_matrix_noise.md | 3 + docs/api/cn/rust/3_qis/7_hamiltonian.md | 4 +- docs/api/cn/rust/3_qis/8_evolution.md | 2 +- docs/api/cn/rust/3_qis/9_metrics_entropy.md | 6 +- docs/api/en/c/3_qis/1_statevector.md | 2 + docs/api/en/c/3_qis/2_density_matrix.md | 6 + docs/api/en/c/3_qis/3_density_matrix_noise.md | 2 + docs/api/en/python/3_qis/1_statevector.md | 2 + docs/api/en/python/3_qis/2_density_matrix.md | 2 + .../en/python/3_qis/3_density_matrix_noise.md | 2 + docs/api/en/python/3_qis/7_hamiltonian.md | 2 + docs/api/en/python/3_qis/9_metrics_entropy.md | 4 +- docs/api/en/rust/3_qis/1_statevector.md | 3 + docs/api/en/rust/3_qis/2_density_matrix.md | 5 + .../en/rust/3_qis/3_density_matrix_noise.md | 3 + docs/api/en/rust/3_qis/7_hamiltonian.md | 4 +- docs/api/en/rust/3_qis/8_evolution.md | 2 +- docs/api/en/rust/3_qis/9_metrics_entropy.md | 6 +- 62 files changed, 1642 insertions(+), 221 deletions(-) create mode 100644 crates/binding-c/tests/qis_constructor_regressions.rs create mode 100644 crates/binding-c/tests/qis_p2_regressions.rs create mode 100644 crates/binding-python/tests/qis/test_p1_regressions.py create mode 100644 crates/binding-python/tests/qis/test_p2_regressions.py create mode 100644 crates/cqlib-core/tests/p1_regressions.rs create mode 100644 crates/cqlib-core/tests/p2_regressions.rs diff --git a/crates/binding-c/src/qis/density_matrix.rs b/crates/binding-c/src/qis/density_matrix.rs index e8119644..3075939d 100644 --- a/crates/binding-c/src/qis/density_matrix.rs +++ b/crates/binding-c/src/qis/density_matrix.rs @@ -94,28 +94,26 @@ fn apply_three( .map_or_else(|e| e, |_| 0) } +/// Returns a newly allocated simulator, or NULL if its dimension is +/// unrepresentable or its state allocation fails. #[unsafe(no_mangle)] pub extern "C" fn density_matrix_new(num_qubits: usize) -> *mut CDensityMatrix { - Box::into_raw(Box::new(CDensityMatrix { - inner: DensityMatrix::new(num_qubits), - })) + match DensityMatrix::try_new(num_qubits) { + Ok(inner) => Box::into_raw(Box::new(CDensityMatrix { inner })), + Err(_) => std::ptr::null_mut(), + } } /// Create a density matrix filled entirely with zeros. This is not a valid /// physical state (trace = 0); it is useful as an accumulator during /// operations like Kraus channel application. Returns NULL when the matrix -/// size overflows the addressable range. +/// size overflows the addressable range or state allocation fails. #[unsafe(no_mangle)] pub extern "C" fn density_matrix_zeros(num_qubits: usize) -> *mut CDensityMatrix { - let Some(bits) = num_qubits.checked_mul(2) else { - return std::ptr::null_mut(); - }; - if bits >= usize::BITS as usize { - return std::ptr::null_mut(); + match DensityMatrix::try_zeros(num_qubits) { + Ok(inner) => Box::into_raw(Box::new(CDensityMatrix { inner })), + Err(_) => std::ptr::null_mut(), } - Box::into_raw(Box::new(CDensityMatrix { - inner: DensityMatrix::zeros(num_qubits), - })) } #[unsafe(no_mangle)] @@ -818,11 +816,14 @@ pub extern "C" fn density_matrix_validate_physical(ptr: *const CDensityMatrix, t } /// Create the maximally mixed state I / 2^N of `num_qubits` qubits. +/// Returns a newly allocated simulator, or NULL if its dimension is +/// unrepresentable or its state allocation fails. #[unsafe(no_mangle)] pub extern "C" fn density_matrix_maximally_mixed(num_qubits: usize) -> *mut CDensityMatrix { - Box::into_raw(Box::new(CDensityMatrix { - inner: DensityMatrix::maximally_mixed(num_qubits), - })) + match DensityMatrix::try_maximally_mixed(num_qubits) { + Ok(inner) => Box::into_raw(Box::new(CDensityMatrix { inner })), + Err(_) => std::ptr::null_mut(), + } } /// Apply a standard gate identified by name (e.g. "H", "CX", "RZZ", "FSIM"). diff --git a/crates/binding-c/src/qis/density_matrix_noise.rs b/crates/binding-c/src/qis/density_matrix_noise.rs index b521a736..c84d7c08 100644 --- a/crates/binding-c/src/qis/density_matrix_noise.rs +++ b/crates/binding-c/src/qis/density_matrix_noise.rs @@ -124,15 +124,18 @@ fn apply_three( /// Create a new noisy simulator with `num_qubits` initialised to |0...0>. /// Pass NULL `noise_model` for an ideal (noise-free) simulation. +/// Returns a newly allocated simulator, or NULL if its dimension is +/// unrepresentable or its state allocation fails. #[unsafe(no_mangle)] pub extern "C" fn density_matrix_noise_new( num_qubits: usize, noise_model: *const CNoiseModel, ) -> *mut CDensityMatrixNoise { let noise = clone_noise_model(noise_model); - Box::into_raw(Box::new(CDensityMatrixNoise { - inner: DensityMatrixNoise::new(num_qubits, noise), - })) + match DensityMatrixNoise::try_new(num_qubits, noise) { + Ok(inner) => Box::into_raw(Box::new(CDensityMatrixNoise { inner })), + Err(_) => std::ptr::null_mut(), + } } /// Free a noisy simulator. NULL is allowed. diff --git a/crates/binding-c/src/qis/statevector.rs b/crates/binding-c/src/qis/statevector.rs index cdb492ca..de03f952 100644 --- a/crates/binding-c/src/qis/statevector.rs +++ b/crates/binding-c/src/qis/statevector.rs @@ -97,11 +97,14 @@ fn apply_three( } /// Create a new statevector with `num_qubits` initialised to |0...0>. +/// Returns a newly allocated simulator, or NULL if its dimension is +/// unrepresentable or its state allocation fails. #[unsafe(no_mangle)] pub extern "C" fn statevector_new(num_qubits: usize) -> *mut CStatevector { - Box::into_raw(Box::new(CStatevector { - inner: Statevector::new(num_qubits), - })) + match Statevector::try_new(num_qubits) { + Ok(inner) => Box::into_raw(Box::new(CStatevector { inner })), + Err(_) => std::ptr::null_mut(), + } } /// Free a statevector. NULL is allowed. diff --git a/crates/binding-c/tests/qis_constructor_regressions.rs b/crates/binding-c/tests/qis_constructor_regressions.rs new file mode 100644 index 00000000..6fc86bae --- /dev/null +++ b/crates/binding-c/tests/qis_constructor_regressions.rs @@ -0,0 +1,51 @@ +// This code is part of Cqlib. +// +// (C) Copyright China Telecom Quantum Group 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. + +//! Invalid dimensions must return NULL instead of unwinding across the C ABI. +use binding_c::qis::*; + +#[test] +fn invalid_dimensions_return_null() { + for n in [usize::BITS as usize - 4, usize::BITS as usize, usize::MAX] { + assert!(statevector_new(n).is_null()); + } + for n in [ + (usize::BITS as usize - 4).div_ceil(2), + usize::BITS as usize, + usize::MAX, + ] { + assert!(density_matrix_new(n).is_null()); + assert!(density_matrix_zeros(n).is_null()); + assert!(density_matrix_maximally_mixed(n).is_null()); + assert!(density_matrix_noise_new(n, std::ptr::null()).is_null()); + } +} + +#[test] +fn normal_and_zero_qubit_constructors_remain_valid() { + for n in [0, 1, 3] { + let sv = statevector_new(n); + assert!(!sv.is_null()); + statevector_free(sv); + for dm in [ + density_matrix_new(n), + density_matrix_zeros(n), + density_matrix_maximally_mixed(n), + ] { + assert!(!dm.is_null()); + density_matrix_free(dm); + } + let noisy = density_matrix_noise_new(n, std::ptr::null()); + assert!(!noisy.is_null()); + density_matrix_noise_free(noisy); + } +} diff --git a/crates/binding-c/tests/qis_p2_regressions.rs b/crates/binding-c/tests/qis_p2_regressions.rs new file mode 100644 index 00000000..5e532b95 --- /dev/null +++ b/crates/binding-c/tests/qis_p2_regressions.rs @@ -0,0 +1,49 @@ +// This code is part of Cqlib. +// +// (C) Copyright China Telecom Quantum Group 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. + +use binding_c::qis::*; +use num_complex::Complex64; + +#[test] +fn nonfinite_states_return_null() { + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY, f64::MAX] { + for value in [Complex64::new(bad, 0.0), Complex64::new(0.0, bad)] { + let data = [value, Complex64::new(0.0, 0.0)]; + assert!(statevector_from_state(1, data.as_ptr(), data.len()).is_null()); + assert!(density_matrix_from_state(1, data.as_ptr(), data.len()).is_null()); + } + } +} + +#[test] +fn nonhermitian_expectation_reports_error_without_overwriting_output() { + let sv = statevector_new(1); + let dm = density_matrix_new(1); + let h = hamiltonian_new(1); + let p = pauli_string_parse(c"Z".as_ptr()); + assert_eq!(hamiltonian_add_term(h, p, 0.0, 1.0), 0); + let mut out = 123.0; + assert_ne!(statevector_expectation(sv, h, &mut out), 0); + assert_eq!(out, 123.0); + assert_ne!(density_matrix_expectation(dm, h, &mut out), 0); + assert_eq!(out, 123.0); + // The imaginary parts cancel: the observable is now Z. + assert_eq!(hamiltonian_add_term(h, p, 1.0, -1.0), 0); + assert_eq!(statevector_expectation(sv, h, &mut out), 0); + assert_eq!(out, 1.0); + assert_eq!(density_matrix_expectation(dm, h, &mut out), 0); + assert_eq!(out, 1.0); + pauli_string_free(p); + hamiltonian_free(h); + statevector_free(sv); + density_matrix_free(dm); +} diff --git a/crates/binding-python/cqlib/qis/entropy.pyi b/crates/binding-python/cqlib/qis/entropy.pyi index 2c21df18..3dc1d1ae 100644 --- a/crates/binding-python/cqlib/qis/entropy.pyi +++ b/crates/binding-python/cqlib/qis/entropy.pyi @@ -38,8 +38,9 @@ def renyi_entropy(dm: DensityMatrix, alpha: float) -> float: The Rényi entropy is defined as :math:`S_\\alpha(\\rho) = \\frac{1}{1-\\alpha} \\log_2(\\text{Tr}(\\rho^\\alpha))`. Special cases: - - When ``alpha -> 1``: Approaches Von Neumann entropy. + - When ``alpha = 1``: Returns Von Neumann entropy. - When ``alpha = 2``: Collision entropy. + - When ``alpha = +inf``: Min-entropy. Args: dm (DensityMatrix): The density matrix representing the quantum state. @@ -49,7 +50,13 @@ def renyi_entropy(dm: DensityMatrix, alpha: float) -> float: float: The Rényi entropy in bits (base-2 logarithm). Raises: - ValueError: If `alpha <= 0` or if eigendecomposition fails. + ValueError: If alpha is NaN or <= 0, eigendecomposition fails, + or no positive spectrum remains. + + Notes: + Non-diagonal spectra discard eigenvalues at or below + dimension * machine epsilon * largest eigenvalue as numerical zeros. + Diagonal spectra retain all positive entries. """ ... diff --git a/crates/binding-python/cqlib/qis/hamiltonian.pyi b/crates/binding-python/cqlib/qis/hamiltonian.pyi index 3210e3a8..7780e5b5 100644 --- a/crates/binding-python/cqlib/qis/hamiltonian.pyi +++ b/crates/binding-python/cqlib/qis/hamiltonian.pyi @@ -206,6 +206,10 @@ class Hamiltonian: Returns: The real expectation value. + + Raises: + ValueError: If the observable is not Hermitian or input dimensions differ. + ValueError: If original or merged coefficients are non-finite. """ ... @@ -217,6 +221,10 @@ class Hamiltonian: Returns: The real expectation value. + + Raises: + ValueError: If the observable is not Hermitian or input dimensions differ. + ValueError: If original or merged coefficients are non-finite. """ ... @@ -231,6 +239,11 @@ class Hamiltonian: Returns: The real expectation value. + + Raises: + ValueError: If the observable is not Hermitian, input dimensions differ, + or no compatible measurement basis is available. + ValueError: If original or merged coefficients are non-finite. """ ... @@ -332,7 +345,8 @@ class Hamiltonian: """Computes the variance for a statevector. Raises: - ValueError: If the Hamiltonian is not Hermitian or qubit counts differ. + ValueError: If the Hamiltonian is not Hermitian, qubit counts differ, + or original or merged coefficients are non-finite. """ ... diff --git a/crates/binding-python/cqlib/qis/pauli.pyi b/crates/binding-python/cqlib/qis/pauli.pyi index c7f0ee15..dd14c3d8 100644 --- a/crates/binding-python/cqlib/qis/pauli.pyi +++ b/crates/binding-python/cqlib/qis/pauli.pyi @@ -289,7 +289,8 @@ class PauliString: The expectation value as a float. Raises: - ValueError: If state string length doesn't match num_qubits or contains invalid chars. + ValueError: If the Pauli string is not Hermitian, or a state string + has the wrong length or contains invalid characters. """ ... @@ -301,6 +302,9 @@ class PauliString: Returns: The real expectation value. + + Raises: + ValueError: If the observable is not Hermitian or input dimensions differ. """ ... @@ -312,6 +316,9 @@ class PauliString: Returns: The real expectation value. + + Raises: + ValueError: If the observable is not Hermitian or input dimensions differ. """ ... @@ -326,6 +333,10 @@ class PauliString: Returns: The real expectation value. + + Raises: + ValueError: If the observable is not Hermitian, input dimensions differ, + or no compatible measurement basis is available. """ ... diff --git a/crates/binding-python/cqlib/qis/state/density_matrix.pyi b/crates/binding-python/cqlib/qis/state/density_matrix.pyi index 61515eff..994e151b 100644 --- a/crates/binding-python/cqlib/qis/state/density_matrix.pyi +++ b/crates/binding-python/cqlib/qis/state/density_matrix.pyi @@ -59,6 +59,9 @@ class DensityMatrix: Returns: A new DensityMatrix instance in the ground state + Raises: + ValueError: If the dimension is unrepresentable or state allocation fails. + Examples: >>> dm = DensityMatrix(2) # |00⟩⟨00| state """ @@ -572,7 +575,8 @@ class DensityMatrix: The expectation value as a real number Raises: - ValueError: If the qubit counts don't match or the observable type is invalid + ValueError: If qubit counts differ, the observable type is invalid, + the observable is not Hermitian, or its coefficients are non-finite. """ ... diff --git a/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi b/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi index fc37fcd2..f0fadb76 100644 --- a/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi +++ b/crates/binding-python/cqlib/qis/state/density_matrix_noise.pyi @@ -60,6 +60,9 @@ class DensityMatrixNoise: Returns: A new DensityMatrixNoise instance + Raises: + ValueError: If the dimension is unrepresentable or state allocation fails. + Examples: >>> from cqlib.qis.state import DensityMatrixNoise >>> # Simulator without noise (ideal simulation) @@ -304,7 +307,8 @@ class DensityMatrixNoise: The expectation value as a real number Raises: - ValueError: If the qubit counts don't match or the observable type is invalid + ValueError: If qubit counts differ, the observable type is invalid, + the observable is not Hermitian, or its coefficients are non-finite. """ ... diff --git a/crates/binding-python/cqlib/qis/state/statevector.pyi b/crates/binding-python/cqlib/qis/state/statevector.pyi index bd997414..6a8e080e 100644 --- a/crates/binding-python/cqlib/qis/state/statevector.pyi +++ b/crates/binding-python/cqlib/qis/state/statevector.pyi @@ -61,6 +61,9 @@ class Statevector: Returns: A new Statevector instance in the ground state + Raises: + ValueError: If the dimension is unrepresentable or state allocation fails. + Examples: >>> sv = Statevector(2) # |00⟩ state """ @@ -559,7 +562,8 @@ class Statevector: The expectation value as a real number Raises: - ValueError: If the qubit counts don't match or the observable type is invalid + ValueError: If qubit counts differ, the observable type is invalid, + the observable is not Hermitian, or its coefficients are non-finite. """ ... diff --git a/crates/binding-python/src/qis/entropy.rs b/crates/binding-python/src/qis/entropy.rs index 7f7ff362..756f9599 100644 --- a/crates/binding-python/src/qis/entropy.rs +++ b/crates/binding-python/src/qis/entropy.rs @@ -43,13 +43,20 @@ pub fn linear_entropy(dm: &PyDensityMatrix) -> PyResult { /// /// Args: /// dm (DensityMatrix): The density matrix representing the quantum state. -/// alpha (float): The order parameter. Must be positive. +/// alpha (float): Positive order, including 1 (Von Neumann entropy) +/// and positive infinity (min-entropy). /// /// Returns: /// float: The Rényi entropy in bits (base-2 logarithm). /// /// Raises: -/// ValueError: If alpha <= 0 or if eigendecomposition fails. +/// ValueError: If alpha is NaN or <= 0, eigendecomposition fails, +/// or no positive spectrum remains. +/// +/// Notes: +/// Non-diagonal spectra discard eigenvalues at or below +/// dimension * machine epsilon * largest eigenvalue as numerical zeros. +/// Diagonal spectra retain all positive entries. #[pyfunction] pub fn renyi_entropy(dm: &PyDensityMatrix, alpha: f64) -> PyResult { core_entropy::renyi_entropy(&dm.inner, alpha).map_err(qis_error_to_py_err) diff --git a/crates/binding-python/src/qis/hamiltonian.rs b/crates/binding-python/src/qis/hamiltonian.rs index 4ef69645..0addb18b 100644 --- a/crates/binding-python/src/qis/hamiltonian.rs +++ b/crates/binding-python/src/qis/hamiltonian.rs @@ -466,6 +466,10 @@ impl PyHamiltonian { } /// Computes the expectation value for a statevector. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. + /// ValueError: If original or merged coefficients are non-finite. fn expectation_statevector(&self, sv: &PyStatevector) -> PyResult { self.inner .expectation_statevector(&sv.inner) @@ -473,6 +477,10 @@ impl PyHamiltonian { } /// Computes the expectation value for a density matrix. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. + /// ValueError: If original or merged coefficients are non-finite. fn expectation_density_matrix(&self, dm: &PyDensityMatrix) -> PyResult { self.inner .expectation_density_matrix(&dm.inner) @@ -480,6 +488,10 @@ impl PyHamiltonian { } /// Computes the expectation value from measurement probabilities. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. + /// ValueError: If original or merged coefficients are non-finite. fn expectation_probs(&self, measurements: &Bound<'_, PyAny>) -> PyResult { let mut rust_measurements = Vec::new(); for item in measurements.try_iter()? { @@ -510,6 +522,10 @@ impl PyHamiltonian { } /// Computes the variance for a statevector. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. + /// ValueError: If original or merged coefficients are non-finite. fn variance_statevector(&self, sv: &PyStatevector) -> PyResult { self.inner .variance_statevector(&sv.inner) diff --git a/crates/binding-python/src/qis/pauli.rs b/crates/binding-python/src/qis/pauli.rs index 948631bf..e9157964 100644 --- a/crates/binding-python/src/qis/pauli.rs +++ b/crates/binding-python/src/qis/pauli.rs @@ -507,7 +507,8 @@ impl PyPauliString { /// The expectation value as a float. /// /// Raises: - /// ValueError: If state string length doesn't match num_qubits or contains invalid chars. + /// ValueError: If the Pauli string is not Hermitian, or a state string + /// has the wrong length or contains invalid characters. /// /// Examples: /// >>> ps = PauliString.from_str("IZ") # Z on qubit 0 @@ -527,6 +528,9 @@ impl PyPauliString { } /// Computes the expectation value for a statevector. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. fn expectation_statevector(&self, sv: &PyStatevector) -> PyResult { self.inner .expectation_statevector(&sv.inner) @@ -534,6 +538,9 @@ impl PyPauliString { } /// Computes the expectation value for a density matrix. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. fn expectation_density_matrix(&self, dm: &PyDensityMatrix) -> PyResult { self.inner .expectation_density_matrix(&dm.inner) @@ -541,6 +548,9 @@ impl PyPauliString { } /// Computes the expectation value from measurement probabilities. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. fn expectation_probs(&self, measurements: &Bound<'_, PyAny>) -> PyResult { let mut rust_measurements = Vec::new(); for item in measurements.try_iter()? { @@ -571,6 +581,9 @@ impl PyPauliString { } /// Computes the variance for a statevector. + /// + /// Raises: + /// ValueError: If the observable is not Hermitian or input dimensions differ. fn variance_statevector(&self, sv: &PyStatevector) -> PyResult { self.inner .variance_statevector(&sv.inner) diff --git a/crates/binding-python/src/qis/state/density_matrix.rs b/crates/binding-python/src/qis/state/density_matrix.rs index 44e1b174..d9fef301 100644 --- a/crates/binding-python/src/qis/state/density_matrix.rs +++ b/crates/binding-python/src/qis/state/density_matrix.rs @@ -1,4 +1,5 @@ // This code is part of Cqlib. +// Modified: fix P1 numerical boundaries and simulator error handling. // Modified to support reproducible state sampling. // // (C) Copyright China Telecom Quantum Group 2026 @@ -83,14 +84,17 @@ impl PyDensityMatrix { /// Returns: /// A new DensityMatrix instance in the ground state /// + /// Raises: + /// ValueError: If the dimension is unrepresentable or state allocation fails. + /// /// Examples: /// >>> from cqlib.qis import DensityMatrix /// >>> dm = DensityMatrix(2) # |00><00| state #[new] - fn new(num_qubits: usize) -> Self { - Self { - inner: DensityMatrix::new(num_qubits), - } + fn new(num_qubits: usize) -> PyResult { + Ok(Self { + inner: DensityMatrix::try_new(num_qubits).map_err(qis_error_to_py_err)?, + }) } /// Creates a density matrix from an initial statevector (pure state). @@ -239,10 +243,10 @@ impl PyDensityMatrix { /// >>> dm.probabilities() /// [0.5, 0.5] #[staticmethod] - fn maximally_mixed(num_qubits: usize) -> Self { - Self { - inner: DensityMatrix::maximally_mixed(num_qubits), - } + fn maximally_mixed(num_qubits: usize) -> PyResult { + Ok(Self { + inner: DensityMatrix::try_maximally_mixed(num_qubits).map_err(qis_error_to_py_err)?, + }) } /// Applies a circuit to this density matrix in place. @@ -848,7 +852,8 @@ impl PyDensityMatrix { /// The expectation value as a real number /// /// Raises: - /// ValueError: If the qubit counts don't match or the observable type is invalid + /// ValueError: If qubit counts differ, the observable type is invalid, + /// the observable is not Hermitian, or its coefficients are non-finite. fn expectation(&self, py: Python<'_>, observable: &Bound<'_, PyAny>) -> PyResult { if let Ok(h) = observable.extract::() { py.detach(|| self.inner.expectation(&h.inner)) diff --git a/crates/binding-python/src/qis/state/density_matrix_noise.rs b/crates/binding-python/src/qis/state/density_matrix_noise.rs index 84bfe638..8916923b 100644 --- a/crates/binding-python/src/qis/state/density_matrix_noise.rs +++ b/crates/binding-python/src/qis/state/density_matrix_noise.rs @@ -81,6 +81,9 @@ impl PyDensityMatrixNoise { /// Returns: /// A new DensityMatrixNoise instance /// + /// Raises: + /// ValueError: If the dimension is unrepresentable or state allocation fails. + /// /// Examples: /// >>> from cqlib.qis.state import DensityMatrixNoise /// >>> # Simulator without noise (ideal simulation) @@ -90,11 +93,11 @@ impl PyDensityMatrixNoise { /// >>> sim = DensityMatrixNoise(2, NoiseModel()) #[new] #[pyo3(signature = (num_qubits, noise_model=None))] - fn new(num_qubits: usize, noise_model: Option) -> Self { + fn new(num_qubits: usize, noise_model: Option) -> PyResult { let model = noise_model.map(|m| m.inner); - Self { - inner: DensityMatrixNoise::new(num_qubits, model), - } + Ok(Self { + inner: DensityMatrixNoise::try_new(num_qubits, model).map_err(qis_error_to_py_err)?, + }) } /// Simulates a circuit, applying noise after each operation. @@ -434,6 +437,10 @@ impl PyDensityMatrixNoise { } /// Computes the expectation value of an observable. + /// + /// Raises: + /// ValueError: If qubit counts differ, the observable type is invalid, + /// the observable is not Hermitian, or its coefficients are non-finite. fn expectation(&self, py: Python<'_>, observable: &Bound<'_, PyAny>) -> PyResult { if let Ok(h) = observable.extract::() { py.detach(|| self.inner.expectation(&h.inner)) diff --git a/crates/binding-python/src/qis/state/statevector.rs b/crates/binding-python/src/qis/state/statevector.rs index 506bb3e7..702769db 100644 --- a/crates/binding-python/src/qis/state/statevector.rs +++ b/crates/binding-python/src/qis/state/statevector.rs @@ -83,14 +83,17 @@ impl PyStatevector { /// Returns: /// A new Statevector instance in the ground state /// + /// Raises: + /// ValueError: If the dimension is unrepresentable or state allocation fails. + /// /// Examples: /// >>> from cqlib.qis import Statevector /// >>> sv = Statevector(2) # |00⟩ state #[new] - fn new(num_qubits: usize) -> Self { - Self { - inner: Statevector::new(num_qubits), - } + fn new(num_qubits: usize) -> PyResult { + Ok(Self { + inner: Statevector::try_new(num_qubits).map_err(qis_error_to_py_err)?, + }) } /// Creates a statevector from initial amplitudes. @@ -194,7 +197,8 @@ impl PyStatevector { /// The expectation value as a real number /// /// Raises: - /// ValueError: If the qubit counts don't match or the observable type is invalid + /// ValueError: If qubit counts differ, the observable type is invalid, + /// the observable is not Hermitian, or its coefficients are non-finite. fn expectation(&self, py: Python<'_>, observable: &Bound<'_, PyAny>) -> PyResult { if let Ok(h) = observable.extract::() { py.detach(|| self.inner.expectation(&h.inner)) diff --git a/crates/binding-python/tests/qis/test_p1_regressions.py b/crates/binding-python/tests/qis/test_p1_regressions.py new file mode 100644 index 00000000..9ad84773 --- /dev/null +++ b/crates/binding-python/tests/qis/test_p1_regressions.py @@ -0,0 +1,124 @@ +# This code is part of Cqlib. +# +# (C) Copyright China Telecom Quantum Group 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. + +"""Public API regressions for P1 numerical and simulation bugs.""" + +import struct + +import numpy as np +import pytest +from cqlib.circuit import Circuit, Qubit, StandardGate +from cqlib.device import NoiseModel, SingleQubitNoise +from cqlib.qis import ( + Statevector, + DensityMatrix, + DensityMatrixNoise, + Hamiltonian, + PauliString, + TrotterMode, + entropy, +) + + +@pytest.mark.parametrize( + "constructor", + [Statevector, DensityMatrix, DensityMatrixNoise, DensityMatrix.maximally_mixed], +) +def test_invalid_dimensions_raise_value_error(constructor): + bits = struct.calcsize("P") * 8 + for n in [bits - 4, bits, 2**bits - 1]: + with pytest.raises(ValueError, match="dimension"): + constructor(n) + for n in [0, 1, 2]: + assert np.sum(constructor(n).probabilities()) == pytest.approx(1.0) + + +def test_pauli_probability_boundary(): + noise = SingleQubitNoise.pauli(0.3, 0.6, 0.1) + assert noise.is_valid() + assert all(np.isfinite(k).all() for k in noise.to_kraus()) + model = NoiseModel() + model.add_single_qubit_error(StandardGate.X, 0, noise) + sim = DensityMatrixNoise(1, model) + sim.apply_x(0) + np.testing.assert_allclose(sim.probabilities(), [0.9, 0.1], atol=1e-12) + + +def test_gate_noise_uses_sparse_labels(): + model = NoiseModel() + model.add_single_qubit_error(StandardGate.X, 5, SingleQubitNoise.bit_flip(1.0)) + circuit = Circuit([Qubit(5)]) + circuit.x(Qubit(5)) + sim = DensityMatrixNoise.from_circuit(circuit, model) + np.testing.assert_allclose(sim.probabilities(), [1.0, 0.0], atol=1e-12) + sim.apply_x(0) + np.testing.assert_allclose(sim.probabilities(), [1.0, 0.0], atol=1e-12) + + +def test_formation_entropy_roundoff(): + state = np.array( + [ + -0.5917562523120286 - 0.14345273890626087j, + 0.12739228627789673 + 0.3361800930257551j, + -0.06600502193827436 + 0.35339664207352467j, + -0.5571108451004554 + 0.2457267928442915j, + ] + ) + assert entropy.entanglement_of_formation( + DensityMatrix.from_state(2, state) + ) == pytest.approx(1.0, abs=1e-12) + + +def test_renyi_extreme_orders(): + dm = DensityMatrix.maximally_mixed(1) + for alpha in [1.0 - 1e-12, 1.0 + 1e-12, 2000.0, float("inf")]: + assert entropy.renyi_entropy(dm, alpha) == pytest.approx(1.0, abs=1e-10) + with pytest.raises(ValueError): + entropy.renyi_entropy(dm, float("nan")) + small = DensityMatrix.from_density_matrix( + 1, np.diag([1 - 1e-14, 1e-14]).astype(complex).flatten() + ) + assert entropy.renyi_entropy(small, 0.1) == pytest.approx( + 0.06257881098503021, abs=1e-12 + ) + + +@pytest.mark.parametrize("method", ["to_trotter_circuit", "to_evolution_circuit"]) +def test_evolution_retains_small_terms(method): + h = Hamiltonian(1) + h.add_term(PauliString.from_str("Z"), 1e-11) + circuit = getattr(h, method)(1e11, 3, TrotterMode.first_order()) + # Probabilities alone cannot expose this regression: compare the full matrix. + from cqlib.circuit import circuit_to_matrix + + np.testing.assert_allclose( + circuit_to_matrix(circuit), np.diag(np.exp([-1j, 1j])), atol=1e-12 + ) + + +@pytest.mark.parametrize("state_type", [Statevector, DensityMatrix, DensityMatrixNoise]) +def test_oversized_circuit_raises_value_error(state_type): + with pytest.raises(ValueError, match="dimension"): + state_type.from_circuit(Circuit(struct.calcsize("P") * 8)) + + +@pytest.mark.parametrize( + "constructor", + [ + Statevector.from_state, + DensityMatrix.from_state, + DensityMatrix.from_density_matrix, + ], +) +def test_oversized_array_constructors_raise_value_error(constructor): + with pytest.raises(ValueError, match="dimension"): + constructor(struct.calcsize("P") * 8, []) diff --git a/crates/binding-python/tests/qis/test_p2_regressions.py b/crates/binding-python/tests/qis/test_p2_regressions.py new file mode 100644 index 00000000..b809b5e8 --- /dev/null +++ b/crates/binding-python/tests/qis/test_p2_regressions.py @@ -0,0 +1,240 @@ +# This code is part of Cqlib. +# +# (C) Copyright China Telecom Quantum Group 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. + +import itertools +import numpy as np +import pytest +from cqlib.circuit import ( + Circuit, + Qubit, + StandardGate, + Instruction, + Parameter, + circuit_to_matrix, +) +from cqlib.device import NoiseModel, SingleQubitNoise +from cqlib.qis import ( + Statevector, + DensityMatrix, + DensityMatrixNoise, + Hamiltonian, + PauliString, + entropy, +) +from cqlib.error_mitigation.zne import ZNEMitigation + + +@pytest.mark.parametrize( + "constructor", [Statevector.from_state, DensityMatrix.from_state] +) +@pytest.mark.parametrize( + "bad", [float("nan"), float("inf"), -float("inf"), np.finfo(float).max] +) +def test_nonfinite_norm_rejected(constructor, bad): + for amp in [complex(bad, 0), complex(0, bad)]: + with pytest.raises(ValueError, match="normalized"): + constructor(1, [amp, 0j]) + assert constructor(0, [1j]).probabilities() == pytest.approx([1]) + assert constructor(1, [0.6, 0.8j]).probabilities() == pytest.approx([0.36, 0.64]) + + +@pytest.mark.parametrize("bad", [float("nan"), float("inf"), -float("inf")]) +def test_nonfinite_density_matrix_rejected(bad): + for pos in range(4): + for value in [complex(bad, 0), complex(0, bad)]: + data = np.eye(2, dtype=complex).flatten() / 2 + data[pos] = value + with pytest.raises(ValueError): + DensityMatrix.from_density_matrix(1, data) + + +@pytest.mark.parametrize("label", ["iZ", "-iI", "iX"]) +def test_real_expectations_reject_nonhermitian(label): + p = PauliString.from_str(label) + h = Hamiltonian(1) + h.add_term(p, 1) + measurements = [(p, {"0": 1.0})] + for observable in [p, h]: + for state in [Statevector(1), DensityMatrix(1), DensityMatrixNoise(1)]: + with pytest.raises(ValueError, match="Hermitian"): + state.expectation(observable) + with pytest.raises(ValueError, match="Hermitian"): + observable.expectation_probs(measurements) + with pytest.raises(ValueError, match="Hermitian"): + p.expectation({"0": 1.0}) + + +def test_merged_hermitian_and_small_coefficients(): + for terms, expected in [ + ([("iZ", -1j)], 1), + ([("Z", 1 + 2j), ("Z", -2j)], 1), + ([("Z", 1e-11)], 1e-11), + ]: + h = Hamiltonian(1) + for label, coeff in terms: + h.add_term(PauliString.from_str(label), coeff) + assert h.expectation_statevector(Statevector(1)) == pytest.approx( + expected, abs=1e-14 + ) + assert h.expectation_density_matrix(DensityMatrix(1)) == pytest.approx( + expected, abs=1e-14 + ) + assert h.expectation_probs( + [(PauliString.from_str("Z"), {"0": 1.0})] + ) == pytest.approx(expected, abs=1e-14) + + +def test_pauli_expectations_against_independent_dense_matrices(): + rng = np.random.default_rng(7451) + matrices = { + "I": np.eye(2), + "X": np.array([[0, 1], [1, 0]]), + "Y": np.array([[0, -1j], [1j, 0]]), + "Z": np.diag([1, -1]), + } + for n in [1, 2, 3]: + psi = rng.normal(size=2**n) + 1j * rng.normal(size=2**n) + psi /= np.linalg.norm(psi) + sv = Statevector.from_state(n, psi) + rho = 0.7 * np.outer(psi, psi.conj()) + 0.3 * np.eye(2**n) / 2**n + dm = DensityMatrix.from_density_matrix(n, rho.flatten()) + for labels in itertools.product("IXYZ", repeat=n): + matrix = np.array([[1]]) + for label in labels: + matrix = np.kron(matrix, matrices[label]) + for sign in [1, -1]: + p = PauliString.from_str(("-" if sign == -1 else "") + "".join(labels)) + assert sv.expectation(p) == pytest.approx( + sign * np.vdot(psi, matrix @ psi).real, abs=1e-12 + ) + assert dm.expectation(p) == pytest.approx( + sign * np.trace(rho @ matrix).real, abs=1e-12 + ) + + +@pytest.mark.parametrize("p", [0, 0.25, 1]) +def test_identity_noise_sparse_labels(p): + model = NoiseModel() + model.add_single_qubit_error(StandardGate.I, 5, SingleQubitNoise.bit_flip(p)) + for count in [1, 2, 3]: + c = Circuit([Qubit(12), Qubit(5)]) + for _ in range(count): + c.i(Qubit(5)) + flipped = (1 - (1 - 2 * p) ** count) / 2 + assert DensityMatrixNoise.from_circuit( + c, model + ).probabilities() == pytest.approx([1 - flipped, 0, flipped, 0]) + + +@pytest.mark.parametrize("symbolic", [False, True]) +@pytest.mark.parametrize("with_gates", [False, True]) +def test_zne_preserves_phase_and_parameter_binding(symbolic, with_gates): + c = Circuit([Qubit(5)]) + c.set_global_phase(Parameter("phi") if symbolic else 0.37) + if with_gates: + c.h(Qubit(5)) + c.t(Qubit(5)) + expected = circuit_to_matrix(c.assign_parameters({"phi": 0.37})) + for selection in [None, [], [Instruction.from_standard_gate(StandardGate.H)]]: + for folded in ZNEMitigation(c, [0, 1, 2]).fold_circuits(selection): + np.testing.assert_allclose( + circuit_to_matrix(folded.assign_parameters({"phi": 0.37})), + expected, + atol=1e-12, + ) + + +def test_renyi_known_nondiagonal_spectrum(): + # Rotated mixed state: eigenvalues stay known without a reference eigensolver. + spectrum = np.array([0.8, 0.2]) + u = np.array([[1, 1j], [1j, 1]]) / np.sqrt(2) + rho = u @ np.diag(spectrum) @ u.conj().T + dm = DensityMatrix.from_density_matrix(1, rho.flatten()) + shannon = -np.sum(spectrum * np.log2(spectrum)) + for alpha in [0.1, 0.5, 1 - 1e-12, 1, 1 + 1e-12, 2, 20, float("inf")]: + if abs(alpha - 1) < 1e-8: + expected = shannon + elif np.isinf(alpha): + expected = -np.log2(max(spectrum)) + else: + expected = np.log2(np.sum(spectrum**alpha)) / (1 - alpha) + assert entropy.renyi_entropy(dm, alpha) == pytest.approx(expected, abs=1e-10) + + +@pytest.mark.parametrize("bad", [complex(float("nan"), 0), complex(0, float("inf"))]) +def test_nonfinite_observable_coefficients(bad): + h = Hamiltonian(1) + h.add_term(PauliString.from_str("Z"), bad) + for call in [ + lambda: h.expectation_statevector(Statevector(1)), + lambda: h.expectation_density_matrix(DensityMatrix(1)), + lambda: h.expectation_probs([(PauliString.from_str("Z"), {"0": 1.0})]), + ]: + with pytest.raises(ValueError, match="finite"): + call() + + +@pytest.mark.parametrize("alpha", [0.01, 0.1, 0.5, 1, 2, float("inf")]) +def test_renyi_nondiagonal_pure_state_is_zero(alpha): + psi = np.array([1, 2, 3, 4], dtype=complex) / np.sqrt(30) + dm = DensityMatrix.from_state(2, psi) + assert entropy.renyi_entropy(dm, alpha) == pytest.approx(0, abs=1e-10) + + +@pytest.mark.parametrize("alpha", [0.01, 0.1, 0.5, 1, 2, float("inf")]) +def test_renyi_rank_deficient_spectrum_is_invariant_under_rotation(alpha): + spectrum = np.array([0.8, 0.2, 0, 0]) + indices = np.arange(4) + u = np.exp(2j * np.pi * np.outer(indices, indices) / 4) / 2 + expected = ( + -np.sum(spectrum[:2] * np.log2(spectrum[:2])) + if alpha == 1 + else -np.log2(0.8) + if np.isinf(alpha) + else np.log2(np.sum(spectrum[:2] ** alpha)) / (1 - alpha) + ) + for rho in [np.diag(spectrum), u @ np.diag(spectrum) @ u.conj().T]: + dm = DensityMatrix.from_density_matrix(2, rho.astype(complex).flatten()) + assert entropy.renyi_entropy(dm, alpha) == pytest.approx(expected, abs=1e-10) + + +@pytest.mark.parametrize("bad", [float("nan"), float("inf"), -float("inf")]) +def test_variance_rejects_nonfinite_coefficients(bad): + for coeff in [complex(bad, 0), complex(0, bad)]: + h = Hamiltonian(1) + h.add_term(PauliString.from_str("Z"), coeff) + with pytest.raises(ValueError, match="finite"): + h.variance_statevector(Statevector(1)) + + +def test_variance_validates_merged_coefficients_and_preserves_small_terms(): + h = Hamiltonian(1) + for _ in range(2): + h.add_term(PauliString.from_str("Z"), np.finfo(float).max) + with pytest.raises(ValueError, match="finite"): + h.variance_statevector(Statevector(1)) + for terms, expected in [ + ([("iX", -1j)], 1), + ([("X", 1 + 2j), ("X", -2j)], 1), + ([("X", 1e-11)], 1e-22), + ]: + h = Hamiltonian(1) + for label, coeff in terms: + h.add_term(PauliString.from_str(label), coeff) + assert h.variance_statevector(Statevector(1)) == pytest.approx( + expected, rel=1e-12, abs=1e-34 + ) + invalid = Hamiltonian.from_pauli(PauliString.from_str("iZ")) + with pytest.raises(ValueError, match="Hermitian"): + invalid.variance_statevector(Statevector(1)) + with pytest.raises(ValueError, match="Qubit count mismatch"): + invalid.variance_statevector(Statevector(2)) diff --git a/crates/cqlib-core/src/device/noise.rs b/crates/cqlib-core/src/device/noise.rs index 9bcd9528..903367f4 100644 --- a/crates/cqlib-core/src/device/noise.rs +++ b/crates/cqlib-core/src/device/noise.rs @@ -164,7 +164,7 @@ impl SingleQubitNoise { ] } Self::Pauli { px, py, pz } => { - let pi = (1.0 - px - py - pz).sqrt(); + let pi = (1.0 - (px + py + pz)).sqrt(); vec![ Pauli::I.to_matrix() * pi, Pauli::X.to_matrix() * px.sqrt(), diff --git a/crates/cqlib-core/src/error_mitigation/zne_mitigation.rs b/crates/cqlib-core/src/error_mitigation/zne_mitigation.rs index 33ec9ce4..5223ca3e 100644 --- a/crates/cqlib-core/src/error_mitigation/zne_mitigation.rs +++ b/crates/cqlib-core/src/error_mitigation/zne_mitigation.rs @@ -348,6 +348,7 @@ impl ZNEMitigation { fn fold_all(&self, level: usize) -> Result { let mut folded = Circuit::from_qubits(self.circuit.qubits())?; + folded.set_global_phase(self.circuit.global_phase()); let inverse = self.circuit.inverse()?; self.append_circuit_ops(&mut folded, &self.circuit)?; @@ -366,6 +367,7 @@ impl ZNEMitigation { ) -> Result { let gate_names: HashSet = gate_set.iter().map(Instruction::name).collect(); let mut folded = Circuit::from_qubits(self.circuit.qubits())?; + folded.set_global_phase(self.circuit.global_phase()); for op in self.circuit.operations() { self.append_operation(&mut folded, &self.circuit, op)?; diff --git a/crates/cqlib-core/src/qis/entropy.rs b/crates/cqlib-core/src/qis/entropy.rs index 60b23f31..162c7cc0 100644 --- a/crates/cqlib-core/src/qis/entropy.rs +++ b/crates/cqlib-core/src/qis/entropy.rs @@ -92,12 +92,14 @@ pub fn linear_entropy(dm: &DensityMatrix) -> Result { /// # Arguments /// /// * `dm` - The density matrix representing the quantum state. -/// * `alpha` - The order parameter. Must be positive and not equal to 1.0 -/// (use Von Neumann entropy directly for α = 1). +/// * `alpha` - Positive order parameter, including 1.0 and positive infinity. +/// NaN is rejected. At α = 1 this returns the Von Neumann entropy; +/// positive infinity returns the min-entropy. /// /// # Returns /// -/// The Rényi entropy in bits (base-2 logarithm), or an error if α ≤ 0. +/// The Rényi entropy in bits (base-2 logarithm), or an error if α ≤ 0, +/// α is NaN, eigendecomposition fails, or no positive spectrum remains. /// /// # Examples /// @@ -115,8 +117,13 @@ pub fn linear_entropy(dm: &DensityMatrix) -> Result { /// /// For α very close to 1.0 (within `f64::EPSILON`), this function automatically /// falls back to computing the Von Neumann entropy for numerical stability. +/// Nearby orders use `expm1` and `log1p` to avoid cancellation. For non-diagonal +/// matrices, eigenvalues at or below `dimension * f64::EPSILON * max_eigenvalue` +/// are treated as numerical zeros. True eigenvalues below this resolution +/// cannot reliably be distinguished from eigensolver roundoff. Diagonal +/// matrices retain all positive entries without this cutoff. pub fn renyi_entropy(dm: &DensityMatrix, alpha: f64) -> Result { - if alpha <= 0.0 { + if alpha.is_nan() || alpha <= 0.0 { return Err(QisError::InvalidParameterValue( "Rényi entropy order alpha must be positive".to_string(), )); @@ -128,20 +135,57 @@ pub fn renyi_entropy(dm: &DensityMatrix, alpha: f64) -> Result { } let mat = super::metrics::density_matrix_to_faer(dm); - let eigenvalues: Vec = mat.self_adjoint_eigenvalues(Side::Lower).map_err(|e| { - QisError::UnsupportedOperation(format!("Eigendecomposition failed: {:?}", e)) - })?; - - // Compute sum of eigenvalues^alpha, handling numerical noise - let mut sum_power = 0.0; - for &eigval in &eigenvalues { - if eigval > 1e-12 { - sum_power += eigval.powf(alpha); - } + // A diagonal spectrum is already available exactly. Eigensolvers can lose + // relative accuracy in tiny eigenvalues, which matter at small alpha. + let diagonal = (0..mat.nrows()) + .all(|i| (0..mat.ncols()).all(|j| i == j || mat[(i, j)] == faer::c64::new(0.0, 0.0))); + let eigenvalues: Vec = if diagonal { + (0..mat.nrows()).map(|i| mat[(i, i)].re).collect() + } else { + mat.self_adjoint_eigenvalues(Side::Lower).map_err(|e| { + QisError::UnsupportedOperation(format!("Eigendecomposition failed: {:?}", e)) + })? + }; + + let max_eigenvalue = eigenvalues.iter().copied().fold(0.0, f64::max); + if !max_eigenvalue.is_finite() || max_eigenvalue == 0.0 { + return Err(QisError::NotNormalized); } - - // Rényi entropy formula: 1/(1-α) * log2(sum(λ_i^α)) - let entropy = sum_power.log2() / (1.0 - alpha); + // Small alpha amplifies positive roundoff in the null space of low-rank + // states. Scale the cutoff with eigensolver precision instead of using a + // fixed probability threshold, while keeping exact diagonal spectra. + let cutoff = if diagonal { + 0.0 + } else { + mat.nrows() as f64 * f64::EPSILON * max_eigenvalue + }; + let eigenvalues: Vec = eigenvalues.into_iter().filter(|&v| v > cutoff).collect(); + if alpha == f64::INFINITY { + return Ok((-max_eigenvalue.log2()).max(0.0)); + } + let delta = alpha - 1.0; + let entropy = if delta.abs() < 1e-5 { + // expm1/log1p avoid cancellation near alpha=1. Normalize the spectral + // weights here so trace roundoff is not amplified by division by delta. + let trace: f64 = eigenvalues.iter().sum(); + let change: f64 = eigenvalues + .iter() + .map(|&v| { + let p = v / trace; + p * (delta * p.ln()).exp_m1() + }) + .sum(); + -change.ln_1p() / delta / std::f64::consts::LN_2 + } else { + let max_log = max_eigenvalue.ln(); + let scaled_sum: f64 = eigenvalues + .iter() + .map(|v| (alpha * (v.ln() - max_log)).exp()) + .sum(); + // Divide before multiplying to avoid overflow even for alpha=f64::MAX. + ((alpha / (1.0 - alpha)) * max_log + scaled_sum.ln() / (1.0 - alpha)) + / std::f64::consts::LN_2 + }; Ok(entropy.max(0.0)) } @@ -442,7 +486,7 @@ pub fn concurrence(dm: &DensityMatrix) -> Result { /// This measure represents the minimum number of Bell pairs required to asymptotically /// create the state ρ using LOCC (Local Operations and Classical Communication). pub fn entanglement_of_formation(dm: &DensityMatrix) -> Result { - let c = concurrence(dm)?; + let c = concurrence(dm)?.clamp(0.0, 1.0); if c < 1e-15 { return Ok(0.0); diff --git a/crates/cqlib-core/src/qis/evolution.rs b/crates/cqlib-core/src/qis/evolution.rs index 2079e981..e81dc3a9 100644 --- a/crates/cqlib-core/src/qis/evolution.rs +++ b/crates/cqlib-core/src/qis/evolution.rs @@ -522,7 +522,7 @@ impl Hamiltonian { } let mut simplified_h = self.clone(); - simplified_h.simplify(); + simplified_h.simplify_with_tolerance(0.0); for (_, coeff) in &simplified_h.terms { if coeff.im.abs() > 1e-10 { @@ -642,7 +642,7 @@ impl Hamiltonian { } let mut simplified_h = self.clone(); - simplified_h.simplify(); + simplified_h.simplify_with_tolerance(0.0); for (_, coeff) in &simplified_h.terms { if coeff.im.abs() > 1e-10 { diff --git a/crates/cqlib-core/src/qis/hamiltonian.rs b/crates/cqlib-core/src/qis/hamiltonian.rs index 71a1819d..015ecfbd 100644 --- a/crates/cqlib-core/src/qis/hamiltonian.rs +++ b/crates/cqlib-core/src/qis/hamiltonian.rs @@ -190,6 +190,11 @@ impl Hamiltonian { /// /// This method is crucial for optimizing performance before executing quantum simulations. pub fn simplify(&mut self) { + self.simplify_with_tolerance(1e-10); + } + + // Evolution must merge terms without pruning small, nonzero coefficients. + pub(crate) fn simplify_with_tolerance(&mut self, tolerance: f64) { if self.terms.is_empty() { return; } @@ -235,7 +240,9 @@ impl Hamiltonian { } else { // New operator found, save the previous one if its coefficient is significant // Filter out terms with a coefficient near zero - if current_coeff.norm() > 1e-10 { + if (tolerance == 0.0 && current_coeff != Complex64::new(0.0, 0.0)) + || current_coeff.norm() > tolerance + { merged.push((current_op, current_coeff)); } current_op = op; @@ -243,7 +250,9 @@ impl Hamiltonian { } } // Push the last accumulated term - if current_coeff.norm() > 1e-10 { + if (tolerance == 0.0 && current_coeff != Complex64::new(0.0, 0.0)) + || current_coeff.norm() > tolerance + { merged.push((current_op, current_coeff)); } } diff --git a/crates/cqlib-core/src/qis/observable.rs b/crates/cqlib-core/src/qis/observable.rs index 26b8712d..ac0e4958 100644 --- a/crates/cqlib-core/src/qis/observable.rs +++ b/crates/cqlib-core/src/qis/observable.rs @@ -76,6 +76,8 @@ const VARIANCE_TOLERANCE: f64 = 1e-10; /// /// Methods may return [`QisError::QubitMismatch`] if the observable /// and the quantum state have incompatible qubit counts. +/// Real expectation methods return [`QisError::NotHermitian`] for non-Hermitian +/// observables; Hamiltonian coefficients must also be finite. pub trait Observable { /// Computes the expectation value for a statevector: $\langle\psi|O|\psi\rangle$. /// @@ -183,6 +185,30 @@ fn finalize_variance(variance: f64) -> Result { } } +// Validate the operator itself, including cancellations between repeated terms. +// Do not prune small coefficients while checking a real-valued API's contract. +fn ensure_hermitian_hamiltonian(h: &Hamiltonian) -> Result<(), QisError> { + let finite = |c: &Complex64| c.re.is_finite() && c.im.is_finite(); + if h.terms.iter().any(|(_, c)| !finite(c)) { + return Err(QisError::InvalidParameterValue( + "observable coefficients must be finite".into(), + )); + } + let mut merged = h.clone(); + merged.simplify_with_tolerance(0.0); + for (_, coeff) in &merged.terms { + if !finite(coeff) { + return Err(QisError::InvalidParameterValue( + "merged observable coefficients must be finite".into(), + )); + } + if coeff.im.abs() > VARIANCE_TOLERANCE { + return Err(QisError::NotHermitian); + } + } + Ok(()) +} + impl Observable for Hamiltonian { fn expectation_statevector(&self, sv: &Statevector) -> Result { if sv.num_qubits != self.num_qubits { @@ -192,6 +218,7 @@ impl Observable for Hamiltonian { }); } + ensure_hermitian_hamiltonian(self)?; let mut expected_value = 0.0; // Iterate over each term in the Hamiltonian @@ -240,6 +267,7 @@ impl Observable for Hamiltonian { }); } + ensure_hermitian_hamiltonian(self)?; let mut expected_value = 0.0; let n = self.num_qubits; let dim = 1 << n; // 2^n @@ -286,6 +314,7 @@ impl Observable for Hamiltonian { &self, measurements: &[(PauliString, HashMap)], ) -> Result { + ensure_hermitian_hamiltonian(self)?; let mut expected_value = 0.0; for (term_pauli, coeff) in &self.terms { @@ -368,8 +397,9 @@ impl Observable for Hamiltonian { }); } + ensure_hermitian_hamiltonian(self)?; let mut simplified = self.clone(); - simplified.simplify(); + simplified.simplify_with_tolerance(0.0); for (_, coeff) in &mut simplified.terms { if coeff.im.abs() > VARIANCE_TOLERANCE { @@ -401,6 +431,7 @@ impl Observable for PauliString { }); } + self.ensure_hermitian()?; let x_mask = self.x_mask(); let z_mask = self.z_mask(); let y_phase = self.y_phase(); @@ -433,6 +464,7 @@ impl Observable for PauliString { }); } + self.ensure_hermitian()?; let n = self.num_qubits; let dim = 1 << n; let x_mask = self.x_mask(); @@ -461,6 +493,7 @@ impl Observable for PauliString { &self, measurements: &[(PauliString, HashMap)], ) -> Result { + self.ensure_hermitian()?; let term_x_mask = self.x_mask(); let term_z_mask = self.z_mask(); let p_active = term_x_mask | term_z_mask; diff --git a/crates/cqlib-core/src/qis/pauli.rs b/crates/cqlib-core/src/qis/pauli.rs index 543efb97..12a8a5c1 100644 --- a/crates/cqlib-core/src/qis/pauli.rs +++ b/crates/cqlib-core/src/qis/pauli.rs @@ -859,10 +859,12 @@ impl PauliString { /// let exp = ps.expectation(&probs).unwrap(); /// assert!((exp).abs() < 1e-10); /// ``` + /// Returns `QisError::NotHermitian` for an imaginary Pauli phase. pub fn expectation( &self, probs: &std::collections::HashMap, ) -> Result { + self.ensure_hermitian()?; // Check if this Pauli string contains X or Y operators // X is represented by x[i]=1, z[i]=0 // Y is represented by x[i]=1, z[i]=1 @@ -914,6 +916,13 @@ impl PauliString { Ok((exp_value * global_phase).re) } + pub(crate) fn ensure_hermitian(&self) -> Result<(), crate::qis::error::QisError> { + match self.phase { + Phase::Plus | Phase::Minus => Ok(()), + Phase::I | Phase::MinusI => Err(crate::qis::error::QisError::NotHermitian), + } + } + /// Computes the phase exponent for multiplying two Pauli operators. /// /// Returns the exponent `n` such that `P1 * P2 = i^n * P3`. diff --git a/crates/cqlib-core/src/qis/state/aligned_buffer.rs b/crates/cqlib-core/src/qis/state/aligned_buffer.rs index 7784083e..608f5d52 100644 --- a/crates/cqlib-core/src/qis/state/aligned_buffer.rs +++ b/crates/cqlib-core/src/qis/state/aligned_buffer.rs @@ -61,6 +61,15 @@ impl AlignedBuffer { AlignedBuffer { ptr, len, layout } } + /// Fallible allocation used by checked state constructors. + pub(super) fn try_new_zeroed(len: usize) -> Option { + let size = len.checked_mul(size_of::())?; + let layout = Layout::from_size_align(size.max(1), 64.max(align_of::())).ok()?; + // SAFETY: layout is nonzero and validated; null is returned as failure. + let ptr = NonNull::new(unsafe { alloc_zeroed(layout) } as *mut T)?; + Some(Self { ptr, len, layout }) + } + /// Returns the raw pointer to the first element (64-byte aligned). #[allow(dead_code)] #[inline(always)] diff --git a/crates/cqlib-core/src/qis/state/density_matrix.rs b/crates/cqlib-core/src/qis/state/density_matrix.rs index b3560659..fbec87a1 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix.rs @@ -197,17 +197,35 @@ impl DensityMatrix { /// /// # Arguments /// * `num_qubits` - Number of qubits in the system. + /// + /// # Panics + /// + /// Panics if the dimension is unrepresentable or allocation fails. + /// Use [`Self::try_new`] to receive an error instead. pub fn new(num_qubits: usize) -> Self { - let size = 1 << (2 * num_qubits); - let mut data = vec![Complex64::new(0.0, 0.0); size]; - if size > 0 { - data[0] = Complex64::new(1.0, 0.0); - } - Self { + Self::try_new(num_qubits).expect("cannot allocate density matrix") + } + + /// Creates the ground state, returning an error for an unrepresentable + /// dimension or a failed allocation. + /// + /// # Errors + /// + /// Returns [`QisError::InvalidParameterValue`] for an unrepresentable + /// dimension, or [`QisError::UnsupportedOperation`] if allocation fails. + pub fn try_new(num_qubits: usize) -> Result { + let size = super::checked_state_len(num_qubits, true)?; + let mut data = Vec::new(); + data.try_reserve_exact(size).map_err(|_| { + QisError::UnsupportedOperation("density matrix allocation failed".into()) + })?; + data.resize(size, Complex64::new(0.0, 0.0)); + data[0] = Complex64::new(1.0, 0.0); + Ok(Self { data, num_qubits, qubit_map: None, - } + }) } /// Creates the maximally mixed state $I / 2^N$. @@ -218,18 +236,29 @@ impl DensityMatrix { /// /// Like every dense density-matrix constructor, this allocates $O(4^N)$ /// complex values. + /// + /// # Panics + /// + /// Panics if the dimension is unrepresentable or allocation fails. + /// Use [`Self::try_maximally_mixed`] to receive an error instead. pub fn maximally_mixed(num_qubits: usize) -> Self { - let dim = 1 << num_qubits; - let mut data = vec![Complex64::new(0.0, 0.0); dim * dim]; + Self::try_maximally_mixed(num_qubits).expect("cannot allocate density matrix") + } + + /// Fallible variant of [`Self::maximally_mixed`]. + /// + /// # Errors + /// + /// Returns [`QisError::InvalidParameterValue`] for an unrepresentable + /// dimension, or [`QisError::UnsupportedOperation`] if allocation fails. + pub fn try_maximally_mixed(num_qubits: usize) -> Result { + let mut state = Self::try_new(num_qubits)?; + let dim = 1usize << num_qubits; let probability = Complex64::new(1.0 / dim as f64, 0.0); for index in 0..dim { - data[index * dim + index] = probability; - } - Self { - data, - num_qubits, - qubit_map: None, + state.data[index * dim + index] = probability; } + Ok(state) } /// Creates a density matrix from an initial statevector (pure state). @@ -244,17 +273,17 @@ impl DensityMatrix { /// Returns [`QisError::InvalidStateDimension`] if the `initial_state` length is /// incorrect, or [`QisError::NotNormalized`] if it is not normalized. pub fn from_state(num_qubits: usize, initial_state: Vec) -> Result { - let dim = 1 << num_qubits; + super::checked_state_len(num_qubits, true)?; + let dim = 1usize << num_qubits; if initial_state.len() != dim { return Err(QisError::InvalidStateDimension(initial_state.len())); } let norm: f64 = initial_state.iter().map(|c| c.norm_sqr()).sum(); - if (norm - 1.0).abs() >= 1e-10 { + if !norm.is_finite() || (norm - 1.0).abs() >= 1e-10 { return Err(QisError::NotNormalized); } - let size = 1 << (2 * num_qubits); - let mut data = vec![Complex64::new(0.0, 0.0); size]; + let mut data = Self::try_new(num_qubits)?.data; let kernel = |(i, chunk): (usize, &mut [Complex64])| { let alpha_i = initial_state[i]; @@ -293,7 +322,7 @@ impl DensityMatrix { num_qubits: usize, dm_state: Vec, ) -> Result { - let size = 1 << (2 * num_qubits); + let size = super::checked_state_len(num_qubits, true)?; if dm_state.len() != size { return Err(QisError::InvalidStateDimension(dm_state.len())); } @@ -337,7 +366,7 @@ impl DensityMatrix { /// * `Ok(DensityMatrix)` - The resulting density matrix after execution. /// * `Err(QisError)` - If the circuit contains unsupported operations. pub fn from_circuit(circuit: &Circuit) -> Result { - let mut dm = DensityMatrix::new(circuit.num_qubits()); + let mut dm = DensityMatrix::try_new(circuit.num_qubits())?; dm.apply_circuit(circuit)?; Ok(dm) } @@ -573,13 +602,25 @@ impl DensityMatrix { /// /// # Arguments /// * `num_qubits` - Number of qubits in the system. + /// + /// # Panics + /// + /// Panics if the dimension is unrepresentable or allocation fails. + /// Use [`Self::try_zeros`] to receive an error instead. pub fn zeros(num_qubits: usize) -> Self { - let size = 1 << (2 * num_qubits); - Self { - data: vec![Complex64::new(0.0, 0.0); size], - num_qubits, - qubit_map: None, - } + Self::try_zeros(num_qubits).expect("cannot allocate density matrix") + } + + /// Fallible variant of [`Self::zeros`]. + /// + /// # Errors + /// + /// Returns [`QisError::InvalidParameterValue`] for an unrepresentable + /// dimension, or [`QisError::UnsupportedOperation`] if allocation fails. + pub fn try_zeros(num_qubits: usize) -> Result { + let mut state = Self::try_new(num_qubits)?; + state.data[0] = Complex64::new(0.0, 0.0); + Ok(state) } /// Computes the trace of the density matrix. @@ -611,6 +652,13 @@ impl DensityMatrix { /// # Returns /// `true` if the matrix is Hermitian within the specified tolerance. pub fn is_hermitian(&self, tol: f64) -> bool { + if self + .data + .iter() + .any(|v| !v.re.is_finite() || !v.im.is_finite()) + { + return false; + } let dim = 1 << self.num_qubits; for i in 0..dim { for j in 0..i { diff --git a/crates/cqlib-core/src/qis/state/density_matrix_noise.rs b/crates/cqlib-core/src/qis/state/density_matrix_noise.rs index 7e0d41b1..d76fe5c0 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix_noise.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix_noise.rs @@ -92,11 +92,26 @@ impl DensityMatrixNoise { /// // Simulator with empty noise model /// let sim = DensityMatrixNoise::new(2, Some(cqlib_core::device::NoiseModel::new())); /// ``` + /// + /// # Panics + /// + /// Panics if the dimension is unrepresentable or allocation fails. + /// Use [`Self::try_new`] to receive an error instead. pub fn new(num_qubits: usize, noise_model: Option) -> Self { - Self { - state: DensityMatrix::new(num_qubits), + Self::try_new(num_qubits, noise_model).expect("cannot allocate noisy density matrix") + } + + /// Fallible ground-state constructor, preserving the supplied noise model. + /// + /// # Errors + /// + /// Returns [`QisError::InvalidParameterValue`] for an unrepresentable + /// dimension, or [`QisError::UnsupportedOperation`] if allocation fails. + pub fn try_new(num_qubits: usize, noise_model: Option) -> Result { + Ok(Self { + state: DensityMatrix::try_new(num_qubits)?, noise_model, - } + }) } /// Simulates a circuit, applying noise after each operation. @@ -139,7 +154,7 @@ impl DensityMatrixNoise { noise_model: Option, ) -> Result { let mut dms = Self { - state: DensityMatrix::new(circuit.num_qubits()), + state: DensityMatrix::try_new(circuit.num_qubits())?, noise_model, }; dms.apply_circuit(circuit)?; @@ -199,114 +214,130 @@ impl DensityMatrixNoise { .map(|p| p.evaluate(&None).ok()) .collect(); - for op in circuit.operations() { - let params: Vec = op - .params - .iter() - .map(|p| match p { - CircuitParam::Fixed(v) => Ok(*v), - CircuitParam::Index(idx) => parameter_values - .get(*idx as usize) - .copied() - .flatten() - .ok_or(QisError::CircuitError( - crate::circuit::CircuitError::SymbolicParameterError, - )), - }) - .collect::, QisError>>()?; - - let qs: Result, QisError> = op - .qubits - .iter() - .map(|q| { - qubit_map.get(q).copied().ok_or_else(|| { - QisError::CircuitError(crate::circuit::CircuitError::QubitNotFound(q.id())) + // Expose labels to gate-noise lookup during execution, restoring the + // previous binding on error (already executed gates are not rolled back). + let previous_map = sim.state.qubit_map.replace(qubit_map); + let result = (|| { + for op in circuit.operations() { + let params: Vec = op + .params + .iter() + .map(|p| match p { + CircuitParam::Fixed(v) => Ok(*v), + CircuitParam::Index(idx) => parameter_values + .get(*idx as usize) + .copied() + .flatten() + .ok_or(QisError::CircuitError( + crate::circuit::CircuitError::SymbolicParameterError, + )), }) - }) - .collect(); - let qs = qs?; + .collect::, QisError>>()?; + + let qs: Result, QisError> = op + .qubits + .iter() + .map(|q| { + sim.state + .qubit_map + .as_ref() + .unwrap() + .get(q) + .copied() + .ok_or_else(|| { + QisError::CircuitError(crate::circuit::CircuitError::QubitNotFound( + q.id(), + )) + }) + }) + .collect(); + let qs = qs?; - match &op.instruction { - Instruction::Standard(gate) => { - sim.apply_standard_gate_noise(*gate, &qs, ¶ms)?; - } - Instruction::McGate(mc_gate) => { - let num_controls = mc_gate.num_ctrl_qubits(); - let base_gate = mc_gate.base_gate(); - - if num_controls == 1 { - let c = qs[0]; - let t = qs[1]; - match base_gate { - StandardGate::X => sim.apply_cx(c, t)?, - StandardGate::Y => sim.apply_cy(c, t)?, - StandardGate::Z => sim.apply_cz(c, t)?, - StandardGate::RX => sim.apply_crx(c, t, params[0])?, - StandardGate::RY => sim.apply_cry(c, t, params[0])?, - StandardGate::RZ => sim.apply_crz(c, t, params[0])?, - _ => { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError( - crate::circuit::CircuitError::NoMatrixRepresentation, - ) - })?; - sim.apply_unitary_gate(&qs, &matrix)?; - sim.apply_noise(*base_gate, &qs)?; + match &op.instruction { + Instruction::Standard(gate) => { + sim.apply_standard_gate_noise(*gate, &qs, ¶ms)?; + } + Instruction::McGate(mc_gate) => { + let num_controls = mc_gate.num_ctrl_qubits(); + let base_gate = mc_gate.base_gate(); + + if num_controls == 1 { + let c = qs[0]; + let t = qs[1]; + match base_gate { + StandardGate::X => sim.apply_cx(c, t)?, + StandardGate::Y => sim.apply_cy(c, t)?, + StandardGate::Z => sim.apply_cz(c, t)?, + StandardGate::RX => sim.apply_crx(c, t, params[0])?, + StandardGate::RY => sim.apply_cry(c, t, params[0])?, + StandardGate::RZ => sim.apply_crz(c, t, params[0])?, + _ => { + let matrix = mc_gate.matrix(¶ms).map_err(|_| { + QisError::CircuitError( + crate::circuit::CircuitError::NoMatrixRepresentation, + ) + })?; + sim.apply_unitary_gate(&qs, &matrix)?; + sim.apply_noise(*base_gate, &qs)?; + } } + } else if num_controls == 2 && *base_gate == StandardGate::X { + sim.apply_ccx(qs[0], qs[1], qs[2])?; + } else { + let matrix = mc_gate.matrix(¶ms).map_err(|_| { + QisError::CircuitError( + crate::circuit::CircuitError::NoMatrixRepresentation, + ) + })?; + sim.apply_unitary_gate(&qs, &matrix)?; + sim.apply_noise(*base_gate, &qs)?; } - } else if num_controls == 2 && *base_gate == StandardGate::X { - sim.apply_ccx(qs[0], qs[1], qs[2])?; - } else { - let matrix = mc_gate.matrix(¶ms).map_err(|_| { - QisError::CircuitError( + } + Instruction::UnitaryGate(u_gate) => { + if let Some(matrix) = u_gate.matrix() { + sim.apply_unitary_gate(&qs, matrix)?; + sim.apply_noise(StandardGate::U, &qs)?; + } else { + return Err(QisError::CircuitError( crate::circuit::CircuitError::NoMatrixRepresentation, - ) - })?; - sim.apply_unitary_gate(&qs, &matrix)?; - sim.apply_noise(*base_gate, &qs)?; + )); + } } - } - Instruction::UnitaryGate(u_gate) => { - if let Some(matrix) = u_gate.matrix() { - sim.apply_unitary_gate(&qs, matrix)?; - sim.apply_noise(StandardGate::U, &qs)?; - } else { + Instruction::CircuitGate(_) => { return Err(QisError::CircuitError( - crate::circuit::CircuitError::NoMatrixRepresentation, + crate::circuit::CircuitError::InvalidOperation( + "CircuitGate should have been decomposed".to_string(), + ), )); } - } - Instruction::CircuitGate(_) => { - return Err(QisError::CircuitError( - crate::circuit::CircuitError::InvalidOperation( - "CircuitGate should have been decomposed".to_string(), - ), - )); - } - Instruction::Directive(directive) => match directive { - Directive::Barrier => {} - Directive::Measure => { - return Err(QisError::UnsupportedOperation( + Instruction::Directive(directive) => match directive { + Directive::Barrier => {} + Directive::Measure => { + return Err(QisError::UnsupportedOperation( "legacy Measure directive is not supported by noisy density matrix circuit evolution; use Circuit::measure as an output declaration".to_string(), )); - } - Directive::Reset => { - let qubit = qs[0]; - sim.reset(qubit)?; - } - }, - Instruction::Delay => continue, - Instruction::ClassicalControl(_) => { - return Err(QisError::UnsupportedOperation( + } + Directive::Reset => { + let qubit = qs[0]; + sim.reset(qubit)?; + } + }, + Instruction::Delay => continue, + Instruction::ClassicalControl(_) => { + return Err(QisError::UnsupportedOperation( "classical control flow is not supported in noisy density matrix simulation" .to_string(), )); + } + Instruction::ClassicalData(_) => {} } - Instruction::ClassicalData(_) => {} } + Ok(()) + })(); + if result.is_err() { + sim.state.qubit_map = previous_map; } - sim.state.qubit_map = Some(qubit_map); - Ok(()) + result } /// Applies a standard gate and then applies any matching gate noise. @@ -342,7 +373,15 @@ impl DensityMatrixNoise { ))); } match gate { - StandardGate::I => {} + StandardGate::I => { + if qs[0] >= self.state.num_qubits { + return Err(QisError::IndexOutOfBounds { + index: qs[0], + max: self.state.num_qubits.saturating_sub(1), + }); + } + self.apply_noise(StandardGate::I, qs)?; + } StandardGate::X => self.apply_x(qs[0])?, StandardGate::Y => self.apply_y(qs[0])?, StandardGate::Z => self.apply_z(qs[0])?, @@ -399,9 +438,19 @@ impl DensityMatrixNoise { /// converts them to Kraus operators, and applies them to the state. /// Supports single-qubit, two-qubit, and three-qubit gates. fn apply_noise(&mut self, gate: StandardGate, qubits: &[usize]) -> Result<(), QisError> { + let labelled_qubit = |position: usize| { + self.state + .qubit_map + .as_ref() + .and_then(|map| { + map.iter() + .find_map(|(&qubit, &index)| (index == position).then_some(qubit)) + }) + .unwrap_or_else(|| Qubit::new(position as u32)) + }; if let Some(noise_model) = &self.noise_model { if qubits.len() == 1 { - let q0 = Qubit::new(qubits[0] as u32); + let q0 = labelled_qubit(qubits[0]); let key = OperationKey::new_single(gate, q0); if let Some(errors) = noise_model.get_single_qubit_errors(&key) { for error in errors { @@ -411,8 +460,8 @@ impl DensityMatrixNoise { } } } else if qubits.len() == 2 { - let q0 = Qubit::new(qubits[0] as u32); - let q1 = Qubit::new(qubits[1] as u32); + let q0 = labelled_qubit(qubits[0]); + let q1 = labelled_qubit(qubits[1]); if let Ok(key) = OperationKey::new_double(gate, q0, q1) && let Some(errors) = noise_model.get_two_qubit_errors(&key) { @@ -423,9 +472,9 @@ impl DensityMatrixNoise { } } } else if qubits.len() == 3 { - let q0 = Qubit::new(qubits[0] as u32); - let q1 = Qubit::new(qubits[1] as u32); - let q2 = Qubit::new(qubits[2] as u32); + let q0 = labelled_qubit(qubits[0]); + let q1 = labelled_qubit(qubits[1]); + let q2 = labelled_qubit(qubits[2]); if let Ok(_key) = OperationKey::new_triple(gate, q0, q1, q2) { // Current noise model struct doesn't have `get_three_qubit_errors`. // Ready for future extension without silently failing or panicking. diff --git a/crates/cqlib-core/src/qis/state/density_matrix_test.rs b/crates/cqlib-core/src/qis/state/density_matrix_test.rs index 34bc2572..91595c4b 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix_test.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix_test.rs @@ -1226,3 +1226,17 @@ fn test_psd_inf_rejected() { }; assert!(!dm.is_positive_semidefinite_approx(1e-10)); } + +#[test] +fn hermiticity_rejects_nonfinite_entries() { + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + for index in 0..4 { + for value in [Complex64::new(bad, 0.0), Complex64::new(0.0, bad)] { + let mut dm = DensityMatrix::new(1); + dm.data[index] = value; + assert!(!dm.is_hermitian(1e-10)); + assert!(dm.validate_physical(1e-10).is_err()); + } + } + } +} diff --git a/crates/cqlib-core/src/qis/state/mod.rs b/crates/cqlib-core/src/qis/state/mod.rs index f5216562..6c23a9a2 100644 --- a/crates/cqlib-core/src/qis/state/mod.rs +++ b/crates/cqlib-core/src/qis/state/mod.rs @@ -172,3 +172,19 @@ pub use statevector::Statevector; #[cfg(test)] #[path = "state_test.rs"] mod state_test; + +// Validate both the exponent and allocation layout before shifting or allocating. +fn checked_state_len(num_qubits: usize, density: bool) -> Result { + let invalid = || { + QisError::InvalidParameterValue(format!( + "state dimension is too large for {num_qubits} qubits" + )) + }; + let exponent = num_qubits + .checked_mul(if density { 2 } else { 1 }) + .ok_or_else(invalid)?; + let exponent = u32::try_from(exponent).map_err(|_| invalid())?; + let len = 1usize.checked_shl(exponent).ok_or_else(invalid)?; + std::alloc::Layout::array::(len).map_err(|_| invalid())?; + Ok(len) +} diff --git a/crates/cqlib-core/src/qis/state/statevector.rs b/crates/cqlib-core/src/qis/state/statevector.rs index 307972d7..9b94a36a 100644 --- a/crates/cqlib-core/src/qis/state/statevector.rs +++ b/crates/cqlib-core/src/qis/state/statevector.rs @@ -351,16 +351,34 @@ impl Statevector { /// assert_eq!(sv.num_qubits, 2); /// assert_eq!(sv.data().len(), 4); // 2^2 amplitudes /// ``` + /// + /// # Panics + /// + /// Panics if the dimension is unrepresentable or allocation fails. + /// Use [`Self::try_new`] to receive an error instead. pub fn new(num_qubits: usize) -> Self { - let size = 1 << num_qubits; + Self::try_new(num_qubits).expect("cannot allocate statevector") + } + + /// Creates the ground state, returning an error for an unrepresentable + /// dimension or a failed allocation. + /// + /// # Errors + /// + /// Returns [`QisError::InvalidParameterValue`] for an unrepresentable + /// dimension, or [`QisError::UnsupportedOperation`] if allocation fails. + pub fn try_new(num_qubits: usize) -> Result { + let size = super::checked_state_len(num_qubits, false)?; // 64-byte aligned allocation — zero-initialised (all amplitudes = 0+0i). - let mut data = AlignedBuffer::::new_zeroed(size); + let mut data = AlignedBuffer::::try_new_zeroed(size).ok_or_else(|| { + QisError::UnsupportedOperation("statevector allocation failed".into()) + })?; data[0] = Complex64::new(1.0, 0.0); // |0...0⟩ - Statevector { + Ok(Statevector { data, num_qubits, qubit_map: None, - } + }) } /// Creates a statevector from initial amplitudes with normalization check. @@ -389,7 +407,7 @@ impl Statevector { /// let sv = Statevector::from_state(2, amps); /// ``` pub fn from_state(num_qubits: usize, initial_state: Vec) -> Result { - let size = 1 << num_qubits; + let size = super::checked_state_len(num_qubits, false)?; if initial_state.len() != size { return Err(QisError::InvalidStateDimension(initial_state.len())); } @@ -397,12 +415,14 @@ impl Statevector { // summing many amplitudes accumulates floating-point roundoff. let norm: f64 = initial_state.iter().map(|c| c.norm_sqr()).sum(); let tolerance = 1e-10 * (size.max(1) as f64); - if (norm - 1.0).abs() >= tolerance { + if !norm.is_finite() || (norm - 1.0).abs() >= tolerance { return Err(QisError::NotNormalized); } // Copy caller-provided Vec into 64-byte aligned buffer. - let mut data = AlignedBuffer::::new_zeroed(size); + let mut data = AlignedBuffer::::try_new_zeroed(size).ok_or_else(|| { + QisError::UnsupportedOperation("statevector allocation failed".into()) + })?; data.as_mut_slice().copy_from_slice(&initial_state); Ok(Statevector { data, @@ -444,7 +464,7 @@ impl Statevector { /// let sv = Statevector::from_circuit(&circuit).unwrap(); /// ``` pub fn from_circuit(circuit: &Circuit) -> Result { - let mut sv = Statevector::new(circuit.num_qubits()); + let mut sv = Statevector::try_new(circuit.num_qubits())?; sv.apply_circuit(circuit)?; Ok(sv) } diff --git a/crates/cqlib-core/tests/p1_regressions.rs b/crates/cqlib-core/tests/p1_regressions.rs new file mode 100644 index 00000000..045c64ee --- /dev/null +++ b/crates/cqlib-core/tests/p1_regressions.rs @@ -0,0 +1,232 @@ +// This code is part of Cqlib. +// +// (C) Copyright China Telecom Quantum Group 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. + +//! Regression tests for numerical boundaries and labelled noisy simulation. +use cqlib_core::circuit::{Circuit, Qubit, StandardGate, circuit_to_matrix}; +use cqlib_core::device::{ + NoiseModel, + noise::{SingleQubitNoise, TwoQubitNoise}, +}; +use cqlib_core::qis::entropy::{entanglement_of_formation, renyi_entropy}; +use cqlib_core::qis::evolution::TrotterMode; +use cqlib_core::qis::{DensityMatrix, DensityMatrixNoise, Hamiltonian}; +use num_complex::Complex64; + +#[test] +fn pauli_probability_boundary_is_trace_preserving() { + let noise = SingleQubitNoise::Pauli { + px: 0.3, + py: 0.6, + pz: 0.1, + }; + assert!(noise.is_valid()); + let ops = noise.to_kraus(); + let mut completeness = ndarray::Array2::::zeros((2, 2)); + for op in &ops { + assert!(op.iter().all(|v| v.re.is_finite() && v.im.is_finite())); + completeness += &op.t().mapv(|v| v.conj()).dot(op); + } + for i in 0..2 { + for j in 0..2 { + assert!( + (completeness[[i, j]] - Complex64::from(if i == j { 1.0 } else { 0.0 })).norm() + < 1e-12 + ); + } + } + let mut model = NoiseModel::new(); + model + .add_single_qubit_error(StandardGate::X, Qubit::new(0), noise) + .unwrap(); + let mut sim = DensityMatrixNoise::new(1, Some(model)); + sim.apply_x(0).unwrap(); + assert!((sim.probabilities()[0] - 0.9).abs() < 1e-12); + assert!((sim.probabilities().iter().sum::() - 1.0).abs() < 1e-12); + assert!( + !SingleQubitNoise::Pauli { + px: 0.4, + py: 0.6, + pz: 0.1 + } + .is_valid() + ); +} + +#[test] +fn gate_noise_uses_original_labels_and_operand_order() { + for ids in [[0, 1], [5, 12], [12, 5]] { + let qs = ids.map(Qubit::new); + let mut circuit = Circuit::from_qubits(qs.to_vec()).unwrap(); + circuit.x(qs[0]).unwrap(); + let mut model = NoiseModel::new(); + model + .add_single_qubit_error(StandardGate::X, qs[0], SingleQubitNoise::BitFlip(1.0)) + .unwrap(); + let mut sim = DensityMatrixNoise::from_circuit(&circuit, Some(model)).unwrap(); + assert!((sim.probabilities()[0] - 1.0).abs() < 1e-12); + sim.apply_x(0).unwrap(); + assert!((sim.probabilities()[0] - 1.0).abs() < 1e-12); + for reverse in [false, true] { + let mut circuit = Circuit::from_qubits(qs.to_vec()).unwrap(); + let [control, target] = if reverse { [qs[1], qs[0]] } else { qs }; + circuit.cx(control, target).unwrap(); + let mut model = NoiseModel::new(); + model + .add_two_qubit_error( + StandardGate::CX, + control, + target, + TwoQubitNoise::Independent { + q0_noise: SingleQubitNoise::BitFlip(1.0), + q1_noise: SingleQubitNoise::BitFlip(0.0), + }, + ) + .unwrap(); + let sim = DensityMatrixNoise::from_circuit(&circuit, Some(model)).unwrap(); + let expected = if reverse { 2 } else { 1 }; + assert!((sim.probabilities()[expected] - 1.0).abs() < 1e-12); + } + } +} + +#[test] +fn formation_entropy_at_concurrence_roundoff_boundary() { + let state = vec![ + Complex64::new(-0.5917562523120286, -0.14345273890626087), + Complex64::new(0.12739228627789673, 0.3361800930257551), + Complex64::new(-0.06600502193827436, 0.35339664207352467), + Complex64::new(-0.5571108451004554, 0.2457267928442915), + ]; + let dm = DensityMatrix::from_state(2, state).unwrap(); + assert!((entanglement_of_formation(&dm).unwrap() - 1.0).abs() < 1e-12); +} + +#[test] +fn renyi_extreme_orders_and_small_positive_eigenvalues() { + for n in 1..=3 { + let dm = DensityMatrix::maximally_mixed(n); + for alpha in [ + 0.1, + 1.0 - 1e-12, + 1.0, + 1.0 + 1e-12, + 2000.0, + f64::MAX, + f64::INFINITY, + ] { + assert!( + (renyi_entropy(&dm, alpha).unwrap() - n as f64).abs() < 1e-10, + "n={n}, alpha={alpha}" + ); + assert!(renyi_entropy(&DensityMatrix::new(n), alpha).unwrap().abs() < 1e-10); + } + } + let dm = DensityMatrix::from_density_matrix_state( + 1, + vec![(1.0 - 1e-14).into(), 0.0.into(), 0.0.into(), 1e-14.into()], + ) + .unwrap(); + assert!( + (renyi_entropy(&dm, 0.1).unwrap() - 0.06257881098503021).abs() < 1e-12, + "got {}", + renyi_entropy(&dm, 0.1).unwrap() + ); + assert!(renyi_entropy(&dm, f64::NAN).is_err()); +} + +#[test] +fn evolution_retains_small_coefficients_before_time_scaling() { + for coefficients in [vec![1e-11], vec![2e-11, -1e-11], vec![1e-11, -1e-11]] { + let mut h = Hamiltonian::new(1); + for &c in &coefficients { + h.add_term("Z".parse().unwrap(), c.into()).unwrap(); + } + let angle = coefficients.iter().sum::() * 1e11; + for mode in [TrotterMode::FirstOrder, TrotterMode::SecondOrder] { + for circuit in [ + h.to_trotter_circuit(1e11, 3, mode).unwrap(), + h.to_evolution_circuit(1e11, 3, mode).unwrap(), + ] { + let m = circuit_to_matrix(&circuit, None).unwrap(); + assert!((m[[0, 0]] - Complex64::from_polar(1.0, -angle)).norm() < 1e-10); + assert!((m[[1, 1]] - Complex64::from_polar(1.0, angle)).norm() < 1e-10); + assert!(m[[0, 1]].norm() < 1e-10 && m[[1, 0]].norm() < 1e-10); + } + } + h.simplify(); + assert!(h.terms.is_empty()); + } +} + +#[test] +fn checked_constructors_reject_dimensions_without_allocating() { + use cqlib_core::qis::Statevector; + for n in [usize::BITS as usize - 4, usize::BITS as usize, usize::MAX] { + assert!(Statevector::try_new(n).is_err()); + assert!(Statevector::from_state(n, vec![]).is_err()); + } + for n in [ + (usize::BITS as usize - 4).div_ceil(2), + usize::BITS as usize, + usize::MAX, + ] { + assert!(DensityMatrix::try_new(n).is_err()); + assert!(DensityMatrixNoise::try_new(n, None).is_err()); + assert!(DensityMatrix::from_state(n, vec![]).is_err()); + assert!(DensityMatrix::from_density_matrix_state(n, vec![]).is_err()); + assert!(DensityMatrix::try_maximally_mixed(n).is_err()); + } + for n in [0, 1, 3] { + assert_eq!( + Statevector::try_new(n).unwrap().data(), + Statevector::new(n).data() + ); + assert_eq!( + DensityMatrix::try_new(n).unwrap().data(), + DensityMatrix::new(n).data() + ); + assert!( + (DensityMatrixNoise::try_new(n, None) + .unwrap() + .probabilities()[0] + - 1.0) + .abs() + < 1e-12 + ); + } +} + +#[test] +fn failed_noisy_circuit_preserves_label_binding() { + use cqlib_core::circuit::Parameter; + let q = Qubit::new(5); + let mut bad = Circuit::from_qubits(vec![q]).unwrap(); + bad.x(q).unwrap(); + bad.rz(q, Parameter::symbol("unbound")).unwrap(); + let mut model = NoiseModel::new(); + model + .add_single_qubit_error(StandardGate::X, q, SingleQubitNoise::BitFlip(1.0)) + .unwrap(); + let mut sim = DensityMatrixNoise::new(1, Some(model.clone())); + assert!(sim.apply_circuit(&bad).is_err()); + // Completed gates used label 5, but the failed application did not bind it. + assert!((sim.probabilities()[0] - 1.0).abs() < 1e-12); + sim.apply_circuit(&Circuit::new(1)).unwrap(); + + let mut good = Circuit::from_qubits(vec![q]).unwrap(); + good.x(q).unwrap(); + let mut bound = DensityMatrixNoise::from_circuit(&good, Some(model)).unwrap(); + assert!(bound.apply_circuit(&bad).is_err()); + bound.apply_x(0).unwrap(); + assert!((bound.probabilities()[0] - 1.0).abs() < 1e-12); + assert!(bound.apply_circuit(&Circuit::new(1)).is_err()); +} diff --git a/crates/cqlib-core/tests/p2_regressions.rs b/crates/cqlib-core/tests/p2_regressions.rs new file mode 100644 index 00000000..ffa6faaa --- /dev/null +++ b/crates/cqlib-core/tests/p2_regressions.rs @@ -0,0 +1,299 @@ +// This code is part of Cqlib. +// +// (C) Copyright China Telecom Quantum Group 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. + +//! Regressions for noisy identity gates, observable validation and circuit phase. +use cqlib_core::circuit::{Circuit, Instruction, Parameter, Qubit, StandardGate}; +use cqlib_core::device::{NoiseModel, noise::SingleQubitNoise}; +use cqlib_core::error_mitigation::ZNEMitigation; +use cqlib_core::qis::entropy::renyi_entropy; +use cqlib_core::qis::error::QisError; +use cqlib_core::qis::{ + DensityMatrix, DensityMatrixNoise, Hamiltonian, Observable, PauliString, Statevector, +}; +use num_complex::Complex64 as C; +use std::collections::HashMap; + +#[test] +fn identity_noise_respects_labels_and_repetition() { + for ids in [[0, 1], [12, 5]] { + for p in [0.0, 0.25, 1.0] { + for count in [1, 2, 3] { + let qs = ids.map(Qubit::new); + let mut c = Circuit::from_qubits(qs.to_vec()).unwrap(); + for _ in 0..count { + c.i(qs[1]).unwrap(); + } + let mut model = NoiseModel::new(); + model + .add_single_qubit_error(StandardGate::I, qs[1], SingleQubitNoise::BitFlip(p)) + .unwrap(); + let sim = DensityMatrixNoise::from_circuit(&c, Some(model)).unwrap(); + let flipped = (1.0 - (1.0 - 2.0 * p).powi(count)) / 2.0; + for (actual, expected) in + sim.probabilities() + .iter() + .zip([1.0 - flipped, 0.0, flipped, 0.0]) + { + assert!((actual - expected).abs() < 1e-12); + } + assert_eq!( + DensityMatrixNoise::from_circuit(&c, None) + .unwrap() + .probabilities(), + vec![1.0, 0.0, 0.0, 0.0] + ); + } + } + } + let mut sim = DensityMatrixNoise::new(1, None); + assert!( + sim.apply_standard_gate_noise(StandardGate::I, &[1], &[]) + .is_err() + ); +} + +#[test] +fn nonfinite_pure_states_are_rejected() { + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY, f64::MAX] { + for amp in [C::new(bad, 0.0), C::new(0.0, bad)] { + let state = vec![amp, C::new(0.0, 0.0)]; + assert!(matches!( + Statevector::from_state(1, state.clone()), + Err(QisError::NotNormalized) + )); + assert!(matches!( + DensityMatrix::from_state(1, state), + Err(QisError::NotNormalized) + )); + } + } + for n in [0, 1] { + let mut v = vec![C::new(0.0, 0.0); 1 << n]; + v[0] = C::new(0.0, 1.0); + assert!(Statevector::from_state(n, v.clone()).is_ok()); + assert!(DensityMatrix::from_state(n, v).is_ok()); + } +} + +#[test] +fn nonhermitian_observables_are_rejected_even_at_zero_expectation() { + let sv = Statevector::new(1); + let dm = DensityMatrix::new(1); + let probs = HashMap::from([("0".to_string(), 1.0)]); + for label in ["iZ", "-iI", "iX"] { + let p: PauliString = label.parse().unwrap(); + let h = Hamiltonian::from_pauli(p.clone()); + let measurements = vec![(p.clone(), probs.clone())]; + for obs in [&p as &dyn Observable, &h as &dyn Observable] { + assert!(matches!( + obs.expectation_statevector(&sv), + Err(QisError::NotHermitian) + )); + assert!(matches!( + obs.expectation_density_matrix(&dm), + Err(QisError::NotHermitian) + )); + assert!(matches!( + obs.expectation_probs(&measurements), + Err(QisError::NotHermitian) + )); + } + assert!(matches!(p.expectation(&probs), Err(QisError::NotHermitian))); + } +} + +#[test] +fn hermiticity_is_checked_after_phase_and_duplicate_merging() { + let sv = Statevector::new(1); + let dm = DensityMatrix::new(1); + let probs = vec![( + "Z".parse().unwrap(), + HashMap::from([("0".to_string(), 1.0)]), + )]; + for terms in [ + vec![("iZ", C::new(0.0, -1.0))], + vec![("Z", C::new(1.0, 2.0)), ("Z", C::new(0.0, -2.0))], + vec![("Z", C::new(1e-11, 0.0))], + ] { + let h = Hamiltonian::from_list( + terms + .iter() + .map(|(p, c)| (p.parse().unwrap(), *c)) + .collect(), + ) + .unwrap(); + let expected = if terms[0].1.re == 1e-11 { 1e-11 } else { 1.0 }; + for actual in [ + h.expectation_statevector(&sv).unwrap(), + h.expectation_density_matrix(&dm).unwrap(), + h.expectation_probs(&probs).unwrap(), + ] { + assert!((actual - expected).abs() < 1e-14); + } + } + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + let h = Hamiltonian::from_list(vec![("Z".parse().unwrap(), C::new(bad, 0.0))]).unwrap(); + assert!(matches!( + h.expectation_statevector(&sv), + Err(QisError::InvalidParameterValue(_)) + )); + assert!(matches!( + h.expectation_density_matrix(&dm), + Err(QisError::InvalidParameterValue(_)) + )); + assert!(matches!( + h.expectation_probs(&probs), + Err(QisError::InvalidParameterValue(_)) + )); + } + let h = Hamiltonian::from_pauli("iZZ".parse().unwrap()); + assert!(matches!( + h.expectation_statevector(&sv), + Err(QisError::QubitMismatch { .. }) + )); +} + +#[test] +fn folding_preserves_full_matrix_including_global_phase() { + for gates in [false, true] { + let mut c = Circuit::from_qubits(vec![Qubit::new(5)]).unwrap(); + c.set_global_phase(Parameter::from(0.37)); + if gates { + c.h(Qubit::new(5)).unwrap(); + c.t(Qubit::new(5)).unwrap(); + } + let expected = c.to_matrix(None).unwrap(); + let zne = ZNEMitigation::new(c, vec![0, 1, 2]); + let selected = [Instruction::Standard(StandardGate::H)]; + for selection in [None, Some(selected.as_slice()), Some([].as_slice())] { + for folded in zne.fold_circuits(selection).unwrap() { + let actual = folded.to_matrix(None).unwrap(); + assert!( + actual + .iter() + .zip(expected.iter()) + .all(|(a, b)| (*a - *b).norm() < 1e-12) + ); + } + } + } +} + +#[test] +fn hermiticity_tolerance_and_merged_overflow() { + for (imaginary, valid) in [(0.5e-10, true), (2e-10, false)] { + let h = + Hamiltonian::from_list(vec![("Z".parse().unwrap(), C::new(1.0, imaginary))]).unwrap(); + assert_eq!( + h.expectation_statevector(&Statevector::new(1)).is_ok(), + valid + ); + assert_eq!( + h.expectation_density_matrix(&DensityMatrix::new(1)).is_ok(), + valid + ); + } + let h = Hamiltonian::from_list(vec![("Z".parse().unwrap(), C::new(f64::MAX, 0.0)); 2]).unwrap(); + assert!(matches!( + h.expectation_statevector(&Statevector::new(1)), + Err(QisError::InvalidParameterValue(_)) + )); +} + +#[test] +fn renyi_entropy_of_nondiagonal_pure_state_is_zero() { + let state = [1.0, 2.0, 3.0, 4.0] + .map(|v| C::new(v / 30.0_f64.sqrt(), 0.0)) + .to_vec(); + let dm = DensityMatrix::from_state(2, state).unwrap(); + for alpha in [0.01, 0.1, 0.5, 1.0, 2.0, f64::INFINITY] { + let actual = renyi_entropy(&dm, alpha).unwrap(); + assert!(actual.abs() < 1e-10, "alpha={alpha}, entropy={actual}"); + } +} + +#[test] +fn renyi_entropy_preserves_rank_deficient_spectrum_under_unitary_rotation() { + // The first two columns of the four-dimensional Fourier matrix are + // orthonormal, so this state has spectrum [0.8, 0.2, 0, 0]. + let column = [ + C::new(1.0, 0.0), + C::new(0.0, 1.0), + C::new(-1.0, 0.0), + C::new(0.0, -1.0), + ]; + let mut rho = Vec::new(); + let mut diagonal = vec![C::new(0.0, 0.0); 16]; + diagonal[0] = 0.8.into(); + diagonal[5] = 0.2.into(); + for i in 0..4 { + for j in 0..4 { + rho.push((C::new(0.8, 0.0) + 0.2 * column[i] * column[j].conj()) / 4.0); + } + } + for data in [diagonal, rho] { + let dm = DensityMatrix::from_density_matrix_state(2, data).unwrap(); + for alpha in [0.01, 0.1, 0.5, 1.0, 2.0, f64::INFINITY] { + let expected = if alpha == 1.0 { + -0.8_f64 * 0.8_f64.log2() - 0.2_f64 * 0.2_f64.log2() + } else if alpha == f64::INFINITY { + -0.8_f64.log2() + } else { + (0.8_f64.powf(alpha) + 0.2_f64.powf(alpha)).log2() / (1.0 - alpha) + }; + assert!((renyi_entropy(&dm, alpha).unwrap() - expected).abs() < 1e-10); + } + } +} + +#[test] +fn variance_validates_coefficients_before_merging_and_retains_small_terms() { + let sv = Statevector::new(1); + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + for coeff in [C::new(bad, 0.0), C::new(0.0, bad)] { + let h = Hamiltonian::from_list(vec![("Z".parse().unwrap(), coeff)]).unwrap(); + assert!(matches!( + h.variance_statevector(&sv), + Err(QisError::InvalidParameterValue(_)) + )); + } + } + let h = Hamiltonian::from_list(vec![("Z".parse().unwrap(), C::new(f64::MAX, 0.0)); 2]).unwrap(); + assert!(matches!( + h.variance_statevector(&sv), + Err(QisError::InvalidParameterValue(_)) + )); + for terms in [ + vec![("iX", C::new(0.0, -1.0))], + vec![("X", C::new(1.0, 2.0)), ("X", C::new(0.0, -2.0))], + ] { + let h = Hamiltonian::from_list( + terms + .into_iter() + .map(|(p, c)| (p.parse().unwrap(), c)) + .collect(), + ) + .unwrap(); + assert!((h.variance_statevector(&sv).unwrap() - 1.0).abs() < 1e-12); + } + let h = Hamiltonian::from_list(vec![("X".parse().unwrap(), C::new(1e-11, 0.0))]).unwrap(); + assert!((h.variance_statevector(&sv).unwrap() - 1e-22).abs() < 1e-34); + let invalid = Hamiltonian::from_pauli("iZ".parse().unwrap()); + assert!(matches!( + invalid.variance_statevector(&sv), + Err(QisError::NotHermitian) + )); + assert!(matches!( + invalid.variance_statevector(&Statevector::new(2)), + Err(QisError::QubitMismatch { .. }) + )); +} diff --git a/docs/api/cn/c/3_qis/1_statevector.md b/docs/api/cn/c/3_qis/1_statevector.md index fface258..dbca7a73 100644 --- a/docs/api/cn/c/3_qis/1_statevector.md +++ b/docs/api/cn/c/3_qis/1_statevector.md @@ -8,6 +8,8 @@ ### statevector_new(num_qubits) +状态维度无法表示或状态内存分配失败时返回 `NULL`。支持零比特状态。 + ```c struct CStatevector *statevector_new(uintptr_t num_qubits); ``` diff --git a/docs/api/cn/c/3_qis/2_density_matrix.md b/docs/api/cn/c/3_qis/2_density_matrix.md index b0037fcf..c57a12e5 100644 --- a/docs/api/cn/c/3_qis/2_density_matrix.md +++ b/docs/api/cn/c/3_qis/2_density_matrix.md @@ -8,6 +8,8 @@ ### density_matrix_new(num_qubits) +状态维度无法表示或状态内存分配失败时返回 `NULL`。支持零比特状态。 + ```c struct CDensityMatrix *density_matrix_new(uintptr_t num_qubits); ``` @@ -20,6 +22,8 @@ struct CDensityMatrix *density_matrix_new(uintptr_t num_qubits); ### density_matrix_maximally_mixed(num_qubits) +状态维度无法表示或状态内存分配失败时返回 `NULL`。支持零比特状态。 + ```c struct CDensityMatrix *density_matrix_maximally_mixed(uintptr_t num_qubits); ``` @@ -32,6 +36,8 @@ struct CDensityMatrix *density_matrix_maximally_mixed(uintptr_t num_qubits); ### density_matrix_zeros(num_qubits) +状态维度无法表示或状态内存分配失败时返回 `NULL`。支持零比特状态。 + ```c struct CDensityMatrix *density_matrix_zeros(uintptr_t num_qubits); ``` diff --git a/docs/api/cn/c/3_qis/3_density_matrix_noise.md b/docs/api/cn/c/3_qis/3_density_matrix_noise.md index 03ddac33..843d85fc 100644 --- a/docs/api/cn/c/3_qis/3_density_matrix_noise.md +++ b/docs/api/cn/c/3_qis/3_density_matrix_noise.md @@ -8,6 +8,8 @@ ### density_matrix_noise_new(num_qubits, noise_model) +状态维度无法表示或状态内存分配失败时返回 `NULL`。支持零比特状态。 + ```c struct CDensityMatrixNoise *density_matrix_noise_new(uintptr_t num_qubits, const struct CNoiseModel *noise_model); diff --git a/docs/api/cn/python/3_qis/1_statevector.md b/docs/api/cn/python/3_qis/1_statevector.md index 15445836..1e9f91fd 100644 --- a/docs/api/cn/python/3_qis/1_statevector.md +++ b/docs/api/cn/python/3_qis/1_statevector.md @@ -20,6 +20,8 @@ from cqlib.qis.state import Statevector ### Statevector(num_qubits) +状态维度无法表示或状态内存分配失败时抛出 `ValueError`。支持零比特状态。 + 创建一个处于 |0…0⟩ 的状态。 参数： diff --git a/docs/api/cn/python/3_qis/2_density_matrix.md b/docs/api/cn/python/3_qis/2_density_matrix.md index 1b053966..8588a82c 100644 --- a/docs/api/cn/python/3_qis/2_density_matrix.md +++ b/docs/api/cn/python/3_qis/2_density_matrix.md @@ -20,6 +20,8 @@ from cqlib.qis.state import DensityMatrix ### DensityMatrix(num_qubits) +状态维度无法表示或状态内存分配失败时抛出 `ValueError`。支持零比特状态。 + 创建一个处于 |0…0⟩⟨0…0| 的纯态。 参数： diff --git a/docs/api/cn/python/3_qis/3_density_matrix_noise.md b/docs/api/cn/python/3_qis/3_density_matrix_noise.md index a2ec8681..07cd7924 100644 --- a/docs/api/cn/python/3_qis/3_density_matrix_noise.md +++ b/docs/api/cn/python/3_qis/3_density_matrix_noise.md @@ -21,6 +21,8 @@ from cqlib.qis.state import DensityMatrixNoise ### DensityMatrixNoise(num_qubits, noise_model=None) +状态维度无法表示或状态内存分配失败时抛出 `ValueError`。支持零比特状态。 + 创建一个模拟器。省略噪声模型时等价于理想模拟。 参数： diff --git a/docs/api/cn/python/3_qis/7_hamiltonian.md b/docs/api/cn/python/3_qis/7_hamiltonian.md index 8e5627ee..d5f5b525 100644 --- a/docs/api/cn/python/3_qis/7_hamiltonian.md +++ b/docs/api/cn/python/3_qis/7_hamiltonian.md @@ -113,6 +113,8 @@ class Observable(Protocol): ### 异常情况 +实数期望值和方差接口要求吸收 Pauli 相位、合并重复项后的算符为厄米算符；原始及合并后的系数必须有限，否则抛出 `ValueError`。计算保留非零小系数。 + - `ValueError`：两个算符比特数不一致（`+`、`+=`）、项的比特数与算符不一致（`from_list`、`add_term`）、系数类型不受支持、观测量与量子态比特数不一致（`expectation_*`、`variance_statevector`）、测量概率中找不到可用的测量基（`expectation_probs`）。 - `ValueError`：`to_trotter_circuit` 与 `to_evolution_circuit` 在 `steps` 为 `0`、算符为空、系数非 Hermitian、或 Pauli 字符串相位非 Hermitian 时抛出。 diff --git a/docs/api/cn/python/3_qis/9_metrics_entropy.md b/docs/api/cn/python/3_qis/9_metrics_entropy.md index 70ed16c9..aa976528 100644 --- a/docs/api/cn/python/3_qis/9_metrics_entropy.md +++ b/docs/api/cn/python/3_qis/9_metrics_entropy.md @@ -370,6 +370,8 @@ assert abs(entropy.linear_entropy(dm_pure) - 0.0) < 1e-10 ### `renyi_entropy(dm, alpha) -> float` +阶数必须为正且不能为 NaN。`alpha = 1` 返回冯·诺依曼熵，正无穷返回最小熵；接近 1 时使用稳定计算。非对角矩阵中不超过 `矩阵维度 × f64::EPSILON × 最大特征值` 的特征值视为数值零；对角矩阵保留全部正特征值。 + 计算阶数为 `alpha` 的 Rényi 熵。 参数： @@ -383,7 +385,7 @@ assert abs(entropy.linear_entropy(dm_pure) - 0.0) < 1e-10 异常情况： -- `ValueError`：`alpha` 不大于零，或特征分解失败。 +- `ValueError`：`alpha` 为 NaN 或不大于零，或特征分解失败。 示例： @@ -528,7 +530,7 @@ assert abs(entropy.entanglement_of_formation(dm_sep) - 0.0) < 1e-10 | 异常 | 触发场景 | | --- | --- | -| `ValueError` | 参与比较的两个量子态比特数不一致；`renyi_entropy` 的 `alpha` 不大于零；`entanglement_entropy_pure` 的子系统为空、包含全部比特或含重复下标；`concurrence` 与 `entanglement_of_formation` 的态不是 2 比特；特征分解失败。 | +| `ValueError` | 参与比较的两个量子态比特数不一致；`renyi_entropy` 的 `alpha` 为 NaN 或不大于零；`entanglement_entropy_pure` 的子系统为空、包含全部比特或含重复下标；`concurrence` 与 `entanglement_of_formation` 的态不是 2 比特；特征分解失败。 | | `IndexError` | `partial_transpose`、`logarithmic_negativity`、`negativity`、`entanglement_entropy_pure` 的比特下标越界。 | 注意 `cqlib.qis.entropy` 是子模块名，`cqlib.qis.metrics` 中另有一个同名函数 `entropy`，二者互不相干。 diff --git a/docs/api/cn/rust/3_qis/1_statevector.md b/docs/api/cn/rust/3_qis/1_statevector.md index 33f14470..ed357327 100644 --- a/docs/api/cn/rust/3_qis/1_statevector.md +++ b/docs/api/cn/rust/3_qis/1_statevector.md @@ -32,7 +32,10 @@ pub struct Statevector { ### 构造与振幅访问 +`new` 在维度无法表示或状态内存分配失败时会 panic。需要处理错误时使用对应的 `try_*` 构造函数。支持零比特状态。 + - `fn new(num_qubits: usize) -> Self`：构造 `|0...0>` 态，首分量置 1，其余为 0。 +- `fn try_new(num_qubits: usize) -> Result`：构造基态；维度无法表示时返回 `QisError::InvalidParameterValue`，状态内存分配失败时返回 `QisError::UnsupportedOperation`。 - `fn from_state(num_qubits: usize, initial_state: Vec) -> Result`：用给定振幅构造，要求长度等于 `2^num_qubits` 且已经归一化。 - `fn data(&self) -> &[Complex64]`：按共享切片读取全部振幅。 - `fn data_mut(&mut self) -> &mut [Complex64]`：按可变切片读取全部振幅；改写可能破坏归一化。 diff --git a/docs/api/cn/rust/3_qis/2_density_matrix.md b/docs/api/cn/rust/3_qis/2_density_matrix.md index a28a52dc..456f1902 100644 --- a/docs/api/cn/rust/3_qis/2_density_matrix.md +++ b/docs/api/cn/rust/3_qis/2_density_matrix.md @@ -32,9 +32,14 @@ pub struct DensityMatrix { ### 构造与数据访问 +`new`、`zeros` 和 `maximally_mixed` 在维度无法表示或状态内存分配失败时会 panic。需要处理错误时使用对应的 `try_*` 构造函数。支持零比特状态。 + - `fn new(num_qubits: usize) -> Self`：构造纯态 `|0...0><0...0|`。 +- `fn try_new(num_qubits: usize) -> Result`：构造基态；维度无法表示时返回 `QisError::InvalidParameterValue`，状态内存分配失败时返回 `QisError::UnsupportedOperation`。 - `fn maximally_mixed(num_qubits: usize) -> Self`：构造最大混合态 `I / 2^N`，各计算基态等概率且无相干项。 +- `fn try_maximally_mixed(num_qubits: usize) -> Result`：对应的可失败构造，错误类型与 `try_new` 相同。 - `fn zeros(num_qubits: usize) -> Self`：构造全零矩阵；它不是合法物理态（迹为 0），用于作为量子信道等操作的累加起点。 +- `fn try_zeros(num_qubits: usize) -> Result`：对应的可失败构造，错误类型与 `try_new` 相同。 - `fn from_state(num_qubits: usize, initial_state: Vec) -> Result`：由纯态振幅做外积 `ρ = |ψ><ψ|` 构造。 - `fn from_density_matrix_state(num_qubits: usize, dm_state: Vec) -> Result`：直接由展平的 `2^N × 2^N` 矩阵构造，并校验 Hermitian 性、半正定性与单位迹。 - `fn data(&self) -> &[Complex64]`：按共享切片读取展平的矩阵元素。 diff --git a/docs/api/cn/rust/3_qis/3_density_matrix_noise.md b/docs/api/cn/rust/3_qis/3_density_matrix_noise.md index 05af1f38..47cc8b4b 100644 --- a/docs/api/cn/rust/3_qis/3_density_matrix_noise.md +++ b/docs/api/cn/rust/3_qis/3_density_matrix_noise.md @@ -35,7 +35,10 @@ pub struct DensityMatrixNoise { ### 构造与线路对接 +`new` 在维度无法表示或状态内存分配失败时会 panic。需要处理错误时使用对应的 `try_*` 构造函数。支持零比特状态。 + - `fn new(num_qubits: usize, noise_model: Option) -> Self`：按比特数与噪声模型构造，初始态为 `|0...0><0...0|`。 +- `fn try_new(num_qubits: usize, noise_model: Option) -> Result`：构造基态；维度无法表示时返回 `QisError::InvalidParameterValue`，状态内存分配失败时返回 `QisError::UnsupportedOperation`。 - `fn from_circuit(circuit: &Circuit, noise_model: Option) -> Result`：执行线路并返回模拟器；线路会先被分解为基础门，噪声在每个门之后施加。 - `fn apply_circuit(&mut self, circuit: &Circuit) -> Result<(), QisError>`：把线路原地作用到当前模拟器上，噪声按已配置的模型施加。 diff --git a/docs/api/cn/rust/3_qis/7_hamiltonian.md b/docs/api/cn/rust/3_qis/7_hamiltonian.md index a8a5cf38..14757c31 100644 --- a/docs/api/cn/rust/3_qis/7_hamiltonian.md +++ b/docs/api/cn/rust/3_qis/7_hamiltonian.md @@ -1,6 +1,6 @@ # Hamiltonian 与 Observable -`cqlib_core::qis::hamiltonian` +`cqlib_core::qis::hamiltonian` `cqlib_core::qis::observable` 本页介绍由 Pauli 串与复系数构成的观测量 `Hamiltonian`，以及计算期望值的统一接口 `Observable` trait。 @@ -175,6 +175,8 @@ pub trait Observable { - `fn variance_statevector(&self, sv: &Statevector) -> Result`：计算纯态方差 $\text{Var}(O) = \langle O^2 \rangle - \langle O \rangle^2$。trait 中的默认实现返回 `QisError::UnsupportedOperation`，支持方差计算的观测量需要覆写该方法。 +实数期望值和方差接口在吸收相位、合并重复项后检查厄米性，不删除非零小系数。原始或合并后系数非有限时返回 `QisError::InvalidParameterValue`，非厄米算符返回 `QisError::NotHermitian`。 + ### 实现者 | 类型 | 说明 | diff --git a/docs/api/cn/rust/3_qis/8_evolution.md b/docs/api/cn/rust/3_qis/8_evolution.md index 73e420db..2870f773 100644 --- a/docs/api/cn/rust/3_qis/8_evolution.md +++ b/docs/api/cn/rust/3_qis/8_evolution.md @@ -193,7 +193,7 @@ assert!(circuit.pauli_evolution(&pauli, 1.0, &qubits).is_err()); - `steps` (`usize`)：Trotter 步数 $n$，必须大于 0。 - `mode` (`TrotterMode`)：分解模式。 -两者都先复制一份 Hamilton 量并化简，再检查系数的虚部：虚部绝对值超过 $10^{-10}$ 时返回 `QisError::NotHermitian`，因为只有 Hermitian 算符才能生成酉演化。 +两者都先复制一份 Hamilton 量，吸收 Pauli 相位并合并重复项，保留非零小系数，再检查系数的虚部：虚部绝对值超过 $10^{-10}$ 时返回 `QisError::NotHermitian`，因为只有 Hermitian 算符才能生成酉演化。 ### `fn to_trotter_circuit(&self, time: f64, steps: usize, mode: TrotterMode) -> Result` diff --git a/docs/api/cn/rust/3_qis/9_metrics_entropy.md b/docs/api/cn/rust/3_qis/9_metrics_entropy.md index 26598d7c..7a41b162 100644 --- a/docs/api/cn/rust/3_qis/9_metrics_entropy.md +++ b/docs/api/cn/rust/3_qis/9_metrics_entropy.md @@ -1,6 +1,6 @@ # 度量与熵 -`cqlib_core::qis::metrics` +`cqlib_core::qis::metrics` `cqlib_core::qis::entropy` 本页介绍量子态之间的度量与信息量函数。`metrics` 集合纯度、保真度、迹距离、冯·诺依曼熵与偏转置相关的量；`entropy` 集合线性熵、Rényi 熵，以及针对纠缠的度量。两模块的函数都以态矢量或密度矩阵为输入，返回 `Result`（`density_matrix_to_faer` 除外）。 @@ -143,6 +143,8 @@ assert!((trace_distance_mixed(&dm0, &dm1).unwrap() - 1.0).abs() < 1e-10); ### `renyi_entropy(dm: &DensityMatrix, alpha: f64) -> Result` +阶数必须为正且不能为 NaN。`alpha = 1` 返回冯·诺依曼熵，正无穷返回最小熵；接近 1 时使用稳定计算。非对角矩阵中不超过 `矩阵维度 × f64::EPSILON × 最大特征值` 的特征值视为数值零；对角矩阵保留全部正特征值。 + Rényi 熵 $S_\alpha(\rho) = \frac{1}{1-\alpha}\log_2\left(\sum_i \lambda_i^\alpha\right)$，$\lambda_i$ 为密度矩阵的特征值。$\alpha$ 不大于 0 时返回 `QisError::InvalidParameterValue`；$\alpha$ 与 1 的差小于 `f64::EPSILON` 时改用冯·诺依曼熵，避免除零。结果截断为非负值。 示例： @@ -283,7 +285,7 @@ assert!(partial_transpose(&dm, &[2]).is_err()); | 错误 | 触发场景 | | --- | --- | | `QisError::QubitMismatch` | 参与运算的两个态比特数不一致。 | -| `QisError::InvalidParameterValue` | `renyi_entropy()` 的 `alpha` 不大于 0。 | +| `QisError::InvalidParameterValue` | `renyi_entropy()` 的 `alpha` 为 NaN 或不大于 0。 | | `QisError::InvalidSubsystem` | `entanglement_entropy_pure()` 的子系统为空、包含全部比特或存在重复下标。 | | `QisError::IndexOutOfBounds` | `partial_transpose()` 的目标下标越界，或 `entanglement_entropy_pure()` 的子系统下标越界。 | | `QisError::UnsupportedDimension` | `concurrence()` 与 `entanglement_of_formation()` 的输入不是双比特。 | diff --git a/docs/api/en/c/3_qis/1_statevector.md b/docs/api/en/c/3_qis/1_statevector.md index 9bc3e97c..c5ee415a 100644 --- a/docs/api/en/c/3_qis/1_statevector.md +++ b/docs/api/en/c/3_qis/1_statevector.md @@ -8,6 +8,8 @@ The `CStatevector*` handle is the pure-state simulator: the quantum state is rep ### statevector_new(num_qubits) +Returns `NULL` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + ```c struct CStatevector *statevector_new(uintptr_t num_qubits); ``` diff --git a/docs/api/en/c/3_qis/2_density_matrix.md b/docs/api/en/c/3_qis/2_density_matrix.md index d66548f6..94dc372f 100644 --- a/docs/api/en/c/3_qis/2_density_matrix.md +++ b/docs/api/en/c/3_qis/2_density_matrix.md @@ -8,6 +8,8 @@ The `CDensityMatrix*` handle is the mixed-state simulator: the quantum state is ### density_matrix_new(num_qubits) +Returns `NULL` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + ```c struct CDensityMatrix *density_matrix_new(uintptr_t num_qubits); ``` @@ -20,6 +22,8 @@ Returns: a newly allocated `CDensityMatrix*`; free it with `density_matrix_free` ### density_matrix_maximally_mixed(num_qubits) +Returns `NULL` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + ```c struct CDensityMatrix *density_matrix_maximally_mixed(uintptr_t num_qubits); ``` @@ -32,6 +36,8 @@ Returns: a newly allocated `CDensityMatrix*`; free it with `density_matrix_free` ### density_matrix_zeros(num_qubits) +Returns `NULL` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + ```c struct CDensityMatrix *density_matrix_zeros(uintptr_t num_qubits); ``` diff --git a/docs/api/en/c/3_qis/3_density_matrix_noise.md b/docs/api/en/c/3_qis/3_density_matrix_noise.md index 74b85ee7..b542c22e 100644 --- a/docs/api/en/c/3_qis/3_density_matrix_noise.md +++ b/docs/api/en/c/3_qis/3_density_matrix_noise.md @@ -8,6 +8,8 @@ The `CDensityMatrixNoise*` handle is the noisy mixed-state simulator: it layers ### density_matrix_noise_new(num_qubits, noise_model) +Returns `NULL` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + ```c struct CDensityMatrixNoise *density_matrix_noise_new(uintptr_t num_qubits, const struct CNoiseModel *noise_model); diff --git a/docs/api/en/python/3_qis/1_statevector.md b/docs/api/en/python/3_qis/1_statevector.md index 3a31f359..691df57d 100644 --- a/docs/api/en/python/3_qis/1_statevector.md +++ b/docs/api/en/python/3_qis/1_statevector.md @@ -20,6 +20,8 @@ Amplitudes are stored contiguously by index in memory, and the amplitude of the ### Statevector(num_qubits) +Raises `ValueError` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + Create a state in |0…0⟩. Parameters: diff --git a/docs/api/en/python/3_qis/2_density_matrix.md b/docs/api/en/python/3_qis/2_density_matrix.md index c1d50ebe..7f9641a2 100644 --- a/docs/api/en/python/3_qis/2_density_matrix.md +++ b/docs/api/en/python/3_qis/2_density_matrix.md @@ -20,6 +20,8 @@ Data is stored contiguously in row-major order. Qubit i corresponds to bit i of ### DensityMatrix(num_qubits) +Raises `ValueError` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + Create a pure state in |0…0⟩⟨0…0|. Parameters: diff --git a/docs/api/en/python/3_qis/3_density_matrix_noise.md b/docs/api/en/python/3_qis/3_density_matrix_noise.md index a2983041..1a318e58 100644 --- a/docs/api/en/python/3_qis/3_density_matrix_noise.md +++ b/docs/api/en/python/3_qis/3_density_matrix_noise.md @@ -21,6 +21,8 @@ Data is stored contiguously in row-major order. Qubit i corresponds to bit i of ### DensityMatrixNoise(num_qubits, noise_model=None) +Raises `ValueError` if the state dimension is unrepresentable or state allocation fails. Zero-qubit states are supported. + Create a simulator. Omitting the noise model is equivalent to an ideal simulation. Parameters: diff --git a/docs/api/en/python/3_qis/7_hamiltonian.md b/docs/api/en/python/3_qis/7_hamiltonian.md index 9f27e7b2..8afa6585 100644 --- a/docs/api/en/python/3_qis/7_hamiltonian.md +++ b/docs/api/en/python/3_qis/7_hamiltonian.md @@ -113,6 +113,8 @@ Raises: ### Raises +Real expectation and variance methods require a Hermitian operator after Pauli phases and duplicate terms are merged. Original and merged coefficients must be finite; violations raise `ValueError`. Small nonzero coefficients are retained. + - `ValueError`: the two operators have inconsistent qubit counts (`+`, `+=`), the qubit count of a term is inconsistent with the operator (`from_list`, `add_term`), the coefficient type is not supported, the observable and the quantum state have inconsistent qubit counts (`expectation_*`, `variance_statevector`), or no usable measurement basis is found in the measurement probabilities (`expectation_probs`). - `ValueError`: `to_trotter_circuit` and `to_evolution_circuit` raise when `steps` is `0`, the operator is empty, a coefficient is not Hermitian, or a Pauli string phase is not Hermitian. diff --git a/docs/api/en/python/3_qis/9_metrics_entropy.md b/docs/api/en/python/3_qis/9_metrics_entropy.md index d7cf7246..56c3be53 100644 --- a/docs/api/en/python/3_qis/9_metrics_entropy.md +++ b/docs/api/en/python/3_qis/9_metrics_entropy.md @@ -370,6 +370,8 @@ assert abs(entropy.linear_entropy(dm_pure) - 0.0) < 1e-10 ### `renyi_entropy(dm, alpha) -> float` +The order must be positive and not NaN. At `alpha = 1` the result is the von Neumann entropy; positive infinity returns min-entropy. Nearby orders use stable evaluation. Non-diagonal spectra treat eigenvalues at or below `dimension × f64::EPSILON × largest eigenvalue` as numerical zeros; diagonal spectra retain all positive entries. + Compute the Rényi entropy of order `alpha`. Parameters: @@ -383,7 +385,7 @@ Returns: Raises: -- `ValueError`: `alpha` is not greater than zero, or eigendecomposition failure. +- `ValueError`: `alpha` is NaN or not greater than zero, or eigendecomposition failure. Example: diff --git a/docs/api/en/rust/3_qis/1_statevector.md b/docs/api/en/rust/3_qis/1_statevector.md index bc6a0a7a..42a89a49 100644 --- a/docs/api/en/rust/3_qis/1_statevector.md +++ b/docs/api/en/rust/3_qis/1_statevector.md @@ -32,7 +32,10 @@ Methods: ### Construction and amplitude access +`new` panic when the dimension is unrepresentable or state allocation fails. Use the corresponding `try_*` constructor to handle these errors. Zero-qubit states are supported. + - `fn new(num_qubits: usize) -> Self`: construct the `|0...0>` state with the first component set to 1 and the rest 0. +- `fn try_new(num_qubits: usize) -> Result`: construct the ground state, returning `QisError::InvalidParameterValue` for an unrepresentable dimension or `QisError::UnsupportedOperation` when state allocation fails. - `fn from_state(num_qubits: usize, initial_state: Vec) -> Result`: construct from the given amplitudes; the length must equal `2^num_qubits` and the state must already be normalized. - `fn data(&self) -> &[Complex64]`: read all amplitudes as a shared slice. - `fn data_mut(&mut self) -> &mut [Complex64]`: read all amplitudes as a mutable slice; modification may break normalization. diff --git a/docs/api/en/rust/3_qis/2_density_matrix.md b/docs/api/en/rust/3_qis/2_density_matrix.md index 8a1a4602..22bb4467 100644 --- a/docs/api/en/rust/3_qis/2_density_matrix.md +++ b/docs/api/en/rust/3_qis/2_density_matrix.md @@ -32,9 +32,14 @@ Methods: ### Construction and data access +`new`, `zeros` and `maximally_mixed` panic when the dimension is unrepresentable or state allocation fails. Use the corresponding `try_*` constructor to handle these errors. Zero-qubit states are supported. + - `fn new(num_qubits: usize) -> Self`: construct the pure state `|0...0><0...0|`. +- `fn try_new(num_qubits: usize) -> Result`: construct the ground state, returning `QisError::InvalidParameterValue` for an unrepresentable dimension or `QisError::UnsupportedOperation` when state allocation fails. - `fn maximally_mixed(num_qubits: usize) -> Self`: construct the maximally mixed state `I / 2^N`, where every computational basis state is equally probable and there are no coherence terms. +- `fn try_maximally_mixed(num_qubits: usize) -> Result`: fallible variant with the same errors as `try_new`. - `fn zeros(num_qubits: usize) -> Self`: construct an all-zero matrix; it is not a valid physical state (trace 0) and is used as the accumulation starting point for operations such as quantum channels. +- `fn try_zeros(num_qubits: usize) -> Result`: fallible variant with the same errors as `try_new`. - `fn from_state(num_qubits: usize, initial_state: Vec) -> Result`: construct by the outer product of pure state amplitudes `ρ = |ψ><ψ|`. - `fn from_density_matrix_state(num_qubits: usize, dm_state: Vec) -> Result`: construct directly from a flattened `2^N × 2^N` matrix, validating Hermiticity, positive semidefiniteness and unit trace. - `fn data(&self) -> &[Complex64]`: read the flattened matrix elements as a shared slice. diff --git a/docs/api/en/rust/3_qis/3_density_matrix_noise.md b/docs/api/en/rust/3_qis/3_density_matrix_noise.md index aade08d4..fbf788ba 100644 --- a/docs/api/en/rust/3_qis/3_density_matrix_noise.md +++ b/docs/api/en/rust/3_qis/3_density_matrix_noise.md @@ -35,7 +35,10 @@ Methods: ### Construction and circuit integration +`new` panic when the dimension is unrepresentable or state allocation fails. Use the corresponding `try_*` constructor to handle these errors. Zero-qubit states are supported. + - `fn new(num_qubits: usize, noise_model: Option) -> Self`: construct from a qubit count and a noise model; the initial state is `|0...0><0...0|`. +- `fn try_new(num_qubits: usize, noise_model: Option) -> Result`: construct the ground state, returning `QisError::InvalidParameterValue` for an unrepresentable dimension or `QisError::UnsupportedOperation` when state allocation fails. - `fn from_circuit(circuit: &Circuit, noise_model: Option) -> Result`: execute a circuit and return the simulator; the circuit is first decomposed into basis gates, and noise is applied after every gate. - `fn apply_circuit(&mut self, circuit: &Circuit) -> Result<(), QisError>`: apply a circuit to the current simulator in place; noise is applied according to the configured model. diff --git a/docs/api/en/rust/3_qis/7_hamiltonian.md b/docs/api/en/rust/3_qis/7_hamiltonian.md index a5f1ae93..44318a14 100644 --- a/docs/api/en/rust/3_qis/7_hamiltonian.md +++ b/docs/api/en/rust/3_qis/7_hamiltonian.md @@ -1,6 +1,6 @@ # Hamiltonian and Observable -`cqlib_core::qis::hamiltonian` +`cqlib_core::qis::hamiltonian` `cqlib_core::qis::observable` This page covers the observable `Hamiltonian`, formed by Pauli strings and complex coefficients, and the unified interface `Observable` trait for computing expectation values. @@ -175,6 +175,8 @@ pub trait Observable { - `fn variance_statevector(&self, sv: &Statevector) -> Result`: compute the pure-state variance $\text{Var}(O) = \langle O^2 \rangle - \langle O \rangle^2$. The default implementation in the trait returns `QisError::UnsupportedOperation`; an observable that supports variance computation must override this method. +Real expectation and variance methods validate Hermiticity after merging phases and repeated terms, without pruning small nonzero coefficients. Non-finite original or merged coefficients return `QisError::InvalidParameterValue`; non-Hermitian operators return `QisError::NotHermitian`. + ### Implementors | Type | Description | diff --git a/docs/api/en/rust/3_qis/8_evolution.md b/docs/api/en/rust/3_qis/8_evolution.md index b453802b..2b018d86 100644 --- a/docs/api/en/rust/3_qis/8_evolution.md +++ b/docs/api/en/rust/3_qis/8_evolution.md @@ -193,7 +193,7 @@ Both entry points accept the same parameters: - `steps` (`usize`): the number of Trotter steps $n$, which must be greater than 0. - `mode` (`TrotterMode`): the decomposition mode. -Both first copy the Hamiltonian and simplify it, then check the imaginary part of the coefficients: when the absolute value of the imaginary part exceeds $10^{-10}$, `QisError::NotHermitian` is returned, because only Hermitian operators can generate unitary evolution. +Both first copy the Hamiltonian, absorb Pauli phases and merge repeated terms without pruning small nonzero coefficients, then check the imaginary part of the coefficients: when the absolute value of the imaginary part exceeds $10^{-10}$, `QisError::NotHermitian` is returned, because only Hermitian operators can generate unitary evolution. ### `fn to_trotter_circuit(&self, time: f64, steps: usize, mode: TrotterMode) -> Result` diff --git a/docs/api/en/rust/3_qis/9_metrics_entropy.md b/docs/api/en/rust/3_qis/9_metrics_entropy.md index cda03779..9d59322a 100644 --- a/docs/api/en/rust/3_qis/9_metrics_entropy.md +++ b/docs/api/en/rust/3_qis/9_metrics_entropy.md @@ -1,6 +1,6 @@ # Metrics and Entropy -`cqlib_core::qis::metrics` +`cqlib_core::qis::metrics` `cqlib_core::qis::entropy` This page covers the metrics and information measures between quantum states. `metrics` collects purity, fidelity, trace distance, von Neumann entropy and partial transpose related quantities; `entropy` collects linear entropy, Rényi entropy, and the measures aimed at entanglement. The functions of both modules take a statevector or a density matrix as input and return `Result` (except `density_matrix_to_faer`). @@ -143,6 +143,8 @@ The linear entropy $S_L(\rho) = 1 - \text{Tr}(\rho^2)$. It is `0.0` for a pure s ### `renyi_entropy(dm: &DensityMatrix, alpha: f64) -> Result` +The order must be positive and not NaN. At `alpha = 1` the result is the von Neumann entropy; positive infinity returns min-entropy. Nearby orders use stable evaluation. Non-diagonal spectra treat eigenvalues at or below `dimension × f64::EPSILON × largest eigenvalue` as numerical zeros; diagonal spectra retain all positive entries. + The Rényi entropy $S_\alpha(\rho) = \frac{1}{1-\alpha}\log_2\left(\sum_i \lambda_i^\alpha\right)$, where $\lambda_i$ are the eigenvalues of the density matrix. When $\alpha$ is not greater than 0, `QisError::InvalidParameterValue` is returned; when the difference between $\alpha$ and 1 is smaller than `f64::EPSILON`, the von Neumann entropy is used instead, avoiding division by zero. The result is truncated to a non-negative value. Example: @@ -283,7 +285,7 @@ Convert a density matrix into the matrix type of the eigendecomposition library, | Error | When it occurs | | --- | --- | | `QisError::QubitMismatch` | The two states involved in the operation have inconsistent qubit counts. | -| `QisError::InvalidParameterValue` | `alpha` of `renyi_entropy()` is not greater than 0. | +| `QisError::InvalidParameterValue` | `alpha` of `renyi_entropy()` is NaN or not greater than 0. | | `QisError::InvalidSubsystem` | The subsystem of `entanglement_entropy_pure()` is empty, contains all qubits, or has duplicate indices. | | `QisError::IndexOutOfBounds` | A target index of `partial_transpose()` is out of range, or a subsystem index of `entanglement_entropy_pure()` is out of range. | | `QisError::UnsupportedDimension` | The input of `concurrence()` and `entanglement_of_formation()` is not two-qubit. | From 3ce2049292f40c29c150c9acad7e366fec4e466b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=81=A5=E5=81=A5?= Date: Thu, 8 Oct 2026 11:50:46 +0800 Subject: [PATCH 2/4] feat(binding-python): support Hamiltonian copy protocols Implement copy.copy and copy.deepcopy using the existing owned-value clone. Include type declarations, English and Chinese documentation, and tests for empty operators, independent terms, and alias preservation in nested deepcopy. Validation: all 58 Hamiltonian Python tests passed against an extension built from the staged source snapshot. --- .../binding-python/cqlib/qis/hamiltonian.pyi | 8 +++++ crates/binding-python/src/qis/hamiltonian.rs | 10 +++++++ .../tests/qis/test_hamiltonian.py | 29 +++++++++++++++++++ docs/api/cn/python/3_qis/7_hamiltonian.md | 1 + docs/api/en/python/3_qis/7_hamiltonian.md | 1 + 5 files changed, 49 insertions(+) diff --git a/crates/binding-python/cqlib/qis/hamiltonian.pyi b/crates/binding-python/cqlib/qis/hamiltonian.pyi index 7780e5b5..3a7ee341 100644 --- a/crates/binding-python/cqlib/qis/hamiltonian.pyi +++ b/crates/binding-python/cqlib/qis/hamiltonian.pyi @@ -357,3 +357,11 @@ class Hamiltonian: A new Hamiltonian instance with the same terms. """ ... + + def __copy__(self) -> "Hamiltonian": + """Return an independent copy via copy.copy().""" + ... + + def __deepcopy__(self, memo: dict) -> "Hamiltonian": + """Return an independent copy via copy.deepcopy().""" + ... diff --git a/crates/binding-python/src/qis/hamiltonian.rs b/crates/binding-python/src/qis/hamiltonian.rs index 0addb18b..b81ca7f6 100644 --- a/crates/binding-python/src/qis/hamiltonian.rs +++ b/crates/binding-python/src/qis/hamiltonian.rs @@ -465,6 +465,16 @@ impl PyHamiltonian { } } + /// Implements Python's shallow-copy protocol with an independent value. + fn __copy__(&self) -> Self { + self.copy() + } + + /// Implements Python's deep-copy protocol; all Rust fields are owned values. + fn __deepcopy__(&self, _memo: &Bound<'_, PyAny>) -> Self { + self.copy() + } + /// Computes the expectation value for a statevector. /// /// Raises: diff --git a/crates/binding-python/tests/qis/test_hamiltonian.py b/crates/binding-python/tests/qis/test_hamiltonian.py index e9065722..969870fb 100644 --- a/crates/binding-python/tests/qis/test_hamiltonian.py +++ b/crates/binding-python/tests/qis/test_hamiltonian.py @@ -721,3 +721,32 @@ def test_hamiltonian_from_list_rejects_mixed_qubit_counts(): (PauliString.from_str("ZZ"), 1.0), ] ) + + +@pytest.mark.parametrize("deep", [False, True]) +@pytest.mark.parametrize("empty", [False, True]) +def test_python_copy_protocols(deep, empty): + from copy import copy, deepcopy + + original = Hamiltonian(2) + if not empty: + original.add_term(PauliString.from_str("ZZ"), 0.5 + 0.25j) + duplicate = (deepcopy if deep else copy)(original) + assert duplicate is not original + assert duplicate == original + assert duplicate.num_qubits == 2 + duplicate.add_term(PauliString.from_str("XX"), 0.3) + assert duplicate.num_terms == original.num_terms + 1 + + +def test_nested_deepcopy_preserves_aliases_without_sharing_source(): + from copy import deepcopy + + original = Hamiltonian(1) + original.add_term(PauliString.from_str("Z"), 1.0) + snapshot = deepcopy({"first": original, "nested": [original]}) + assert snapshot["first"] is snapshot["nested"][0] + assert snapshot["first"] is not original + snapshot["first"].scale(2.0) + assert original.terms[0][1] == 1.0 + assert snapshot["first"].terms[0][1] == 2.0 diff --git a/docs/api/cn/python/3_qis/7_hamiltonian.md b/docs/api/cn/python/3_qis/7_hamiltonian.md index d5f5b525..6cb11f52 100644 --- a/docs/api/cn/python/3_qis/7_hamiltonian.md +++ b/docs/api/cn/python/3_qis/7_hamiltonian.md @@ -99,6 +99,7 @@ class Observable(Protocol): - `to_trotter_circuit(time, steps, mode) -> Circuit`：生成 Trotter-Suzuki 分解给出的时间演化线路，用一串 Pauli 旋转逼近 $U(t) = e^{-iHt}$。 - `to_evolution_circuit(time, steps, mode) -> Circuit`：生成时间演化线路。所有项对易时走精确的单遍分解；存在非对易项时回退到指定的 Trotter 模式与步数。 - `copy() -> Hamiltonian`：返回副本。 +- 标准库 `copy.copy(h)` 和 `copy.deepcopy(h)` 同样返回项数据独立的副本。 - `expectation_statevector(sv) -> float`：由态矢量计算期望值。 - `expectation_density_matrix(dm) -> float`：由密度矩阵计算期望值。 - `expectation_probs(measurements) -> float`：由测量概率计算期望值。 diff --git a/docs/api/en/python/3_qis/7_hamiltonian.md b/docs/api/en/python/3_qis/7_hamiltonian.md index 8afa6585..98adcc34 100644 --- a/docs/api/en/python/3_qis/7_hamiltonian.md +++ b/docs/api/en/python/3_qis/7_hamiltonian.md @@ -99,6 +99,7 @@ Raises: - `to_trotter_circuit(time, steps, mode) -> Circuit`: generate the time evolution circuit given by the Trotter-Suzuki decomposition, approximating $U(t) = e^{-iHt}$ by a sequence of Pauli rotations. - `to_evolution_circuit(time, steps, mode) -> Circuit`: generate the time evolution circuit. When all terms commute, an exact single-pass decomposition is used; when non-commuting terms exist, it falls back to the specified Trotter mode and number of steps. - `copy() -> Hamiltonian`: return a copy. +- Standard-library `copy.copy(h)` and `copy.deepcopy(h)` also return independent copies of the terms. - `expectation_statevector(sv) -> float`: compute the expectation value from a statevector. - `expectation_density_matrix(dm) -> float`: compute the expectation value from a density matrix. - `expectation_probs(measurements) -> float`: compute the expectation value from measurement probabilities. From 694607aaa4bb68dc3c70cc38e190566352c973bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=81=A5=E5=81=A5?= Date: Thu, 8 Oct 2026 11:53:45 +0800 Subject: [PATCH 3/4] docs: clarify CQLib conventions and algorithm references Describe endian ordering, release workflow organization, gate-basis tests, and phase conventions directly in terms of CQLib behavior. Cite the original SABRE paper and retain the existing commutative-cancellation implementation reference. Validation: all 99 compiler rewrite tests passed. Runtime behavior is unchanged. --- crates/cqlib-core/src/circuit/circuit_to_matrix.rs | 4 ++-- .../src/compile/knowledge/builtin_rules_contract_test.rs | 4 ++-- .../src/compile/transform/commutative_cancellation/sets.rs | 6 ++++-- crates/cqlib-core/src/compile/transform/layout/sabre.rs | 5 ++++- .../src/compile/transform/rewrite/rewrite_test.rs | 6 +++--- docs/releasing.md | 5 ++--- 6 files changed, 17 insertions(+), 13 deletions(-) diff --git a/crates/cqlib-core/src/circuit/circuit_to_matrix.rs b/crates/cqlib-core/src/circuit/circuit_to_matrix.rs index ed6a0694..9dc5a5d4 100644 --- a/crates/cqlib-core/src/circuit/circuit_to_matrix.rs +++ b/crates/cqlib-core/src/circuit/circuit_to_matrix.rs @@ -26,7 +26,7 @@ //! //! ## Memory Layout & Convention //! -//! - **State Vector Ordering**: Little-Endian (similar to Qiskit). Qubit 0 corresponds to the Least Significant Bit (LSB). +//! - **State Vector Ordering**: Little-Endian. Qubit 0 corresponds to the Least Significant Bit (LSB). //! State $|q_{n-1} \dots q_1 q_0\rangle$. //! - **Parallelization**: Large matrix multiplications (large state spaces) are automatically parallelized using `rayon`. @@ -123,7 +123,7 @@ pub fn circuit_to_matrix( // Map logical qubit index to physical bit position (0 to N-1). // Little-Endian mapping: first qubit in target_order corresponds to LSB (bit 0) // q0 -> bit 0, q1 -> bit 1, ... - // This matches Qiskit's convention. + // Qubit 0 occupies the least significant bit. let qubit_bit_map: HashMap = target_order .iter() .enumerate() diff --git a/crates/cqlib-core/src/compile/knowledge/builtin_rules_contract_test.rs b/crates/cqlib-core/src/compile/knowledge/builtin_rules_contract_test.rs index 619f7a78..e02a04f6 100644 --- a/crates/cqlib-core/src/compile/knowledge/builtin_rules_contract_test.rs +++ b/crates/cqlib-core/src/compile/knowledge/builtin_rules_contract_test.rs @@ -903,8 +903,8 @@ fn decompose_swap_x2p_cz_triple_product_differs_from_swap_by_exactly_pi() { "bare (X2P⊗X2P)·CZ triple product must differ from SWAP by a global phase" ); - // The qiskit SX-convention tail -π/2 is wrong for the phase-free X2P - // convention used here: it yields exp(iπ/2)·SWAP = i·SWAP. + // With X2P = Rx(π/2), a phase tail of -π/2 yields + // exp(iπ/2)·SWAP = i·SWAP instead of exactly SWAP. let mut wrong_phase = repeat_three_times(|c| { c.x2p(q0).unwrap(); c.x2p(q1).unwrap(); diff --git a/crates/cqlib-core/src/compile/transform/commutative_cancellation/sets.rs b/crates/cqlib-core/src/compile/transform/commutative_cancellation/sets.rs index ee090d27..f576d2cd 100644 --- a/crates/cqlib-core/src/compile/transform/commutative_cancellation/sets.rs +++ b/crates/cqlib-core/src/compile/transform/commutative_cancellation/sets.rs @@ -12,8 +12,7 @@ //! Commutation-set grouping and cancellation-key analysis. //! -//! The analysis mirrors qiskit's `CommutativeCancellation` structure: -//! every wire's operation sequence is greedily partitioned into contiguous +//! Each wire's operation sequence is greedily partitioned into contiguous //! sets of pairwise exactly-commuting operations, and self-inverse gates are //! bucketed by a key of `(gate, canonical qargs, per-wire set indices)`. //! Two operations share a key only when everything between them on each of @@ -23,6 +22,9 @@ //! //! Sets are tracked per wire independently; a two-qubit gate therefore //! records one set index per qubit and the key keeps both as a tuple. +//! +//! Implementation reference for per-wire grouping and cancellation keys: +//! Qiskit's `CommutativeCancellation`. use crate::circuit::{Instruction, Operation, Parameter, Qubit, StandardGate}; use crate::compile::commutation::{Commutation, CommutationChecker}; diff --git a/crates/cqlib-core/src/compile/transform/layout/sabre.rs b/crates/cqlib-core/src/compile/transform/layout/sabre.rs index aa7790e1..d5da8aca 100644 --- a/crates/cqlib-core/src/compile/transform/layout/sabre.rs +++ b/crates/cqlib-core/src/compile/transform/layout/sabre.rs @@ -12,6 +12,9 @@ //! SABRE initial-layout adapter. //! +//! Algorithm reference: Gushu Li, Yufei Ding, and Yuan Xie, +//! [Tackling the Qubit Mapping Problem for NISQ-Era Quantum Devices](https://arxiv.org/abs/1809.02573). +//! //! This module owns layout-candidate generation, objective scoring, and //! [`LayoutResult`] construction. The SABRE core remains in //! [`crate::compile::sabre`] and is used here only through crate-internal @@ -418,7 +421,7 @@ pub(crate) struct PreparedSabreParetoSearch { } impl PreparedSabreParetoSearch { - /// Builds the Qiskit-style four-refinement search profile. A supplied + /// Builds the multi-start four-refinement search profile. A supplied /// layout is normalized and kept fixed; generated layouts deliberately do /// not use the dense-interaction refinement shortcut employed by candidate /// zero. diff --git a/crates/cqlib-core/src/compile/transform/rewrite/rewrite_test.rs b/crates/cqlib-core/src/compile/transform/rewrite/rewrite_test.rs index a4f9a915..c9e5b9e9 100644 --- a/crates/cqlib-core/src/compile/transform/rewrite/rewrite_test.rs +++ b/crates/cqlib-core/src/compile/transform/rewrite/rewrite_test.rs @@ -533,7 +533,7 @@ fn lowering_reaches_target_basis_through_multiple_steps() { } #[test] -fn lowers_ccx_directly_to_multiple_qiskit_basis_sets() { +fn lowers_ccx_directly_to_multiple_native_basis_sets() { let q0 = Qubit::new(0); let q1 = Qubit::new(1); let q2 = Qubit::new(2); @@ -582,7 +582,7 @@ fn lowers_ccx_directly_to_multiple_qiskit_basis_sets() { } #[test] -fn qiskit_rzz_ccx_rule_uses_five_entanglers() { +fn rzz_ccx_rule_uses_five_entanglers() { let q0 = Qubit::new(0); let q1 = Qubit::new(1); let q2 = Qubit::new(2); @@ -612,7 +612,7 @@ fn qiskit_rzz_ccx_rule_uses_five_entanglers() { } #[test] -fn qiskit_cz_ccx_rule_uses_native_x2p_x_rz_cz_template() { +fn cz_ccx_rule_uses_native_x2p_x_rz_cz_template() { let q0 = Qubit::new(0); let q1 = Qubit::new(1); let q2 = Qubit::new(2); diff --git a/docs/releasing.md b/docs/releasing.md index f7b9b3d7..afe01cbd 100644 --- a/docs/releasing.md +++ b/docs/releasing.md @@ -2,9 +2,8 @@ ## File Organization -The structure follows Qiskit's separation of entry-point workflows, reusable -workflows, and tools scripts, without adopting its automated publishing steps -or modifying project source code, existing tests, or Cargo configuration. +CQLib separates entry-point workflows, reusable workflows, and tool scripts. +Publishing is triggered manually through release tags. | File | Responsibility | | --- | --- | From f08479f037bc8ddaf570b81c79622ea2326eb0a2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=81=A5=E5=81=A5?= Date: Thu, 8 Oct 2026 13:41:57 +0800 Subject: [PATCH 4/4] test(qis): make formation entropy regressions portable Use a square-root machine-epsilon tolerance for rotated pure-state formation entropy, while checking finite results and physical bounds. Null-space eigenvalue roundoff is amplified by square roots and varies between eigensolver backends. Extract the existing concurrence-to-entropy calculation into a private helper and directly test clamping above one, below zero, and an interior value. Preserve the public API and numerical algorithm. Validation on macOS: 36 entropy unit tests, 16 core regression tests, and 1205 Python QIS tests passed against a freshly built extension. Formatting, Ruff, and Clippy checks passed. Windows and Linux CI verification remains pending. --- .../tests/qis/test_p1_regressions.py | 10 ++++++---- crates/cqlib-core/src/qis/entropy.rs | 10 +++++++--- crates/cqlib-core/src/qis/entropy_test.rs | 15 +++++++++++++++ crates/cqlib-core/tests/p1_regressions.rs | 15 +++++++++++++-- 4 files changed, 41 insertions(+), 9 deletions(-) diff --git a/crates/binding-python/tests/qis/test_p1_regressions.py b/crates/binding-python/tests/qis/test_p1_regressions.py index 9ad84773..f8e1f2dd 100644 --- a/crates/binding-python/tests/qis/test_p1_regressions.py +++ b/crates/binding-python/tests/qis/test_p1_regressions.py @@ -64,7 +64,7 @@ def test_gate_noise_uses_sparse_labels(): np.testing.assert_allclose(sim.probabilities(), [1.0, 0.0], atol=1e-12) -def test_formation_entropy_roundoff(): +def test_formation_entropy_of_rotated_maximally_entangled_state(): state = np.array( [ -0.5917562523120286 - 0.14345273890626087j, @@ -73,9 +73,11 @@ def test_formation_entropy_roundoff(): -0.5571108451004554 + 0.2457267928442915j, ] ) - assert entropy.entanglement_of_formation( - DensityMatrix.from_state(2, state) - ) == pytest.approx(1.0, abs=1e-12) + result = entropy.entanglement_of_formation(DensityMatrix.from_state(2, state)) + assert np.isfinite(result) and 0 <= result <= 1 + # Square roots amplify null-space eigenvalue noise to O(sqrt(machine epsilon)). + tolerance = 8 * np.sqrt(np.finfo(float).eps) + assert result == pytest.approx(1.0, rel=0, abs=tolerance) def test_renyi_extreme_orders(): diff --git a/crates/cqlib-core/src/qis/entropy.rs b/crates/cqlib-core/src/qis/entropy.rs index 162c7cc0..ec0791e1 100644 --- a/crates/cqlib-core/src/qis/entropy.rs +++ b/crates/cqlib-core/src/qis/entropy.rs @@ -486,10 +486,14 @@ pub fn concurrence(dm: &DensityMatrix) -> Result { /// This measure represents the minimum number of Bell pairs required to asymptotically /// create the state ρ using LOCC (Local Operations and Classical Communication). pub fn entanglement_of_formation(dm: &DensityMatrix) -> Result { - let c = concurrence(dm)?.clamp(0.0, 1.0); + Ok(entanglement_of_formation_from_concurrence(concurrence(dm)?)) +} + +fn entanglement_of_formation_from_concurrence(c: f64) -> f64 { + let c = c.clamp(0.0, 1.0); if c < 1e-15 { - return Ok(0.0); + return 0.0; } // Compute x = (1 + sqrt(1 - C^2)) / 2 @@ -506,7 +510,7 @@ pub fn entanglement_of_formation(dm: &DensityMatrix) -> Result { binary_entropy -= (1.0 - x) * (1.0 - x).log2(); } - Ok(binary_entropy) + binary_entropy } /// Helper function to compute the matrix square root of a positive semi-definite matrix. diff --git a/crates/cqlib-core/src/qis/entropy_test.rs b/crates/cqlib-core/src/qis/entropy_test.rs index e0cf1ef4..1c8e6f24 100644 --- a/crates/cqlib-core/src/qis/entropy_test.rs +++ b/crates/cqlib-core/src/qis/entropy_test.rs @@ -666,6 +666,21 @@ fn test_eof_bell_states() { ); } +#[test] +fn test_eof_clamps_concurrence_roundoff_before_square_root() { + for c in [1.0, 1.0 + f64::EPSILON, 1.0 + 8.0 * f64::EPSILON] { + let eof = super::entanglement_of_formation_from_concurrence(c); + assert_eq!(eof, 1.0, "concurrence={c}, entropy={eof}"); + } + for c in [0.0, -f64::EPSILON] { + assert_eq!(super::entanglement_of_formation_from_concurrence(c), 0.0); + } + // An interior value still uses the binary entropy formula: C=0.6 gives x=0.9. + let expected = -0.9_f64 * 0.9_f64.log2() - 0.1_f64 * 0.1_f64.log2(); + let actual = super::entanglement_of_formation_from_concurrence(0.6); + assert!((actual - expected).abs() < 1e-14); +} + #[test] fn test_eof_separable_states() { // Separable states have EOF = 0 diff --git a/crates/cqlib-core/tests/p1_regressions.rs b/crates/cqlib-core/tests/p1_regressions.rs index 045c64ee..252fc083 100644 --- a/crates/cqlib-core/tests/p1_regressions.rs +++ b/crates/cqlib-core/tests/p1_regressions.rs @@ -99,7 +99,7 @@ fn gate_noise_uses_original_labels_and_operand_order() { } #[test] -fn formation_entropy_at_concurrence_roundoff_boundary() { +fn formation_entropy_of_rotated_maximally_entangled_state() { let state = vec![ Complex64::new(-0.5917562523120286, -0.14345273890626087), Complex64::new(0.12739228627789673, 0.3361800930257551), @@ -107,7 +107,18 @@ fn formation_entropy_at_concurrence_roundoff_boundary() { Complex64::new(-0.5571108451004554, 0.2457267928442915), ]; let dm = DensityMatrix::from_state(2, state).unwrap(); - assert!((entanglement_of_formation(&dm).unwrap() - 1.0).abs() < 1e-12); + let eof = entanglement_of_formation(&dm).unwrap(); + assert!( + eof.is_finite() && (0.0..=1.0).contains(&eof), + "entropy={eof}" + ); + // Square roots amplify O(epsilon) noise in zero eigenvalues to O(sqrt(epsilon)). + // Eigensolver backends differ in the sign of this noise on rank-one states. + let tolerance = 8.0 * f64::EPSILON.sqrt(); + assert!( + (eof - 1.0).abs() < tolerance, + "entropy={eof}, tolerance={tolerance}" + ); } #[test]