Fix QIS numerical edge cases and support Hamiltonian copy protocols

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").

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.

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.

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);
+    }
+}

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);
+}

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.
     """
     ...
 

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()."""
+        ...

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.
         """
         ...
 

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.
         """
         ...
 

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.
         """
         ...
 

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.
         """
         ...