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..3a7ee341 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. """ ... @@ -343,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/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..b81ca7f6 100644 --- a/crates/binding-python/src/qis/hamiltonian.rs +++ b/crates/binding-python/src/qis/hamiltonian.rs @@ -465,7 +465,21 @@ 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: + /// 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 +487,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 +498,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 +532,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_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/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..f8e1f2dd --- /dev/null +++ b/crates/binding-python/tests/qis/test_p1_regressions.py @@ -0,0 +1,126 @@ +# 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_of_rotated_maximally_entangled_state(): + state = np.array( + [ + -0.5917562523120286 - 0.14345273890626087j, + 0.12739228627789673 + 0.3361800930257551j, + -0.06600502193827436 + 0.35339664207352467j, + -0.5571108451004554 + 0.2457267928442915j, + ] + ) + 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(): + 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/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/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..ec0791e1 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,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)?; + 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 @@ -462,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/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..252fc083 --- /dev/null +++ b/crates/cqlib-core/tests/p1_regressions.rs @@ -0,0 +1,243 @@ +// 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_of_rotated_maximally_entangled_state() { + 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(); + 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] +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..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`：由测量概率计算期望值。 @@ -113,6 +114,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..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. @@ -113,6 +114,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. | 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 | | --- | --- |