diff --git a/crates/binding-python/cqlib/circuit/ansatz/hamiltonian_evolution.pyi b/crates/binding-python/cqlib/circuit/ansatz/hamiltonian_evolution.pyi index c68ff21c..a64c9cc4 100644 --- a/crates/binding-python/cqlib/circuit/ansatz/hamiltonian_evolution.pyi +++ b/crates/binding-python/cqlib/circuit/ansatz/hamiltonian_evolution.pyi @@ -331,9 +331,10 @@ class PauliEvolutionAnsatz: def __init__(self, hamiltonian: Hamiltonian) -> None: """Creates a new ``PauliEvolutionAnsatz``. - The Hamiltonian is automatically **simplified** before processing: - Pauli phases are absorbed into coefficients, duplicate terms are merged, - and near-zero terms are removed. The default strategy is + The Hamiltonian is normalized before processing: Pauli phases are absorbed, + duplicate terms are merged, and imaginary residuals of magnitude <= 1e-10 + are projected to zero. Larger imaginary parts are rejected; every nonzero + real coefficient is retained. The default strategy is :meth:`EvolutionStrategy.auto` with ``steps=1``. Args: diff --git a/crates/binding-python/cqlib/qis/__init__.py b/crates/binding-python/cqlib/qis/__init__.py index 37eff524..bcadae72 100644 --- a/crates/binding-python/cqlib/qis/__init__.py +++ b/crates/binding-python/cqlib/qis/__init__.py @@ -65,6 +65,15 @@ # Hamiltonian module from .hamiltonian import Hamiltonian as Hamiltonian +from . import fermion as fermion +from .fermion import ( + FermionicLadder as FermionicLadder, + FermionicOp as FermionicOp, + FermionToQubitMapper as FermionToQubitMapper, + OccupationStateMapper as OccupationStateMapper, + JordanWignerMapper as JordanWignerMapper, +) + # Evolution module from .evolution import TrotterMode as TrotterMode @@ -103,6 +112,12 @@ def num_qubits(self) -> int: ... "Pauli", "PauliString", "Hamiltonian", + "fermion", + "FermionicLadder", + "FermionicOp", + "FermionToQubitMapper", + "OccupationStateMapper", + "JordanWignerMapper", "TrotterMode", "state", "Statevector", diff --git a/crates/binding-python/cqlib/qis/__init__.pyi b/crates/binding-python/cqlib/qis/__init__.pyi index fb4dbe49..e3df9927 100644 --- a/crates/binding-python/cqlib/qis/__init__.pyi +++ b/crates/binding-python/cqlib/qis/__init__.pyi @@ -33,6 +33,15 @@ from .pauli import PauliString as PauliString # Hamiltonian module from .hamiltonian import Hamiltonian as Hamiltonian +from . import fermion as fermion +from .fermion import ( + FermionicLadder as FermionicLadder, + FermionicOp as FermionicOp, + FermionToQubitMapper as FermionToQubitMapper, + OccupationStateMapper as OccupationStateMapper, + JordanWignerMapper as JordanWignerMapper, +) + # Evolution module from .evolution import TrotterMode as TrotterMode @@ -83,6 +92,12 @@ __all__ = [ "Pauli", "PauliString", "Hamiltonian", + "fermion", + "FermionicLadder", + "FermionicOp", + "FermionToQubitMapper", + "OccupationStateMapper", + "JordanWignerMapper", "TrotterMode", "state", "Statevector", diff --git a/crates/binding-python/cqlib/qis/fermion.py b/crates/binding-python/cqlib/qis/fermion.py new file mode 100644 index 00000000..ffabc726 --- /dev/null +++ b/crates/binding-python/cqlib/qis/fermion.py @@ -0,0 +1,65 @@ +# 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. + +"""Structured fermionic algebra and fermion-to-qubit mappings. + +Modes are numbered from zero. A product A * B applies B first. Jordan-Wigner +maps mode p to qubit p, with a Z chain on modes below p; qubit 0 is the least +significant bit. Core imposes no chemical or spin-orbital ordering. + +The excitation generator is G = i(T - T†), where targets supply creators in +input order and sources supply annihilators in reverse input order. Evolving +the mapped generator implements exp(-i theta G) = exp(theta (T - T†)). +""" + +from __future__ import annotations + +from typing import Protocol, Sequence, runtime_checkable + +from .._native import qis as _qis_module +from .hamiltonian import Hamiltonian + +FermionicLadder = _qis_module.fermion.FermionicLadder +FermionicOp = _qis_module.fermion.FermionicOp +JordanWignerMapper = _qis_module.fermion.JordanWignerMapper + + +@runtime_checkable +class FermionToQubitMapper(Protocol): + """Structural contract for operator mappings, including non-Hermitian inputs. + + Implementations need not inherit this protocol. Runtime checks establish + member presence only; callers must still honor signatures and semantics. + """ + + def map_operator(self, op: FermionicOp) -> Hamiltonian: ... + def num_qubits(self, num_modes: int) -> int: ... + + +@runtime_checkable +class OccupationStateMapper(Protocol): + """Separate contract for computational-basis occupation encodings. + + occupation[p] denotes mode p and returned bits[p] addresses qubit p. + Not every operator mapping has a computational-basis occupation encoding. + """ + + def map_occupation(self, occupation: Sequence[int]) -> list[bool]: ... + + +__all__ = [ + "FermionicLadder", + "FermionicOp", + "FermionToQubitMapper", + "OccupationStateMapper", + "JordanWignerMapper", +] diff --git a/crates/binding-python/cqlib/qis/fermion.pyi b/crates/binding-python/cqlib/qis/fermion.pyi new file mode 100644 index 00000000..422893eb --- /dev/null +++ b/crates/binding-python/cqlib/qis/fermion.pyi @@ -0,0 +1,149 @@ +# 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. + +"""Typed interfaces for structured fermionic algebra and sparse mappings.""" + +from typing import ( + Iterable, + Literal, + Protocol, + Sequence, + TypeAlias, + final, + runtime_checkable, +) +from .hamiltonian import Hamiltonian + +_Coefficient: TypeAlias = complex | float | int + +@final +class FermionicLadder: + """Immutable ladder, constructed with create or annihilate, not directly.""" + @staticmethod + def create(mode: int) -> FermionicLadder: ... + @staticmethod + def annihilate(mode: int) -> FermionicLadder: ... + @property + def mode(self) -> int: ... + @property + def kind(self) -> Literal["create", "annihilate"]: ... + def adjoint(self) -> FermionicLadder: ... + def __eq__(self, other: object, /) -> bool: ... + def __hash__(self) -> int: ... + def __repr__(self) -> str: ... + +@final +class FermionicOp: + """Mutable complex sum of ordered products. B acts first in A * B. + + Identical ordered products merge exactly; no implicit normal ordering or + numerical truncation occurs. Mode counts are retained for zero and identity. + """ + + __hash__: None # type: ignore[assignment] + def __new__(cls, num_modes: int) -> FermionicOp: ... + @staticmethod + def identity(num_modes: int) -> FermionicOp: ... + @staticmethod + def from_term( + num_modes: int, + term: Iterable[FermionicLadder], + coefficient: _Coefficient = 1.0, + ) -> FermionicOp: ... + @staticmethod + def from_terms( + num_modes: int, + terms: Iterable[tuple[Iterable[FermionicLadder], _Coefficient]], + ) -> FermionicOp: ... + @staticmethod + def creation(num_modes: int, mode: int) -> FermionicOp: ... + @staticmethod + def annihilation(num_modes: int, mode: int) -> FermionicOp: ... + @staticmethod + def number(num_modes: int, mode: int) -> FermionicOp: ... + @staticmethod + def total_number(num_modes: int) -> FermionicOp: ... + @staticmethod + def excitation_generator( + num_modes: int, sources: Sequence[int], targets: Sequence[int] + ) -> FermionicOp: + """G = i(T - T†); targets in order, source annihilators in reverse order. + + Lists must be nonempty, equally sized, unique and disjoint. Input order + is preserved, with no implicit sorting or factorial normalization. + """ + ... + @property + def num_modes(self) -> int: ... + @property + def num_terms(self) -> int: ... + @property + def terms(self) -> list[tuple[tuple[FermionicLadder, ...], complex]]: + """Independent snapshot in deterministic structural order.""" + ... + def is_zero(self) -> bool: + """Whether the stored sum is empty; use equiv to test algebraic zero.""" + ... + def add_term( + self, term: Iterable[FermionicLadder], coefficient: _Coefficient = 1.0 + ) -> None: + """Modify in place; errors leave this object unchanged.""" + ... + def adjoint(self) -> FermionicOp: ... + def normal_order(self) -> FermionicOp: + """Return a new expansion: creators first, ascending modes in each group.""" + ... + def scaled(self, factor: _Coefficient) -> FermionicOp: ... + def chop(self, *, atol: float = 0.0) -> None: + """Remove terms with complex magnitude <= atol, in place.""" + ... + def equiv(self, other: FermionicOp, *, atol: float = 0.0) -> bool: ... + def is_hermitian(self, *, atol: float = 0.0) -> bool: ... + def __add__(self, other: FermionicOp, /) -> FermionicOp: ... + def __sub__(self, other: FermionicOp, /) -> FermionicOp: ... + def __mul__(self, other: FermionicOp | _Coefficient, /) -> FermionicOp: ... + def __rmul__(self, other: _Coefficient, /) -> FermionicOp: ... + def __neg__(self) -> FermionicOp: ... + def __eq__(self, other: object, /) -> bool: + """Structural equality; use equiv for algebraic equivalence.""" + ... + def copy(self) -> FermionicOp: ... + def __copy__(self) -> FermionicOp: ... + def __deepcopy__(self, memo: dict[int, object]) -> FermionicOp: ... + def __repr__(self) -> str: ... + def __str__(self) -> str: ... + +@runtime_checkable +class FermionToQubitMapper(Protocol): + def map_operator(self, op: FermionicOp) -> Hamiltonian: ... + def num_qubits(self, num_modes: int) -> int: ... + +@runtime_checkable +class OccupationStateMapper(Protocol): + def map_occupation(self, occupation: Sequence[int]) -> list[bool]: ... + +@final +class JordanWignerMapper: + """Mode p maps to qubit p; Z chain uses modes below p. Qubit 0 is LSB.""" + def __init__(self) -> None: ... + def map_operator(self, op: FermionicOp) -> Hamiltonian: ... + def num_qubits(self, num_modes: int) -> int: ... + def map_occupation(self, occupation: Sequence[int]) -> list[bool]: ... + def __repr__(self) -> str: ... + +__all__ = [ + "FermionicLadder", + "FermionicOp", + "FermionToQubitMapper", + "OccupationStateMapper", + "JordanWignerMapper", +] diff --git a/crates/binding-python/cqlib/qis/hamiltonian.py b/crates/binding-python/cqlib/qis/hamiltonian.py index c5f6b95f..20a1a58d 100644 --- a/crates/binding-python/cqlib/qis/hamiltonian.py +++ b/crates/binding-python/cqlib/qis/hamiltonian.py @@ -12,8 +12,12 @@ """Hamiltonian observable wrapper.""" +from typing import TypeAlias + from .._native import qis as _qis_module +Coefficient: TypeAlias = float | int | complex | tuple[float, float] + Hamiltonian = _qis_module.Hamiltonian __all__ = ["Hamiltonian"] diff --git a/crates/binding-python/cqlib/qis/hamiltonian.pyi b/crates/binding-python/cqlib/qis/hamiltonian.pyi index 3a7ee341..61d43dd6 100644 --- a/crates/binding-python/cqlib/qis/hamiltonian.pyi +++ b/crates/binding-python/cqlib/qis/hamiltonian.pyi @@ -10,7 +10,7 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. -from typing import List, Tuple, Union, final, Sequence, Dict +from typing import List, Tuple, Union, final, Sequence, Dict, Iterable import numpy as np @@ -180,6 +180,52 @@ class Hamiltonian: """Returns ``True`` if all Pauli terms mutually commute.""" ... + @staticmethod + def from_terms( + num_qubits: int, + terms: Iterable[tuple[PauliString, complex | float | int]], + ) -> Hamiltonian: + """Validate, absorb phases and merge duplicates, removing exact zeros only.""" + ... + + def normalized(self) -> Hamiltonian: + """Return a validated, merged copy without numerical truncation.""" + ... + + def normalized_hermitian(self) -> Hamiltonian: + """Normalize and project imaginary residuals <= 1e-10; reject larger ones.""" + ... + + def add(self, other: Hamiltonian) -> Hamiltonian: + """Return a validated, merged sum; legacy + and += retain concatenation.""" + ... + + def scaled(self, factor: complex | float | int) -> Hamiltonian: ... + def adjoint(self) -> Hamiltonian: ... + def chop(self, *, atol: float = 0.0) -> None: + """Normalize and truncate by complex magnitude, in place and atomically.""" + ... + + def equiv(self, other: Hamiltonian, *, atol: float = 0.0) -> bool: ... + def is_hermitian(self, *, atol: float = 0.0) -> bool: + """Compare with adjoint; coefficient criterion is 2*abs(imag) <= atol.""" + ... + + def tensor(self, other: Hamiltonian) -> Hamiltonian: + """Return self tensor other, placing other on low qubits.""" + ... + + def remap_qubits(self, mapping: Sequence[int], target_width: int) -> Hamiltonian: + """Embed mapping[source] = destination, with identity on unused targets.""" + ... + + def __sub__(self, other: Hamiltonian, /) -> Hamiltonian: ... + def __mul__(self, other: Hamiltonian | complex | float | int, /) -> Hamiltonian: + """Operator product AB (B acts first), or scalar multiplication.""" + ... + + def __rmul__(self, other: complex | float | int, /) -> Hamiltonian: ... + def __neg__(self) -> Hamiltonian: ... def to_matrix(self) -> np.ndarray: """Computes the dense matrix representation as a NumPy array. @@ -247,7 +293,7 @@ class Hamiltonian: """ ... - def __add__(self, other: "Hamiltonian") -> "Hamiltonian": + def __add__(self, other: "Hamiltonian", /) -> "Hamiltonian": """Adds two Hamiltonians together. Note: This performs a simple lazy concatenation of the term lists. @@ -273,7 +319,7 @@ class Hamiltonian: """ ... - def __iadd__(self, other: "Hamiltonian") -> "Hamiltonian": + def __iadd__(self, other: "Hamiltonian", /) -> "Hamiltonian": """In-place addition of another Hamiltonian. Args: @@ -287,7 +333,7 @@ class Hamiltonian: """ ... - def __eq__(self, other: object) -> bool: + def __eq__(self, other: object, /) -> bool: """Checks if two Hamiltonians are equal.""" ... @@ -365,3 +411,5 @@ class Hamiltonian: def __deepcopy__(self, memo: dict) -> "Hamiltonian": """Return an independent copy via copy.deepcopy().""" ... + +__all__ = ["Hamiltonian"] diff --git a/crates/binding-python/src/circuit/ansatz/hamiltonian_evolution.rs b/crates/binding-python/src/circuit/ansatz/hamiltonian_evolution.rs index ab7e3c6b..dc439f35 100644 --- a/crates/binding-python/src/circuit/ansatz/hamiltonian_evolution.rs +++ b/crates/binding-python/src/circuit/ansatz/hamiltonian_evolution.rs @@ -361,8 +361,9 @@ impl From for PauliEvolutionAnsatz { impl PyPauliEvolutionAnsatz { /// Creates a new ``PauliEvolutionAnsatz`` from a Hamiltonian. /// - /// The Hamiltonian is automatically **simplified** (phases absorbed, duplicate - /// terms merged, near-zero terms removed) before any further processing. + /// The Hamiltonian is normalized (phases absorbed and duplicate terms merged). + /// Imaginary residuals with magnitude <= 1e-10 are projected to zero; larger + /// ones are rejected. Every nonzero real coefficient is retained. /// /// Default strategy is :meth:`EvolutionStrategy.auto` with ``steps=1``. /// diff --git a/crates/binding-python/src/qis/fermion.rs b/crates/binding-python/src/qis/fermion.rs new file mode 100644 index 00000000..83a519cd --- /dev/null +++ b/crates/binding-python/src/qis/fermion.rs @@ -0,0 +1,506 @@ +// 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. + +//! Structured fermionic operators and sparse mappings; algebra stays in Core. + +use super::hamiltonian::PyHamiltonian; +use super::qis_error_to_py_err; +use cqlib_core::qis::{ + FermionToQubitMapper, FermionicLadder, FermionicOp, JordanWignerMapper, OccupationStateMapper, +}; +use num_complex::Complex64; +use pyo3::basic::CompareOp; +use pyo3::exceptions::{PyIndexError, PyValueError}; +use pyo3::prelude::*; +use pyo3::types::{PyComplex, PyList, PyTuple}; +use std::collections::hash_map::DefaultHasher; +use std::hash::{Hash, Hasher}; + +// Respect Python's integer protocol, including type-level method lookup and +// validation of __index__ return values, rather than invoking the method directly. +fn integer_index<'py>(value: &Bound<'py, PyAny>) -> PyResult> { + // SAFETY: value is a live Python object under the GIL. PyNumber_Index returns + // a new owned reference, or NULL with an exception; Bound handles both cases. + unsafe { Bound::from_owned_ptr_or_err(value.py(), pyo3::ffi::PyNumber_Index(value.as_ptr())) } +} + +pub(super) fn index(value: &Bound<'_, PyAny>, is_count: bool) -> PyResult { + let integer = integer_index(value)?; + if integer.rich_compare(0, CompareOp::Lt)?.is_truthy()? { + return Err(if is_count { + PyValueError::new_err("register size must be nonnegative") + } else { + PyIndexError::new_err("mode indices must be nonnegative") + }); + } + integer.extract() +} + +fn indices(values: &Bound<'_, PyAny>) -> PyResult> { + values.try_iter()?.map(|v| index(&v?, false)).collect() +} + +fn ladders(term: &Bound<'_, PyAny>) -> PyResult> { + term.try_iter()? + .map(|v| Ok(v?.extract::>()?.inner)) + .collect() +} + +/// An immutable creation or annihilation operator on a zero-based mode. +#[pyclass( + name = "FermionicLadder", + module = "cqlib.qis.fermion", + frozen, + skip_from_py_object +)] +#[derive(Clone)] +pub struct PyFermionicLadder { + inner: FermionicLadder, +} + +#[pymethods] +impl PyFermionicLadder { + #[staticmethod] + /// Create an immutable creation ladder on a zero-based mode. + fn create(mode: &Bound<'_, PyAny>) -> PyResult { + Ok(Self { + inner: FermionicLadder::Create(index(mode, false)?), + }) + } + + #[staticmethod] + /// Create an immutable annihilation ladder on a zero-based mode. + fn annihilate(mode: &Bound<'_, PyAny>) -> PyResult { + Ok(Self { + inner: FermionicLadder::Annihilate(index(mode, false)?), + }) + } + + #[getter] + /// Zero-based mode index; the containing operator validates its range. + fn mode(&self) -> usize { + self.inner.mode() + } + + #[getter] + /// Ladder kind, either "create" or "annihilate". + fn kind(&self) -> &'static str { + match self.inner { + FermionicLadder::Create(_) => "create", + FermionicLadder::Annihilate(_) => "annihilate", + } + } + + /// Return a ladder on the same mode with creation and annihilation exchanged. + fn adjoint(&self) -> Self { + Self { + inner: self.inner.adjoint(), + } + } + + fn __eq__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + let py = other.py(); + match other.extract::>() { + Ok(rhs) => Ok((self.inner == rhs.inner) + .into_pyobject(py)? + .to_owned() + .into_any() + .unbind()), + Err(_) => Ok(py.NotImplemented()), + } + } + + fn __hash__(&self) -> isize { + let mut hasher = DefaultHasher::new(); + self.inner.hash(&mut hasher); + let hash = hasher.finish() as isize; + if hash == -1 { -2 } else { hash } + } + + fn __repr__(&self) -> String { + format!("FermionicLadder.{}({})", self.kind(), self.mode()) + } +} + +/// A finite complex sum of ordered fermionic products on a declared mode space. +/// +/// In A * B, B acts first. Products are not automatically normal ordered. +/// Only exact zero coefficients are removed; use chop for explicit truncation. +#[pyclass( + name = "FermionicOp", + module = "cqlib.qis.fermion", + skip_from_py_object +)] +#[derive(Clone)] +pub struct PyFermionicOp { + inner: FermionicOp, +} + +impl From for PyFermionicOp { + fn from(inner: FermionicOp) -> Self { + Self { inner } + } +} + +#[pymethods] +impl PyFermionicOp { + /// Construct the zero operator; num_modes is retained even for zero modes. + #[new] + fn new(num_modes: &Bound<'_, PyAny>) -> PyResult { + Ok(FermionicOp::new(index(num_modes, true)?).into()) + } + + #[staticmethod] + /// Return the identity on the declared mode space, including zero modes. + fn identity(num_modes: &Bound<'_, PyAny>) -> PyResult { + Ok(FermionicOp::identity(index(num_modes, true)?).into()) + } + + /// An empty product is identity, whereas an empty sum is zero. + #[staticmethod] + #[pyo3(signature = (num_modes, term, coefficient=Complex64::new(1.0, 0.0)), text_signature = "(num_modes, term, coefficient=1.0)")] + fn from_term( + num_modes: &Bound<'_, PyAny>, + term: &Bound<'_, PyAny>, + coefficient: Complex64, + ) -> PyResult { + FermionicOp::from_term(index(num_modes, true)?, ladders(term)?, coefficient) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Accept an iterable of (ordered ladders, coefficient) pairs and merge duplicates. + #[staticmethod] + fn from_terms(num_modes: &Bound<'_, PyAny>, terms: &Bound<'_, PyAny>) -> PyResult { + let mut result = FermionicOp::new(index(num_modes, true)?); + for item in terms.try_iter()? { + let item = item?; + let pair = item.cast::()?; + if pair.len() != 2 { + return Err(PyValueError::new_err( + "each term must be a pair (ladders, coefficient)", + )); + } + result + .add_term(ladders(&pair.get_item(0)?)?, pair.get_item(1)?.extract()?) + .map_err(qis_error_to_py_err)?; + } + Ok(result.into()) + } + + #[staticmethod] + /// Return a creation operator; invalid mode indices raise IndexError. + fn creation(num_modes: &Bound<'_, PyAny>, mode: &Bound<'_, PyAny>) -> PyResult { + FermionicOp::creation(index(num_modes, true)?, index(mode, false)?) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + #[staticmethod] + /// Return an annihilation operator; invalid mode indices raise IndexError. + fn annihilation(num_modes: &Bound<'_, PyAny>, mode: &Bound<'_, PyAny>) -> PyResult { + FermionicOp::annihilation(index(num_modes, true)?, index(mode, false)?) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + #[staticmethod] + /// Return the mode number operator a†_p a_p. + fn number(num_modes: &Bound<'_, PyAny>, mode: &Bound<'_, PyAny>) -> PyResult { + FermionicOp::number(index(num_modes, true)?, index(mode, false)?) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + #[staticmethod] + /// Return total particle number; zero modes yield the zero operator. + fn total_number(num_modes: &Bound<'_, PyAny>) -> PyResult { + Ok(FermionicOp::total_number(index(num_modes, true)?).into()) + } + + /// G = i(T - T†), with target creators followed by reversed source annihilators. + /// Input order is preserved; no spin-orbital ordering or factorial is imposed. + #[staticmethod] + fn excitation_generator( + num_modes: &Bound<'_, PyAny>, + sources: &Bound<'_, PyAny>, + targets: &Bound<'_, PyAny>, + ) -> PyResult { + FermionicOp::excitation_generator( + index(num_modes, true)?, + &indices(sources)?, + &indices(targets)?, + ) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + #[getter] + /// Declared number of modes, retained even for zero and identity operators. + fn num_modes(&self) -> usize { + self.inner.num_modes() + } + + #[getter] + /// Number of stored ordered products, before algebraic normal ordering. + fn num_terms(&self) -> usize { + self.inner.num_terms() + } + + /// Return an independent snapshot in deterministic structural order. + #[getter] + fn terms<'py>(&self, py: Python<'py>) -> PyResult> { + let result = PyList::empty(py); + for (term, coefficient) in self.inner.terms() { + let term = PyTuple::new( + py, + term.iter() + .map(|&inner| Py::new(py, PyFermionicLadder { inner })) + .collect::>>()?, + )?; + result.append((term, coefficient))?; + } + Ok(result) + } + + /// Test whether storage is empty; use equiv for algebraic zero. + fn is_zero(&self) -> bool { + self.inner.is_zero() + } + + /// Add a term in place and return None; errors leave this operator unchanged. + #[pyo3(signature = (term, coefficient=Complex64::new(1.0, 0.0)), text_signature = "($self, term, coefficient=1.0)")] + fn add_term(&mut self, term: &Bound<'_, PyAny>, coefficient: Complex64) -> PyResult<()> { + self.inner + .add_term(ladders(term)?, coefficient) + .map_err(qis_error_to_py_err) + } + + /// Return a conjugate transpose, reversing products and exchanging ladder kinds. + fn adjoint(&self) -> Self { + self.inner.adjoint().into() + } + + /// Return a normal-ordered expansion, including contractions and fermionic signs. + /// Creators precede annihilators; each group uses ascending mode indices. + fn normal_order(&self) -> PyResult { + self.inner + .normal_order() + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Return a scaled copy; the factor and resulting coefficients must be finite. + fn scaled(&self, factor: Complex64) -> PyResult { + self.inner + .scaled(factor) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Explicitly remove coefficients whose complex magnitude is at most atol, in place. + #[pyo3(signature = (*, atol=0.0))] + fn chop(&mut self, atol: f64) -> PyResult<()> { + self.inner.chop(atol).map_err(qis_error_to_py_err) + } + + /// Compare algebraically after normal ordering the difference. + #[pyo3(signature = (other, *, atol=0.0))] + fn equiv(&self, other: &Self, atol: f64) -> PyResult { + self.inner + .equiv(&other.inner, atol) + .map_err(qis_error_to_py_err) + } + + #[pyo3(signature = (*, atol=0.0))] + /// Compare algebraically with the adjoint; atol must be finite and nonnegative. + fn is_hermitian(&self, atol: f64) -> PyResult { + self.inner.is_hermitian(atol).map_err(qis_error_to_py_err) + } + + fn __add__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + match other.extract::>() { + Ok(rhs) => Ok(Py::new( + other.py(), + Self::from( + self.inner + .try_add(&rhs.inner) + .map_err(qis_error_to_py_err)?, + ), + )? + .into_any()), + Err(_) => Ok(other.py().NotImplemented()), + } + } + + fn __sub__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + match other.extract::>() { + Ok(rhs) => Ok(Py::new( + other.py(), + Self::from( + self.inner + .try_sub(&rhs.inner) + .map_err(qis_error_to_py_err)?, + ), + )? + .into_any()), + Err(_) => Ok(other.py().NotImplemented()), + } + } + + fn __mul__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + let result = if let Ok(rhs) = other.extract::>() { + self.inner.try_mul(&rhs.inner) + } else if let Some(factor) = scalar_operand(other)? { + self.inner.scaled(factor) + } else { + return Ok(other.py().NotImplemented()); + }; + Ok(Py::new(other.py(), Self::from(result.map_err(qis_error_to_py_err)?))?.into_any()) + } + + fn __rmul__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + match scalar_operand(other)? { + Some(factor) => Ok(Py::new(other.py(), self.scaled(factor)?)?.into_any()), + None => Ok(other.py().NotImplemented()), + } + } + + fn __neg__(&self) -> PyResult { + self.scaled(Complex64::new(-1.0, 0.0)) + } + + fn __eq__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + match other.extract::>() { + Ok(rhs) => Ok((self.inner == rhs.inner) + .into_pyobject(other.py())? + .to_owned() + .into_any() + .unbind()), + Err(_) => Ok(other.py().NotImplemented()), + } + } + + /// Return an independent copy of this mutable operator. + fn copy(&self) -> Self { + self.clone() + } + fn __copy__(&self) -> Self { + self.copy() + } + fn __deepcopy__(&self, _memo: &Bound<'_, PyAny>) -> Self { + self.copy() + } + + fn __repr__(&self) -> String { + format!( + "FermionicOp(num_modes={}, num_terms={})", + self.num_modes(), + self.num_terms() + ) + } + + fn __str__(&self) -> String { + if self.is_zero() { + return "0".into(); + } + self.inner + .terms() + .map(|(term, coefficient)| { + let product = if term.is_empty() { + "I".into() + } else { + term.iter() + .map(|ladder| match ladder { + FermionicLadder::Create(mode) => format!("a†_{mode}"), + FermionicLadder::Annihilate(mode) => format!("a_{mode}"), + }) + .collect::>() + .join(" ") + }; + format!("({coefficient}) {product}") + }) + .collect::>() + .join(" + ") + } +} + +// Decline foreign types for Python's reflected arithmetic, but propagate numeric errors. +pub(super) fn scalar_operand(value: &Bound<'_, PyAny>) -> PyResult> { + if value.is_instance_of::() + || value.hasattr("__complex__")? + || value.hasattr("__float__")? + || value.hasattr("__index__")? + { + value.extract().map(Some) + } else { + Ok(None) + } +} + +/// Jordan-Wigner encoding: mode p maps to qubit p, with Z factors on modes below p. +#[pyclass(name = "JordanWignerMapper", module = "cqlib.qis.fermion", frozen)] +pub struct PyJordanWignerMapper; + +#[pymethods] +impl PyJordanWignerMapper { + #[new] + fn new() -> Self { + Self + } + + /// Map to a sparse Hamiltonian; non-Hermitian inputs are valid here. + fn map_operator(&self, op: &PyFermionicOp) -> PyResult { + JordanWignerMapper + .map_operator(&op.inner) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Return the register width; Jordan-Wigner uses one qubit per mode. + fn num_qubits(&self, num_modes: &Bound<'_, PyAny>) -> PyResult { + Ok(JordanWignerMapper.num_qubits(index(num_modes, true)?)) + } + + /// bits[p] addresses qubit p (qubit 0 is the least significant bit). + /// Accept only integer or boolean zero/one occupations. Empty input is vacuum. + fn map_occupation(&self, occupation: &Bound<'_, PyAny>) -> PyResult> { + let mut occupations = Vec::new(); + for value in occupation.try_iter()? { + let value = value?; + let integer = integer_index(&value)?; + let bit = if integer.rich_compare(0, CompareOp::Eq)?.is_truthy()? { + 0 + } else if integer.rich_compare(1, CompareOp::Eq)?.is_truthy()? { + 1 + } else { + return Err(PyValueError::new_err("occupations must be zero or one")); + }; + occupations.push(bit); + } + JordanWignerMapper + .map_occupation(&occupations) + .map_err(qis_error_to_py_err) + } + + fn __repr__(&self) -> &'static str { + "JordanWignerMapper()" + } +} + +pub fn register_fermion_module(parent: &Bound<'_, PyModule>) -> PyResult<()> { + let module = PyModule::new(parent.py(), "fermion")?; + module.add_class::()?; + module.add_class::()?; + module.add_class::()?; + parent.add_submodule(&module) +} diff --git a/crates/binding-python/src/qis/hamiltonian.rs b/crates/binding-python/src/qis/hamiltonian.rs index b81ca7f6..afae3ee4 100644 --- a/crates/binding-python/src/qis/hamiltonian.rs +++ b/crates/binding-python/src/qis/hamiltonian.rs @@ -188,6 +188,166 @@ impl PyHamiltonian { Ok(Self { inner }) } + /// Construct a validated Pauli sum, absorbing phases and merging duplicates. + /// Only exact zeros are removed. Empty input retains the declared width. + #[staticmethod] + fn from_terms(num_qubits: &Bound<'_, PyAny>, terms: &Bound<'_, PyAny>) -> PyResult { + let width = super::fermion::index(num_qubits, true)?; + let mut parsed = Vec::new(); + for item in terms.try_iter()? { + let item = item?; + let pair = item.cast::()?; + if pair.len() != 2 { + return Err(PyValueError::new_err( + "each term must be (PauliString, coefficient)", + )); + } + let pauli = pair + .get_item(0)? + .extract::>()? + .inner + .clone(); + let coefficient = pair.get_item(1)?.extract::()?; + parsed.push((pauli, coefficient)); + } + Hamiltonian::from_terms(width, parsed) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Return a validated sum with phases absorbed and duplicates merged. + /// Retain all nonzero coefficients; legacy simplify() also truncates small terms. + fn normalized(&self) -> PyResult { + self.inner + .normalized() + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Normalize and project imaginary roundoff <= 1e-10 to zero. + /// Larger imaginary coefficients raise ValueError; all nonzero real terms remain. + fn normalized_hermitian(&self) -> PyResult { + self.inner + .normalized_hermitian() + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Return a validated, merged sum without changing either operand. + /// Legacy + and += concatenate stored terms instead. + fn add(&self, other: &Self) -> PyResult { + self.inner + .try_add(&other.inner) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Return a validated scaled copy, retaining all nonzero coefficients. + /// Unlike legacy scale(), validate the factor and resulting coefficients. + fn scaled(&self, factor: Complex64) -> PyResult { + self.inner + .scaled(factor) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Return the conjugate transpose after absorbing Pauli phases. + fn adjoint(&self) -> PyResult { + self.inner + .adjoint() + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Normalize and remove coefficients with complex magnitude <= atol, in place. + /// Return None; errors leave this Hamiltonian unchanged. + #[pyo3(signature = (*, atol=0.0))] + fn chop(&mut self, atol: f64) -> PyResult<()> { + self.inner.chop(atol).map_err(qis_error_to_py_err) + } + + /// Compare normalized coefficient differences, requiring equal register widths. + /// atol is a coefficient tolerance, not a matrix operator-norm bound. + #[pyo3(signature = (other, *, atol=0.0))] + fn equiv(&self, other: &Self, atol: f64) -> PyResult { + self.inner + .equiv(&other.inner, atol) + .map_err(qis_error_to_py_err) + } + + /// Compare with the adjoint; normalized coefficients satisfy 2*abs(imag) <= atol. + /// This check does not project imaginary parts to zero. + #[pyo3(signature = (*, atol=0.0))] + fn is_hermitian(&self, atol: f64) -> PyResult { + self.inner.is_hermitian(atol).map_err(qis_error_to_py_err) + } + + /// Return self tensor other: other occupies low qubits, self occupies high qubits. + fn tensor(&self, other: &Self) -> PyResult { + self.inner + .tensor(&other.inner) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Embed qubits using mapping[source] = destination; unused targets are identity. + /// Destinations must be unique and in range. This is not a fermionic permutation. + fn remap_qubits( + &self, + mapping: &Bound<'_, PyAny>, + target_width: &Bound<'_, PyAny>, + ) -> PyResult { + let destinations = mapping + .try_iter()? + .map(|v| super::fermion::index(&v?, false)) + .collect::>>()?; + self.inner + .remap_qubits(&destinations, super::fermion::index(target_width, true)?) + .map(Into::into) + .map_err(qis_error_to_py_err) + } + + /// Return a validated difference on the same register, without truncation. + fn __sub__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + match other.extract::>() { + Ok(rhs) => Ok(Py::new( + other.py(), + Self::from( + self.inner + .try_sub(&rhs.inner) + .map_err(qis_error_to_py_err)?, + ), + )? + .into_any()), + Err(_) => Ok(other.py().NotImplemented()), + } + } + + /// Return AB (B acts first) or a scalar multiple; validate and merge the result. + fn __mul__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + let result = if let Ok(rhs) = other.extract::>() { + self.inner.try_mul(&rhs.inner) + } else if let Some(factor) = super::fermion::scalar_operand(other)? { + self.inner.scaled(factor) + } else { + return Ok(other.py().NotImplemented()); + }; + Ok(Py::new(other.py(), Self::from(result.map_err(qis_error_to_py_err)?))?.into_any()) + } + + /// Return a validated scalar multiple. + fn __rmul__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + match super::fermion::scalar_operand(other)? { + Some(factor) => Ok(Py::new(other.py(), self.scaled(factor)?)?.into_any()), + None => Ok(other.py().NotImplemented()), + } + } + + /// Return a validated negation without modifying this operator. + fn __neg__(&self) -> PyResult { + self.scaled(Complex64::new(-1.0, 0.0)) + } + /// Adds a new Pauli string term with a given coefficient to the Hamiltonian. /// /// Args: @@ -195,7 +355,7 @@ impl PyHamiltonian { /// coeff: The coefficient for this term. Can be float, int, complex, or tuple (real, imag). /// /// Returns: - /// self (for method chaining) + /// None. This method modifies the Hamiltonian in place. /// /// Raises: /// ValueError: If the number of qubits in the operator does not match the Hamiltonian. @@ -343,15 +503,18 @@ impl PyHamiltonian { /// >>> h2 = Hamiltonian(2) /// >>> h2.add_term(PauliString.from_str("XX"), 0.3) /// >>> h3 = h1 + h2 # Contains both ZZ and XX terms - fn __add__(&self, other: &PyHamiltonian) -> PyResult { - if self.inner.num_qubits != other.inner.num_qubits { + fn __add__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + let Ok(rhs) = other.extract::>() else { + return Ok(other.py().NotImplemented()); + }; + if self.inner.num_qubits != rhs.inner.num_qubits { return Err(PyValueError::new_err(format!( "Cannot add Hamiltonians with different numbers of qubits: {} vs {}", - self.inner.num_qubits, other.inner.num_qubits + self.inner.num_qubits, rhs.inner.num_qubits ))); } - let result = self.inner.clone() + other.inner.clone(); - Ok(Self { inner: result }) + let result = self.inner.clone() + rhs.inner.clone(); + Ok(Py::new(other.py(), Self { inner: result })?.into_any()) } /// In-place addition of another Hamiltonian. @@ -369,8 +532,15 @@ impl PyHamiltonian { Ok(()) } - fn __eq__(&self, other: &PyHamiltonian) -> bool { - self.inner == other.inner + fn __eq__(&self, other: &Bound<'_, PyAny>) -> PyResult> { + match other.extract::>() { + Ok(rhs) => Ok((self.inner == rhs.inner) + .into_pyobject(other.py())? + .to_owned() + .into_any() + .unbind()), + Err(_) => Ok(other.py().NotImplemented()), + } } fn __str__(&self) -> String { diff --git a/crates/binding-python/src/qis/mod.rs b/crates/binding-python/src/qis/mod.rs index 305ef130..daba0b6c 100644 --- a/crates/binding-python/src/qis/mod.rs +++ b/crates/binding-python/src/qis/mod.rs @@ -18,6 +18,7 @@ use pyo3::prelude::*; pub mod entropy; pub mod evolution; +pub mod fermion; pub mod hamiltonian; pub mod metrics; pub mod pauli; @@ -27,7 +28,9 @@ pub mod state; /// Converts a QisError to the appropriate Python exception. pub(crate) fn qis_error_to_py_err(err: QisError) -> PyErr { match err { - QisError::IndexOutOfBounds { .. } => PyIndexError::new_err(err.to_string()), + QisError::IndexOutOfBounds { .. } | QisError::ModeOutOfBounds { .. } => { + PyIndexError::new_err(err.to_string()) + } _ => PyValueError::new_err(err.to_string()), } } @@ -46,6 +49,7 @@ pub fn register_qis_module(parent: &Bound<'_, PyModule>) -> PyResult<()> { qis_module.add_class::()?; qis_module.add_class::()?; qis_module.add_class::()?; + fermion::register_fermion_module(&qis_module)?; // Register state submodule state::register_state_module(&qis_module)?; diff --git a/crates/binding-python/tests/parameter/test_arithmetic.py b/crates/binding-python/tests/parameter/test_arithmetic.py index 94b6a89a..3a769c9c 100644 --- a/crates/binding-python/tests/parameter/test_arithmetic.py +++ b/crates/binding-python/tests/parameter/test_arithmetic.py @@ -26,6 +26,28 @@ from cqlib.circuit import Parameter +def test_affine_coefficients_and_cancellation_residues(): + import math + + theta = Parameter("theta") + for coefficient in (1e-15, math.nextafter(1.0, 2.0)): + scaled = (theta * coefficient).simplify() + assert scaled.evaluate({"theta": 1.0}) == coefficient + assert scaled.replace("theta", Parameter(1.0)).evaluate() == coefficient + coefficient = math.nextafter(1.0, 2.0) + residue = (theta * coefficient - theta).simplify() + assert residue.evaluate({"theta": 1.0}) == coefficient - 1.0 + assert (Parameter(1e-15) - Parameter(0.0)).simplify().evaluate() == 1e-15 + + +def test_subnormal_parameter_division_preserves_representable_quotient(): + smallest = float.fromhex("0x0.0000000000001p-1022") + theta = Parameter("theta") + scaled = theta * smallest + assert scaled.simplify().evaluate({"theta": 1.0}) == smallest + assert (scaled / smallest).simplify().evaluate({"theta": 1.0}) == 1.0 + + class TestAddition: """Tests for addition operation.""" diff --git a/crates/binding-python/tests/qis/test_fermion.py b/crates/binding-python/tests/qis/test_fermion.py new file mode 100644 index 00000000..125f38fd --- /dev/null +++ b/crates/binding-python/tests/qis/test_fermion.py @@ -0,0 +1,378 @@ +# 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. + +"""Fermionic Python contracts and independent Fock-space comparisons.""" + +import copy +import itertools + +import numpy as np +import pytest + +from cqlib.qis import ( + FermionicLadder as Ladder, + FermionicOp as Op, + FermionToQubitMapper, + Hamiltonian, + JordanWignerMapper, + OccupationStateMapper, + Statevector, +) +from cqlib.circuit.ansatz import BasisEncoding, EvolutionStrategy, PauliEvolutionAnsatz + + +def fock_ladder(num_modes, mode, create): + """Independent ladder action on little-endian occupation basis states.""" + matrix = np.zeros((1 << num_modes, 1 << num_modes), dtype=complex) + for column in range(1 << num_modes): + occupied = bool(column & (1 << mode)) + if occupied == create: + continue + lower = column & ((1 << mode) - 1) + sign = (-1) ** bin(lower).count("1") + matrix[column ^ (1 << mode), column] = sign + return matrix + + +def fock_matrix(op): + dim = 1 << op.num_modes + result = np.zeros((dim, dim), dtype=complex) + for term, coefficient in op.terms: + matrix = np.eye(dim, dtype=complex) + for ladder in term: + matrix = matrix @ fock_ladder( + op.num_modes, ladder.mode, ladder.kind == "create" + ) + result += coefficient * matrix + return result + + +def test_exports_and_ladder_contract(): + from cqlib.qis.fermion import FermionicOp + + assert FermionicOp is Op + assert Op.__module__ == "cqlib.qis.fermion" + create = Ladder.create(2) + assert create.mode == 2 and create.kind == "create" + assert create.adjoint() == Ladder.annihilate(2) + assert create.adjoint().adjoint() == create + assert len({create, Ladder.create(2)}) == 1 + assert create != 2 + assert repr(create) == "FermionicLadder.create(2)" + with pytest.raises(AttributeError): + create.mode = 1 + with pytest.raises(TypeError): + Ladder() + + +def test_construction_snapshots_and_exact_zero(): + create = Ladder.create(1) + op = Op.from_terms(3, iter([([create], 2), ((create,), -1), ((), 1j)])) + assert op.num_modes == 3 and op.num_terms == 2 + assert all(isinstance(c, complex) for _, c in op.terms) + snapshot = op.terms + assert all(isinstance(term, tuple) for term, _ in snapshot) + snapshot.clear() + assert op.num_terms == 2 + assert Op.from_terms(3, []).is_zero() + assert Op.from_term(3, []).equiv(Op.identity(3)) + tiny = Op.from_term(3, [create], 1e-15) + assert not tiny.is_zero() + assert not tiny.equiv(Op(3)) + assert tiny.chop(atol=1e-15) is None + assert tiny.is_zero() + with pytest.raises(AttributeError): + op.num_modes = 1 + with pytest.raises(TypeError): + hash(op) + + +def test_arithmetic_adjoint_and_normal_order(): + a = Op.annihilation(2, 0) + b = Op.annihilation(2, 1) + identity = Op.identity(2) + assert (a * a.adjoint() + a.adjoint() * a).equiv(identity) + assert (a * b + b * a).equiv(Op(2)) + assert (a * b.adjoint() + b.adjoint() * a).equiv(Op(2)) + contraction = a * a.adjoint() + assert contraction != contraction.normal_order() + assert contraction.equiv(identity - Op.number(2, 0)) + assert not (a * a).is_zero() # Raw product is stored until normal ordering. + assert (a * a).normal_order().is_zero() + assert (a * b).adjoint() == b.adjoint() * a.adjoint() + assert (2j * a) == a * 2j == a.scaled(2j) + assert -a == a * -1 + assert sum([a, b], start=Op(2)) == a + b + assert Op.total_number(2).is_hermitian() + alias = a + a += b + assert alias == Op.annihilation(2, 0) + assert a == alias + b + + +def test_copy_and_mutation_failure_atomicity(): + op = Op.number(2, 0) + for duplicate in (op.copy(), copy.copy(op), copy.deepcopy(op)): + assert duplicate == op and duplicate is not op + duplicate.add_term([], 2) + assert duplicate != op + before = op.copy() + for term, coefficient, error in ( + ([Ladder.create(2)], 0, IndexError), + (["bad"], 1, TypeError), + ([], complex(float("nan"), 0), ValueError), + ): + with pytest.raises(error): + op.add_term(term, coefficient) + assert op == before + assert "num_modes=2" in repr(op) + assert "a†_0 a_0" in str(op) + assert str(Op(0)) == "0" and "I" in str(Op.identity(0)) + + +@pytest.mark.parametrize("value", [-1, 2**100]) +def test_index_errors(value): + error = IndexError if value < 0 else OverflowError + with pytest.raises(error): + Ladder.create(value) + with pytest.raises(error): + Op.creation(2, value) + with pytest.raises(ValueError if value < 0 else OverflowError): + Op(value) + with pytest.raises(ValueError if value < 0 else OverflowError): + JordanWignerMapper().num_qubits(value) + + +@pytest.mark.parametrize("value", [1.2, "1", None]) +def test_type_errors(value): + with pytest.raises(TypeError): + Ladder.create(value) + with pytest.raises(TypeError): + Op(value) + if not isinstance(value, float): + with pytest.raises(TypeError): + Op.from_term(2, [], value) + + +@pytest.mark.parametrize("atol", [-1, float("nan"), float("inf")]) +def test_invalid_tolerances(atol): + op = Op.number(2, 0) + before = op.copy() + for operation in ( + lambda: op.chop(atol=atol), + lambda: op.equiv(before, atol=atol), + lambda: op.is_hermitian(atol=atol), + ): + with pytest.raises(ValueError): + operation() + assert op == before + + +def test_arithmetic_errors_and_reflected_dispatch(): + a = Op.number(2, 0) + b = Op.number(3, 0) + for operation in (lambda: a + b, lambda: a - b, lambda: a * b, lambda: a.equiv(b)): + with pytest.raises(ValueError): + operation() + for operation in (lambda: a + 1, lambda: a * "x", lambda: "x" * a): + with pytest.raises(TypeError): + operation() + assert a.__add__(object()) is NotImplemented + assert a.__sub__(object()) is NotImplemented + assert a.__mul__(object()) is NotImplemented + assert a.__rmul__(object()) is NotImplemented + assert a.__eq__(object()) is NotImplemented + + class Reflected: + def __rmul__(self, other): + return "reflected" + + assert a * Reflected() == "reflected" + for value in (float("nan"), float("inf"), complex(0, float("inf"))): + with pytest.raises(ValueError): + Op(2) * value + with pytest.raises(ValueError): + Op.from_term(2, [], value) + large = Op.from_term(2, [], 1e308) + before = large.copy() + with pytest.raises(ValueError): + large.add_term([], 1e308) + assert large == before + + +@pytest.mark.parametrize( + "sources,targets", [([], []), ([0], [1, 2]), ([0, 0], [1, 2]), ([0], [0])] +) +def test_excitation_validation(sources, targets): + with pytest.raises(ValueError): + Op.excitation_generator(4, sources, targets) + + +def test_mapping_against_independent_fock_products(): + mapper = JordanWignerMapper() + ladders = [ + factory(p) for factory in (Ladder.create, Ladder.annihilate) for p in range(2) + ] + # Exhaustive ordered products exercise JW phases, contractions and fermionic signs. + for length in range(4): + for term in itertools.product(ladders, repeat=length): + op = Op.from_term(2, term, 0.3 + 0.2j) + expected = fock_matrix(op) + np.testing.assert_allclose( + mapper.map_operator(op).to_matrix(), expected, atol=1e-15 + ) + np.testing.assert_allclose( + fock_matrix(op.normal_order()), expected, atol=1e-15 + ) + a = Op.annihilation(2, 0) + np.testing.assert_allclose( + mapper.map_operator(a.adjoint()).to_matrix(), fock_matrix(a).conj().T + ) + combo = Op.from_terms( + 2, [((), 1j), ([Ladder.create(1)], 0.7), ([Ladder.create(1)], 0.3)] + ) + np.testing.assert_allclose( + mapper.map_operator(combo).to_matrix(), fock_matrix(combo) + ) + for op in (Op.number(2, 1), Op.total_number(2)): + np.testing.assert_allclose(mapper.map_operator(op).to_matrix(), fock_matrix(op)) + + +def test_zero_identity_and_occupations(): + mapper = JordanWignerMapper() + for n in (0, 1, 3): + zero = mapper.map_operator(Op(n)) + assert ( + isinstance(zero, Hamiltonian) + and zero.num_qubits == n + and zero.num_terms == 0 + ) + np.testing.assert_allclose(zero.to_matrix(), np.zeros((1 << n, 1 << n))) + np.testing.assert_allclose( + mapper.map_operator(Op.identity(n)).to_matrix(), np.eye(1 << n) + ) + assert mapper.map_occupation([]) == [] + assert mapper.map_occupation([1, False, True]) == [True, False, True] + for value in (-1, 2, 256, 2**100): + with pytest.raises(ValueError): + mapper.map_occupation([value]) + with pytest.raises(TypeError): + mapper.map_occupation([1.0]) + circuit = BasisEncoding(mapper.map_occupation([1, 0, 1])).build_circuit("hf") + state = Statevector.from_circuit(circuit) + np.testing.assert_allclose(state.data, np.eye(8)[5]) + + +def test_integer_protocol_validation(): + class ValidIndex: + def __index__(self): + return 1 + + class InvalidIndex: + def __index__(self): + return 1.0 + + class InstanceOnlyIndex: + def __init__(self): + self.__index__ = lambda: 1 + + mapper = JordanWignerMapper() + assert Ladder.create(ValidIndex()) == Ladder.create(1) + assert mapper.map_occupation([ValidIndex()]) == [True] + assert mapper.num_qubits(ValidIndex()) == 1 + assert Op(np.int64(2)).num_modes == 2 + assert mapper.map_occupation(np.array([0, 1], dtype=np.int64)) == [False, True] + for value in (InvalidIndex(), InstanceOnlyIndex()): + for operation in ( + lambda: Ladder.create(value), + lambda: Op(value), + lambda: mapper.num_qubits(value), + lambda: mapper.map_occupation([value]), + ): + with pytest.raises(TypeError): + operation() + + +def test_protocols_and_python_extension(): + mapper = JordanWignerMapper() + assert isinstance(mapper, FermionToQubitMapper) + assert isinstance(mapper, OccupationStateMapper) + + class CustomMapper: + def map_operator(self, op): + return mapper.map_operator(op) + + def num_qubits(self, num_modes): + return num_modes + + custom = CustomMapper() + assert isinstance(custom, FermionToQubitMapper) + assert not isinstance(custom, OccupationStateMapper) + np.testing.assert_allclose( + custom.map_operator(Op.number(2, 0)).to_matrix(), fock_matrix(Op.number(2, 0)) + ) + + +def test_mapping_preserves_representable_subnormal_and_near_unit_coefficients(): + smallest = float.fromhex("0x0.0000000000001p-1022") + mapper = JordanWignerMapper() + scalar = Op.from_term(0, [], smallest) + assert scalar.normal_order().terms[0][1].real == smallest + assert mapper.map_operator(scalar).terms[0][1].real == smallest + creation = Op.creation(1, 0).scaled(smallest * 2) + assert all( + abs(coefficient) == smallest + for _, coefficient in mapper.map_operator(creation).terms + ) + assert mapper.map_operator(creation).to_matrix()[1, 0].real == smallest * 2 + almost_one = float(np.nextafter(1.0, 2.0)) + number = Op.number(1, 0).scaled(almost_one) + assert number.terms[0][1].real == almost_one + assert mapper.map_operator(number).to_matrix()[1, 1].real == almost_one + + +def test_excitation_evolution_with_independent_parameters(): + mapper = JordanWignerMapper() + generators = [ + Op.excitation_generator(4, [0, 1], [2, 3]), + Op.excitation_generator(4, [0], [2]), + ] + assert generators[0].is_hermitian() + assert Op.excitation_generator(4, [1, 0], [2, 3]).equiv(-generators[0]) + circuit = BasisEncoding(mapper.map_occupation([1, 1, 0, 0])).build_circuit("hf") + for name, generator in zip(("theta_double", "theta_single"), generators): + block = ( + PauliEvolutionAnsatz(mapper.map_operator(generator)) + .with_strategy(EvolutionStrategy.exact()) + .with_time_param_name(name) + .build_circuit("ignored") + ) + circuit.compose(block) + assert set(circuit.symbols) == {"theta_double", "theta_single"} + angles = {"theta_double": 0.23, "theta_single": -0.37} + partial = circuit.assign_parameters({"theta_double": angles["theta_double"]}) + assert partial.symbols == ["theta_single"] + bound = partial.assign_parameters({"theta_single": angles["theta_single"]}) + actual = Statevector.from_circuit(bound) + expected = np.eye(16, dtype=complex)[3] + for generator, theta in zip(generators, angles.values()): + matrix = fock_matrix(generator) + eigenvalues, eigenvectors = np.linalg.eigh(matrix) + expected = ( + (eigenvectors * np.exp(-1j * theta * eigenvalues)) + @ eigenvectors.conj().T + @ expected + ) + np.testing.assert_allclose(actual.data, expected, atol=1e-14) + assert mapper.map_operator(Op.total_number(4)).expectation_statevector( + actual + ) == pytest.approx(2) diff --git a/crates/binding-python/tests/qis/test_hamiltonian.py b/crates/binding-python/tests/qis/test_hamiltonian.py index 969870fb..1abcf710 100644 --- a/crates/binding-python/tests/qis/test_hamiltonian.py +++ b/crates/binding-python/tests/qis/test_hamiltonian.py @@ -16,6 +16,172 @@ from cqlib.qis import Hamiltonian, PauliString, Phase, TrotterMode +def test_checked_sparse_algebra_and_legacy_compatibility(): + x = Hamiltonian.from_terms(1, [(PauliString.from_str("X"), 1)]) + z = Hamiltonian.from_pauli(PauliString.from_str("Z")) + np.testing.assert_array_equal((x * z).to_matrix(), x.to_matrix() @ z.to_matrix()) + np.testing.assert_array_equal( + (x * z).adjoint().to_matrix(), (x.to_matrix() @ z.to_matrix()).conj().T + ) + assert (x - x).num_terms == 0 + assert (2j * x).equiv(x * 2j) + assert (-x).equiv(x.scaled(-1)) + assert x.add(x).equiv(x * 2) + assert x.add(x).num_terms == 1 + # Existing addition and mutation contracts remain compatible. + legacy = x + x + assert legacy.num_terms == 2 + legacy += x + assert legacy.num_terms == 3 + assert legacy.normalized().num_terms == 1 + assert legacy.add_term(PauliString.from_str("Z"), 1) is None + assert legacy.scale(2) is None + + +def test_checked_normalization_and_truncation(): + p = PauliString.from_str("X") + h = Hamiltonian.from_terms(1, [(p, 1e-15), (PauliString.from_str("+iX"), 2e-15)]) + assert h.num_terms == 1 + assert h.terms[0][1] == complex(1e-15, 2e-15) + assert h.normalized().terms == h.terms + assert not h.equiv(Hamiltonian(1)) + truncated = h.copy() + assert truncated.chop(atol=3e-15) is None + assert truncated.num_terms == 0 and truncated.num_qubits == 1 + legacy = h.copy() + legacy.simplify() + assert legacy.num_terms == 0 # Legacy behavior is deliberately retained. + assert Hamiltonian.from_terms(5, []).num_qubits == 5 + assert Hamiltonian.from_terms(0, []).to_matrix().shape == (1, 1) + + +def test_checked_tensor_and_qubit_embedding(): + x = Hamiltonian.from_pauli(PauliString.from_str("X")) + z = Hamiltonian.from_pauli(PauliString.from_str("Z")) + np.testing.assert_array_equal( + x.tensor(z).to_matrix(), np.kron(x.to_matrix(), z.to_matrix()) + ) + embedded = x.remap_qubits([2], target_width=3) + assert embedded.equiv(Hamiltonian.from_pauli(PauliString.from_str("XII"))) + assert x.tensor(Hamiltonian(2)).num_qubits == 3 + wide = x.remap_qubits([69], target_width=70) + assert wide.num_qubits == 70 and wide.num_terms == 1 + assert wide.equiv(wide.adjoint()) + for mapping, width in (([], 1), ([1], 1), ([0, 1], 2)): + with pytest.raises(ValueError): + x.remap_qubits(mapping, width) + with pytest.raises(ValueError): + x.tensor(z).remap_qubits([0, 0], 2) + with pytest.raises(IndexError): + x.remap_qubits([-1], 1) + with pytest.raises(TypeError): + x.remap_qubits([0.5], 1) + + +def test_checked_validation_and_failure_atomicity(): + x = Hamiltonian.from_pauli(PauliString.from_str("X")) + with pytest.raises(ValueError): + Hamiltonian.from_terms(2, [(PauliString.from_str("X"), 0)]) + with pytest.raises(ValueError): + Hamiltonian.from_terms(1, [(PauliString.from_str("X"), float("nan"))]) + with pytest.raises(TypeError): + Hamiltonian.from_terms(1, [(PauliString.from_str("X"), "invalid")]) + with pytest.raises(ValueError): + Hamiltonian.from_terms(-1, []) + for operation in ( + lambda: x.add(Hamiltonian(2)), + lambda: x - Hamiltonian(2), + lambda: x * Hamiltonian(2), + lambda: x.equiv(Hamiltonian(2)), + lambda: Hamiltonian(1).scaled(float("inf")), + ): + with pytest.raises(ValueError): + operation() + assert x.__mul__(object()) is NotImplemented + assert x.__sub__(object()) is NotImplemented + assert x.__rmul__(object()) is NotImplemented + assert x.__add__(object()) is NotImplemented + assert x.__eq__(object()) is NotImplemented + assert x != object() + for atol in (-1, float("nan"), float("inf")): + before = x.copy() + with pytest.raises(ValueError): + x.chop(atol=atol) + assert x == before + with pytest.raises(ValueError): + x.equiv(before, atol=atol) + with pytest.raises(ValueError): + x.is_hermitian(atol=atol) + overflow = x.scaled(1e308) + with pytest.raises(ValueError): + overflow.add(overflow) + with pytest.raises(ValueError): + overflow * 2 + invalid = x.copy() + invalid.scale(float("inf")) # Legacy scale retains its historical behavior. + before = invalid.copy() + with pytest.raises(ValueError): + invalid.chop() + assert invalid.num_qubits == before.num_qubits + assert invalid.num_terms == before.num_terms + for (_, actual), (_, expected) in zip(invalid.terms, before.terms): + assert actual.real == expected.real + assert math.isnan(actual.imag) and math.isnan(expected.imag) + + +def test_hermitian_projection_preserves_small_real_coefficients(): + h = Hamiltonian.from_terms(1, [(PauliString.from_str("X"), 1e-15 + 1e-11j)]) + assert not h.is_hermitian() + assert h.is_hermitian(atol=2e-11) + projected = h.normalized_hermitian() + assert projected.terms[0][1] == complex(1e-15) + assert h.terms[0][1].imag == 1e-11 # Projection returns a copy. + purely_imaginary = Hamiltonian.from_terms(1, [(PauliString.from_str("X"), 1e-11j)]) + assert purely_imaginary.normalized_hermitian().num_terms == 0 + with pytest.raises(ValueError): + Hamiltonian.from_terms( + 1, [(PauliString.from_str("X"), 2e-10j)] + ).normalized_hermitian() + + +def test_sparse_subnormal_and_near_unit_coefficients(): + smallest = float.fromhex("0x0.0000000000001p-1022") + scalar = Hamiltonian.from_terms(0, [(PauliString(0), smallest)]) + assert scalar.normalized().terms[0][1].real == smallest + assert scalar.scaled(2).terms[0][1].real == smallest * 2 + assert scalar.to_matrix()[0, 0].real == smallest + almost_one = float(np.nextafter(1.0, 2.0)) + x = Hamiltonian.from_pauli(PauliString.from_str("X")) + scaled = almost_one * x + assert scaled.terms[0][1].real == almost_one + assert scaled.adjoint().terms[0][1].real == almost_one + residue = scaled - x + assert residue.terms[0][1].real == almost_one - 1.0 + + +def test_small_symbolic_evolution_preserves_rotation_and_global_phase(): + from cqlib.circuit.ansatz import EvolutionStrategy, PauliEvolutionAnsatz + from cqlib.qis import Statevector + + small = 1e-12 + for label in ("X", "I"): + h = Hamiltonian.from_terms(1, [(PauliString.from_str(label), small)]) + circuit = ( + PauliEvolutionAnsatz(h) + .with_strategy(EvolutionStrategy.exact()) + .with_time_param_name("theta") + .build_circuit("small") + ) + with pytest.raises(ValueError): + Statevector.from_circuit( + circuit + ) # Identity-only symbolic phase must also be bound. + bound = circuit.assign_parameters({"theta": 1.0}) + state = Statevector.from_circuit(bound).data + index = 1 if label == "X" else 0 + assert state[index].imag == pytest.approx(-math.sin(small), rel=1e-14, abs=0) + + def test_hamiltonian_initialization(): h = Hamiltonian(2) assert h.num_qubits == 2 diff --git a/crates/cqlib-core/README.md b/crates/cqlib-core/README.md index 4a7e3911..6498e6b9 100644 --- a/crates/cqlib-core/README.md +++ b/crates/cqlib-core/README.md @@ -77,7 +77,8 @@ assert!(bound.parameters().is_empty()); - **Simulation** — Statevector, density-matrix, noisy density-matrix, and stabilizer simulation - **Quantum information** — Pauli operators, Hamiltonians, observables, - evolution, metrics, and entropy functions + evolution, fermionic algebra, excitation generators, Jordan–Wigner mapping, + metrics, and entropy functions - **Error mitigation** — Zero-noise extrapolation and virtual distillation - **Interchange formats** — OpenQASM 2, OpenQASM 3, and QCIS - **Visualization** — Text and SVG circuit rendering plus state and result plots diff --git a/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution.rs b/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution.rs index e9096ac2..27f2b4a8 100644 --- a/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution.rs +++ b/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution.rs @@ -181,7 +181,7 @@ pub struct EvolutionInfo { /// /// # Requirements /// -/// - The Hamiltonian must be **Hermitian**: after calling `simplify()`, all +/// - The Hamiltonian must be **Hermitian**: after phase normalization and merging, all /// coefficients must have negligible imaginary part ($|\text{Im}(c_k)| < 10^{-10}$). /// - The Hamiltonian must be **non-empty**. /// - For [`EvolutionStrategy::Exact`], all terms must mutually commute. @@ -232,7 +232,7 @@ pub struct PauliEvolutionAnsatz { /// The simplified Hamiltonian (normalized, de-duplicated). /// /// **Invariant**: always in simplified form (phases absorbed, duplicates merged, - /// near-zero terms removed). Any future mutating API must re-run `simplify()` and + /// exact-zero terms removed). Any future mutating API must re-run normalization and /// re-validate Hermiticity before updating this field. hamiltonian: Hamiltonian, /// Cached result of `hamiltonian.all_terms_commute()`, computed once in `new()`. @@ -252,20 +252,20 @@ pub struct PauliEvolutionAnsatz { impl PauliEvolutionAnsatz { /// Creates a new ansatz from a Hamiltonian. /// - /// Internally calls [`Hamiltonian::simplify`] to normalize the Hamiltonian - /// (merge duplicate terms, absorb phases into coefficients, remove near-zero terms) - /// before any further processing. + /// Validates and normalizes the Hamiltonian (merges duplicate terms, absorbs + /// phases into coefficients and removes exact zeros). Small nonzero coefficients + /// are retained, matching the numeric Hamiltonian evolution APIs. Imaginary + /// roundoff residuals of magnitude at most `1e-10` are projected to zero. /// /// The default strategy is [`EvolutionStrategy::Auto`] with `steps = 1`. /// /// # Errors /// /// Returns [`QisError::InvalidParameterValue`] if the Hamiltonian is empty after - /// simplification, or [`QisError::NotHermitian`] if any coefficient has a non-zero - /// imaginary part after simplification. + /// real projection, or [`QisError::NotHermitian`] if any merged coefficient + /// has an imaginary part of magnitude greater than `1e-10`. pub fn new(hamiltonian: Hamiltonian) -> Result { - let mut h = hamiltonian; - h.simplify(); + let h = hamiltonian.normalized_hermitian()?; if h.terms.is_empty() { return Err(QisError::InvalidParameterValue( @@ -274,15 +274,6 @@ impl PauliEvolutionAnsatz { )); } - // Hermitian check: all coefficients must be real after simplify(). - // simplify() absorbs any PauliString phase into the coefficient, so a - // term like (0.5) * (-iY) becomes (-0.5i) * Y — which fails this check. - for (_pauli, coeff) in &h.terms { - if coeff.im.abs() > 1e-10 { - return Err(QisError::NotHermitian); - } - } - // Cache commutativity once; reused by validate(), evolution_info(), build_circuit(). let all_terms_commute = h.all_terms_commute(); diff --git a/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution_test.rs b/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution_test.rs index 024615a4..1b09d36b 100644 --- a/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution_test.rs +++ b/crates/cqlib-core/src/circuit/ansatz/hamiltonian_evolution_test.rs @@ -26,7 +26,7 @@ use ndarray::Array2; use num_complex::Complex64; use crate::circuit::ansatz::Ansatz; -use crate::circuit::circuit_to_matrix; +use crate::circuit::circuit_to_matrix::circuit_to_matrix; use crate::qis::evolution::TrotterMode; use crate::qis::hamiltonian::Hamiltonian; @@ -86,6 +86,120 @@ fn circuit_matrix_at_t( circuit_to_matrix(&bound, None).unwrap() } +#[test] +fn small_coefficients_survive_symbolic_evolution_and_trotter_scaling() { + for coefficient in [1e-16, 1e-12, -1e-16] { + for pauli in ["Z", "I", "X", "Y"] { + let h = Hamiltonian::from_list(vec![( + pauli.parse().unwrap(), + Complex64::new(coefficient, 0.0), + )]) + .unwrap(); + let time = 1.0 / coefficient; + let expected = matrix_exp_iht(&hamiltonian_matrix(&h), time); + for strategy in [ + EvolutionStrategy::Exact, + EvolutionStrategy::Auto { steps: 512 }, + EvolutionStrategy::Trotter { + mode: TrotterMode::FirstOrder, + steps: 512, + }, + EvolutionStrategy::Trotter { + mode: TrotterMode::SecondOrder, + steps: 512, + }, + EvolutionStrategy::Trotter { + mode: TrotterMode::Randomized(7), + steps: 512, + }, + ] { + let ansatz = PauliEvolutionAnsatz::new(h.clone()) + .unwrap() + .with_strategy(strategy.clone()); + let circuit = ansatz.build_circuit("small").unwrap(); + assert_eq!(circuit.symbols().len(), ansatz.num_parameters()); + assert!(circuit.symbols().contains("small_t")); + // An incomplete binding must preserve the time dependence too. + let partial = circuit.assign_parameters(&Some(HashMap::new())).unwrap(); + assert!(partial.symbols().contains("small_t")); + let bound = partial + .assign_parameters(&Some(HashMap::from([("small_t", time)]))) + .unwrap(); + assert!(frob_dist(&circuit_to_matrix(&bound, None).unwrap(), &expected) < 1e-10); + let actual = circuit_matrix_at_t(&ansatz, "small", "small_t", time); + assert!( + frob_dist(&actual, &expected) < 1e-10, + "{pauli}, coefficient={coefficient}, strategy={strategy:?}" + ); + } + } + } +} + +#[test] +fn imaginary_residuals_are_projected_before_ansatz_zero_check() { + for pauli in ["I", "Z"] { + let h = Hamiltonian::from_list(vec![(pauli.parse().unwrap(), Complex64::new(0.0, 1e-12))]) + .unwrap(); + assert!(matches!( + PauliEvolutionAnsatz::new(h), + Err(crate::qis::QisError::InvalidParameterValue(_)) + )); + } + + // Real parts cancel, leaving only an accepted imaginary roundoff residual. + let h = Hamiltonian::from_list(vec![ + ("Z".parse().unwrap(), Complex64::new(1.0, 1e-12)), + ("Z".parse().unwrap(), Complex64::new(-1.0, 0.0)), + ]) + .unwrap(); + assert!(PauliEvolutionAnsatz::new(h).is_err()); + + // Projection must retain even a very small nonzero real generator. + let h = + Hamiltonian::from_list(vec![("Z".parse().unwrap(), Complex64::new(1e-16, 1e-12))]).unwrap(); + let ansatz = PauliEvolutionAnsatz::new(h).unwrap(); + assert_eq!(ansatz.hamiltonian.terms[0].1, Complex64::new(1e-16, 0.0)); + let actual = circuit_matrix_at_t(&ansatz, "residual", "residual_t", 1e16); + let expected = matrix_exp_iht(&hamiltonian_matrix(&ansatz.hamiltonian), 1e16); + assert!(frob_dist(&actual, &expected) < 1e-12); +} + +#[test] +fn mapped_small_number_operator_retains_time_and_global_phase() { + use crate::qis::fermion::{FermionToQubitMapper, FermionicOp, JordanWignerMapper}; + + let op = FermionicOp::number(1, 0) + .unwrap() + .scaled(Complex64::new(1e-16, 0.0)) + .unwrap(); + let h = JordanWignerMapper::new().map_operator(&op).unwrap(); + let ansatz = PauliEvolutionAnsatz::new(h).unwrap(); + let actual = circuit_matrix_at_t(&ansatz, "number", "number_t", 1e16); + // exp(-i n): vacuum has phase 1 and the occupied state has phase exp(-i). + let mut expected = Array2::::eye(2); + expected[[1, 1]] = Complex64::new(0.0, -1.0).exp(); + assert!(frob_dist(&actual, &expected) < 1e-12); +} + +#[test] +fn small_noncommuting_auto_evolution_matches_numeric_trotter() { + let h = Hamiltonian::from_list(vec![ + ("X".parse().unwrap(), Complex64::new(1e-12, 0.0)), + ("Z".parse().unwrap(), Complex64::new(-1e-12, 0.0)), + ]) + .unwrap(); + let numeric = h + .to_trotter_circuit(1e12, 512, TrotterMode::FirstOrder) + .unwrap(); + let expected = circuit_to_matrix(&numeric, None).unwrap(); + let ansatz = PauliEvolutionAnsatz::new(h) + .unwrap() + .with_strategy(EvolutionStrategy::Auto { steps: 512 }); + let actual = circuit_matrix_at_t(&ansatz, "auto", "auto_t", 1e12); + assert!(frob_dist(&actual, &expected) < 1e-12); +} + #[test] fn test_new_simplifies_hamiltonian() { // Add two identical terms: should merge into one with double coefficient diff --git a/crates/cqlib-core/src/circuit/circuit_param.rs b/crates/cqlib-core/src/circuit/circuit_param.rs index 15721c06..01a8e308 100644 --- a/crates/cqlib-core/src/circuit/circuit_param.rs +++ b/crates/cqlib-core/src/circuit/circuit_param.rs @@ -102,8 +102,10 @@ impl ParameterValue { /// /// Fixed values use floating-point multiplication. Symbolic values use /// expression multiplication and remain `Param`, even if the resulting - /// expression is constant (including multiplication by zero). This method - /// does not explicitly evaluate, simplify, or canonicalize expressions. + /// expression is constant (including multiplication by zero). In affine + /// expressions with unresolved symbols, coefficients are preserved without + /// tolerance pruning. + /// This method does not explicitly evaluate or canonicalize expressions. /// Unlike converting a `Parameter` with `Into`, it does not fold constants /// into `Fixed` values. /// diff --git a/crates/cqlib-core/src/circuit/circuit_param_test.rs b/crates/cqlib-core/src/circuit/circuit_param_test.rs index 45690540..8a99fe28 100644 --- a/crates/cqlib-core/src/circuit/circuit_param_test.rs +++ b/crates/cqlib-core/src/circuit/circuit_param_test.rs @@ -59,6 +59,38 @@ fn scaled_constant_expressions_remain_param() { } } +#[test] +fn small_symbolic_scales_survive_canonicalization_and_substitution() { + let input = ParameterValue::Param(Parameter::symbol("theta")); + for factor in [1e-16, -1e-16, 1e-12 / 512.0, f64::from_bits(1)] { + let ParameterValue::Param(scaled) = input.scaled(factor) else { + panic!("symbolic input must remain Param"); + }; + let canonical = scaled.canonicalized().unwrap(); + assert!(canonical.get_symbols().contains("theta")); + assert_eq!( + canonical + .evaluate(&Some(HashMap::from([("theta", 1.0)]))) + .unwrap(), + factor + ); + + let renamed = canonical.replace("theta", Parameter::symbol("time")); + let rescaled = ParameterValue::Param(renamed).scaled(-0.5); + let ParameterValue::Param(rescaled) = rescaled else { + unreachable!() + }; + let canonical = rescaled.canonicalized().unwrap(); + let expected = factor * -0.5; + assert_eq!( + canonical + .evaluate(&Some(HashMap::from([("time", 1e16)]))) + .unwrap(), + expected * 1e16 + ); + } +} + #[test] fn scaled_rejects_non_finite_factors_for_both_variants() { for input in [ diff --git a/crates/cqlib-core/src/circuit/parameter.rs b/crates/cqlib-core/src/circuit/parameter.rs index 10c0ef31..cc497cc9 100644 --- a/crates/cqlib-core/src/circuit/parameter.rs +++ b/crates/cqlib-core/src/circuit/parameter.rs @@ -40,7 +40,237 @@ use std::collections::{HashMap, HashSet}; use std::fmt; use std::ops::{Add, Div, Mul, Neg, Sub}; use std::str::FromStr; -use symb_anafis::{Expr, Simplify, parse, symb}; +use symb_anafis::{Expr, Simplify, parse, symb, visitor::ExprView}; + +// Integral coefficients in the exactly representable integer range cannot +// trigger the engine's near-zero/near-one pruning. Ordinary sums of symbols +// therefore avoid coefficient collection and expression rebuilding entirely. +fn needs_affine_protection(expr: &Expr) -> bool { + match expr.view() { + ExprView::Number(value) => value.fract() != 0.0 || value.abs() > (1u64 << 53) as f64, + ExprView::Sum(children) | ExprView::Product(children) => { + children.iter().any(|child| needs_affine_protection(child)) + } + ExprView::Div(numerator, denominator) => { + needs_affine_protection(numerator) + || denominator.as_number().is_some_and(|value| { + value != 0.0 && needs_affine_protection(&Expr::number(1.0 / value)) + }) + } + _ => false, + } +} + +type AffineCoefficients = HashMap, f64>; + +// Coefficients are collected once per input, using hashed symbol keys. Building +// a replacement AST is deferred until the engine actually changes a coefficient. +fn affine_coefficients(inputs: &[(&Expr, f64)]) -> Option { + fn collect(expr: &Expr, scale: f64, terms: &mut AffineCoefficients) -> Option<()> { + match expr.view() { + ExprView::Number(value) => *terms.entry(None).or_default() += scale * value, + ExprView::Symbol(_) => *terms.entry(Some(expr.clone())).or_default() += scale, + ExprView::Sum(children) => { + for child in children.iter() { + collect(child, scale, terms)?; + } + } + ExprView::Product(children) => { + let mut scale = scale; + let mut symbolic = None; + for child in children.iter() { + if let Some(value) = child.as_number() { + scale *= value; + } else if symbolic.replace(child).is_some() { + return None; + } + } + if let Some(child) = symbolic { + collect(child, scale, terms)?; + } else { + *terms.entry(None).or_default() += scale; + } + } + ExprView::Div(numerator, denominator) => { + let divisor = denominator.as_number()?; + if divisor == 0.0 { + return None; + } + collect(numerator, scale / divisor, terms)?; + } + _ => return None, + } + Some(()) + } + + let mut terms = HashMap::new(); + for (expr, scale) in inputs { + collect(expr, *scale, &mut terms)?; + } + if terms.values().any(|coefficient| !coefficient.is_finite()) { + return None; + } + terms.retain(|_, coefficient| *coefficient != 0.0); + Some(terms) +} + +fn weighted_symbol(symbol: Expr, coefficient: f64) -> Expr { + if coefficient == 0.0 { + return Expr::number(0.0); + } + if coefficient == 1.0 { + return symbol; + } + let product = Expr::product(vec![Expr::number(coefficient), symbol.clone()]); + if let ExprView::Product(factors) = product.view() + && factors + .iter() + .any(|factor| factor.as_number() == Some(coefficient)) + { + return product; + } + // Public map() reconstructs products without approximate numeric pruning. + // If construction rounded a near-unit coefficient, insert its exact value + // into a non-unit product template after construction. + Expr::product(vec![Expr::number(2.0), symbol]).map(|factor| { + if factor.as_number() == Some(2.0) { + Expr::number(coefficient) + } else { + factor.clone() + } + }) +} + +fn affine_expr(terms: &AffineCoefficients) -> Expr { + // Construct a sum of distinct unit-weight symbols before inserting weights, + // preventing automatic polynomial conversion from deleting small terms. + let unit_sum = Expr::sum( + terms + .keys() + .map(|symbol| symbol.clone().unwrap_or_else(|| Expr::number(1.0))) + .collect(), + ); + unit_sum.map(|node| { + if node.as_number() == Some(1.0) { + Expr::number(terms.get(&None).copied().unwrap_or(1.0)) + } else if matches!(node.view(), ExprView::Symbol(_)) { + weighted_symbol( + node.clone(), + terms.get(&Some(node.clone())).copied().unwrap_or(1.0), + ) + } else { + node.clone() + } + }) +} + +fn preserve_affine(expected: Option, fallback: Expr) -> Expr { + match expected { + Some(expected) if expected.keys().all(Option::is_none) => { + Expr::number(expected.get(&None).copied().unwrap_or(0.0)) + } + Some(expected) if affine_coefficients(&[(&fallback, 1.0)]).as_ref() != Some(&expected) => { + affine_expr(&expected) + } + _ => fallback, + } +} + +fn add_parameter_expr(lhs: Expr, rhs: Expr) -> Expr { + if let (Some(a), Some(b)) = (lhs.as_number(), rhs.as_number()) { + return Expr::number(a + b); + } + if !needs_affine_protection(&lhs) && !needs_affine_protection(&rhs) { + return lhs + rhs; + } + let affine = affine_coefficients(&[(&lhs, 1.0), (&rhs, 1.0)]); + preserve_affine(affine, lhs + rhs) +} + +fn sub_parameter_expr(lhs: Expr, rhs: Expr) -> Expr { + if let (Some(a), Some(b)) = (lhs.as_number(), rhs.as_number()) { + return Expr::number(a - b); + } + if !needs_affine_protection(&lhs) && !needs_affine_protection(&rhs) { + return lhs - rhs; + } + let affine = affine_coefficients(&[(&lhs, 1.0), (&rhs, -1.0)]); + preserve_affine(affine, lhs - rhs) +} + +// The common rotation-angle case needs neither a coefficient map nor a second +// symbolic operation. Handle one weighted symbol directly, including rescaling. +fn scaled_symbol(expr: &Expr, scale: f64) -> Option { + let (symbol, coefficient) = match expr.view() { + ExprView::Symbol(_) => (expr.clone(), scale), + ExprView::Product(children) if children.len() == 2 => { + let (number, symbol) = if let Some(number) = children[0].as_number() { + (number, &children[1]) + } else { + (children[1].as_number()?, &children[0]) + }; + if !matches!(symbol.view(), ExprView::Symbol(_)) { + return None; + } + ((**symbol).clone(), number * scale) + } + _ => return None, + }; + coefficient + .is_finite() + .then(|| weighted_symbol(symbol, coefficient)) +} + +fn mul_parameter_expr(lhs: Expr, rhs: Expr) -> Expr { + if let (Some(a), Some(b)) = (lhs.as_number(), rhs.as_number()) { + return Expr::number(a * b); + } + let affine = if let Some(scale) = lhs.as_number() { + if let Some(expr) = scaled_symbol(&rhs, scale) { + return expr; + } + if !needs_affine_protection(&lhs) && !needs_affine_protection(&rhs) { + return lhs * rhs; + } + affine_coefficients(&[(&rhs, scale)]) + } else if let Some(scale) = rhs.as_number() { + if let Some(expr) = scaled_symbol(&lhs, scale) { + return expr; + } + if !needs_affine_protection(&lhs) && !needs_affine_protection(&rhs) { + return lhs * rhs; + } + affine_coefficients(&[(&lhs, scale)]) + } else { + None + }; + preserve_affine(affine, lhs * rhs) +} + +fn div_parameter_expr(lhs: Expr, rhs: Expr) -> Expr { + if let (Some(a), Some(b)) = (lhs.as_number(), rhs.as_number()) + && b != 0.0 + { + return Expr::number(a / b); + } + let affine = if let Some(divisor) = rhs.as_number().filter(|value| *value != 0.0) { + affine_coefficients(&[(&lhs, 1.0)]).and_then(|mut terms| { + // Divide coefficients directly: the reciprocal of a subnormal + // divisor can overflow even when the quotient is representable. + for coefficient in terms.values_mut() { + *coefficient /= divisor; + } + if terms.values().any(|coefficient| !coefficient.is_finite()) { + return None; + } + terms.retain(|_, coefficient| *coefficient != 0.0); + Some(terms) + }) + } else { + None + }; + preserve_affine(affine, lhs / rhs) +} /// A symbolic or numeric parameter used to represent gate angles and other /// continuous values in a quantum circuit. @@ -261,6 +491,11 @@ impl Parameter { /// /// Simplification is performed in a domain-safe mode to avoid invalid /// transformations (e.g. taking the log of a potentially negative number). + /// For affine sums, scalar products and division by numeric constants, + /// nonzero `f64` coefficients and cancellation results are preserved without + /// tolerance pruning. Numeric literals are retained exactly. Floating-point + /// rounding, overflow and underflow still apply. Parsing and general + /// nonlinear simplification use the symbolic engine's own numeric policy. /// /// # Examples /// @@ -273,8 +508,16 @@ impl Parameter { /// assert_eq!(simplified.to_string(), "θ"); /// ``` pub fn simplify(&self) -> Result { + if self.expr.as_number().is_some() { + return Ok(self.clone()); + } + let affine = needs_affine_protection(&self.expr) + .then(|| affine_coefficients(&[(&self.expr, 1.0)])) + .flatten(); let expr = Simplify::new().domain_safe(true).simplify(&self.expr)?; - Ok(Self { expr }) + Ok(Self { + expr: preserve_affine(affine, expr), + }) } /// Returns the numeric value of this parameter without symbolic @@ -292,10 +535,11 @@ impl Parameter { /// Returns the canonical storage form of this parameter expression. /// - /// The parameter is first simplified. If the simplified expression has no - /// free symbols, it is evaluated to a concrete finite value and rebuilt as - /// a numeric parameter. This also normalizes `-0.0` to `0.0`, giving circuit - /// operation storage a stable representation for equivalent constants. + /// A constant expression is evaluated before symbolic simplification. If + /// simplification eliminates all symbols, the result is also evaluated and + /// rebuilt as a finite numeric parameter. This also normalizes `-0.0` to + /// `0.0`, giving circuit operation storage a stable representation for + /// equivalent constants. /// /// Symbolic expressions are returned in their simplified form. This method /// does not intern symbolic parameters into any circuit; use @@ -308,6 +552,11 @@ impl Parameter { /// Returns [`ParameterError`] if simplification fails, a constant expression /// cannot be evaluated, or the evaluated constant is not finite. pub fn canonicalized(&self) -> Result { + // Evaluate constant expressions before any approximate symbolic rewrite. + if self.get_symbols().is_empty() { + let value = self.evaluate(&None)?; + return Ok(Parameter::from(if value == 0.0 { 0.0 } else { value })); + } let simplified = self.simplify()?; if simplified.get_symbols().is_empty() { @@ -812,7 +1061,8 @@ impl From for Parameter { /// type `$t` in both directions (e.g. `Parameter + f64` and `f64 + Parameter`). /// /// The implementation converts `$t` to `Parameter` via `Parameter::from`, extracts the -/// underlying `Expr`, and delegates to the corresponding `Expr` operator. +/// underlying `Expr`, and preserves affine coefficients before falling back to +/// the corresponding symbolic engine operator. macro_rules! impl_ops_for_type { ($t:ty) => { // Parameter + T @@ -822,8 +1072,8 @@ macro_rules! impl_ops_for_type { Self { // 1. `Parameter::from(rhs)` converts the primitive to `Parameter`. // 2. `.expr` unwraps the underlying `Expr` node. - // 3. `self.expr + ...` invokes `symb_anafis`'s `Add` impl. - expr: self.expr + Parameter::from(rhs).expr, + // 3. Preserve affine coefficients before general symbolic addition. + expr: add_parameter_expr(self.expr, Parameter::from(rhs).expr), } } } @@ -833,7 +1083,7 @@ macro_rules! impl_ops_for_type { type Output = Parameter; fn add(self, rhs: Parameter) -> Self::Output { Parameter { - expr: Parameter::from(self).expr + rhs.expr, + expr: add_parameter_expr(Parameter::from(self).expr, rhs.expr), } } } @@ -843,7 +1093,7 @@ macro_rules! impl_ops_for_type { type Output = Parameter; fn sub(self, rhs: $t) -> Self::Output { Self { - expr: self.expr - Parameter::from(rhs).expr, + expr: sub_parameter_expr(self.expr, Parameter::from(rhs).expr), } } } @@ -853,7 +1103,7 @@ macro_rules! impl_ops_for_type { type Output = Parameter; fn sub(self, rhs: Parameter) -> Self::Output { Parameter { - expr: Parameter::from(self).expr - rhs.expr, + expr: sub_parameter_expr(Parameter::from(self).expr, rhs.expr), } } } @@ -863,7 +1113,7 @@ macro_rules! impl_ops_for_type { type Output = Parameter; fn mul(self, rhs: $t) -> Self::Output { Self { - expr: self.expr * Parameter::from(rhs).expr, + expr: mul_parameter_expr(self.expr, Parameter::from(rhs).expr), } } } @@ -873,7 +1123,7 @@ macro_rules! impl_ops_for_type { type Output = Parameter; fn mul(self, rhs: Parameter) -> Self::Output { Parameter { - expr: Parameter::from(self).expr * rhs.expr, + expr: mul_parameter_expr(Parameter::from(self).expr, rhs.expr), } } } @@ -883,7 +1133,7 @@ macro_rules! impl_ops_for_type { type Output = Parameter; fn div(self, rhs: $t) -> Self::Output { Self { - expr: self.expr / Parameter::from(rhs).expr, + expr: div_parameter_expr(self.expr, Parameter::from(rhs).expr), } } } @@ -893,7 +1143,7 @@ macro_rules! impl_ops_for_type { type Output = Parameter; fn div(self, rhs: Parameter) -> Self::Output { Parameter { - expr: Parameter::from(self).expr / rhs.expr, + expr: div_parameter_expr(Parameter::from(self).expr, rhs.expr), } } } @@ -919,15 +1169,15 @@ impl_ops_for_type!(u32); /// Cloning `Expr` is inexpensive because `symb_anafis` uses reference-counted /// node sharing internally. macro_rules! impl_binary_op_ref { - ($($trait:ident, $method:ident),* $(,)?) => { + ($($trait:ident, $method:ident, $expr_method:ident);* $(;)?) => { $( // 1. Parameter op Parameter — both sides consumed. impl $trait for Parameter { type Output = Parameter; fn $method(self, rhs: Parameter) -> Self::Output { Parameter { - // Delegate directly to the underlying `Expr` operator. - expr: self.expr.$method(rhs.expr), + // Preserve affine coefficients before general symbolic arithmetic. + expr: $expr_method(self.expr, rhs.expr), } } } @@ -939,7 +1189,7 @@ macro_rules! impl_binary_op_ref { Parameter { // `Expr::clone` is cheap — the tree uses Arc-based // node sharing inside `symb_anafis`. - expr: self.expr.clone().$method(rhs.expr.clone()), + expr: $expr_method(self.expr.clone(), rhs.expr.clone()), } } } @@ -949,7 +1199,7 @@ macro_rules! impl_binary_op_ref { type Output = Parameter; fn $method(self, rhs: &'a Parameter) -> Self::Output { Parameter { - expr: self.expr.$method(rhs.expr.clone()), + expr: $expr_method(self.expr, rhs.expr.clone()), } } } @@ -959,7 +1209,7 @@ macro_rules! impl_binary_op_ref { type Output = Parameter; fn $method(self, rhs: Parameter) -> Self::Output { Parameter { - expr: self.expr.clone().$method(rhs.expr), + expr: $expr_method(self.expr.clone(), rhs.expr), } } } @@ -968,10 +1218,10 @@ macro_rules! impl_binary_op_ref { } impl_binary_op_ref! { - Add, add, - Sub, sub, - Mul, mul, - Div, div, + Add, add, add_parameter_expr; + Sub, sub, sub_parameter_expr; + Mul, mul, mul_parameter_expr; + Div, div, div_parameter_expr; } // Implement Neg for owned Parameter @@ -979,7 +1229,7 @@ impl Neg for Parameter { type Output = Parameter; fn neg(self) -> Self::Output { Parameter { - expr: self.expr.neg(), + expr: mul_parameter_expr(Expr::number(-1.0), self.expr), } } } @@ -989,7 +1239,7 @@ impl Neg for &Parameter { type Output = Parameter; fn neg(self) -> Self::Output { Parameter { - expr: self.expr.clone().neg(), + expr: mul_parameter_expr(Expr::number(-1.0), self.expr.clone()), } } } diff --git a/crates/cqlib-core/src/circuit/parameter_test.rs b/crates/cqlib-core/src/circuit/parameter_test.rs index 8fe576f2..e26af669 100644 --- a/crates/cqlib-core/src/circuit/parameter_test.rs +++ b/crates/cqlib-core/src/circuit/parameter_test.rs @@ -11,6 +11,8 @@ // that they have been altered from the originals. use super::Parameter; +use crate::circuit::circuit_impl::Circuit; +use crate::circuit::circuit_to_matrix::circuit_to_matrix; use std::collections::{HashMap, HashSet}; use std::f64::consts::{E, PI}; @@ -1136,6 +1138,111 @@ fn test_canonicalized_folds_constant_expression() { assert_eq!(canonical.evaluate(&None).unwrap(), 0.0); } +#[test] +fn affine_cancellation_preserves_small_constants_without_false_symbols() { + let t = Parameter::symbol("t"); + let u = Parameter::symbol("u"); + let q: Parameter = t * 1e-16 + u * 2e-16; + for constant in [3e-16, -3e-16, f64::from_bits(1), 0.0] { + let p = q.clone() + constant; + let result = (p - &q).canonicalized().unwrap(); + assert!(result.get_symbols().is_empty(), "{constant:?}: {result}"); + assert_eq!(result.evaluate(&None).unwrap(), constant); + assert_eq!( + result.simplify().unwrap().evaluate(&None).unwrap(), + constant + ); + assert_eq!(result.canonicalized().unwrap(), result); + assert_eq!( + (result.clone() + 0.0_f64).evaluate(&None).unwrap(), + constant + ); + if constant.abs() == 3e-16 { + let angle: Parameter = result * 1e16; + let mut circuit = Circuit::new(1); + circuit.rz(circuit.qubits()[0], angle).unwrap(); + assert!(circuit.symbols().is_empty()); + let matrix = circuit_to_matrix(&circuit, None).unwrap(); + let expected = num_complex::Complex64::from_polar(1.0, -constant * 1e16 / 2.0); + assert!((matrix[[0, 0]] - expected).norm() < 1e-14); + } + } +} + +#[test] +fn affine_scaling_preserves_near_unit_and_subnormal_weights() { + let t = Parameter::symbol("t"); + let bindings = Some(HashMap::from([("t", 1.0)])); + for coefficient in [ + f64::from_bits(1), + 1e-16, + -1e-16, + 0.375, + 1.0 - f64::EPSILON, + 1.0 + f64::EPSILON, + ] { + let weighted = t.clone() * coefficient; + let canonical = weighted.canonicalized().unwrap(); + assert!(canonical.get_symbols().contains("t")); + assert_eq!(canonical.evaluate(&bindings).unwrap(), coefficient); + let rescaled: Parameter = canonical * 2.0; + assert_eq!( + rescaled + .canonicalized() + .unwrap() + .evaluate(&bindings) + .unwrap(), + coefficient * 2.0 + ); + } +} + +#[test] +fn affine_constant_division_preserves_small_coefficients() { + let t = Parameter::symbol("t"); + let divisor = Parameter::from(1e16); + for (expr, expected) in [ + (t.clone() / 1e16, 1.0), + (&t / &divisor, 1.0), + ((t.clone() * 1e-16) / 2.0, 0.5), + ((t.clone() / 1e16) / 2.0, 0.5), + (t.clone() / -1e16, -1.0), + ] { + let canonical = expr.canonicalized().unwrap(); + assert!(canonical.get_symbols().contains("t")); + assert_eq!( + canonical + .evaluate(&Some(HashMap::from([("t", 1e16)]))) + .unwrap(), + expected + ); + } + let tiny = f64::from_bits(1); + assert_eq!((Parameter::from(tiny) / tiny).evaluate(&None).unwrap(), 1.0); + let quotient: Parameter = (t.clone() * tiny) / tiny; + assert_eq!( + quotient + .canonicalized() + .unwrap() + .evaluate(&Some(HashMap::from([("t", 1.0)]))) + .unwrap(), + 1.0 + ); + assert_eq!( + Parameter::from(tiny) + .canonicalized() + .unwrap() + .evaluate(&None) + .unwrap(), + tiny + ); + let cancelled: Parameter = t.clone() / 1e16 - t / 1e16; + assert_eq!( + cancelled.canonicalized().unwrap().evaluate(&None).unwrap(), + 0.0 + ); +} + #[test] fn test_canonicalized_normalizes_negative_zero() { let canonical = Parameter::from(-0.0).canonicalized().unwrap(); diff --git a/crates/cqlib-core/src/compile/transform/canonicalize/canonicalize_test.rs b/crates/cqlib-core/src/compile/transform/canonicalize/canonicalize_test.rs index 5488b9de..566e5fb1 100644 --- a/crates/cqlib-core/src/compile/transform/canonicalize/canonicalize_test.rs +++ b/crates/cqlib-core/src/compile/transform/canonicalize/canonicalize_test.rs @@ -11,11 +11,13 @@ // that they have been altered from the originals. use super::{CanonicalizeConfig, Canonicalizer, canonicalize_circuit}; +use crate::circuit::circuit_impl::Circuit; +use crate::circuit::circuit_to_matrix::circuit_to_matrix; use crate::circuit::gate::FrozenCircuit; use crate::circuit::{ - Circuit, CircuitError, CircuitGate, CircuitParam, ClassicalControlOp, ClassicalDataOp, - ClassicalExpr, ClassicalType, Directive, Instruction, MCGate, Parameter, ParameterValue, Qubit, - StandardGate, UnitaryGate, ValueInstruction, ValueOperation, circuit_to_matrix, + CircuitError, CircuitGate, CircuitParam, ClassicalControlOp, ClassicalDataOp, ClassicalExpr, + ClassicalType, Directive, Instruction, MCGate, Parameter, ParameterValue, Qubit, StandardGate, + UnitaryGate, ValueInstruction, ValueOperation, }; use crate::compile::test_utils::generated_small_matrix_circuit; use crate::compile::transform::{RewriteEdits, TransformOutcome, TransformerTestExt}; @@ -1043,19 +1045,43 @@ fn canonicalize_phase_fast_path_does_not_accept_negative_zero() { } #[test] -fn canonical_transform_normalizes_global_phase_accumulator_residue() { +fn canonical_transform_preserves_nonzero_global_phase_residue() { let mut circuit = Circuit::new(1); - circuit.set_global_phase(Parameter::from(-3.330_669_073_875_469_6e-16)); + let phase = Parameter::from(-3.330_669_073_875_469_6e-16); + circuit.set_global_phase(phase.clone()); circuit.h(Qubit::new(0)).unwrap(); let first = Canonicalizer::production().run(&circuit).unwrap(); - assert!(first.changed); - assert_eq!(first.circuit.global_phase(), Parameter::from(0.0)); + assert!(!first.changed); + assert_eq!(first.circuit.global_phase(), phase); let second = Canonicalizer::production().run(&first.circuit).unwrap(); assert!(!second.changed); } +#[test] +fn canonical_transform_preserves_small_nonzero_rotation_angles() { + let mut circuit = Circuit::new(1); + circuit.rz(Qubit::new(0), Parameter::from(3e-16)).unwrap(); + let result = Canonicalizer::production().run(&circuit).unwrap(); + assert_eq!(result.circuit.operations().len(), 1); + assert_eq!( + result + .circuit + .resolve_parameter(&result.circuit.operations()[0].params[0]) + .unwrap() + .evaluate(&None) + .unwrap(), + 3e-16 + ); + assert!( + !Canonicalizer::production() + .run(&result.circuit) + .unwrap() + .changed + ); +} + #[test] fn canonicalization_is_idempotent_for_mixed_production_input() { let mut circuit = Circuit::new(3); diff --git a/crates/cqlib-core/src/compile/transform/canonicalize/canonicalizer.rs b/crates/cqlib-core/src/compile/transform/canonicalize/canonicalizer.rs index df075031..5a0ee49c 100644 --- a/crates/cqlib-core/src/compile/transform/canonicalize/canonicalizer.rs +++ b/crates/cqlib-core/src/compile/transform/canonicalize/canonicalizer.rs @@ -267,15 +267,6 @@ impl<'a> CanonicalizeRound<'a> { fn run(mut self) -> Result<(Circuit, Vec>), CompilerError> { self.top_phase = self.canonicalize_parameter(self.top_phase.clone())?; - if !self.source.operations().is_empty() && !self.top_phase.is_canonical_positive_zero() { - // The previous operation loop normalized `current + 0` once for - // every source operation, including operations with no phase - // contribution. Preserve that representation for tiny numeric - // residues with one up-front normalization instead of repeating - // the same symbolic work for every operation. - self.top_phase = - self.canonicalize_parameter(self.top_phase.clone() + Parameter::from(0.0))?; - } let mut top_level = Vec::with_capacity(self.source.operations().len()); let mut provenance = Vec::with_capacity(self.source.operations().len()); diff --git a/crates/cqlib-core/src/compile/transform/canonicalize/verify.rs b/crates/cqlib-core/src/compile/transform/canonicalize/verify.rs index 337b3cf4..a15fdee9 100644 --- a/crates/cqlib-core/src/compile/transform/canonicalize/verify.rs +++ b/crates/cqlib-core/src/compile/transform/canonicalize/verify.rs @@ -72,23 +72,6 @@ pub fn verify_circuit(circuit: &Circuit, mode: VerifyMode<'_>) -> Result<(), Com verify_parameter_finite(&global_phase, "global phase")?; if matches!(mode, VerifyMode::Output { .. }) { verify_parameter_canonical(&global_phase, circuit.global_phase_param(), "global phase")?; - // A canonicalization round accumulates a zero phase contribution for - // every retained root operation. The expression engine can normalize - // a tiny numeric residue while constructing that sum even when the - // residue alone is already a canonical numeric parameter. Reject that - // representation here so the verified fast path is a true fixed point. - if !circuit.operations().is_empty() && !global_phase.is_canonical_positive_zero() { - let accumulated = (global_phase.clone() + Parameter::from(0.0)) - .canonicalized() - .map_err(|error| { - CompilerError::Circuit(crate::circuit::CircuitError::InvalidParameter(error)) - })?; - if accumulated != global_phase { - return Err(CompilerError::InvariantViolation( - "global phase is not in accumulator-normal form".to_string(), - )); - } - } } verify_operations(circuit, circuit.operations(), "root", mode)?; diff --git a/crates/cqlib-core/src/device/noise.rs b/crates/cqlib-core/src/device/noise.rs index 903367f4..7b861e45 100644 --- a/crates/cqlib-core/src/device/noise.rs +++ b/crates/cqlib-core/src/device/noise.rs @@ -558,4 +558,36 @@ mod tests { result => panic!("unexpected result: {result:?}"), } } + + #[test] + fn pauli_probability_boundary_has_complete_finite_kraus_operators() { + 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 + ); + } + } + assert!( + !SingleQubitNoise::Pauli { + px: 0.4, + py: 0.6, + pz: 0.1 + } + .is_valid() + ); + } } diff --git a/crates/cqlib-core/src/error_mitigation/zne_mitigation_test.rs b/crates/cqlib-core/src/error_mitigation/zne_mitigation_test.rs index 282cd262..4ded2ee0 100644 --- a/crates/cqlib-core/src/error_mitigation/zne_mitigation_test.rs +++ b/crates/cqlib-core/src/error_mitigation/zne_mitigation_test.rs @@ -347,3 +347,30 @@ fn test_extrapolate_api_exponential() { .unwrap(); assert!((extrapolated - a).abs() < 1e-10); } + +#[test] +fn folding_preserves_full_matrix_including_global_phase() { + use crate::circuit::Parameter; + 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) + ); + } + } + } +} diff --git a/crates/cqlib-core/src/qis/entropy_test.rs b/crates/cqlib-core/src/qis/entropy_test.rs index 1c8e6f24..d6d1be9c 100644 --- a/crates/cqlib-core/src/qis/entropy_test.rs +++ b/crates/cqlib-core/src/qis/entropy_test.rs @@ -761,3 +761,105 @@ fn test_eof_concurrence_relationship() { assert!(c_sep.abs() < EPSILON); assert!(eof_sep.abs() < EPSILON); } + +#[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 renyi_entropy_of_nondiagonal_pure_state_is_zero() { + let state = [1.0, 2.0, 3.0, 4.0] + .map(|v| Complex64::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 = [ + Complex64::new(1.0, 0.0), + Complex64::new(0.0, 1.0), + Complex64::new(-1.0, 0.0), + Complex64::new(0.0, -1.0), + ]; + let mut rho = Vec::new(); + let mut diagonal = vec![Complex64::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((Complex64::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); + } + } +} diff --git a/crates/cqlib-core/src/qis/error.rs b/crates/cqlib-core/src/qis/error.rs index 808e764c..6cbf3aff 100644 --- a/crates/cqlib-core/src/qis/error.rs +++ b/crates/cqlib-core/src/qis/error.rs @@ -11,6 +11,7 @@ // that they have been altered from the originals. use crate::circuit::CircuitError; +use num_complex::Complex64; use thiserror::Error; /// Comprehensive error type for the `qis` (Quantum Information Science) module. @@ -19,6 +20,21 @@ use thiserror::Error; /// density matrix simulation, observables, Hamiltonians, and Pauli operator manipulations. #[derive(Debug, Error)] pub enum QisError { + /// A fermionic mode lies outside the declared mode space. + #[error("Fermionic mode {mode} is out of bounds for {num_modes} modes")] + ModeOutOfBounds { mode: usize, num_modes: usize }, + + /// Fermionic operators must declare the same mode count. + #[error("Mode count mismatch: expected {expected}, got {actual}")] + ModeMismatch { expected: usize, actual: usize }, + + /// Operator coefficients and arithmetic results must be finite. + #[error("Operator coefficient must be finite, got {coefficient}")] + NonFiniteCoefficient { coefficient: Complex64 }, + + /// A fermionic occupation must be either zero or one. + #[error("Invalid occupation at mode {index}: expected 0 or 1, got {value}")] + InvalidOccupation { index: usize, value: u8 }, /// Errors originating from quantum circuit operations, such as unsupported gates /// or symbolic parameters. #[error(transparent)] @@ -90,6 +106,14 @@ pub enum QisError { NonCliffordGate(String), } +pub(crate) fn validate_coefficient(coefficient: Complex64) -> Result<(), QisError> { + if coefficient.re.is_finite() && coefficient.im.is_finite() { + Ok(()) + } else { + Err(QisError::NonFiniteCoefficient { coefficient }) + } +} + /// Error type for parsing PauliString from a string representation. #[derive(Debug, Clone, PartialEq, Eq, Error)] pub enum PauliStringParseError { diff --git a/crates/cqlib-core/src/qis/evolution.rs b/crates/cqlib-core/src/qis/evolution.rs index e81dc3a9..dc4bde31 100644 --- a/crates/cqlib-core/src/qis/evolution.rs +++ b/crates/cqlib-core/src/qis/evolution.rs @@ -449,6 +449,16 @@ pub(crate) fn trotter_second_order_core( Ok(()) } +// Numeric evolution historically reports coefficient overflow as unsupported. +fn normalized_evolution(h: &Hamiltonian) -> Result { + h.normalized_hermitian().map_err(|error| match error { + nonfinite @ QisError::NonFiniteCoefficient { .. } => { + QisError::UnsupportedOperation(nonfinite.to_string()) + } + other => other, + }) +} + impl Hamiltonian { /// Converts the Hamiltonian to a Trotterized time evolution circuit. /// @@ -459,6 +469,8 @@ impl Hamiltonian { /// the Hamiltonian terms commute. For an auto-selecting method that uses the /// mathematically exact single-pass decomposition when all terms commute, see /// [`to_evolution_circuit`](Hamiltonian::to_evolution_circuit). + /// Nonempty sums that cancel to zero produce an identity circuit, retaining + /// the declared qubit count. An empty input sum is rejected. /// /// # Arguments /// @@ -521,15 +533,7 @@ impl Hamiltonian { )); } - let mut simplified_h = self.clone(); - simplified_h.simplify_with_tolerance(0.0); - - for (_, coeff) in &simplified_h.terms { - if coeff.im.abs() > 1e-10 { - return Err(QisError::NotHermitian); - } - } - + let simplified_h = normalized_evolution(self)?; let mut circuit = Circuit::new(simplified_h.num_qubits); let qubits: Vec = circuit.qubits(); let t = ParameterValue::Fixed(time); @@ -591,6 +595,8 @@ impl Hamiltonian { /// This mirrors the behavior of [`PauliEvolutionAnsatz`] with /// [`EvolutionStrategy::Auto`], closing the capability gap between the numeric /// and parametric APIs. + /// Nonempty sums that cancel to zero produce an identity circuit, retaining + /// the declared qubit count. An empty input sum is rejected. /// /// # Arguments /// @@ -641,15 +647,7 @@ impl Hamiltonian { )); } - let mut simplified_h = self.clone(); - simplified_h.simplify_with_tolerance(0.0); - - for (_, coeff) in &simplified_h.terms { - if coeff.im.abs() > 1e-10 { - return Err(QisError::NotHermitian); - } - } - + let simplified_h = normalized_evolution(self)?; let mut circuit = Circuit::new(simplified_h.num_qubits); let qubits: Vec = circuit.qubits(); diff --git a/crates/cqlib-core/src/qis/evolution_test.rs b/crates/cqlib-core/src/qis/evolution_test.rs index 31312e9b..9ace77e2 100644 --- a/crates/cqlib-core/src/qis/evolution_test.rs +++ b/crates/cqlib-core/src/qis/evolution_test.rs @@ -732,3 +732,27 @@ fn test_evolution_circuit_exact_matrix_correctness() { assert_abs_diff_eq!(matrix[[0, 1]].re, 0.0, epsilon = 1e-12); assert_abs_diff_eq!(matrix[[1, 0]].re, 0.0, epsilon = 1e-12); } + +#[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()); + } +} diff --git a/crates/cqlib-core/src/qis/fermion/mapping/jordan_wigner.rs b/crates/cqlib-core/src/qis/fermion/mapping/jordan_wigner.rs new file mode 100644 index 00000000..db60aaf9 --- /dev/null +++ b/crates/cqlib-core/src/qis/fermion/mapping/jordan_wigner.rs @@ -0,0 +1,106 @@ +// 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. + +//! Jordan–Wigner mapping using sparse Pauli arithmetic. + +use super::{FermionToQubitMapper, OccupationStateMapper}; +use crate::qis::fermion::{FermionicLadder, FermionicOp}; +use crate::qis::{Hamiltonian, Pauli, PauliString, QisError}; +use num_complex::Complex64; +use std::collections::BTreeMap; + +/// Jordan–Wigner encoding with mode p mapped to qubit p. +/// +/// The Z chain covers modes `0..p`: creation is `Z_chain (X_p - iY_p)/2`, +/// annihilation is `Z_chain (X_p + iY_p)/2`. Qubit zero is the least-significant +/// basis bit, and the rightmost character of a displayed Pauli string. +/// Mapping never constructs a dense matrix or implicitly normal orders its input. +/// A product of k ladders can expand to 2^k terms before sparse merging. +#[derive(Debug, Clone, Copy, Default, PartialEq, Eq)] +pub struct JordanWignerMapper; + +impl JordanWignerMapper { + /// Creates a stateless mapper. + pub fn new() -> Self { + Self + } +} + +impl FermionToQubitMapper for JordanWignerMapper { + fn num_qubits(&self, num_modes: usize) -> usize { + num_modes + } + + fn map_operator(&self, op: &FermionicOp) -> Result { + let n = self.num_qubits(op.num_modes()); + let mut result = Hamiltonian::new(n); + // Build only the components needed by this input, and reuse them across terms. + let mut components = BTreeMap::new(); + for (term, coefficient) in op.terms() { + let mut product = Hamiltonian::new(n); + // Merge the unit-weight product before applying the input coefficient, + // so representable small results are not lost in ladder intermediates. + product.add_term(PauliString::new(n), Complex64::new(1.0, 0.0))?; + for ladder in term { + let p = ladder.mode(); + let (x, y) = components.entry(p).or_insert_with(|| { + let mut x = PauliString::new(n); + let mut y = PauliString::new(n); + for j in 0..p { + x.set_pauli(j, Pauli::Z); + y.set_pauli(j, Pauli::Z); + } + x.set_pauli(p, Pauli::X); + y.set_pauli(p, Pauli::Y); + (x, y) + }); + let imaginary = match ladder { + FermionicLadder::Create(_) => -0.5, + FermionicLadder::Annihilate(_) => 0.5, + }; + let mut next = Hamiltonian::new(n); + for (pauli, c) in &product.terms { + next.add_term(pauli * &*x, *c * 0.5)?; + next.add_term(pauli * &*y, *c * Complex64::new(0.0, imaginary))?; + } + product = next.normalized()?; + if product.terms.is_empty() { + break; + } + } + result.terms.extend( + product + .terms + .into_iter() + .map(|(pauli, c)| (pauli, c * coefficient)), + ); + } + result.normalized() + } +} + +impl OccupationStateMapper for JordanWignerMapper { + fn map_occupation(&self, occupation: &[u8]) -> Result, QisError> { + occupation + .iter() + .enumerate() + .map(|(index, value)| match value { + 0 => Ok(false), + 1 => Ok(true), + value => Err(QisError::InvalidOccupation { + index, + value: *value, + }), + }) + .collect() + } +} diff --git a/crates/cqlib-core/src/qis/fermion/mapping/mod.rs b/crates/cqlib-core/src/qis/fermion/mapping/mod.rs new file mode 100644 index 00000000..b7fed5dd --- /dev/null +++ b/crates/cqlib-core/src/qis/fermion/mapping/mod.rs @@ -0,0 +1,37 @@ +// 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. + +//! Contracts for fermion-to-qubit operators and occupation-state encodings. + +pub mod jordan_wigner; +pub use jordan_wigner::JordanWignerMapper; + +use super::FermionicOp; +use crate::qis::{Hamiltonian, QisError}; + +/// Maps finite fermionic operators, including non-Hermitian operators, to Pauli sums. +/// Circuit synthesis and Hermiticity requirements belong to downstream consumers. +pub trait FermionToQubitMapper { + /// Maps an operator while preserving its declared mode space. + fn map_operator(&self, op: &FermionicOp) -> Result; + /// Returns the qubit count required to encode the given mode space. + fn num_qubits(&self, num_modes: usize) -> usize; +} + +/// Encodes occupations into computational-basis bits, where `bits[p]` addresses qp. +/// This is separate from operator mapping because some encodings do not map +/// occupations to individual computational basis states. +pub trait OccupationStateMapper { + /// Maps a vector of zero/one occupations; its length defines the mode count. + /// An empty vector encodes the zero-mode vacuum as an empty bit vector. + fn map_occupation(&self, occupation: &[u8]) -> Result, QisError>; +} diff --git a/crates/cqlib-core/src/qis/fermion/mod.rs b/crates/cqlib-core/src/qis/fermion/mod.rs new file mode 100644 index 00000000..324fcdbe --- /dev/null +++ b/crates/cqlib-core/src/qis/fermion/mod.rs @@ -0,0 +1,74 @@ +// 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. + +//! Fermionic algebra, qubit mappings and occupation encodings. +//! +//! Modes are numbered from zero. Ordered products act right-to-left. The Fock +//! basis is `|n_{N-1} ... n_0> = (a†_0)^n_0 ... (a†_{N-1})^n_{N-1} |vac>`. +//! Qubit zero is the least-significant bit; state indices are `sum(n_p * 2^p)`. +//! Core does not assign spin or chemical-orbital meanings to these indices. +//! Zero and identity operators retain their declared mode count, including zero +//! modes. Circuit consumers can have stricter requirements: basis encoding needs +//! nonempty bits, and `PauliEvolutionAnsatz` needs a nonzero Hermitian operator. +//! Numeric evolution rejects empty input sums but returns identity evolution +//! when a nonempty sum cancels to zero. +//! +//! # Particle-number observable +//! ``` +//! use cqlib_core::qis::{FermionicOp, FermionToQubitMapper, JordanWignerMapper, +//! Observable, Statevector}; +//! let mapper = JordanWignerMapper::new(); +//! let number = mapper.map_operator(&FermionicOp::total_number(3))?; +//! let mut state = Statevector::new(3); +//! state.apply_x(0); +//! state.apply_x(2); +//! assert_eq!(number.expectation_statevector(&state)?, 2.0); +//! # Ok::<(), cqlib_core::qis::QisError>(()) +//! ``` +//! +//! # Parameterized excitation evolution +//! ``` +//! use cqlib_core::qis::{FermionicOp, +//! FermionToQubitMapper, JordanWignerMapper}; +//! use cqlib_core::circuit::ansatz::{Ansatz, PauliEvolutionAnsatz}; +//! let h = FermionicOp::excitation_generator(4, &[0, 1], &[2, 3])?; +//! // e^(-i H theta) = e^(theta (T - T†)). +//! let mapped = JordanWignerMapper::new().map_operator(&h)?; +//! let circuit = PauliEvolutionAnsatz::new(mapped)?.build_circuit("excitation")?; +//! assert_eq!(circuit.num_qubits(), 4); +//! # Ok::<(), cqlib_core::qis::QisError>(()) +//! ``` +//! +//! Each excitation can have an independent parameter: build one evolution +//! block per generator with a distinct `with_time_param_name`, then append +//! the blocks using [`crate::circuit::Circuit::compose`]. A single ansatz built +//! from a sum of generators instead has one shared time parameter. +//! +//! # Occupation-state encoding +//! ``` +//! use cqlib_core::qis::{JordanWignerMapper, OccupationStateMapper, Statevector}; +//! use cqlib_core::circuit::ansatz::{Ansatz, BasisEncoding}; +//! let bits = JordanWignerMapper::new().map_occupation(&[1, 0, 0])?; +//! let circuit = BasisEncoding::new(bits).build_circuit("reference")?; +//! let state = Statevector::from_circuit(&circuit)?; +//! assert_eq!(state.probabilities()[1], 1.0); // displayed as |001> +//! # Ok::<(), cqlib_core::qis::QisError>(()) +//! ``` + +pub mod mapping; +pub mod operator; + +pub use mapping::{FermionToQubitMapper, JordanWignerMapper, OccupationStateMapper}; +pub use operator::{FermionicLadder, FermionicOp}; + +#[cfg(test)] +mod tests; diff --git a/crates/cqlib-core/src/qis/fermion/operator.rs b/crates/cqlib-core/src/qis/fermion/operator.rs new file mode 100644 index 00000000..4b234d96 --- /dev/null +++ b/crates/cqlib-core/src/qis/fermion/operator.rs @@ -0,0 +1,516 @@ +// 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. + +//! Structured fermionic operators and canonical anticommutation algebra. + +use crate::qis::QisError; +use crate::qis::error::validate_coefficient; +use num_complex::Complex64; +use std::collections::{BTreeMap, BTreeSet}; +use std::ops::{Add, Mul, Neg, Sub}; + +/// A creation or annihilation operator on a zero-based fermionic mode. +/// +/// A sequence `[A, B]` denotes the mathematical product `A B`: `B` acts first. +#[derive(Debug, Clone, Copy, PartialEq, Eq, PartialOrd, Ord, Hash)] +pub enum FermionicLadder { + /// Creates a particle in the given mode. + Create(usize), + /// Removes a particle from the given mode. + Annihilate(usize), +} + +impl FermionicLadder { + /// Returns the mode on which this operator acts. + pub fn mode(self) -> usize { + match self { + Self::Create(mode) | Self::Annihilate(mode) => mode, + } + } + + /// Exchanges creation and annihilation. + pub fn adjoint(self) -> Self { + match self { + Self::Create(mode) => Self::Annihilate(mode), + Self::Annihilate(mode) => Self::Create(mode), + } + } +} + +/// A finite complex linear combination of ordered fermionic products. +/// +/// The declared mode count is retained even for zero and identity operators. +/// Equal ordered products are merged automatically; only exact zero coefficients +/// are removed. Products are not automatically normal ordered. Use [`Self::chop`] +/// for explicit numerical truncation and [`Self::equiv`] for algebraic equivalence. +/// Structural equality also compares the declared mode count. +/// +/// Fallible arithmetic operators return `Result`, for example `(&a + &b)?`. +#[derive(Debug, Clone, PartialEq)] +pub struct FermionicOp { + num_modes: usize, + terms: BTreeMap, Complex64>, +} + +pub(crate) fn validate_tolerance(atol: f64) -> Result<(), QisError> { + if !atol.is_finite() || atol < 0.0 { + return Err(QisError::InvalidParameterValue( + "tolerance must be finite and nonnegative".into(), + )); + } + Ok(()) +} + +impl FermionicOp { + /// Creates the zero operator on the declared mode space, including zero modes. + pub fn new(num_modes: usize) -> Self { + Self { + num_modes, + terms: BTreeMap::new(), + } + } + + /// Creates the identity. An empty product is the identity, not the zero operator. + pub fn identity(num_modes: usize) -> Self { + let mut result = Self::new(num_modes); + result.terms.insert(Vec::new(), Complex64::new(1.0, 0.0)); + result + } + + /// Creates a weighted ordered product. Even zero-weight terms are validated. + pub fn from_term( + num_modes: usize, + term: impl IntoIterator, + coefficient: Complex64, + ) -> Result { + let mut result = Self::new(num_modes); + result.add_term(term, coefficient)?; + Ok(result) + } + + /// Creates a sum of products, merging duplicates without numerical truncation. + pub fn from_terms( + num_modes: usize, + terms: impl IntoIterator, Complex64)>, + ) -> Result { + let mut result = Self::new(num_modes); + for (term, coefficient) in terms { + result.add_term(term, coefficient)?; + } + Ok(result) + } + + /// Adds a weighted product atomically; an error leaves this operator unchanged. + pub fn add_term( + &mut self, + term: impl IntoIterator, + coefficient: Complex64, + ) -> Result<(), QisError> { + validate_coefficient(coefficient)?; + let term: Vec<_> = term.into_iter().collect(); + for ladder in &term { + let mode = ladder.mode(); + if mode >= self.num_modes { + return Err(QisError::ModeOutOfBounds { + mode, + num_modes: self.num_modes, + }); + } + } + let merged = self.terms.get(&term).copied().unwrap_or_default() + coefficient; + validate_coefficient(merged)?; + if merged == Complex64::new(0.0, 0.0) { + self.terms.remove(&term); + } else { + self.terms.insert(term, merged); + } + Ok(()) + } + + /// Creates a single creation operator. + pub fn creation(num_modes: usize, mode: usize) -> Result { + Self::from_term(num_modes, [FermionicLadder::Create(mode)], 1.0.into()) + } + + /// Creates a single annihilation operator. + pub fn annihilation(num_modes: usize, mode: usize) -> Result { + Self::from_term(num_modes, [FermionicLadder::Annihilate(mode)], 1.0.into()) + } + + /// Creates the number operator `a†_p a_p`. + pub fn number(num_modes: usize, mode: usize) -> Result { + Self::from_term( + num_modes, + [ + FermionicLadder::Create(mode), + FermionicLadder::Annihilate(mode), + ], + 1.0.into(), + ) + } + + /// Creates the total particle-number operator, or zero for zero modes. + pub fn total_number(num_modes: usize) -> Self { + let terms = (0..num_modes) + .map(|mode| { + ( + vec![ + FermionicLadder::Create(mode), + FermionicLadder::Annihilate(mode), + ], + 1.0.into(), + ) + }) + .collect(); + Self { num_modes, terms } + } + + /// Creates the unit-weight Hermitian excitation generator `G = i(T - T†)`. + /// + /// For `sources = [i_0, ..., i_{r-1}]` and `targets = [a_0, ..., a_{r-1}]`, + /// the ordered excitation is + /// `T = a†_{a_0} ... a†_{a_{r-1}} a_{i_{r-1}} ... a_{i_0}`. + /// The rightmost operator acts first. Thus `[i, j] -> [a, b]` means + /// `T = a†_a a†_b a_j a_i`. + /// + /// Input order is preserved, without implicit sorting or normal ordering. + /// An odd permutation of either list negates `G`; applying odd permutations + /// to both lists preserves it. Explicit [`Self::normal_order`] preserves + /// these fermionic signs. No factorial or other normalization is applied. + /// + /// Mapping `G` and evolving it with [`crate::circuit::ansatz::PauliEvolutionAnsatz`] + /// implements `exp(-i θ G) = exp(θ (T - T†))`. The evolution parameter belongs + /// to the ansatz, not this operator. The generator conserves total particle + /// number. Reference occupations and spin interpretations are caller choices. + /// + /// # Errors + /// + /// Returns [`QisError::InvalidParameterValue`] if the lists are empty, have + /// unequal lengths, contain repeated modes, or overlap. Returns + /// [`QisError::ModeOutOfBounds`] if an index is outside `0..num_modes`. + /// + /// # Example + /// ``` + /// use cqlib_core::qis::FermionicOp; + /// let single = FermionicOp::excitation_generator(4, &[0], &[2])?; + /// let double = FermionicOp::excitation_generator(4, &[0, 1], &[2, 3])?; + /// assert!(single.is_hermitian(0.0)?); + /// assert!(double.is_hermitian(0.0)?); + /// # Ok::<(), cqlib_core::qis::QisError>(()) + /// ``` + pub fn excitation_generator( + num_modes: usize, + sources: &[usize], + targets: &[usize], + ) -> Result { + if sources.len() != targets.len() { + return Err(QisError::InvalidParameterValue( + "excitation sources and targets must have equal lengths".into(), + )); + } + if sources.is_empty() { + return Err(QisError::InvalidParameterValue( + "excitation sources and targets must be nonempty".into(), + )); + } + + let mut source_modes = BTreeSet::new(); + let mut target_modes = BTreeSet::new(); + for (label, modes, seen) in [ + ("sources", sources, &mut source_modes), + ("targets", targets, &mut target_modes), + ] { + for &mode in modes { + if mode >= num_modes { + return Err(QisError::ModeOutOfBounds { mode, num_modes }); + } + if !seen.insert(mode) { + return Err(QisError::InvalidParameterValue(format!( + "excitation {label} contain repeated mode {mode}" + ))); + } + } + } + if let Some(mode) = source_modes.intersection(&target_modes).next() { + return Err(QisError::InvalidParameterValue(format!( + "excitation sources and targets overlap at mode {mode}" + ))); + } + + let term = targets.iter().copied().map(FermionicLadder::Create).chain( + sources + .iter() + .rev() + .copied() + .map(FermionicLadder::Annihilate), + ); + let excitation = Self::from_term(num_modes, term, 1.0.into())?; + excitation + .try_sub(&excitation.adjoint())? + .scaled(Complex64::new(0.0, 1.0)) + } + + /// Returns the declared mode count. + pub fn num_modes(&self) -> usize { + self.num_modes + } + + /// Returns the number of stored ordered products. + pub fn num_terms(&self) -> usize { + self.terms.len() + } + + /// Iterates over products and their coefficients in deterministic structural order. + pub fn terms(&self) -> impl ExactSizeIterator + '_ { + self.terms + .iter() + .map(|(term, coefficient)| (term.as_slice(), *coefficient)) + } + + /// Reports whether the stored sum is zero. Algebraic cancellations require normal ordering. + pub fn is_zero(&self) -> bool { + self.terms.is_empty() + } + + fn check_modes(&self, other: &Self) -> Result<(), QisError> { + if self.num_modes != other.num_modes { + return Err(QisError::ModeMismatch { + expected: self.num_modes, + actual: other.num_modes, + }); + } + Ok(()) + } + + /// Adds operators on the same mode space, without modifying either operand. + pub fn try_add(&self, other: &Self) -> Result { + self.check_modes(other)?; + let mut result = self.clone(); + for (term, coefficient) in other.terms() { + result.add_term(term.iter().copied(), coefficient)?; + } + Ok(result) + } + + /// Subtracts operators on the same mode space. + pub fn try_sub(&self, other: &Self) -> Result { + self.check_modes(other)?; + let mut result = self.clone(); + for (term, coefficient) in other.terms() { + result.add_term(term.iter().copied(), -coefficient)?; + } + Ok(result) + } + + /// Forms `self * other`; the right operand acts first. Does not normal order. + pub fn try_mul(&self, other: &Self) -> Result { + self.check_modes(other)?; + let mut result = Self::new(self.num_modes); + for (left, a) in self.terms() { + for (right, b) in other.terms() { + result.add_term(left.iter().chain(right).copied(), a * b)?; + } + } + Ok(result) + } + + /// Scales all coefficients. The factor is validated even for the zero operator. + pub fn scaled(&self, factor: Complex64) -> Result { + validate_coefficient(factor)?; + Self::from_terms( + self.num_modes, + self.terms + .iter() + .map(|(term, c)| (term.clone(), *c * factor)), + ) + } + + /// Takes the conjugate transpose, reversing each product without reordering it. + pub fn adjoint(&self) -> Self { + Self { + num_modes: self.num_modes, + terms: self + .terms + .iter() + .map(|(term, coefficient)| { + ( + term.iter().rev().map(|ladder| ladder.adjoint()).collect(), + coefficient.conj(), + ) + }) + .collect(), + } + } + + /// Expands using the canonical anticommutation relations. + /// + /// Creators precede annihilators; each group uses ascending mode indices. + /// This is an equivalent operator expansion, including contraction terms: + /// `a_p a†_p = I - a†_p a_p`. Small nonzero coefficients are retained. + /// An explicit work stack avoids recursion in long products. + /// Sorting a branch of k ladders takes O(k²) comparisons; contractions can + /// produce additional branches. + pub fn normal_order(&self) -> Result { + let mut result = Self::new(self.num_modes); + let mut pending: Vec<_> = self + .terms + .iter() + .map(|(term, c)| (term.clone(), *c, 0usize)) + .collect(); + while let Some((mut term, mut coefficient, mut index)) = pending.pop() { + let mut complete = true; + while index + 1 < term.len() { + let left = term[index]; + let right = term[index + 1]; + if left == right { + complete = false; + break; + } + if left > right { + // Swapping distinct fermionic ladders contributes a minus sign. + // Equal-mode annihilator/creator pairs also contribute a contraction. + if let (FermionicLadder::Annihilate(p), FermionicLadder::Create(q)) = + (left, right) + && p == q + { + let mut contracted = term.clone(); + contracted.drain(index..index + 2); + // The prefix is already sorted. Removing this pair only + // requires rechecking the new boundary with that prefix. + pending.push((contracted, coefficient, index.saturating_sub(1))); + } + term.swap(index, index + 1); + coefficient = -coefficient; + // Revisit the changed boundary, rather than scanning from zero. + index = index.saturating_sub(1); + } else { + index += 1; + } + } + if complete { + result.add_term(term, coefficient)?; + } + } + Ok(result) + } + + /// Explicitly removes coefficients whose complex magnitude is at most `atol`. + /// Invalid tolerances leave the operator unchanged. + pub fn chop(&mut self, atol: f64) -> Result<(), QisError> { + validate_tolerance(atol)?; + self.terms + .retain(|_, coefficient| coefficient.norm() > atol); + Ok(()) + } + + /// Compares operators algebraically after normal ordering their difference. + /// Different mode counts are an error; `atol` must be finite and nonnegative. + pub fn equiv(&self, other: &Self, atol: f64) -> Result { + validate_tolerance(atol)?; + let difference = self.try_sub(other)?.normal_order()?; + Ok(difference + .terms + .values() + .all(|coefficient| coefficient.norm() <= atol)) + } + + /// Checks Hermiticity algebraically, allowing non-normal-ordered inputs. + pub fn is_hermitian(&self, atol: f64) -> Result { + self.equiv(&self.adjoint(), atol) + } +} + +macro_rules! binary_operator { + ($trait:ident, $method:ident, $checked:ident) => { + impl $trait<&FermionicOp> for &FermionicOp { + type Output = Result; + fn $method(self, rhs: &FermionicOp) -> Self::Output { + self.$checked(rhs) + } + } + impl $trait for FermionicOp { + type Output = Result; + fn $method(self, rhs: FermionicOp) -> Self::Output { + self.$checked(&rhs) + } + } + impl $trait<&FermionicOp> for FermionicOp { + type Output = Result; + fn $method(self, rhs: &FermionicOp) -> Self::Output { + self.$checked(rhs) + } + } + impl $trait for &FermionicOp { + type Output = Result; + fn $method(self, rhs: FermionicOp) -> Self::Output { + self.$checked(&rhs) + } + } + }; +} + +binary_operator!(Add, add, try_add); +binary_operator!(Sub, sub, try_sub); +binary_operator!(Mul, mul, try_mul); + +macro_rules! scalar_operator { + ($scalar:ty) => { + impl Mul<$scalar> for &FermionicOp { + type Output = Result; + fn mul(self, rhs: $scalar) -> Self::Output { + self.scaled(rhs.into()) + } + } + impl Mul<$scalar> for FermionicOp { + type Output = Result; + fn mul(self, rhs: $scalar) -> Self::Output { + self.scaled(rhs.into()) + } + } + impl Mul<&FermionicOp> for $scalar { + type Output = Result; + fn mul(self, rhs: &FermionicOp) -> Self::Output { + rhs.scaled(self.into()) + } + } + impl Mul for $scalar { + type Output = Result; + fn mul(self, rhs: FermionicOp) -> Self::Output { + rhs.scaled(self.into()) + } + } + }; +} + +scalar_operator!(f64); +scalar_operator!(Complex64); + +impl Neg for &FermionicOp { + type Output = FermionicOp; + fn neg(self) -> Self::Output { + let mut result = self.clone(); + for coefficient in result.terms.values_mut() { + *coefficient = -*coefficient; + } + result + } +} + +impl Neg for FermionicOp { + type Output = Self; + fn neg(mut self) -> Self::Output { + for coefficient in self.terms.values_mut() { + *coefficient = -*coefficient; + } + self + } +} diff --git a/crates/cqlib-core/src/qis/fermion/tests.rs b/crates/cqlib-core/src/qis/fermion/tests.rs new file mode 100644 index 00000000..5a90d144 --- /dev/null +++ b/crates/cqlib-core/src/qis/fermion/tests.rs @@ -0,0 +1,1119 @@ +// 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 super::*; +use crate::circuit::ansatz::{Ansatz, BasisEncoding, EvolutionStrategy, PauliEvolutionAnsatz}; +use crate::circuit::circuit_to_matrix::circuit_to_matrix; +use crate::qis::{ + DensityMatrix, Hamiltonian, Observable, PauliString, Phase, QisError, Statevector, TrotterMode, +}; +use FermionicLadder::{Annihilate, Create}; +use ndarray::Array2; +use num_complex::Complex64; +use std::collections::{BTreeMap, HashMap}; + +// Independent Fock-basis action: no Pauli matrices or mapping formulas are used. +fn fock_matrix(op: &FermionicOp) -> Array2 { + let dim = 1 << op.num_modes(); + let mut matrix = Array2::zeros((dim, dim)); + for (word, coefficient) in op.terms() { + for source in 0..dim { + let mut state = source; + let mut amplitude = coefficient; + for ladder in word.iter().rev() { + let p = ladder.mode(); + let occupied = state & (1 << p) != 0; + if matches!(ladder, Create(_)) == occupied { + amplitude = Complex64::new(0.0, 0.0); + break; + } + if (state & ((1 << p) - 1)).count_ones() % 2 == 1 { + amplitude = -amplitude; + } + state ^= 1 << p; + } + matrix[[state, source]] += amplitude; + } + } + matrix +} + +fn assert_matrix_close(left: &Array2, right: &Array2) { + assert_eq!(left.dim(), right.dim()); + for (a, b) in left.iter().zip(right) { + assert!((*a - *b).norm() < 1e-11, "{a} != {b}"); + } +} + +// Independent occupation transitions for T: remove sources in input order, +// then insert targets in reverse order. This does not inspect FermionicOp words +// or use its algebra, normal ordering or JW mapping. +fn excitation_fock_matrix(n: usize, sources: &[usize], targets: &[usize]) -> Array2 { + let dim = 1usize << n; + let mut transition = Array2::::zeros((dim, dim)); + for initial in 0..dim { + let mut state = initial; + let mut sign = 1.0; + let mut allowed = true; + for (&mode, create) in sources + .iter() + .map(|p| (p, false)) + .chain(targets.iter().rev().map(|p| (p, true))) + { + if (state & (1 << mode) != 0) == create { + allowed = false; + break; + } + if (state & ((1 << mode) - 1)).count_ones() % 2 != 0 { + sign = -sign; + } + state ^= 1 << mode; + } + if allowed { + transition[[state, initial]] = Complex64::new(sign, 0.0); + } + } + (&transition - &transition.t().mapv(|c| c.conj())).mapv(|c| Complex64::new(0.0, 1.0) * c) +} + +#[test] +fn excitation_generator_validates_all_index_constraints() { + for (sources, targets) in [ + (vec![], vec![]), + (vec![], vec![1]), + (vec![0, 1], vec![2]), + (vec![0, 0], vec![2, 3]), + (vec![0, 1], vec![2, 2]), + (vec![0], vec![0]), + (vec![0, 1], vec![1, 2]), + ] { + assert!(matches!( + FermionicOp::excitation_generator(4, &sources, &targets), + Err(QisError::InvalidParameterValue(_)) + )); + } + for (n, sources, targets, invalid) in [ + (0, vec![0], vec![1], 0), + (4, vec![4], vec![1], 4), + (4, vec![0], vec![4], 4), + (4, vec![0], vec![usize::MAX], usize::MAX), + ] { + assert!(matches!( + FermionicOp::excitation_generator(n, &sources, &targets), + Err(QisError::ModeOutOfBounds { mode, num_modes }) + if mode == invalid && num_modes == n + )); + } +} + +#[test] +fn excitation_generator_preserves_input_order_and_permutation_signs() { + let g = FermionicOp::excitation_generator(5, &[1, 0], &[4, 2]).unwrap(); + let expected = FermionicOp::from_terms( + 5, + [ + ( + vec![Create(4), Create(2), Annihilate(0), Annihilate(1)], + Complex64::new(0.0, 1.0), + ), + ( + vec![Create(1), Create(0), Annihilate(2), Annihilate(4)], + Complex64::new(0.0, -1.0), + ), + ], + ) + .unwrap(); + assert_eq!(g, expected); + assert_eq!(g.num_modes(), 5); + assert_eq!(g.num_terms(), 2); + assert_eq!(g.adjoint(), g); + let ordered = g.normal_order().unwrap(); + assert_eq!( + FermionicOp::excitation_generator(5, &[0, 1], &[4, 2]) + .unwrap() + .normal_order() + .unwrap(), + -&ordered + ); + assert_eq!( + FermionicOp::excitation_generator(5, &[1, 0], &[2, 4]) + .unwrap() + .normal_order() + .unwrap(), + -&ordered + ); + assert_eq!( + FermionicOp::excitation_generator(5, &[0, 1], &[2, 4]) + .unwrap() + .normal_order() + .unwrap(), + ordered + ); + assert_eq!( + FermionicOp::excitation_generator(5, &[4, 2], &[1, 0]).unwrap(), + -&g + ); + + let triple = FermionicOp::excitation_generator(6, &[2, 0, 1], &[4, 5, 3]).unwrap(); + assert_matrix_close( + &fock_matrix(&triple), + &excitation_fock_matrix(6, &[2, 0, 1], &[4, 5, 3]), + ); +} + +#[test] +fn excitation_generator_single_and_double_have_expected_fock_signs() { + let single = FermionicOp::excitation_generator(3, &[0], &[2]).unwrap(); + let mut expected = Array2::::zeros((8, 8)); + let i = Complex64::new(0.0, 1.0); + expected[[4, 1]] = i; + expected[[1, 4]] = -i; + // A spectator in mode 1 flips the single-excitation sign. + expected[[6, 3]] = -i; + expected[[3, 6]] = i; + assert_eq!(fock_matrix(&single), expected); + + let double = FermionicOp::excitation_generator(4, &[0, 1], &[2, 3]).unwrap(); + let mut expected = Array2::::zeros((16, 16)); + expected[[12, 3]] = i; + expected[[3, 12]] = -i; + assert_eq!(fock_matrix(&double), expected); +} + +#[test] +fn excitation_generators_match_fock_action_and_conserve_particle_number() { + let mapper = JordanWignerMapper::new(); + let check = |n, sources: &[usize], targets: &[usize]| { + let op = FermionicOp::excitation_generator(n, sources, targets).unwrap(); + let expected = excitation_fock_matrix(n, sources, targets); + assert_eq!(fock_matrix(&op), expected); + assert_eq!(expected, expected.t().mapv(|c| c.conj())); + assert!(op.is_hermitian(0.0).unwrap()); + assert_eq!(fock_matrix(&op.normal_order().unwrap()), expected); + assert_matrix_close( + &mapper.map_operator(&op).unwrap().to_matrix().unwrap(), + &expected, + ); + for ((row, col), value) in expected.indexed_iter() { + if value.norm() != 0.0 { + assert_eq!(row.count_ones(), col.count_ones()); + } + } + let number = FermionicOp::total_number(n); + assert!( + number + .try_mul(&op) + .unwrap() + .try_sub(&op.try_mul(&number).unwrap()) + .unwrap() + .normal_order() + .unwrap() + .is_zero() + ); + }; + // Every ordered, disjoint single and double excitation for up to five modes. + for n in 2..=5 { + for source in 0..n { + for target in 0..n { + if source == target { + continue; + } + check(n, &[source], &[target]); + for source2 in 0..n { + if [source, target].contains(&source2) { + continue; + } + for target2 in 0..n { + if [source, target, source2].contains(&target2) { + continue; + } + check(n, &[source, source2], &[target, target2]); + } + } + } + } + } + check(6, &[2, 0, 1], &[4, 5, 3]); +} + +#[test] +fn excitation_generator_evolution_matches_independent_fock_exponential() { + let mapper = JordanWignerMapper::new(); + for (n, sources, targets) in [ + (3, vec![0], vec![2]), + (3, vec![2], vec![0]), + (4, vec![0, 1], vec![2, 3]), + (5, vec![1, 0], vec![4, 2]), + ] { + let op = FermionicOp::excitation_generator(n, &sources, &targets).unwrap(); + let mapped = mapper.map_operator(&op).unwrap(); + assert!(mapped.all_terms_commute()); + let k = + excitation_fock_matrix(n, &sources, &targets).mapv(|g| Complex64::new(0.0, -1.0) * g); + for time in [-0.37, 0.0, 0.61] { + assert_matrix_close( + &bound_evolution(mapped.clone(), time), + &matrix_exp(&k, time), + ); + } + } +} + +#[test] +fn composed_excitations_have_independent_parameters_and_correct_statevectors() { + let mapper = JordanWignerMapper::new(); + let occupation = mapper.map_occupation(&[1, 1, 0, 0]).unwrap(); + let mut circuit = BasisEncoding::new(occupation) + .build_circuit("reference") + .unwrap(); + for (sources, targets, name) in [ + (vec![0, 1], vec![2, 3], "theta_double"), + (vec![0], vec![2], "theta_single"), + ] { + let g = FermionicOp::excitation_generator(4, &sources, &targets).unwrap(); + let block = PauliEvolutionAnsatz::new(mapper.map_operator(&g).unwrap()) + .unwrap() + .with_time_param_name(name) + .with_strategy(EvolutionStrategy::Exact) + .build_circuit("excitation") + .unwrap(); + assert_eq!(block.symbols().len(), 1); + circuit.compose(&block, None).unwrap(); + } + assert_eq!(circuit.symbols().len(), 2); + assert!(circuit.symbols().contains("theta_double")); + assert!(circuit.symbols().contains("theta_single")); + let number = mapper.map_operator(&FermionicOp::total_number(4)).unwrap(); + for (double, single) in [ + (0.0_f64, 0.0_f64), + (0.37, 0.0), + (0.0, -0.21), + (0.37, -0.21), + (-0.43, 0.19), + ] { + // Bind in two stages to exercise independence and partial binding. + let partial = circuit + .assign_parameters(&Some(HashMap::from([("theta_double", double)]))) + .unwrap(); + assert_eq!(partial.symbols().len(), 1); + assert!(partial.symbols().contains("theta_single")); + let bound = partial + .assign_parameters(&Some(HashMap::from([("theta_single", single)]))) + .unwrap(); + assert!(bound.symbols().is_empty()); + let state = Statevector::from_circuit(&bound).unwrap(); + let mut expected = vec![Complex64::new(0.0, 0.0); 16]; + expected[3] = (double.cos() * single.cos()).into(); + expected[6] = (-double.cos() * single.sin()).into(); + expected[12] = (double.sin() * single.cos()).into(); + expected[9] = (-double.sin() * single.sin()).into(); + for (actual, expected) in state.data().iter().zip(expected) { + assert!((*actual - expected).norm() < 1e-11); + } + assert!((number.expectation_statevector(&state).unwrap() - 2.0).abs() < 1e-11); + } +} + +fn matrix_exp(matrix: &Array2, scale: f64) -> Array2 { + let mut term = Array2::::eye(matrix.nrows()); + let mut result = term.clone(); + for order in 1..=48 { + term = term.dot(matrix).mapv(|c| c * scale / order as f64); + result += &term; + } + result +} + +fn bound_evolution(h: Hamiltonian, time: f64) -> Array2 { + let circuit = PauliEvolutionAnsatz::new(h) + .unwrap() + .build_circuit("fermion") + .unwrap(); + let bindings = HashMap::from([("fermion_t", time)]); + let circuit = circuit.assign_parameters(&Some(bindings)).unwrap(); + circuit_to_matrix(&circuit, None).unwrap() +} + +fn pauli_terms(h: &Hamiltonian) -> BTreeMap { + assert!(h.terms.iter().all(|(p, _)| p.phase == Phase::Plus)); + h.terms.iter().map(|(p, c)| (p.to_string(), *c)).collect() +} + +#[test] +fn canonical_anticommutation_relations() { + for p in 0..4 { + for q in 0..4 { + let a = FermionicOp::annihilation(4, p).unwrap(); + let b = FermionicOp::annihilation(4, q).unwrap(); + let bd = b.adjoint(); + let anticommutator = (a.try_mul(&bd).unwrap() + bd.try_mul(&a).unwrap()).unwrap(); + let expected = if p == q { + FermionicOp::identity(4) + } else { + FermionicOp::new(4) + }; + assert_eq!(anticommutator.normal_order().unwrap(), expected); + let annihilators = (a.try_mul(&b).unwrap() + b.try_mul(&a).unwrap()).unwrap(); + assert!(annihilators.normal_order().unwrap().is_zero()); + let creators = + (a.adjoint().try_mul(&bd).unwrap() + bd.try_mul(&a.adjoint()).unwrap()).unwrap(); + assert!(creators.normal_order().unwrap().is_zero()); + } + } +} + +#[test] +fn contractions_and_normal_order_preserve_remaining_branches() { + let a = FermionicOp::annihilation(3, 0).unwrap(); + let ad = a.adjoint(); + let hole = (&a * &ad).unwrap(); + assert_eq!( + hole.normal_order().unwrap(), + (FermionicOp::identity(3) - FermionicOp::number(3, 0).unwrap()).unwrap() + ); + assert_eq!((&hole * &a).unwrap().normal_order().unwrap(), a); + assert_eq!((&ad * &hole).unwrap().normal_order().unwrap(), ad); + let n = FermionicOp::number(3, 0).unwrap(); + assert_eq!((&n * &n).unwrap().normal_order().unwrap(), n); + let input = FermionicOp::from_term( + 3, + [Annihilate(0), Create(0), Create(1), Annihilate(2)], + 1.0.into(), + ) + .unwrap(); + let expected = FermionicOp::from_terms( + 3, + [ + (vec![Create(1), Annihilate(2)], 1.0.into()), + ( + vec![Create(0), Create(1), Annihilate(0), Annihilate(2)], + 1.0.into(), + ), + ], + ) + .unwrap(); + assert_eq!(input.normal_order().unwrap(), expected); + assert_eq!(expected.normal_order().unwrap(), expected); +} + +#[test] +fn long_single_branch_normal_order_preserves_sign_and_nilpotency() { + for n in [399, 400, 799, 800] { + let word: Vec<_> = (0..n).rev().map(Create).collect(); + let input = FermionicOp::from_term(n, word.clone(), 1.0.into()).unwrap(); + let sign = if (n * (n - 1) / 2) % 2 == 0 { + 1.0 + } else { + -1.0 + }; + let expected = FermionicOp::from_term(n, (0..n).map(Create), sign.into()).unwrap(); + assert_eq!(input.normal_order().unwrap(), expected); + assert_eq!(expected.normal_order().unwrap(), expected); + + let mut repeated = word; + repeated.push(Create(n - 1)); + assert!( + FermionicOp::from_term(n, repeated, 1.0.into()) + .unwrap() + .normal_order() + .unwrap() + .is_zero() + ); + } +} + +#[test] +fn multiplication_keeps_order_and_adjoint_reverses_it() { + let a = + FermionicOp::from_term(3, [Create(2), Annihilate(0)], Complex64::new(2.0, 3.0)).unwrap(); + let expected = + FermionicOp::from_term(3, [Create(0), Annihilate(2)], Complex64::new(2.0, -3.0)).unwrap(); + assert_eq!(a.adjoint(), expected); + assert_eq!(a.adjoint().adjoint(), a); + let b = FermionicOp::creation(3, 1).unwrap(); + let product = (&a * &b).unwrap(); + assert_eq!( + product.terms().next().unwrap().0, + [Create(2), Annihilate(0), Create(1)] + ); + assert_eq!(product.adjoint(), (&b.adjoint() * &a.adjoint()).unwrap()); + assert_eq!(Create(1).adjoint(), Annihilate(1)); +} + +#[test] +fn exact_merging_explicit_chopping_and_structural_equality() { + let mut op = FermionicOp::from_terms( + 4, + [ + (vec![Create(1)], 0.5.into()), + (vec![Create(1)], (-0.5).into()), + (vec![Annihilate(3)], Complex64::new(1e-14, 0.0)), + ], + ) + .unwrap(); + assert_eq!(op.num_terms(), 1); + assert_eq!(op.num_modes(), 4); + op.chop(0.0).unwrap(); + assert_eq!(op.num_terms(), 1); + op.chop(1e-14).unwrap(); + assert_eq!(op, FermionicOp::new(4)); + let input = FermionicOp::from_term(2, [Create(1), Create(0)], 1.0.into()).unwrap(); + let ordered = FermionicOp::from_term(2, [Create(0), Create(1)], (-1.0).into()).unwrap(); + assert_ne!(input, ordered); + assert!(input.equiv(&ordered, 0.0).unwrap()); + assert_ne!(FermionicOp::new(1), FermionicOp::new(2)); +} + +#[test] +fn fallible_owned_and_borrowed_arithmetic() { + let a = FermionicOp::creation(2, 0).unwrap(); + let b = FermionicOp::annihilation(2, 1).unwrap(); + let sum = a.try_add(&b).unwrap(); + assert_eq!((&a + &b).unwrap(), sum); + assert_eq!((a.clone() + &b).unwrap(), sum); + assert_eq!((&a + b.clone()).unwrap(), sum); + assert_eq!((a.clone() + b.clone()).unwrap(), sum); + let difference = a.try_sub(&b).unwrap(); + assert_eq!((&a - &b).unwrap(), difference); + assert_eq!((a.clone() - &b).unwrap(), difference); + assert_eq!((&a - b.clone()).unwrap(), difference); + assert_eq!((a.clone() - b.clone()).unwrap(), difference); + let product = a.try_mul(&b).unwrap(); + assert_eq!((&a * &b).unwrap(), product); + assert_eq!((a.clone() * &b).unwrap(), product); + assert_eq!((&a * b.clone()).unwrap(), product); + assert_eq!((a.clone() * b.clone()).unwrap(), product); + for factor in [Complex64::new(2.0, 0.0), Complex64::new(0.0, 1.0)] { + let scaled = a.scaled(factor).unwrap(); + assert_eq!((&a * factor).unwrap(), scaled); + assert_eq!((a.clone() * factor).unwrap(), scaled); + assert_eq!((factor * &a).unwrap(), scaled); + assert_eq!((factor * a.clone()).unwrap(), scaled); + } + assert_eq!((&a * 2.0).unwrap(), (2.0 * &a).unwrap()); + assert_eq!((a.clone() * 2.0).unwrap(), (2.0 * a.clone()).unwrap()); + assert_eq!(-&a, -a.clone()); + assert!((&a + -&a).unwrap().is_zero()); +} + +#[test] +fn invalid_inputs_and_failed_mutation_do_not_corrupt_operators() { + assert!(matches!( + FermionicOp::creation(0, 0), + Err(QisError::ModeOutOfBounds { + mode: 0, + num_modes: 0 + }) + )); + assert!(FermionicOp::annihilation(2, 2).is_err()); + assert!(FermionicOp::number(1, usize::MAX).is_err()); + assert!(FermionicOp::from_term(1, [Create(2)], 0.0.into()).is_err()); + let a = FermionicOp::creation(1, 0).unwrap(); + let b = FermionicOp::new(2); + for result in [&a + &b, &a - &b, &a * &b] { + assert!(matches!( + result, + Err(QisError::ModeMismatch { + expected: 1, + actual: 2 + }) + )); + } + assert!(a.equiv(&b, 0.0).is_err()); + for c in [ + Complex64::new(f64::NAN, 0.0), + Complex64::new(0.0, f64::INFINITY), + ] { + assert!(matches!( + FermionicOp::from_term(1, [], c), + Err(QisError::NonFiniteCoefficient { .. }) + )); + assert!(FermionicOp::new(1).scaled(c).is_err()); + } + let mut large = FermionicOp::from_term(1, [Create(0)], f64::MAX.into()).unwrap(); + let saved = large.clone(); + assert!(large.add_term([Create(0)], f64::MAX.into()).is_err()); + assert_eq!(large, saved); + assert!(large.add_term([Create(1)], 1.0.into()).is_err()); + assert_eq!(large, saved); + assert!(large.scaled(2.0.into()).is_err()); + for tolerance in [-1.0, f64::NAN, f64::INFINITY] { + assert!(large.chop(tolerance).is_err()); + assert_eq!(large, saved); + assert!(large.is_hermitian(tolerance).is_err()); + assert!(large.equiv(&saved, tolerance).is_err()); + } +} + +#[test] +fn hermiticity_checks_use_fermionic_algebra() { + let hole = FermionicOp::from_term(2, [Annihilate(0), Create(0)], 1.0.into()).unwrap(); + assert!(hole.is_hermitian(0.0).unwrap()); + let t = FermionicOp::from_term(2, [Create(1), Annihilate(0)], 1.0.into()).unwrap(); + assert!(!t.is_hermitian(0.0).unwrap()); + let generator = (&t - t.adjoint()) + .unwrap() + .scaled(Complex64::new(0.0, 1.0)) + .unwrap(); + assert!(generator.is_hermitian(0.0).unwrap()); + assert!( + !FermionicOp::identity(0) + .scaled(Complex64::new(0.0, 1.0)) + .unwrap() + .is_hermitian(0.0) + .unwrap() + ); +} + +#[test] +fn ladder_number_hole_and_identity_mappings() { + let mapper = JordanWignerMapper::new(); + let creation = mapper + .map_operator(&FermionicOp::creation(4, 2).unwrap()) + .unwrap(); + assert_eq!( + pauli_terms(&creation), + BTreeMap::from([ + ("+IXZZ".into(), 0.5.into()), + ("+IYZZ".into(), Complex64::new(0.0, -0.5)), + ]) + ); + let annihilation = mapper + .map_operator(&FermionicOp::annihilation(4, 2).unwrap()) + .unwrap(); + assert_eq!( + pauli_terms(&annihilation), + BTreeMap::from([ + ("+IXZZ".into(), 0.5.into()), + ("+IYZZ".into(), Complex64::new(0.0, 0.5)), + ]) + ); + assert_eq!( + pauli_terms( + &mapper + .map_operator(&FermionicOp::number(4, 2).unwrap()) + .unwrap() + ), + BTreeMap::from([ + ("+IIII".into(), 0.5.into()), + ("+IZII".into(), (-0.5).into()), + ]) + ); + let hole = FermionicOp::from_term(4, [Annihilate(2), Create(2)], 1.0.into()).unwrap(); + assert_eq!( + pauli_terms(&mapper.map_operator(&hole).unwrap()), + BTreeMap::from([("+IIII".into(), 0.5.into()), ("+IZII".into(), 0.5.into()),]) + ); + assert_eq!( + pauli_terms(&mapper.map_operator(&FermionicOp::identity(4)).unwrap()), + BTreeMap::from([("+IIII".into(), 1.0.into())]) + ); +} + +#[test] +fn zero_scalar_and_wide_mode_spaces_are_preserved() { + let mapper = JordanWignerMapper::new(); + for n in [0, 1, 4, 129] { + let zero = mapper.map_operator(&FermionicOp::new(n)).unwrap(); + assert_eq!(zero.num_qubits, n); + assert!(zero.terms.is_empty()); + assert_eq!(mapper.num_qubits(n), n); + } + let scalar = FermionicOp::identity(0) + .scaled(Complex64::new(2.0, 3.0)) + .unwrap(); + assert_eq!( + mapper.map_operator(&scalar).unwrap().to_matrix().unwrap()[[0, 0]], + Complex64::new(2.0, 3.0) + ); + assert!(FermionicOp::total_number(0).is_zero()); + let wide = mapper + .map_operator(&FermionicOp::creation(129, 128).unwrap()) + .unwrap(); + assert_eq!(wide.num_qubits, 129); + assert_eq!(wide.terms.len(), 2); + assert!( + wide.terms + .iter() + .all(|(p, _)| (0..128).all(|j| p.get_pauli(j) == crate::qis::Pauli::Z)) + ); + let nilpotent = FermionicOp::from_term(4, [Create(0), Create(0)], 1.0.into()).unwrap(); + let mapped = mapper.map_operator(&nilpotent).unwrap(); + assert!(mapped.terms.is_empty()); + assert_eq!(mapped.num_qubits, 4); +} + +fn check_word(n: usize, word: Vec) { + let op = FermionicOp::from_term(n, word, 1.0.into()).unwrap(); + let normal = op.normal_order().unwrap(); + assert_eq!(normal.normal_order().unwrap(), normal); + let expected = fock_matrix(&op); + assert_matrix_close(&fock_matrix(&normal), &expected); + let mapper = JordanWignerMapper::new(); + assert_matrix_close( + &mapper.map_operator(&op).unwrap().to_matrix().unwrap(), + &expected, + ); + assert_matrix_close( + &mapper.map_operator(&normal).unwrap().to_matrix().unwrap(), + &expected, + ); +} + +#[test] +fn exhaustive_short_products_match_independent_fock_matrices() { + for n in 1usize..=3 { + for length in 0u32..=4 { + for mut encoded in 0..(2 * n).pow(length) { + let word = (0..length) + .map(|_| { + let symbol = encoded % (2 * n); + encoded /= 2 * n; + if symbol < n { + Create(symbol) + } else { + Annihilate(symbol - n) + } + }) + .collect(); + check_word(n, word); + } + } + } +} + +#[test] +fn fixed_seed_longer_products_match_independent_fock_matrices() { + let mut seed = 20261009u64; + for sample in 0..256 { + let word = (0..sample % 9) + .map(|_| { + seed = seed.wrapping_mul(6364136223846793005).wrapping_add(1); + let symbol = (seed >> 32) as usize % 8; + if symbol < 4 { + Create(symbol) + } else { + Annihilate(symbol - 4) + } + }) + .collect(); + check_word(4, word); + } +} + +#[test] +fn mapping_preserves_linear_combinations_products_and_adjoint() { + let mapper = JordanWignerMapper::new(); + let a = FermionicOp::from_terms( + 3, + [ + (vec![Annihilate(1), Create(2)], Complex64::new(1.0, 2.0)), + (vec![], Complex64::new(-2.0, 1.0)), + ], + ) + .unwrap(); + let b = FermionicOp::creation(3, 0).unwrap(); + let ma = mapper.map_operator(&a).unwrap().to_matrix().unwrap(); + let mb = mapper.map_operator(&b).unwrap().to_matrix().unwrap(); + assert_matrix_close(&ma, &fock_matrix(&a)); + assert_matrix_close( + &mapper + .map_operator(&(&a + &b).unwrap()) + .unwrap() + .to_matrix() + .unwrap(), + &(&ma + &mb), + ); + assert_matrix_close( + &mapper + .map_operator(&(&a * &b).unwrap()) + .unwrap() + .to_matrix() + .unwrap(), + &ma.dot(&mb), + ); + assert_matrix_close( + &mapper + .map_operator(&a.adjoint()) + .unwrap() + .to_matrix() + .unwrap(), + &ma.t().mapv(|c| c.conj()), + ); +} + +#[test] +fn occupation_encoding_prepares_basis_states_and_number_observables() { + let mapper = JordanWignerMapper::new(); + let operator_mapper: &dyn FermionToQubitMapper = &mapper; + let state_mapper: &dyn OccupationStateMapper = &mapper; + for n in 1..=4 { + let number = operator_mapper + .map_operator(&FermionicOp::total_number(n)) + .unwrap(); + for state in 0usize..1 << n { + let occupation: Vec<_> = (0..n).map(|p| ((state >> p) & 1) as u8).collect(); + let bits = state_mapper.map_occupation(&occupation).unwrap(); + let circuit = BasisEncoding::new(bits) + .build_circuit("occupation") + .unwrap(); + let sv = Statevector::from_circuit(&circuit).unwrap(); + assert_eq!(sv.probabilities()[state], 1.0); + assert_eq!( + number.expectation_statevector(&sv).unwrap(), + state.count_ones() as f64 + ); + assert_eq!(number.variance_statevector(&sv).unwrap(), 0.0); + let dm = DensityMatrix::from_circuit(&circuit).unwrap(); + assert_eq!( + number.expectation_density_matrix(&dm).unwrap(), + state.count_ones() as f64 + ); + for (p, occupation_value) in occupation.iter().enumerate() { + let np = operator_mapper + .map_operator(&FermionicOp::number(n, p).unwrap()) + .unwrap(); + assert_eq!( + np.expectation_statevector(&sv).unwrap(), + *occupation_value as f64 + ); + } + } + } + assert_eq!(mapper.map_occupation(&[]).unwrap(), Vec::::new()); + for invalid in [2, 255] { + assert!( + matches!(mapper.map_occupation(&[1, invalid]), Err(QisError::InvalidOccupation { index: 1, value }) if value == invalid) + ); + } +} + +#[test] +fn single_and_double_excitation_evolution_matches_fock_exponentials() { + let mapper = JordanWignerMapper::new(); + for (n, word) in [ + (3, vec![Create(2), Annihilate(0)]), + (4, vec![Create(2), Create(3), Annihilate(1), Annihilate(0)]), + ] { + let t = FermionicOp::from_term(n, word, 1.0.into()).unwrap(); + let k = (&t - t.adjoint()).unwrap(); + let h = k.scaled(Complex64::new(0.0, 1.0)).unwrap(); + let mapped = mapper.map_operator(&h).unwrap(); + assert!(mapped.all_terms_commute()); + let time = 0.37; + let expected = matrix_exp(&fock_matrix(&k), time); + assert_matrix_close(&bound_evolution(mapped.clone(), time), &expected); + let numeric = mapped + .to_evolution_circuit(time, 1, TrotterMode::FirstOrder) + .unwrap(); + assert_matrix_close(&circuit_to_matrix(&numeric, None).unwrap(), &expected); + let particle_number = fock_matrix(&FermionicOp::total_number(n)); + let generator = fock_matrix(&k); + assert_matrix_close( + &particle_number.dot(&generator), + &generator.dot(&particle_number), + ); + } +} + +#[test] +fn reversing_double_excitation_annihilators_changes_the_sign() { + let a = FermionicOp::from_term( + 4, + [Create(2), Create(3), Annihilate(1), Annihilate(0)], + 1.0.into(), + ) + .unwrap(); + let b = FermionicOp::from_term( + 4, + [Create(2), Create(3), Annihilate(0), Annihilate(1)], + 1.0.into(), + ) + .unwrap(); + assert!((&a + &b).unwrap().normal_order().unwrap().is_zero()); + assert_eq!(fock_matrix(&a)[[12, 3]], Complex64::new(1.0, 0.0)); + assert_matrix_close(&fock_matrix(&b), &fock_matrix(&a).mapv(|c| -c)); +} + +#[test] +fn small_coefficients_survive_mapping_and_both_evolution_apis() { + let mapper = JordanWignerMapper::new(); + let t = FermionicOp::creation(1, 0).unwrap(); + let k = (&t - t.adjoint()).unwrap(); + let tiny = k.scaled(Complex64::new(0.0, 1e-12)).unwrap(); + let mapped = mapper.map_operator(&tiny).unwrap(); + assert_eq!(mapped.terms.len(), 1); + assert_eq!(mapped.terms[0].1, Complex64::new(1e-12, 0.0)); + assert_eq!( + PauliEvolutionAnsatz::new(mapped.clone()) + .unwrap() + .evolution_info() + .num_terms, + 1 + ); + // Large time makes loss of a small coefficient observable in the matrix. + let numeric = mapped + .to_evolution_circuit(1e12, 1, TrotterMode::FirstOrder) + .unwrap(); + let expected = matrix_exp(&fock_matrix(&k), 1.0); + assert_matrix_close(&circuit_to_matrix(&numeric, None).unwrap(), &expected); + assert_matrix_close(&bound_evolution(mapped, 1e12), &expected); +} + +#[test] +fn identity_evolution_preserves_global_phase_and_nonhermitian_mapping_is_allowed() { + let mapper = JordanWignerMapper::new(); + let h = mapper + .map_operator(&FermionicOp::identity(2).scaled(0.3.into()).unwrap()) + .unwrap(); + let expected = Array2::::eye(4).mapv(|c| c * Complex64::from_polar(1.0, -0.15)); + assert_matrix_close(&bound_evolution(h, 0.5), &expected); + let creation = mapper + .map_operator(&FermionicOp::creation(2, 0).unwrap()) + .unwrap(); + assert!(matches!( + PauliEvolutionAnsatz::new(creation), + Err(QisError::NotHermitian) + )); + assert!(PauliEvolutionAnsatz::new(mapper.map_operator(&FermionicOp::new(2)).unwrap()).is_err()); +} + +#[test] +fn mapping_reports_coefficient_overflow_after_pauli_merging() { + let op = FermionicOp::from_terms( + 1, + [ + (vec![], f64::MAX.into()), + (vec![Create(0), Annihilate(0)], f64::MAX.into()), + (vec![Annihilate(0), Create(0)], f64::MAX.into()), + ], + ) + .unwrap(); + assert!(matches!( + JordanWignerMapper::new().map_operator(&op), + Err(QisError::NonFiniteCoefficient { .. }) + )); +} + +#[test] +fn mapping_applies_subnormal_weights_after_merging_the_unit_product() { + let tiny = f64::from_bits(2); + let half = f64::from_bits(1); + let mapper = JordanWignerMapper::new(); + for coefficient in [ + Complex64::new(tiny, 0.0), + Complex64::new(0.0, -tiny), + Complex64::new(tiny, -tiny), + ] { + for (word, z_sign) in [ + (vec![Create(0), Annihilate(0)], -1.0), + (vec![Annihilate(0), Create(0)], 1.0), + ( + vec![Create(0), Annihilate(0), Create(0), Annihilate(0)], + -1.0, + ), + ] { + let input = FermionicOp::from_term(1, word, coefficient).unwrap(); + let mapped = mapper.map_operator(&input).unwrap(); + let expected_half = Complex64::new( + if coefficient.re == 0.0 { 0.0 } else { half }, + if coefficient.im == 0.0 { 0.0 } else { -half }, + ); + // Compare coefficients exactly: a matrix tolerance would hide subnormals. + assert_eq!( + pauli_terms(&mapped), + BTreeMap::from([ + ("+I".into(), expected_half), + ("+Z".into(), expected_half * z_sign), + ]) + ); + } + let identity = FermionicOp::identity(1).scaled(coefficient).unwrap(); + assert_eq!( + pauli_terms(&mapper.map_operator(&identity).unwrap()), + BTreeMap::from([("+I".into(), coefficient)]) + ); + } +} + +#[test] +fn normal_order_and_arithmetic_report_coefficient_overflow() { + let op = FermionicOp::from_terms( + 1, + [ + (vec![Create(0)], f64::MAX.into()), + (vec![Create(0), Annihilate(0), Create(0)], f64::MAX.into()), + ], + ) + .unwrap(); + assert!(matches!( + op.normal_order(), + Err(QisError::NonFiniteCoefficient { .. }) + )); + let a = FermionicOp::from_term(1, [Create(0)], f64::MAX.into()).unwrap(); + assert!((&a + &a).is_err()); + let twice_identity = FermionicOp::identity(1).scaled(2.0.into()).unwrap(); + assert!((&a * twice_identity).is_err()); + let tiny = FermionicOp::creation(1, 0) + .unwrap() + .scaled(1e-310.into()) + .unwrap(); + let mut retained = tiny.clone(); + retained.chop(0.0).unwrap(); + assert_eq!(retained, tiny); + assert_eq!( + JordanWignerMapper::new() + .map_operator(&tiny) + .unwrap() + .terms + .len(), + 2 + ); +} + +#[test] +fn noncommuting_excitation_sum_uses_the_same_numeric_and_symbolic_product_formula() { + let first = FermionicOp::from_term(3, [Create(1), Annihilate(0)], 1.0.into()).unwrap(); + let second = FermionicOp::from_term(3, [Create(2), Annihilate(1)], 1.0.into()).unwrap(); + let t = (&first + &second).unwrap(); + let k = (&t - t.adjoint()).unwrap(); + let h = JordanWignerMapper::new() + .map_operator(&k.scaled(Complex64::new(0.0, 1.0)).unwrap()) + .unwrap(); + assert!(!h.all_terms_commute()); + assert!( + PauliEvolutionAnsatz::new(h.clone()) + .unwrap() + .with_strategy(EvolutionStrategy::Exact) + .build_circuit("exact") + .is_err() + ); + let time = 0.2; + let numeric = h + .to_evolution_circuit(time, 8, TrotterMode::SecondOrder) + .unwrap(); + let numeric = circuit_to_matrix(&numeric, None).unwrap(); + let circuit = PauliEvolutionAnsatz::new(h) + .unwrap() + .with_strategy(EvolutionStrategy::Trotter { + mode: TrotterMode::SecondOrder, + steps: 8, + }) + .build_circuit("sum") + .unwrap(); + let circuit = circuit + .assign_parameters(&Some(HashMap::from([("sum_t", time)]))) + .unwrap(); + assert_matrix_close(&numeric, &circuit_to_matrix(&circuit, None).unwrap()); + let expected = matrix_exp(&fock_matrix(&k), time); + assert!((&numeric - &expected).iter().all(|c| c.norm() < 1e-3)); +} + +#[test] +fn all_hamiltonian_consumers_validate_original_and_merged_inputs() { + let sv = Statevector::new(1); + let dm = DensityMatrix::new(1); + let x: PauliString = "X".parse().unwrap(); + let mut malformed = x.clone(); + malformed.z.clear(); + let invalid = [ + Hamiltonian { + num_qubits: 1, + terms: vec![(PauliString::new(2), 1.0.into())], + }, + Hamiltonian { + num_qubits: 1, + terms: vec![(malformed, 1.0.into())], + }, + Hamiltonian { + num_qubits: 1, + terms: vec![(x.clone(), f64::NAN.into())], + }, + Hamiltonian { + num_qubits: 1, + terms: vec![ + (x.clone(), f64::INFINITY.into()), + (x.clone(), f64::NEG_INFINITY.into()), + ], + }, + Hamiltonian { + num_qubits: 1, + terms: vec![(x.clone(), f64::MAX.into()), (x.clone(), f64::MAX.into())], + }, + ]; + for h in invalid { + assert!(h.expectation_statevector(&sv).is_err()); + assert!(h.expectation_density_matrix(&dm).is_err()); + assert!(h.expectation_probs(&[]).is_err()); + assert!(h.variance_statevector(&sv).is_err()); + assert!( + h.to_evolution_circuit(1.0, 1, TrotterMode::FirstOrder) + .is_err() + ); + assert!( + h.to_trotter_circuit(1.0, 1, TrotterMode::FirstOrder) + .is_err() + ); + assert!(PauliEvolutionAnsatz::new(h).is_err()); + } + let zero = Hamiltonian::from_list(vec![(x.clone(), 1.0.into()), (x, (-1.0).into())]).unwrap(); + assert_eq!(zero.expectation_statevector(&sv).unwrap(), 0.0); + assert_eq!(zero.expectation_probs(&[]).unwrap(), 0.0); + for mode in [ + TrotterMode::FirstOrder, + TrotterMode::SecondOrder, + TrotterMode::Randomized(7), + ] { + for circuit in [ + zero.to_evolution_circuit(1.0, 3, mode).unwrap(), + zero.to_trotter_circuit(1.0, 3, mode).unwrap(), + ] { + assert_eq!(circuit.num_qubits(), 1); + assert_eq!(circuit_to_matrix(&circuit, None).unwrap(), Array2::eye(2)); + } + } + assert!(PauliEvolutionAnsatz::new(zero).is_err()); +} + +#[test] +fn hamiltonian_consumers_preserve_coefficient_error_categories() { + let sv = Statevector::new(1); + let dm = DensityMatrix::new(1); + let x: PauliString = "X".parse().unwrap(); + for coefficients in [ + vec![f64::NAN], + vec![f64::INFINITY], + vec![f64::MAX, f64::MAX], + ] { + let h = Hamiltonian::from_list( + coefficients + .into_iter() + .map(|c| (x.clone(), c.into())) + .collect(), + ) + .unwrap(); + for result in [ + h.expectation_statevector(&sv), + h.expectation_density_matrix(&dm), + h.expectation_probs(&[]), + h.variance_statevector(&sv), + ] { + assert!(matches!(result, Err(QisError::InvalidParameterValue(_)))); + } + for result in [ + h.to_evolution_circuit(1.0, 1, TrotterMode::FirstOrder), + h.to_trotter_circuit(1.0, 1, TrotterMode::FirstOrder), + ] { + assert!(matches!(result, Err(QisError::UnsupportedOperation(_)))); + } + } +} diff --git a/crates/cqlib-core/src/qis/hamiltonian.rs b/crates/cqlib-core/src/qis/hamiltonian.rs index 015ecfbd..f82d67e2 100644 --- a/crates/cqlib-core/src/qis/hamiltonian.rs +++ b/crates/cqlib-core/src/qis/hamiltonian.rs @@ -13,10 +13,14 @@ //! Hamiltonian representation and operations. //! //! This module provides the [`Hamiltonian`] struct, which represents a quantum -//! observable as a linear combination of Pauli strings (tensor products of Pauli operators) +//! operator as a linear combination of Pauli strings (tensor products of Pauli operators) //! with complex coefficients. It is primarily used for defining system energies, //! observables for expectation value calculations, and operators for time evolution. +//! +//! General algebra permits non-Hermitian operators. Real observables and unitary +//! evolution validate Hermiticity at their consumption boundaries. +use crate::qis::error::validate_coefficient; use crate::qis::pauli::PauliString; use crate::qis::{Phase, QisError}; use ndarray::Array2; @@ -29,6 +33,10 @@ use std::ops::Add; /// A `Hamiltonian` is essentially a sparse representation of a $2^N \times 2^N$ /// matrix, expressed as $H = \sum_k c_k P_k$, where $c_k$ is a complex coefficient /// and $P_k$ is an $N$-qubit Pauli string. +/// +/// Despite its name, this type also represents non-Hermitian operators. Checked +/// algebra returns phase-normalized, merged terms without approximate pruning. +/// [`PartialEq`] compares stored structure; use [`Self::equiv`] for algebraic comparison. #[derive(Debug, Clone, PartialEq)] pub struct Hamiltonian { /// The number of qubits this operator acts on. @@ -51,6 +59,56 @@ impl fmt::Display for Hamiltonian { } impl Hamiltonian { + /// Returns a validated Pauli sum with phases absorbed and exact zeros removed. + /// Retains the declared qubit count, including for zero and scalar operators. + pub fn normalized(&self) -> Result { + Self::from_terms(self.num_qubits, self.terms.iter().cloned()) + } + + fn validate(&self) -> Result<(), QisError> { + for (pauli, coefficient) in &self.terms { + if pauli.num_qubits != self.num_qubits { + return Err(QisError::QubitMismatch { + expected: self.num_qubits, + actual: pauli.num_qubits, + }); + } + for actual in [pauli.x.len(), pauli.z.len()] { + if actual != self.num_qubits { + return Err(QisError::DimensionMismatch { + expected: self.num_qubits, + actual, + }); + } + } + validate_coefficient(*coefficient)?; + } + Ok(()) + } + + /// Shared contract for real-valued observables and Hamiltonian evolution. + /// After merging and absorbing Pauli phases, imaginary parts with magnitude + /// greater than `1e-10` are rejected. Accepted imaginary roundoff residuals + /// are projected to zero, and terms with exactly zero real part are removed. + /// Every nonzero real coefficient is retained, however small. + pub fn normalized_hermitian(&self) -> Result { + let mut normalized = self.normalized()?; + if normalized + .terms + .iter() + .any(|(_, coefficient)| coefficient.im.abs() > 1e-10) + { + return Err(QisError::NotHermitian); + } + for (_, coefficient) in &mut normalized.terms { + coefficient.im = 0.0; + } + normalized + .terms + .retain(|(_, coefficient)| coefficient.re != 0.0); + Ok(normalized) + } + /// Creates a new empty Hamiltonian. /// /// The resulting Hamiltonian represents the zero operator for the given @@ -82,25 +140,49 @@ impl Hamiltonian { /// /// # Errors /// - /// Returns [`QisError::QubitMismatch`] when a term's qubit count differs - /// from [`Hamiltonian::num_qubits`]. + /// Validates term widths, bit-vector lengths and finite coefficients. Returns + /// [`QisError::UnsupportedOperation`] if the dense allocation size overflows + /// or its reservation fails. Matrix-entry overflow is also reported. pub fn to_matrix(&self) -> Result, QisError> { - for (pauli, _) in &self.terms { - if pauli.num_qubits != self.num_qubits { - return Err(QisError::QubitMismatch { - expected: self.num_qubits, - actual: pauli.num_qubits, - }); - } + let normalized = self.normalized()?; + let size_error = || { + QisError::UnsupportedOperation( + "dense Pauli matrix exceeds addressable allocation size".into(), + ) + }; + let shift = u32::try_from(self.num_qubits).map_err(|_| size_error())?; + let dim = 1usize.checked_shl(shift).ok_or_else(size_error)?; + let len = dim.checked_mul(dim).ok_or_else(size_error)?; + let bytes = len + .checked_mul(std::mem::size_of::()) + .ok_or_else(size_error)?; + if bytes > isize::MAX as usize { + return Err(size_error()); } - - let dim = 1usize << self.num_qubits; - let mut matrix = Array2::from_elem((dim, dim), Complex64::new(0.0, 0.0)); - for (pauli, coefficient) in &self.terms { - let term = pauli.to_matrix().mapv(|value| value * *coefficient); - matrix += &term; + let mut data = Vec::new(); + data.try_reserve_exact(len).map_err(|error| { + QisError::UnsupportedOperation(format!("cannot allocate dense Pauli matrix: {error}")) + })?; + data.resize(len, Complex64::new(0.0, 0.0)); + // P|column> flips X bits; Z parity is evaluated before those flips. + for (pauli, coefficient) in &normalized.terms { + let x = pauli.x_mask(); + let z = pauli.z_mask(); + let factor = *coefficient * pauli.y_phase(); + validate_coefficient(factor)?; + for column in 0..dim { + let sign = if (column & z).count_ones() & 1 == 0 { + 1.0 + } else { + -1.0 + }; + let entry = &mut data[(column ^ x) * dim + column]; + *entry += factor * sign; + validate_coefficient(*entry)?; + } } - Ok(matrix) + Array2::from_shape_vec((dim, dim), data) + .map_err(|error| QisError::UnsupportedOperation(error.to_string())) } /// Creates a Hamiltonian from a list of Pauli strings and their corresponding coefficients. @@ -127,6 +209,236 @@ impl From for Hamiltonian { } } +impl Hamiltonian { + /// The identity on the declared register; at zero qubits this is scalar one. + pub fn identity(num_qubits: usize) -> Self { + Self::from_pauli(PauliString::new(num_qubits)) + } + + /// Constructs a validated, phase-normalized sum, merging duplicate terms. + /// Only exact zeros are removed, and empty input retains `num_qubits`. + /// + /// Invalid widths, bit-vector lengths, nonfinite coefficients and intermediate + /// coefficient overflow return an error. No approximate truncation occurs. + pub fn from_terms( + num_qubits: usize, + terms: impl IntoIterator, + ) -> Result { + let mut result = Self { + num_qubits, + terms: terms.into_iter().collect(), + }; + result.validate()?; + result.merge_terms(0.0, true)?; + Ok(result) + } + + fn check_width(&self, other: &Self) -> Result<(), QisError> { + if self.num_qubits != other.num_qubits { + return Err(QisError::QubitMismatch { + expected: self.num_qubits, + actual: other.num_qubits, + }); + } + Ok(()) + } + + fn check_materialized_width(width: usize) -> Result<(), QisError> { + // BitSlice reserves metadata bits, so its length bound is below usize::MAX. + // Zero sums need no bit vectors and may retain larger declared widths. + if width > bitvec::slice::BitSlice::::MAX_BITS { + return Err(QisError::UnsupportedOperation( + "Pauli register width exceeds the BitVec representation limit".into(), + )); + } + Ok(()) + } + + /// Adds two operators on the same register, without modifying either input. + /// Unlike legacy `Add`, this validates inputs and merges the resulting terms. + pub fn try_add(&self, other: &Self) -> Result { + self.check_width(other)?; + let left = self.normalized()?; + let right = other.normalized()?; + Self::from_terms(self.num_qubits, left.terms.into_iter().chain(right.terms)) + } + + /// Subtracts two operators on the same register, retaining every nonzero term. + pub fn try_sub(&self, other: &Self) -> Result { + self.check_width(other)?; + let left = self.normalized()?; + let right = other.normalized()?; + Self::from_terms( + self.num_qubits, + left.terms + .into_iter() + .chain(right.terms.into_iter().map(|(p, c)| (p, -c))), + ) + } + + /// Forms the matrix product `self * other`; the right operand acts first. + /// Both inputs are validated even when either represents zero. Products may + /// be non-Hermitian. Expansion takes O(m k n) work before sorting/merging + /// and may generate m k terms, for term counts m, k and register width n. + /// + /// ``` + /// use cqlib_core::qis::{Hamiltonian, QisError}; + /// use num_complex::Complex64; + /// let x = Hamiltonian::from_pauli("X".parse()?); + /// let z = Hamiltonian::from_pauli("Z".parse()?); + /// let minus_i_y = Hamiltonian::from_terms(1, [ + /// ("Y".parse()?, Complex64::new(0.0, -1.0)), + /// ])?; + /// assert!(x.try_mul(&z)?.equiv(&minus_i_y, 0.0)?); + /// # Ok::<(), QisError>(()) + /// ``` + pub fn try_mul(&self, other: &Self) -> Result { + self.check_width(other)?; + let left = self.normalized()?; + let right = other.normalized()?; + Self::from_terms( + self.num_qubits, + left.terms + .iter() + .flat_map(|(p, a)| right.terms.iter().map(move |(q, b)| (p * q, *a * *b))), + ) + } + + /// Returns a scaled sum, checking the factor even for the zero operator. + pub fn scaled(&self, factor: Complex64) -> Result { + validate_coefficient(factor)?; + let normalized = self.normalized()?; + Self::from_terms( + self.num_qubits, + normalized.terms.into_iter().map(|(p, c)| (p, c * factor)), + ) + } + + /// Returns the conjugate transpose after absorbing Pauli phases. + /// Phase-free Pauli strings are Hermitian, so their coefficients are conjugated. + pub fn adjoint(&self) -> Result { + let mut result = self.normalized()?; + for (_, coefficient) in &mut result.terms { + *coefficient = coefficient.conj(); + } + Ok(result) + } + + fn validate_tolerance(atol: f64) -> Result<(), QisError> { + if !atol.is_finite() || atol < 0.0 { + return Err(QisError::InvalidParameterValue( + "Pauli coefficient tolerance must be finite and nonnegative".into(), + )); + } + Ok(()) + } + + /// Normalizes, then removes coefficients with complex magnitude at most `atol`. + /// On any validation or merging error, the stored operator is unchanged. + pub fn chop(&mut self, atol: f64) -> Result<(), QisError> { + Self::validate_tolerance(atol)?; + let mut result = self.normalized()?; + result.terms.retain(|(_, c)| c.norm() > atol); + *self = result; + Ok(()) + } + + /// Compares normalized coefficients of `self - other` against `atol`. + /// Width mismatch is an error. This coefficient tolerance is not a matrix + /// operator-norm bound, and differs from structural `PartialEq`. + pub fn equiv(&self, other: &Self, atol: f64) -> Result { + Self::validate_tolerance(atol)?; + Ok(self + .try_sub(other)? + .terms + .iter() + .all(|(_, c)| c.norm() <= atol)) + } + + /// Checks `self.equiv(self.adjoint(), atol)` as for `FermionicOp`. + /// On a normalized sum this bounds twice the imaginary coefficient magnitude. + /// It does not project coefficients to real values; `normalized_hermitian` + /// retains its separate consumer contract (imaginary tolerance `1e-10`). + pub fn is_hermitian(&self, atol: f64) -> Result { + Self::validate_tolerance(atol)?; + self.equiv(&self.adjoint()?, atol) + } + + /// Returns `self ⊗ other`: `other` occupies low qubits and `self` high qubits. + /// Register widths are added even for zero operators. No dense matrix is built. + /// A nonzero result with a width beyond BitVec's length limit returns an error. + /// + /// ``` + /// use cqlib_core::qis::{Hamiltonian, QisError}; + /// let x = Hamiltonian::from_pauli("X".parse()?); + /// let z = Hamiltonian::from_pauli("Z".parse()?); + /// assert!(x.tensor(&z)?.equiv(&Hamiltonian::from_pauli("XZ".parse()?), 0.0)?); + /// # Ok::<(), QisError>(()) + /// ``` + pub fn tensor(&self, other: &Self) -> Result { + let width = self + .num_qubits + .checked_add(other.num_qubits) + .ok_or_else(|| { + QisError::InvalidParameterValue("tensor register width overflow".into()) + })?; + let left = self.normalized()?; + let right = other.normalized()?; + if !left.terms.is_empty() && !right.terms.is_empty() { + Self::check_materialized_width(width)?; + } + Self::from_terms( + width, + left.terms.iter().flat_map(|(p, a)| { + right.terms.iter().map(move |(q, b)| { + let mut product = q.clone(); + product.x.extend_from_bitslice(&p.x); + product.z.extend_from_bitslice(&p.z); + product.num_qubits = width; + (product, *a * *b) + }) + }), + ) + } + + /// Relabels qubits using `mapping[source] = destination`, inserting identities + /// on unused target qubits. Requires exactly one distinct destination per + /// source qubit, each less than `target_width`, even for zero/identity input. + /// This is a qubit embedding, not a fermionic mode permutation or a trace. + /// If terms must be constructed, unrepresentable target widths return + /// [`QisError::UnsupportedOperation`] before BitVec allocation. + pub fn remap_qubits(&self, mapping: &[usize], target_width: usize) -> Result { + if mapping.len() != self.num_qubits { + return Err(QisError::InvalidParameterValue( + "qubit mapping length mismatch".into(), + )); + } + let mut destinations = std::collections::HashSet::with_capacity(mapping.len()); + for &destination in mapping { + if destination >= target_width || !destinations.insert(destination) { + return Err(QisError::InvalidParameterValue( + "qubit mapping must contain distinct, in-range destinations".into(), + )); + } + } + let normalized = self.normalized()?; + if !normalized.terms.is_empty() { + Self::check_materialized_width(target_width)?; + } + Self::from_terms( + target_width, + normalized.terms.into_iter().map(|(p, c)| { + let mut mapped = PauliString::new(target_width); + for (source, &destination) in mapping.iter().enumerate() { + mapped.x.set(destination, p.x[source]); + mapped.z.set(destination, p.z[source]); + } + (mapped, c) + }), + ) + } +} + impl TryFrom> for Hamiltonian { type Error = QisError; @@ -195,8 +507,14 @@ impl Hamiltonian { // Evolution must merge terms without pruning small, nonzero coefficients. pub(crate) fn simplify_with_tolerance(&mut self, tolerance: f64) { + // Legacy simplification retains its infallible API and numerical behavior. + self.merge_terms(tolerance, false) + .expect("unchecked merging cannot fail"); + } + + fn merge_terms(&mut self, tolerance: f64, checked: bool) -> Result<(), QisError> { if self.terms.is_empty() { - return; + return Ok(()); } // 1. Normalize the Phase of all PauliStrings to Plus (+1), and multiply the phase into the coefficient. @@ -213,6 +531,9 @@ impl Hamiltonian { // Multiply the phase into the coefficient let new_coeff = coeff * phase_val; + if checked { + validate_coefficient(new_coeff)?; + } simplified_terms.push((pure_op, new_coeff)); } @@ -237,6 +558,9 @@ impl Hamiltonian { if op.x == current_op.x && op.z == current_op.z { // Identical operator found, accumulate the coefficient current_coeff += coeff; + if checked { + validate_coefficient(current_coeff)?; + } } else { // New operator found, save the previous one if its coefficient is significant // Filter out terms with a coefficient near zero @@ -258,6 +582,7 @@ impl Hamiltonian { } self.terms = merged; + Ok(()) } /// Scales all terms in the Hamiltonian by a complex factor. diff --git a/crates/cqlib-core/src/qis/hamiltonian_test.rs b/crates/cqlib-core/src/qis/hamiltonian_test.rs index c3a0417f..be4c2b05 100644 --- a/crates/cqlib-core/src/qis/hamiltonian_test.rs +++ b/crates/cqlib-core/src/qis/hamiltonian_test.rs @@ -15,6 +15,547 @@ use crate::qis::{PauliString, Phase, QisError}; use ndarray::{Array2, arr2}; use num_complex::Complex64; +/// Test-only Kronecker product, independent of production Pauli operations. +fn matrix_tensor(a: &Array2, b: &Array2) -> Array2 { + Array2::from_shape_fn((a.nrows() * b.nrows(), a.ncols() * b.ncols()), |(r, c)| { + a[[r / b.nrows(), c / b.ncols()]] * b[[r % b.nrows(), c % b.ncols()]] + }) +} + +fn matrix_adjoint(a: &Array2) -> Array2 { + a.t().mapv(|c| c.conj()) +} + +/// Literal single-qubit matrices and explicit labels, not production masks/phases. +fn matrix_oracle(terms: &[(&str, Phase, Complex64)]) -> Array2 { + let width = terms[0].0.len(); + let mut sum = Array2::zeros((1 << width, 1 << width)); + let zero = Complex64::new(0.0, 0.0); + let one = Complex64::new(1.0, 0.0); + let i = Complex64::new(0.0, 1.0); + for &(label, phase, coefficient) in terms { + let mut product = arr2(&[[one]]); + for letter in label.chars() { + let local = match letter { + 'I' => arr2(&[[one, zero], [zero, one]]), + 'X' => arr2(&[[zero, one], [one, zero]]), + 'Y' => arr2(&[[zero, -i], [i, zero]]), + 'Z' => arr2(&[[one, zero], [zero, -one]]), + _ => panic!("invalid test label"), + }; + product = matrix_tensor(&product, &local); + } + let phase = match phase { + Phase::Plus => one, + Phase::Minus => -one, + Phase::I => i, + Phase::MinusI => -i, + }; + sum += &product.mapv(|c| c * phase * coefficient); + } + sum +} + +fn raw_sum(terms: &[(&str, Phase, Complex64)]) -> Hamiltonian { + Hamiltonian::from_list( + terms + .iter() + .map(|&(label, phase, coefficient)| { + let mut pauli: PauliString = label.parse().unwrap(); + pauli.phase = phase; + (pauli, coefficient) + }) + .collect(), + ) + .unwrap() +} + +fn assert_matrix_close(actual: &Array2, expected: &Array2) { + assert_eq!(actual.dim(), expected.dim()); + for (a, b) in actual.iter().zip(expected) { + assert!((*a - *b).norm() < 1e-12, "{actual:?} != {expected:?}"); + } +} + +#[test] +fn sparse_arithmetic_matches_independent_matrices() { + let a_terms = [ + ("XYZ", Phase::MinusI, Complex64::new(0.5, 0.25)), + ("XYZ", Phase::Plus, Complex64::new(-0.125, 0.75)), + ("YIX", Phase::Minus, Complex64::new(1.0, -0.5)), + ("III", Phase::I, Complex64::new(-0.25, 0.5)), + ]; + let b_terms = [ + ("ZZY", Phase::I, Complex64::new(0.75, -0.25)), + ("XYY", Phase::Plus, Complex64::new(-0.5, 0.5)), + ]; + let a = raw_sum(&a_terms); + let b = raw_sum(&b_terms); + let saved = (a.clone(), b.clone()); + let am = matrix_oracle(&a_terms); + let bm = matrix_oracle(&b_terms); + assert_matrix_close(&a.to_matrix().unwrap(), &am); + assert_matrix_close(&a.try_add(&b).unwrap().to_matrix().unwrap(), &(&am + &bm)); + assert_matrix_close(&a.try_sub(&b).unwrap().to_matrix().unwrap(), &(&am - &bm)); + assert_matrix_close(&a.try_mul(&b).unwrap().to_matrix().unwrap(), &am.dot(&bm)); + assert_matrix_close( + &a.adjoint().unwrap().to_matrix().unwrap(), + &matrix_adjoint(&am), + ); + let factor = Complex64::new(0.25, -0.5); + assert_matrix_close( + &a.scaled(factor).unwrap().to_matrix().unwrap(), + &am.mapv(|c| c * factor), + ); + assert_eq!((a, b), saved); +} + +#[test] +fn exhaustive_two_qubit_products_include_all_global_phases() { + for a in ['I', 'X', 'Y', 'Z'] { + for b in ['I', 'X', 'Y', 'Z'] { + for c in ['I', 'X', 'Y', 'Z'] { + for d in ['I', 'X', 'Y', 'Z'] { + for phase in [Phase::Plus, Phase::Minus, Phase::I, Phase::MinusI] { + let left_label = format!("{a}{b}"); + let right_label = format!("{c}{d}"); + let left = [(left_label.as_str(), phase, Complex64::new(0.5, 0.25))]; + let right = [( + right_label.as_str(), + Phase::MinusI, + Complex64::new(-0.25, 0.5), + )]; + let product = raw_sum(&left).try_mul(&raw_sum(&right)).unwrap(); + assert_matrix_close( + &product.to_matrix().unwrap(), + &matrix_oracle(&left).dot(&matrix_oracle(&right)), + ); + } + } + } + } + } +} + +#[test] +fn cancellations_and_adjoint_product_order_follow_operator_algebra() { + let x = Hamiltonian::from_pauli("X".parse().unwrap()); + let z = Hamiltonian::from_pauli("Z".parse().unwrap()); + let sum = x.try_add(&z).unwrap(); + assert_eq!( + sum.try_mul(&sum).unwrap(), + Hamiltonian::identity(1).scaled(2.0.into()).unwrap() + ); + let product = x.try_mul(&z).unwrap(); + assert_eq!( + product.terms, + vec![("Y".parse().unwrap(), Complex64::new(0.0, -1.0))] + ); + assert!(!product.is_hermitian(0.0).unwrap()); + assert_eq!( + product.adjoint().unwrap(), + z.adjoint().unwrap().try_mul(&x.adjoint().unwrap()).unwrap() + ); + let commutator = product.try_sub(&z.try_mul(&x).unwrap()).unwrap(); + assert_eq!(commutator.terms[0].1, Complex64::new(0.0, -2.0)); + assert!( + product + .try_add(&z.try_mul(&x).unwrap()) + .unwrap() + .terms + .is_empty() + ); +} + +#[test] +fn checked_construction_and_comparison_keep_legacy_structure_semantics() { + let p: PauliString = "X".parse().unwrap(); + let raw = + Hamiltonian::from_list(vec![(p.clone(), 0.25.into()), (p.clone(), 0.75.into())]).unwrap(); + let canonical = Hamiltonian::from_terms(1, raw.terms.clone()).unwrap(); + assert_eq!(canonical, Hamiltonian::from_pauli(p)); + assert_ne!(raw, canonical); + assert!(raw.equiv(&canonical, 0.0).unwrap()); + assert_eq!((raw.clone() + canonical.clone()).terms.len(), 3); + assert_eq!(raw.try_add(&canonical).unwrap().terms.len(), 1); + assert_eq!(Hamiltonian::from_list(vec![]).unwrap().num_qubits, 0); + assert_eq!(Hamiltonian::from_terms(7, []).unwrap(), Hamiltonian::new(7)); + assert_eq!( + canonical.normalized().unwrap().normalized().unwrap(), + canonical + ); + let reversed = Hamiltonian::from_list(raw.terms.into_iter().rev().collect()).unwrap(); + assert_eq!(reversed.normalized().unwrap(), canonical); +} + +#[test] +fn explicit_chopping_follows_merging_and_preserves_tiny_coefficients() { + let tiny = f64::from_bits(1); + let p: PauliString = "Y".parse().unwrap(); + let mut phased = p.clone(); + phased.phase = Phase::I; + let small = Hamiltonian::from_terms(1, [(phased, tiny.into())]).unwrap(); + assert_eq!(small.terms[0].1, Complex64::new(0.0, tiny)); + for result in [ + small.try_mul(&Hamiltonian::identity(1)).unwrap(), + small.scaled(1.0.into()).unwrap(), + small.try_add(&Hamiltonian::new(1)).unwrap(), + small.adjoint().unwrap().adjoint().unwrap(), + ] { + assert_eq!(result, small); + } + let mut h = Hamiltonian::from_list(vec![(p.clone(), 0.75.into()), (p, 0.75.into())]).unwrap(); + h.chop(1.0).unwrap(); + assert_eq!(h.terms[0].1, Complex64::new(1.5, 0.0)); + h.chop(1.5).unwrap(); + assert_eq!(h, Hamiltonian::new(1)); + let mut small = small; + small.chop(0.0).unwrap(); + assert_eq!(small.terms.len(), 1); + small.chop(tiny).unwrap(); + assert!(small.terms.is_empty()); + let mut complex = + Hamiltonian::from_terms(1, [("X".parse().unwrap(), Complex64::new(3.0, 4.0))]).unwrap(); + complex.chop(4.0).unwrap(); + assert_eq!(complex.terms.len(), 1); + complex.chop(5.0).unwrap(); + assert!(complex.terms.is_empty()); +} + +#[test] +fn hermiticity_uses_adjoint_difference_tolerance_without_projection() { + let h = + Hamiltonian::from_terms(1, [("Z".parse().unwrap(), Complex64::new(1.0, 0.125))]).unwrap(); + assert!(h.is_hermitian(0.25).unwrap()); + assert!(!h.is_hermitian(0.249).unwrap()); + assert_eq!(h.terms[0].1.im, 0.125); + let phased = raw_sum(&[("X", Phase::I, Complex64::new(0.0, -1.0))]); + assert!(phased.is_hermitian(0.0).unwrap()); + assert!(!h.equiv(&Hamiltonian::identity(1), 0.0).unwrap()); +} + +#[test] +fn tensor_and_permutation_embedding_match_independent_matrices() { + let at = [("YX", Phase::MinusI, Complex64::new(0.5, 0.25))]; + let bt = [ + ("Z", Phase::I, Complex64::new(-0.25, 1.0)), + ("X", Phase::Minus, 0.5.into()), + ]; + let a = raw_sum(&at); + let b = raw_sum(&bt); + assert_matrix_close( + &a.tensor(&b).unwrap().to_matrix().unwrap(), + &matrix_tensor(&matrix_oracle(&at), &matrix_oracle(&bt)), + ); + let swapped = a.remap_qubits(&[1, 0], 2).unwrap(); + assert_matrix_close( + &swapped.to_matrix().unwrap(), + &matrix_oracle(&[("XY", Phase::MinusI, Complex64::new(0.5, 0.25))]), + ); + let embedded = a.remap_qubits(&[3, 1], 4).unwrap(); + assert_matrix_close( + &embedded.to_matrix().unwrap(), + &matrix_oracle(&[("XIYI", Phase::MinusI, Complex64::new(0.5, 0.25))]), + ); + let scalar = Hamiltonian::identity(0) + .scaled(Complex64::new(2.0, -0.5)) + .unwrap(); + assert_eq!( + a.tensor(&scalar).unwrap(), + a.scaled(Complex64::new(2.0, -0.5)).unwrap() + ); + assert_eq!( + scalar.tensor(&a).unwrap(), + a.scaled(Complex64::new(2.0, -0.5)).unwrap() + ); +} + +#[test] +fn zero_and_scalar_operators_retain_declared_spaces() { + for width in [0, 1, 128] { + let zero = Hamiltonian::new(width); + let identity = Hamiltonian::identity(width); + assert_eq!(zero.try_mul(&identity).unwrap(), zero); + assert_eq!(identity.try_mul(&zero).unwrap(), zero); + assert_eq!(identity.try_sub(&identity).unwrap(), zero); + assert_eq!(identity.scaled(0.0.into()).unwrap(), zero); + assert_eq!( + zero.tensor(&Hamiltonian::new(2)).unwrap().num_qubits, + width + 2 + ); + let map: Vec<_> = (0..width).map(|q| q + 1).collect(); + assert_eq!( + zero.remap_qubits(&map, width + 2).unwrap(), + Hamiltonian::new(width + 2) + ); + assert_eq!( + identity.remap_qubits(&map, width + 2).unwrap(), + Hamiltonian::identity(width + 2) + ); + } + assert_eq!( + Hamiltonian::identity(0).to_matrix().unwrap(), + arr2(&[[Complex64::new(1.0, 0.0)]]) + ); +} + +#[test] +fn wide_register_arithmetic_never_needs_machine_word_masks() { + use crate::qis::Pauli; + for width in [65, 128, 1000] { + let mut p = PauliString::new(width); + p.set_pauli(width - 1, Pauli::Y); + p.set_pauli(0, Pauli::Z); + let a = Hamiltonian::from_terms(width, [(p.clone(), 0.5.into())]).unwrap(); + assert_eq!( + a.try_mul(&a).unwrap(), + Hamiltonian::identity(width).scaled(0.25.into()).unwrap() + ); + assert!(a.is_hermitian(0.0).unwrap()); + assert!( + a.try_add(&a) + .unwrap() + .equiv(&a.scaled(2.0.into()).unwrap(), 0.0) + .unwrap() + ); + let mapping: Vec<_> = (0..width).rev().collect(); + let reversed = a.remap_qubits(&mapping, width).unwrap(); + assert_eq!(reversed.terms[0].0.get_pauli(0), Pauli::Y); + assert_eq!(reversed.terms[0].0.get_pauli(width - 1), Pauli::Z); + assert_eq!(reversed.remap_qubits(&mapping, width).unwrap(), a); + let tensor = a + .tensor(&Hamiltonian::from_pauli("X".parse().unwrap())) + .unwrap(); + assert_eq!(tensor.num_qubits, width + 1); + assert_eq!(tensor.terms[0].0.get_pauli(0), Pauli::X); + assert_eq!(tensor.terms[0].0.get_pauli(width), Pauli::Y); + } +} + +#[test] +fn invalid_inputs_are_rejected_even_when_other_operand_is_zero() { + let mut short = PauliString::new(1); + short.x.clear(); + let mut long = PauliString::new(1); + long.z.push(false); + for terms in [ + vec![(PauliString::new(2), 1.0.into())], + vec![(short, 1.0.into())], + vec![(long, 1.0.into())], + vec![(PauliString::new(1), f64::NAN.into())], + vec![(PauliString::new(1), f64::INFINITY.into())], + vec![(PauliString::new(1), Complex64::new(0.0, f64::NEG_INFINITY))], + ] { + let h = Hamiltonian { + num_qubits: 1, + terms: terms.clone(), + }; + let zero = Hamiltonian::new(1); + assert!(Hamiltonian::from_terms(1, terms).is_err()); + assert!(h.try_add(&zero).is_err()); + assert!(h.try_sub(&zero).is_err()); + assert!(h.try_mul(&zero).is_err()); + assert!(zero.try_mul(&h).is_err()); + assert!(h.scaled(0.0.into()).is_err()); + assert!(h.adjoint().is_err()); + assert!(h.tensor(&zero).is_err()); + assert!(zero.tensor(&h).is_err()); + assert!(h.remap_qubits(&[0], 1).is_err()); + assert!(h.equiv(&h, 0.0).is_err()); + assert!(h.is_hermitian(0.0).is_err()); + assert!(h.to_matrix().is_err()); + let mut chopped = h.clone(); + assert!(chopped.chop(0.0).is_err()); + // Compare coefficient bits because NaN is not equal to itself. + assert_eq!(chopped.num_qubits, h.num_qubits); + assert_eq!(chopped.terms.len(), h.terms.len()); + for ((actual_p, actual_c), (expected_p, expected_c)) in chopped.terms.iter().zip(&h.terms) { + assert_eq!(actual_p, expected_p); + assert_eq!(actual_c.re.to_bits(), expected_c.re.to_bits()); + assert_eq!(actual_c.im.to_bits(), expected_c.im.to_bits()); + } + } + let a = Hamiltonian::new(1); + let b = Hamiltonian::new(2); + for result in [a.try_add(&b), a.try_sub(&b), a.try_mul(&b)] { + assert!(matches!(result, Err(QisError::QubitMismatch { .. }))); + } + assert!(a.equiv(&b, 0.0).is_err()); +} + +#[test] +fn invalid_tolerances_mappings_and_scalars_do_not_mutate_input() { + let h = Hamiltonian::identity(2); + for atol in [-1.0, f64::NAN, f64::INFINITY] { + let mut chopped = h.clone(); + assert!(chopped.chop(atol).is_err()); + assert_eq!(chopped, h); + assert!(h.equiv(&h, atol).is_err()); + assert!(h.is_hermitian(atol).is_err()); + } + for mapping in [&[0][..], &[0, 0], &[0, 2]] { + assert!(h.remap_qubits(mapping, 2).is_err()); + assert!(Hamiltonian::new(2).remap_qubits(mapping, 2).is_err()); + } + for factor in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + assert!(Hamiltonian::new(2).scaled(factor.into()).is_err()); + } + assert!( + Hamiltonian::new(usize::MAX) + .tensor(&Hamiltonian::new(1)) + .is_err() + ); +} + +#[test] +fn unrepresentable_embeddings_return_errors_without_panicking() { + let limit = bitvec::slice::BitSlice::::MAX_BITS; + for width in [limit + 1, usize::MAX] { + // Empty source mapping is valid, but scalar identity still requires a + // target Pauli string: reject its unrepresentable width before allocation. + let scalar = Hamiltonian::identity(0); + assert!(matches!( + scalar.remap_qubits(&[], width), + Err(QisError::UnsupportedOperation(_)) + )); + assert!(matches!( + Hamiltonian::identity(1).remap_qubits(&[width - 1], width), + Err(QisError::UnsupportedOperation(_)) + )); + // Zero sums do not allocate target bit vectors and retain declared width. + assert_eq!( + Hamiltonian::new(0).remap_qubits(&[], width).unwrap(), + Hamiltonian::new(width) + ); + assert_eq!( + Hamiltonian::new(width).tensor(&scalar).unwrap(), + Hamiltonian::new(width) + ); + } +} + +#[test] +fn intermediate_and_dense_coefficient_overflow_is_reported() { + let max = Hamiltonian::identity(1).scaled(f64::MAX.into()).unwrap(); + let two = Hamiltonian::identity(1).scaled(2.0.into()).unwrap(); + for result in [ + max.try_add(&max), + max.scaled(2.0.into()), + max.try_mul(&two), + max.tensor(&two), + max.try_sub(&max.scaled((-1.0).into()).unwrap()), + ] { + assert!(matches!(result, Err(QisError::NonFiniteCoefficient { .. }))); + } + // A later cancellation cannot hide an earlier nonfinite accumulation. + let raw = Hamiltonian::from_list( + [f64::MAX, f64::MAX, -f64::MAX] + .into_iter() + .map(|c| (PauliString::new(1), c.into())) + .collect(), + ) + .unwrap(); + assert!(matches!( + raw.normalized(), + Err(QisError::NonFiniteCoefficient { .. }) + )); + let mut chopped = raw.clone(); + assert!(chopped.chop(0.0).is_err()); + assert_eq!(chopped, raw); + let dense = Hamiltonian::from_terms( + 1, + [ + ("I".parse().unwrap(), f64::MAX.into()), + ("Z".parse().unwrap(), f64::MAX.into()), + ], + ) + .unwrap(); + assert!(matches!( + dense.to_matrix(), + Err(QisError::NonFiniteCoefficient { .. }) + )); +} + +#[test] +fn dense_conversion_checks_shift_element_and_byte_overflow_before_allocation() { + // On 64-bit: 64 overflows the shift, 32 the element count, 30 the byte range. + for width in [ + usize::BITS as usize, + usize::BITS as usize / 2, + (usize::BITS as usize - 4) / 2, + usize::MAX, + ] { + assert!(matches!( + Hamiltonian::new(width).to_matrix(), + Err(QisError::UnsupportedOperation(_)) + )); + } +} + +#[test] +fn mapped_fermionic_products_and_adjoint_use_public_pauli_algebra() { + use crate::qis::{FermionToQubitMapper, FermionicOp, JordanWignerMapper}; + let mapper = JordanWignerMapper::new(); + for p in 0..3 { + for q in 0..3 { + let a = FermionicOp::creation(3, p).unwrap(); + let b = FermionicOp::annihilation(3, q).unwrap(); + let mapped_a = mapper.map_operator(&a).unwrap(); + let mapped_b = mapper.map_operator(&b).unwrap(); + assert!( + mapped_a + .try_mul(&mapped_b) + .unwrap() + .equiv(&mapper.map_operator(&a.try_mul(&b).unwrap()).unwrap(), 0.0) + .unwrap() + ); + assert!( + mapped_a + .adjoint() + .unwrap() + .equiv(&mapper.map_operator(&a.adjoint()).unwrap(), 0.0) + .unwrap() + ); + } + } +} + +#[test] +fn hermitian_normalization_projects_only_accepted_imaginary_residuals() { + let h = Hamiltonian::from_list(vec![ + ("I".parse().unwrap(), Complex64::new(0.0, 1e-12)), + ("Z".parse().unwrap(), Complex64::new(1e-16, -1e-10)), + ]) + .unwrap(); + let normalized = h.normalized_hermitian().unwrap(); + assert_eq!(normalized.num_qubits, 1); + assert_eq!( + normalized.terms, + vec![("Z".parse().unwrap(), Complex64::new(1e-16, 0.0))] + ); + assert_eq!(normalized.normalized_hermitian().unwrap(), normalized); + // General complex normalization must not apply the real projection. + assert_eq!(h.normalized().unwrap().terms.len(), 2); + + let mut invalid = h.clone(); + invalid.terms[0].1.im = 1.01e-10; + assert!(matches!( + invalid.normalized_hermitian(), + Err(QisError::NotHermitian) + )); + + // Validate the merged imaginary coefficient before removing zero real parts. + let merged_invalid = Hamiltonian::from_list(vec![ + ("Z".parse().unwrap(), Complex64::new(0.0, 0.75e-10)), + ("Z".parse().unwrap(), Complex64::new(0.0, 0.75e-10)), + ]) + .unwrap(); + assert!(matches!( + merged_invalid.normalized_hermitian(), + Err(QisError::NotHermitian) + )); +} + #[test] fn empty_hamiltonian_to_matrix_is_zero() { assert_eq!( diff --git a/crates/cqlib-core/src/qis/mod.rs b/crates/cqlib-core/src/qis/mod.rs index cf196fb3..45906459 100644 --- a/crates/cqlib-core/src/qis/mod.rs +++ b/crates/cqlib-core/src/qis/mod.rs @@ -24,6 +24,7 @@ //! - [`state`]: Quantum state representations (statevector, density matrix, stabilizer) //! - [`pauli`]: Pauli operators and Pauli strings with symplectic encoding //! - [`hamiltonian`]: Hamiltonian construction from Pauli strings +//! - [`fermion`]: Fermionic operator algebra and fermion-to-qubit mappings //! - [`observable`]: Trait for computing expectation values //! - [`entropy`]: Entropy measures and entanglement metrics //! - [`error`]: Error types for QIS operations @@ -79,6 +80,7 @@ pub mod entropy; pub mod error; pub mod evolution; +pub mod fermion; pub mod hamiltonian; pub mod metrics; pub mod observable; @@ -88,6 +90,9 @@ pub mod state; pub use error::{PauliStringParseError, QisError}; pub use evolution::{PauliEvolution, TrotterMode}; +pub use fermion::{ + FermionToQubitMapper, FermionicLadder, FermionicOp, JordanWignerMapper, OccupationStateMapper, +}; pub use hamiltonian::Hamiltonian; pub use observable::Observable; pub use pauli::{Pauli, PauliIter, PauliString, Phase}; diff --git a/crates/cqlib-core/src/qis/observable.rs b/crates/cqlib-core/src/qis/observable.rs index ac0e4958..ce1733e5 100644 --- a/crates/cqlib-core/src/qis/observable.rs +++ b/crates/cqlib-core/src/qis/observable.rs @@ -77,7 +77,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. +/// observables. Nonfinite Hamiltonian coefficients, including merge overflow, +/// return [`QisError::InvalidParameterValue`]. pub trait Observable { /// Computes the expectation value for a statevector: $\langle\psi|O|\psi\rangle$. /// @@ -185,28 +186,14 @@ 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); +// Keep the existing Observable error contract while sharing algebraic validation. +fn normalized_observable(h: &Hamiltonian) -> Result { + h.normalized_hermitian().map_err(|error| match error { + nonfinite @ QisError::NonFiniteCoefficient { .. } => { + QisError::InvalidParameterValue(nonfinite.to_string()) } - } - Ok(()) + other => other, + }) } impl Observable for Hamiltonian { @@ -218,11 +205,11 @@ impl Observable for Hamiltonian { }); } - ensure_hermitian_hamiltonian(self)?; + let normalized = normalized_observable(self)?; let mut expected_value = 0.0; // Iterate over each term in the Hamiltonian - for (pauli_str, coeff) in &self.terms { + for (pauli_str, coeff) in &normalized.terms { let x_mask = pauli_str.x_mask(); let z_mask = pauli_str.z_mask(); @@ -267,12 +254,12 @@ impl Observable for Hamiltonian { }); } - ensure_hermitian_hamiltonian(self)?; + let normalized = normalized_observable(self)?; let mut expected_value = 0.0; let n = self.num_qubits; let dim = 1 << n; // 2^n - for (pauli_str, coeff) in &self.terms { + for (pauli_str, coeff) in &normalized.terms { let x_mask = pauli_str.x_mask(); let z_mask = pauli_str.z_mask(); @@ -314,10 +301,10 @@ impl Observable for Hamiltonian { &self, measurements: &[(PauliString, HashMap)], ) -> Result { - ensure_hermitian_hamiltonian(self)?; + let normalized = normalized_observable(self)?; let mut expected_value = 0.0; - for (term_pauli, coeff) in &self.terms { + for (term_pauli, coeff) in &normalized.terms { let term_x_mask = term_pauli.x_mask(); let term_z_mask = term_pauli.z_mask(); let p_active = term_x_mask | term_z_mask; @@ -397,9 +384,7 @@ impl Observable for Hamiltonian { }); } - ensure_hermitian_hamiltonian(self)?; - let mut simplified = self.clone(); - simplified.simplify_with_tolerance(0.0); + let mut simplified = normalized_observable(self)?; for (_, coeff) in &mut simplified.terms { if coeff.im.abs() > VARIANCE_TOLERANCE { diff --git a/crates/cqlib-core/src/qis/observable_test.rs b/crates/cqlib-core/src/qis/observable_test.rs index 1fc05c34..4b5b4e03 100644 --- a/crates/cqlib-core/src/qis/observable_test.rs +++ b/crates/cqlib-core/src/qis/observable_test.rs @@ -390,3 +390,162 @@ fn test_empty_hamiltonian_variance_is_zero() { fn test_variance_tiny_negative_roundoff_clamps_to_zero() { assert_close(finalize_variance(-VARIANCE_TOLERANCE / 2.0).unwrap(), 0.0); } + +#[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", Complex64::new(0.0, -1.0))], + vec![ + ("Z", Complex64::new(1.0, 2.0)), + ("Z", Complex64::new(0.0, -2.0)), + ], + vec![("Z", Complex64::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(), Complex64::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 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(), Complex64::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(), Complex64::new(f64::MAX, 0.0)); + 2 + ]) + .unwrap(); + assert!(matches!( + h.expectation_statevector(&Statevector::new(1)), + Err(QisError::InvalidParameterValue(_)) + )); +} + +#[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 [Complex64::new(bad, 0.0), Complex64::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(), Complex64::new(f64::MAX, 0.0)); + 2 + ]) + .unwrap(); + assert!(matches!( + h.variance_statevector(&sv), + Err(QisError::InvalidParameterValue(_)) + )); + for terms in [ + vec![("iX", Complex64::new(0.0, -1.0))], + vec![ + ("X", Complex64::new(1.0, 2.0)), + ("X", Complex64::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(), Complex64::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/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs b/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs index 9c0bbfa9..82775269 100644 --- a/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs +++ b/crates/cqlib-core/src/qis/state/density_matrix_noise_test.rs @@ -1459,3 +1459,122 @@ fn test_ccx_with_noise_model_does_not_panic() { "CCX should produce support on |111>, got {probs:?}" ); } + +#[test] +fn pauli_probability_boundary_preserves_simulation_trace() { + let noise = SingleQubitNoise::Pauli { + px: 0.3, + py: 0.6, + pz: 0.1, + }; + 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); +} + +#[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 failed_noisy_circuit_preserves_label_binding() { + use crate::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()); +} + +#[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() + ); +} diff --git a/crates/cqlib-core/src/qis/state/state_test.rs b/crates/cqlib-core/src/qis/state/state_test.rs index 86e40296..2c707bb5 100644 --- a/crates/cqlib-core/src/qis/state/state_test.rs +++ b/crates/cqlib-core/src/qis/state/state_test.rs @@ -1,5 +1,4 @@ // This code is part of Cqlib. -// Modified to support reproducible state sampling. // // (C) Copyright China Telecom Quantum Group 2026 // @@ -10,8 +9,6 @@ // 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. -// -// Modified to validate terminal measurements and clarify or correct reset behavior. //! Simulation consistency tests between Statevector and DensityMatrix @@ -875,3 +872,64 @@ fn full_measurement_uses_the_supplied_random_stream() { ); } } + +#[test] +fn checked_constructors_reject_dimensions_without_allocating() { + 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 nonfinite_pure_states_are_rejected() { + use num_complex::Complex64; + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY, f64::MAX] { + for amp in [Complex64::new(bad, 0.0), Complex64::new(0.0, bad)] { + let state = vec![amp, Complex64::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![Complex64::new(0.0, 0.0); 1 << n]; + v[0] = Complex64::new(0.0, 1.0); + assert!(Statevector::from_state(n, v.clone()).is_ok()); + assert!(DensityMatrix::from_state(n, v).is_ok()); + } +} diff --git a/crates/cqlib-core/src/qis/state/statevector.rs b/crates/cqlib-core/src/qis/state/statevector.rs index 9b94a36a..9b4f85da 100644 --- a/crates/cqlib-core/src/qis/state/statevector.rs +++ b/crates/cqlib-core/src/qis/state/statevector.rs @@ -471,6 +471,7 @@ impl Statevector { /// Applies a quantum circuit to this statevector in-place. /// + /// The circuit's resolved global phase is applied, retaining nonzero residuals. /// Once bound to circuit qubit IDs, subsequent circuits must have the same /// IDs in the same order. The mapping is committed only on success; gates /// completed before an error are not rolled back. @@ -520,6 +521,15 @@ impl Statevector { // Decompose circuit to basic gates let circuit = circuit.decompose()?; super::validate_terminal_measurements(&circuit)?; + let global_phase = circuit + .global_phase() + .evaluate(&None) + .map_err(|_| QisError::CircuitError(CircuitError::SymbolicParameterError))?; + if !global_phase.is_finite() { + return Err(QisError::InvalidParameterValue( + "circuit global phase must be finite".into(), + )); + } let sv = self; let qubit_map = super::prepare_qubit_map(&circuit.qubits(), &sv.qubit_map)?; @@ -644,6 +654,7 @@ impl Statevector { } } + sv.apply_gphase(global_phase)?; sv.qubit_map = Some(qubit_map); Ok(()) } @@ -2688,7 +2699,7 @@ impl Statevector { /// # Arguments /// * `phi` - Phase angle pub fn apply_gphase(&mut self, phi: f64) -> Result<(), QisError> { - if phi.abs() < 1e-15 { + if phi == 0.0 { return Ok(()); } let phase = Complex64::from_polar(1.0, phi); diff --git a/crates/cqlib-core/src/qis/state/statevector_test.rs b/crates/cqlib-core/src/qis/state/statevector_test.rs index de277f7b..28f9d8a3 100644 --- a/crates/cqlib-core/src/qis/state/statevector_test.rs +++ b/crates/cqlib-core/src/qis/state/statevector_test.rs @@ -908,6 +908,35 @@ fn test_global_phase() { ); } +#[test] +fn circuit_evolution_applies_global_phase_including_small_residuals() { + use crate::circuit::Parameter; + for num_qubits in [0, 1] { + for angle in [0.3, 1e-12, 1e-20] { + let mut circuit = Circuit::new(num_qubits); + circuit.set_global_phase(Parameter::from(angle)); + let state = Statevector::from_circuit(&circuit).unwrap(); + let expected = Complex64::from_polar(1.0, angle); + assert_eq!(state.data[0], expected); + let mut incremental = Statevector::new(num_qubits); + incremental.apply_circuit(&circuit).unwrap(); + assert_eq!(incremental.data[0], expected); + } + } +} + +#[test] +fn circuit_evolution_rejects_unbound_global_phase_before_mutating() { + use crate::circuit::Parameter; + let mut circuit = Circuit::new(1); + circuit.x(Qubit::new(0)).unwrap(); + circuit.set_global_phase(Parameter::symbol("phi")); + let mut state = Statevector::new(1); + let before = state.data.clone(); + assert!(state.apply_circuit(&circuit).is_err()); + assert_eq!(&state.data[..], &before[..]); +} + #[test] fn test_relative_phase() { // S gate introduces relative phase diff --git a/crates/cqlib-core/tests/p1_regressions.rs b/crates/cqlib-core/tests/p1_regressions.rs deleted file mode 100644 index 252fc083..00000000 --- a/crates/cqlib-core/tests/p1_regressions.rs +++ /dev/null @@ -1,243 +0,0 @@ -// 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 deleted file mode 100644 index ffa6faaa..00000000 --- a/crates/cqlib-core/tests/p2_regressions.rs +++ /dev/null @@ -1,299 +0,0 @@ -// 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/python/3_qis/0_overview.md b/docs/api/cn/python/3_qis/0_overview.md index 9339a0b1..fb3f4aba 100644 --- a/docs/api/cn/python/3_qis/0_overview.md +++ b/docs/api/cn/python/3_qis/0_overview.md @@ -29,6 +29,8 @@ ### 观测量与时间演化 +[`FermionicOp` 与 `JordanWignerMapper`](10_fermion.md) 提供通用费米代数、激发生成元及费米到 Pauli 和的映射，可直接复用现有演化组件。 + `PauliString` 与 `Hamiltonian` 描述可观测量，两者都可以给出期望值；`Hamiltonian` 还可以通过 Trotter 分解转成演化线路，回到 `cqlib.circuit` 的线路空间中执行。 --- diff --git a/docs/api/cn/python/3_qis/10_fermion.md b/docs/api/cn/python/3_qis/10_fermion.md new file mode 100644 index 00000000..15670693 --- /dev/null +++ b/docs/api/cn/python/3_qis/10_fermion.md @@ -0,0 +1,72 @@ +# 费米算子与 Jordan–Wigner 映射 + +`cqlib.qis.fermion` 提供通用费米代数与稀疏映射。公共类型也可以从 `cqlib.qis` 导入，计算由 Rust Core 完成。 + +```python +from cqlib.qis import FermionicLadder, FermionicOp, JordanWignerMapper + +a = FermionicOp.annihilation(num_modes=4, mode=0) +identity = FermionicOp.identity(4) +assert (a * a.adjoint() + a.adjoint() * a).equiv(identity) + +term = [FermionicLadder.create(2), FermionicLadder.annihilate(0)] +op = FermionicOp.from_term(num_modes=4, term=term, coefficient=0.5j) +hamiltonian = JordanWignerMapper().map_operator(op) +``` + +## 算子与运算 + +`FermionicLadder.create(mode)` 和 `annihilate(mode)` 创建不可变单算子，提供只读 `mode`、`kind`、`adjoint()`、值比较和哈希。它自身不声明模态空间，越界检查发生在构造 `FermionicOp` 时。 + +`FermionicOp(num_modes)` 创建零算子；`identity(num_modes)` 创建恒等算子。`from_term(num_modes, term, coefficient=1.0)` 接受一个有序乘积；`from_terms(num_modes, terms)` 接受 `(term, coefficient)` 二元组的可迭代对象，合并重复项。空乘积是恒等项，空项集合是零算子，二者都保留模态数量，包括零模态。 + +| 接口 | 行为 | +| --- | --- | +| `creation` / `annihilation` / `number` | 参数为 `num_modes, mode`；构造产生、湮灭、模态数算子。 | +| `total_number(num_modes)` | 构造总粒子数算子。 | +| `A + B` / `A - B` / `-A` | 返回新算子；加减要求模态数量一致。 | +| `A * B` | 算子乘积 AB，B 先作用；不自动正规序。 | +| `c * A` / `A * c` / `A.scaled(c)` | 标量乘法，接受有限的整数、浮点数和复数。 | +| `adjoint()` | 共轭系数、反转乘积并交换产生与湮灭。 | +| `normal_order()` | 返回新展开：产生算子在前，每组模态编号升序，保留收缩项与符号。 | +| `add_term(term, coefficient=1.0)` | 原地添加项；失败后对象不变，返回 `None`。 | +| `chop(*, atol=0.0)` | 原地删除复模长不大于容差的项，返回 `None`。 | +| `equiv(other, *, atol=0.0)` | 正规序后判断代数等价，模态数量不一致报错。 | +| `is_hermitian(*, atol=0.0)` | 判断算子与其共轭转置是否代数等价。 | + +只读属性 `num_modes`、`num_terms` 描述存储；`terms` 返回独立的 `list[(tuple[FermionicLadder, ...], complex)]` 快照。`==` 比较存储结构，`is_zero()` 检查存储为空；数学等价与代数零使用 `equiv()` 判断。 + +默认只删除精确零，不隐式截断微小系数。容差必须有限且非负。`copy()`、`copy.copy()`、`copy.deepcopy()` 创建独立副本。算子可变且不可哈希；`A += B` 重新绑定结果，不修改 A 的别名。求和使用 `sum(operators, start=FermionicOp(num_modes))`。 + +## 激发生成元与演化 + +`excitation_generator(num_modes, sources, targets)` 构造 G = i(T − T†)。T 的产生算子按 targets 原顺序排列，湮灭算子按 sources 逆序排列。例如 `[0, 1] -> [2, 3]` 表示 T = a†₂ a†₃ a₁ a₀。列表必须非空、等长、各自无重复且互不重叠；不自动排序，也不引入阶乘归一化。 + +```python +from cqlib.circuit.ansatz import BasisEncoding, EvolutionStrategy, PauliEvolutionAnsatz + +mapper = JordanWignerMapper() +circuit = BasisEncoding(mapper.map_occupation([1, 1, 0, 0])).build_circuit("hf") +for name, sources, targets in [ + ("theta_double", [0, 1], [2, 3]), + ("theta_single", [0], [2]), +]: + generator = FermionicOp.excitation_generator(4, sources, targets) + block = (PauliEvolutionAnsatz(mapper.map_operator(generator)) + .with_strategy(EvolutionStrategy.exact()) + .with_time_param_name(name).build_circuit("excitation")) + circuit.compose(block) +bound = circuit.assign_parameters({"theta_double": 0.2, "theta_single": -0.3}) +``` + +每块演化为 exp(−iθG) = exp(θ(T − T†))。多个独立参数应分别构造、命名并组合演化块；把生成元相加后构造一块演化会共享一个参数。 + +## 映射与扩展 + +`JordanWignerMapper()` 提供 `map_operator(op) -> Hamiltonian`、`num_qubits(num_modes) -> int` 和 `map_occupation(occupation) -> list[bool]`。算子映射支持非厄米、恒等和零算子；线路与期望值消费者自行要求厄米性。 + +模态从 0 编号，模态 p 对应比特 p，Z 链覆盖小于 p 的编号，比特 0 是最低有效位。占据输入必须为整数 0/1 或布尔值；返回 `bits[p]` 对应比特 p，空输入返回空列表。Core 不规定 α/β 或化学轨道顺序。 + +`FermionToQubitMapper` 与 `OccupationStateMapper` 是独立的、可运行时检查的 `typing.Protocol`。开发者可以实现同签名的方法，不必继承协议；运行时检查只确认成员存在，不验证签名和计算正确性。Python 映射器由算法代码直接调用，不自动桥接到 Rust trait。 + +类型错误抛出 `TypeError`；负数或越界模态抛出 `IndexError`；负模态数量、空间不匹配、非有限系数、非法占据值或容差抛出 `ValueError`；模态整数超出 Rust 可表示范围抛出 `OverflowError`。 diff --git a/docs/api/cn/python/3_qis/7_hamiltonian.md b/docs/api/cn/python/3_qis/7_hamiltonian.md index 6cb11f52..fd09bb46 100644 --- a/docs/api/cn/python/3_qis/7_hamiltonian.md +++ b/docs/api/cn/python/3_qis/7_hamiltonian.md @@ -223,3 +223,39 @@ assert trotterized.num_qubits == 1 | `ValueError` | `expectation_probs` 无法从给定的测量基推出待求算符。 | 观测量计算的输入类型见 [Statevector](1_statevector.md) 与 [DensityMatrix](2_density_matrix.md)；Trotter 模式的定义见 [Trotter 演化](8_evolution.md)。 + +## 经过校验的稀疏代数 + +`from_terms(num_qubits, terms)` 明确指定宽度，接受 `(PauliString, 整数/浮点数/复数)` 项。它校验尺寸与有限系数、吸收 Pauli 相位并合并重复项，只删除精确零；空输入保留指定宽度。 + +| 接口 | 行为 | +| --- | --- | +| `normalized()` | 返回校验与合并后的新对象，保留所有非零系数。 | +| `add(other)` | 同宽度的安全加法，返回合并后的新对象。 | +| `A - B`、`A * B`、`-A` | 安全减法、算子乘积及取负；AB 中 B 先作用。 | +| `c * A`、`A * c`、`scaled(c)` | 校验标量及结果系数，返回新对象。 | +| `adjoint()` | 返回共轭转置，包括 Pauli 相位。 | +| `chop(*, atol=0.0)` | 原地规范化，并删除复模长不大于容差的系数；失败后对象不变。 | +| `equiv(other, *, atol=0.0)` | 比较规范化后的系数差，宽度不一致报错；容差不是矩阵算子范数界。 | +| `is_hermitian(*, atol=0.0)` | 与共轭转置比较，即规范化系数满足 2·abs(imag) ≤ atol。 | +| `normalized_hermitian()` | 将模长不大于 1e-10 的虚部残差投影为零，更大虚部报错；保留任意非零实部。 | +| `tensor(other)` | 返回 self ⊗ other，other 占低编号比特。 | +| `remap_qubits(mapping, target_width)` | mapping[source] = destination；目标必须唯一且在范围内，未使用位置补恒等。 | + +这些接口支持非厄米算子，厄米性由演化与期望值消费者校验。所有容差必须有限且非负。比特嵌入是量子比特重标记，不是费米模态置换。 + +```python +from cqlib.qis import Hamiltonian, PauliString + +x = Hamiltonian.from_terms(1, [(PauliString.from_str("X"), 1.0)]) +z = Hamiltonian.from_terms(1, [(PauliString.from_str("Z"), 1.0)]) +assert (x * z).equiv(Hamiltonian.from_terms(1, [(PauliString.from_str("Y"), -1j)])) +assert x.add(x).equiv(2 * x) +assert x.remap_qubits([2], target_width=3).num_qubits == 3 +``` + +### 旧接口兼容 + +`from_list`、`add_term`、`scale`、`+`、`+=` 保留现有行为。其中 `+` 与 `+=` 拼接存储项；`scale` 原地缩放。`simplify()` 使用旧的近零截断策略，会删除微小系数。需要完整保留映射结果时使用 `normalized()`，需要明确截断时使用 `chop(atol=...)`。 + +`add_term()` 和 `scale()` 原地修改并返回 `None`，不支持链式调用。新代数运算遇到宽度不一致、非有限系数或中间系数溢出时抛出 `ValueError`。 diff --git a/docs/api/cn/rust/3_qis/7_hamiltonian.md b/docs/api/cn/rust/3_qis/7_hamiltonian.md index 14757c31..ce86449d 100644 --- a/docs/api/cn/rust/3_qis/7_hamiltonian.md +++ b/docs/api/cn/rust/3_qis/7_hamiltonian.md @@ -3,9 +3,9 @@ `cqlib_core::qis::hamiltonian` `cqlib_core::qis::observable` -本页介绍由 Pauli 串与复系数构成的观测量 `Hamiltonian`，以及计算期望值的统一接口 `Observable` trait。 +本页介绍由 Pauli 串与复系数构成的算符 `Hamiltonian`，以及计算实数期望值的统一接口 `Observable` trait。 -`Hamiltonian` 是 $2^N \times 2^N$ 矩阵的稀疏表示，形式为 $H = \sum_k c_k P_k$，其中 $c_k$ 是复系数，$P_k$ 是 $N$ 比特 Pauli 串。同一个结构既用于物理解释中的哈密顿量，也用于任意需要求期望值的可观测量。 +`Hamiltonian` 是 $2^N \times 2^N$ 矩阵的稀疏表示，形式为 $H = \sum_k c_k P_k$，其中 $c_k$ 是复系数，$P_k$ 是 $N$ 比特 Pauli 串。该结构允许非厄米算符；实数期望值和酉演化入口会额外检查厄米性。 ## 导入 @@ -46,10 +46,59 @@ pub struct Hamiltonian { - `fn simplify(&mut self)`：化简项列表。该操作分两步：把每个 `PauliString` 的内部相位吸收进系数，使所有 Pauli 串归一到 `Phase::Plus`；再按 X、Z 分量向量分组，合并相同算子的系数。绝对值小于 $10^{-10}$ 的系数会被丢弃。 - `fn scale(&mut self, factor: Complex64)`：把所有项的系数乘以同一个复因子。 +### 安全 Pauli 和式代数 + +以下接口返回 `Result`，验证比特数、X/Z 位向量长度、有限系数及系数乘法/累加结果。返回的算符已吸收相位、按 Z 再 X 向量确定性排序并合并重复项，且只删除精确零；不采用 `simplify()` 的默认截断。 + +| 接口 | 语义 | +| --- | --- | +| `identity(num_qubits) -> Self` | 指定寄存器上的单位算符；零比特时为标量 1。 | +| `from_terms(num_qubits, terms) -> Result` | 从 `(PauliString, Complex64)` 的 `IntoIterator` 安全构造；空输入保留声明宽度。 | +| `try_add(&self, other: &Self)` / `try_sub(&self, other: &Self)` | 同宽度加减，返回 `Result`，不修改输入。 | +| `try_mul(&self, other: &Self) -> Result` | 矩阵乘积 $AB$，右操作数先作用；允许非厄米结果。 | +| `scaled(&self, factor: Complex64) -> Result` | 返回缩放结果；零算符也验证因子。 | +| `adjoint(&self) -> Result` | 相位归一化后共轭系数，得到 $A^\dagger$。 | +| `normalized(&self) -> Result` | 验证、吸收相位、合并同类项，仅删除精确零。 | +| `chop(&mut self, atol: f64) -> Result<(), QisError>` | 先归并，再删除复数模长 `≤ atol` 的项；失败时原值不变。 | +| `equiv(&self, other: &Self, atol: f64) -> Result` | 判断归一化后 $A-B$ 的每个系数模长是否 `≤ atol`。 | +| `is_hermitian(&self, atol: f64) -> Result` | 等价于与伴随的 `equiv`，不投影到实数。 | +| `tensor(&self, other: &Self) -> Result` | $A\otimes B$，B 占低比特，A 占高比特。 | +| `remap_qubits(&self, mapping: &[usize], target_width: usize) -> Result` | `mapping[source] = destination`；支持置换和嵌入，未使用位置补 I。 | + +所有容差必须有限且非负。`equiv` 使用逐项系数容差，不是矩阵算符范数界；宽度不同返回 `QubitMismatch`，不会隐式扩展寄存器。`is_hermitian` 的容差衡量 $A-A^\dagger$，归一化后等于虚部幅值的两倍，与 `FermionicOp` 一致。已有 `normalized_hermitian()` 使用独立消费者契约：允许虚部幅值不超过 `1e-10`，将其投影到实数，但保留任何非零实系数。 + +`remap_qubits` 要求映射长度等于源宽度，目标位置不重复且小于 `target_width`，包括零/单位算符。它不表示费米子模式置换或部分迹。零算符保留寄存器宽度，稀疏运算支持超过 64 个量子比特；乘法可展开成两侧项数的乘积，仍可能产生大量项。 + +若张量或嵌入结果需要构造 Pauli 项，宽度超过 BitVec 的可表示长度时返回 `UnsupportedOperation`，不会触发长度断言。零算符不需要位向量，可以保留更大的声明宽度。 + +```rust +use cqlib_core::qis::{Hamiltonian, QisError}; +use num_complex::Complex64; + +let x = Hamiltonian::from_pauli("X".parse()?); +let z = Hamiltonian::from_pauli("Z".parse()?); +let sum = x.try_add(&z)?; +assert!(sum.try_mul(&sum)?.equiv(&Hamiltonian::identity(1).scaled(2.0.into())?, 0.0)?); + +let xz = x.try_mul(&z)?; // XZ = -iY +let minus_i_y = Hamiltonian::from_terms(1, [("Y".parse()?, Complex64::new(0.0, -1.0))])?; +assert!(xz.equiv(&minus_i_y, 0.0)?); +assert!(!xz.is_hermitian(0.0)?); + +// [X,Z]，以及 X ⊗ Z：最右字符对应比特 0。 +let commutator = xz.try_sub(&z.try_mul(&x)?)?; +assert!(commutator.equiv(&minus_i_y.scaled(2.0.into())?, 0.0)?); +assert!(x.tensor(&z)?.equiv(&Hamiltonian::from_pauli("XZ".parse()?), 0.0)?); +assert!(x.remap_qubits(&[2], 3)?.equiv(&Hamiltonian::from_pauli("XII".parse()?), 0.0)?); +# Ok::<(), QisError>(()) +``` + +旧 `Add`、`scale`、`simplify` 和构造器保留原契约；新调用优先使用安全接口。`PartialEq` 比较存储结构，项顺序、重复项或相位分布不同的等价算符可能不相等，应使用 `equiv`。 + ### 查询与转换 - `fn all_terms_commute(&self) -> bool`：判断所有 Pauli 项是否两两对易。为真时 $e^{-iHt}$ 可以按项分解为精确演化。 -- `fn to_matrix(&self) -> Result, QisError>`：返回 $2^N \times 2^N$ 稠密矩阵。任一项的比特数与 `num_qubits` 不一致时返回 `QisError::QubitMismatch`。空算量返回对应维度的零矩阵。 +- `fn to_matrix(&self) -> Result, QisError>`：验证并返回 $2^N \times 2^N$ 稠密矩阵。空算符返回对应维度的零矩阵，零比特输入返回 1×1 矩阵。分配前检查维度、元素数和字节数；大小溢出或内存预留失败返回 `UnsupportedOperation`。非法结构或非有限系数（包括矩阵元素累加溢出）返回相应错误。需要 $O(4^N)$ 内存，仅适用于小系统。 - `fn to_trotter_circuit(&self, time: f64, steps: usize, mode: TrotterMode) -> Result`：见 [Trotter 演化](8_evolution.md)。 - `fn to_evolution_circuit(&self, time: f64, steps: usize, mode: TrotterMode) -> Result`：见 [Trotter 演化](8_evolution.md)。 diff --git a/docs/api/en/python/3_qis/0_overview.md b/docs/api/en/python/3_qis/0_overview.md index e92b10c1..2564dd96 100644 --- a/docs/api/en/python/3_qis/0_overview.md +++ b/docs/api/en/python/3_qis/0_overview.md @@ -29,6 +29,8 @@ All four simulators provide a `from_circuit()` static method, which simulates a ### Observables and time evolution +[`FermionicOp` and `JordanWignerMapper`](10_fermion.md) provide general fermionic algebra, excitation generators and mappings to Pauli sums for the existing evolution components. + `PauliString` and `Hamiltonian` describe observables, and both can give expectation values; `Hamiltonian` can also be converted into an evolution circuit through Trotter decomposition, returning to the circuit space of `cqlib.circuit` for execution. --- diff --git a/docs/api/en/python/3_qis/10_fermion.md b/docs/api/en/python/3_qis/10_fermion.md new file mode 100644 index 00000000..9044bf9e --- /dev/null +++ b/docs/api/en/python/3_qis/10_fermion.md @@ -0,0 +1,70 @@ +# Fermionic Operators and Jordan–Wigner Mapping + +`cqlib.qis.fermion` exposes general fermionic algebra and sparse mappings backed by Rust Core. Public types are also exported from `cqlib.qis`. + +```python +from cqlib.qis import FermionicLadder, FermionicOp, JordanWignerMapper + +a = FermionicOp.annihilation(num_modes=4, mode=0) +assert (a * a.adjoint() + a.adjoint() * a).equiv(FermionicOp.identity(4)) +term = [FermionicLadder.create(2), FermionicLadder.annihilate(0)] +op = FermionicOp.from_term(num_modes=4, term=term, coefficient=0.5j) +hamiltonian = JordanWignerMapper().map_operator(op) +``` + +## Operators and Arithmetic + +`FermionicLadder.create(mode)` and `annihilate(mode)` construct immutable ladders with read-only `mode` and `kind`, `adjoint()`, value equality and hashing. The operator containing a ladder validates its mode against the declared mode space. + +`FermionicOp(num_modes)` constructs zero; `identity(num_modes)` constructs identity. `from_term(num_modes, term, coefficient=1.0)` accepts an ordered product. `from_terms(num_modes, terms)` accepts an iterable of `(term, coefficient)` tuples and merges duplicate products. An empty product is identity; an empty sum is zero. Both retain their declared mode count, including zero modes. + +| Interface | Behavior | +| --- | --- | +| `creation` / `annihilation` / `number` | Arguments are `num_modes, mode`; construct creation, annihilation or mode-number operators. | +| `total_number(num_modes)` | Construct total particle number. | +| `A + B` / `A - B` / `-A` | Return new operators; addition and subtraction require matching mode counts. | +| `A * B` | Operator product AB: B acts first, without implicit normal ordering. | +| `c * A` / `A * c` / `A.scaled(c)` | Scale by a finite integer, float or complex coefficient. | +| `adjoint()` | Conjugate coefficients, reverse products and exchange ladder kinds. | +| `normal_order()` | Return a new expansion: creators first, ascending modes within each group, including contractions and signs. | +| `add_term(term, coefficient=1.0)` | Mutate atomically and return `None`. | +| `chop(*, atol=0.0)` | Remove coefficients with complex magnitude <= tolerance in place; return `None`. | +| `equiv(other, *, atol=0.0)` | Compare algebraically after normal ordering; mismatched mode counts raise an error. | +| `is_hermitian(*, atol=0.0)` | Compare algebraically with the adjoint. | + +Read-only `num_modes` and `num_terms` describe storage. `terms` returns an independent snapshot of `(tuple[FermionicLadder, ...], complex)` pairs. Equality compares stored structure; `is_zero()` tests empty storage. Use `equiv()` for algebraic equality or zero. + +Only exact zeros are removed by default. Tolerances must be finite and nonnegative. `copy()` and Python copy protocols create independent objects. Operators are mutable and unhashable. `A += B` rebinds the result without modifying aliases of A. Use `sum(operators, start=FermionicOp(num_modes))` to sum operators. + +## Excitation Generators and Evolution + +`excitation_generator(num_modes, sources, targets)` constructs G = i(T − T†), with target creators in input order followed by source annihilators in reverse input order. Thus `[0, 1] -> [2, 3]` means T = a†₂ a†₃ a₁ a₀. Lists must be nonempty, equally sized, individually unique and disjoint. No implicit sorting or factorial normalization is applied. + +```python +from cqlib.circuit.ansatz import BasisEncoding, EvolutionStrategy, PauliEvolutionAnsatz + +mapper = JordanWignerMapper() +circuit = BasisEncoding(mapper.map_occupation([1, 1, 0, 0])).build_circuit("hf") +for name, sources, targets in [ + ("theta_double", [0, 1], [2, 3]), + ("theta_single", [0], [2]), +]: + generator = FermionicOp.excitation_generator(4, sources, targets) + block = (PauliEvolutionAnsatz(mapper.map_operator(generator)) + .with_strategy(EvolutionStrategy.exact()) + .with_time_param_name(name).build_circuit("excitation")) + circuit.compose(block) +bound = circuit.assign_parameters({"theta_double": 0.2, "theta_single": -0.3}) +``` + +Each block implements exp(−iθG) = exp(θ(T − T†)). Build and name separate blocks for independent parameters; evolving a sum of generators uses one shared parameter. + +## Mappings and Extension + +`JordanWignerMapper()` provides `map_operator(op) -> Hamiltonian`, `num_qubits(num_modes) -> int` and `map_occupation(occupation) -> list[bool]`. Operator mapping accepts non-Hermitian, identity and zero operators; downstream evolution and expectation consumers impose Hermiticity requirements. + +Mode p maps to qubit p, with a Z chain on lower-numbered modes. Qubit 0 is the least significant bit. Occupations must be integer zero/one or booleans; returned `bits[p]` addresses qubit p. Empty input returns an empty list. Core specifies no chemical or alpha/beta orbital ordering. + +`FermionToQubitMapper` and `OccupationStateMapper` are independent runtime-checkable `typing.Protocol` interfaces. Custom implementations need matching methods without inheriting a binding class. Runtime checks establish member presence only, not signatures or correctness. Algorithms call Python mappers directly; no automatic Rust trait callback bridge is provided. + +Wrong types raise `TypeError`; negative or out-of-range modes raise `IndexError`. Negative mode counts, mismatched spaces, nonfinite coefficients, invalid occupations or tolerances raise `ValueError`. Mode integers outside the Rust integer range raise `OverflowError`. diff --git a/docs/api/en/python/3_qis/7_hamiltonian.md b/docs/api/en/python/3_qis/7_hamiltonian.md index 98adcc34..9286c66c 100644 --- a/docs/api/en/python/3_qis/7_hamiltonian.md +++ b/docs/api/en/python/3_qis/7_hamiltonian.md @@ -223,3 +223,39 @@ assert trotterized.num_qubits == 1 | `ValueError` | `expectation_probs` cannot derive the target operator from the given measurement bases. | The input types of observable computation are given in [Statevector](1_statevector.md) and [DensityMatrix](2_density_matrix.md); the definition of Trotter modes is given in [Trotter evolution](8_evolution.md). + +## Validated Sparse Algebra + +`from_terms(num_qubits, terms)` declares the width explicitly and accepts `(PauliString, int/float/complex)` pairs. It validates widths and finite coefficients, absorbs Pauli phases and merges duplicates, removing only exact zeros. Empty input retains the declared width. + +| Interface | Behavior | +| --- | --- | +| `normalized()` | Return a validated, merged copy retaining every nonzero coefficient. | +| `add(other)` | Return a validated, merged sum on the same register. | +| `A - B`, `A * B`, `-A` | Checked subtraction, operator multiplication and negation; B acts first in AB. | +| `c * A`, `A * c`, `scaled(c)` | Validate the factor and resulting coefficients; return a new object. | +| `adjoint()` | Return the conjugate transpose, accounting for Pauli phases. | +| `chop(*, atol=0.0)` | Normalize and truncate by complex magnitude in place; errors preserve the object. | +| `equiv(other, *, atol=0.0)` | Compare normalized coefficient differences; width mismatch raises an error. Tolerance is not a matrix operator-norm bound. | +| `is_hermitian(*, atol=0.0)` | Compare with the adjoint: normalized coefficients satisfy 2*abs(imag) <= atol. | +| `normalized_hermitian()` | Project imaginary residuals <= 1e-10 to zero; reject larger ones and retain all nonzero real parts. | +| `tensor(other)` | Return self tensor other; other occupies low qubits. | +| `remap_qubits(mapping, target_width)` | Embed mapping[source] = destination; destinations must be unique and in range. Unused targets are identity. | + +General algebra accepts non-Hermitian operators. Evolution and expectation consumers enforce Hermiticity. Tolerances must be finite and nonnegative. Qubit embedding is not a fermionic mode permutation. + +```python +from cqlib.qis import Hamiltonian, PauliString + +x = Hamiltonian.from_terms(1, [(PauliString.from_str("X"), 1.0)]) +z = Hamiltonian.from_terms(1, [(PauliString.from_str("Z"), 1.0)]) +assert (x * z).equiv(Hamiltonian.from_terms(1, [(PauliString.from_str("Y"), -1j)])) +assert x.add(x).equiv(2 * x) +assert x.remap_qubits([2], target_width=3).num_qubits == 3 +``` + +### Legacy Compatibility + +`from_list`, `add_term`, `scale`, `+` and `+=` keep their existing behavior. Addition concatenates stored terms; `scale` mutates in place. Legacy `simplify()` truncates near-zero coefficients. Use `normalized()` to retain small mapped coefficients, or `chop(atol=...)` for explicit truncation. + +`add_term()` and `scale()` return `None` and do not support chaining. New algebra operations raise `ValueError` on mismatched widths, nonfinite coefficients or intermediate coefficient overflow. diff --git a/docs/api/en/rust/3_qis/7_hamiltonian.md b/docs/api/en/rust/3_qis/7_hamiltonian.md index 44318a14..89721379 100644 --- a/docs/api/en/rust/3_qis/7_hamiltonian.md +++ b/docs/api/en/rust/3_qis/7_hamiltonian.md @@ -3,9 +3,9 @@ `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. +This page covers the operator `Hamiltonian`, formed by Pauli strings and complex coefficients, and the unified `Observable` trait for computing real expectation values. -`Hamiltonian` is a sparse representation of a $2^N \times 2^N$ matrix, of the form $H = \sum_k c_k P_k$, where $c_k$ is a complex coefficient and $P_k$ is an $N$-qubit Pauli string. The same structure serves both as the Hamiltonian of a physical interpretation and as an arbitrary observable whose expectation value is needed. +`Hamiltonian` is a sparse representation of a $2^N \times 2^N$ matrix, of the form $H = \sum_k c_k P_k$, where $c_k$ is a complex coefficient and $P_k$ is an $N$-qubit Pauli string. The type permits non-Hermitian operators; real expectation and unitary evolution methods additionally check Hermiticity. ## Import @@ -46,10 +46,59 @@ Fields: - `fn simplify(&mut self)`: simplify the term list. The operation has two steps: absorb the internal phase of each `PauliString` into the coefficient, normalizing all Pauli strings to `Phase::Plus`; then group by the X and Z component vectors and merge the coefficients of identical operators. Coefficients whose absolute value is smaller than $10^{-10}$ are discarded. - `fn scale(&mut self, factor: Complex64)`: multiply the coefficients of all terms by the same complex factor. +### Checked Pauli-sum algebra + +The fallible interfaces below validate register widths, X/Z bit-vector lengths, finite coefficients and coefficient products/accumulations. Returned operators have phases absorbed, terms sorted deterministically by Z then X vectors, duplicates merged and only exact zeros removed. They do not use the default pruning threshold of `simplify()`. + +| Interface | Meaning | +| --- | --- | +| `identity(num_qubits) -> Self` | Identity on the declared register; scalar one at zero qubits. | +| `from_terms(num_qubits, terms) -> Result` | Checked construction from an `IntoIterator` of `(PauliString, Complex64)`; empty input retains the declared width. | +| `try_add(&self, other: &Self)` / `try_sub(&self, other: &Self)` | Same-width addition/subtraction returning `Result`, without modifying inputs. | +| `try_mul(&self, other: &Self) -> Result` | Matrix product $AB$: the right operand acts first. Results may be non-Hermitian. | +| `scaled(&self, factor: Complex64) -> Result` | Returns a scaled operator; checks the factor even for zero input. | +| `adjoint(&self) -> Result` | Conjugates coefficients after phase normalization to obtain $A^\dagger$. | +| `normalized(&self) -> Result` | Validates, absorbs phases, merges terms and removes only exact zeros. | +| `chop(&mut self, atol: f64) -> Result<(), QisError>` | Merges first, then removes coefficients with complex magnitude `≤ atol`; leaves the original unchanged on failure. | +| `equiv(&self, other: &Self, atol: f64) -> Result` | Checks whether every normalized coefficient of $A-B$ has magnitude `≤ atol`. | +| `is_hermitian(&self, atol: f64) -> Result` | Equivalent to `equiv` against the adjoint; does not project coefficients to real values. | +| `tensor(&self, other: &Self) -> Result` | $A\otimes B$: B occupies low qubits, A high qubits. | +| `remap_qubits(&self, mapping: &[usize], target_width: usize) -> Result` | `mapping[source] = destination`; permutations and embeddings insert I at unused positions. | + +Tolerances must be finite and nonnegative. `equiv` uses a per-coefficient tolerance, not a matrix operator-norm bound. Different widths return `QubitMismatch`; registers are never implicitly expanded. `is_hermitian` measures $A-A^\dagger$, bounding twice the imaginary coefficient magnitude after normalization, as for `FermionicOp`. The existing `normalized_hermitian()` has a separate consumer contract: imaginary magnitudes up to `1e-10` are projected to real values, while every nonzero real coefficient is retained. + +`remap_qubits` requires a mapping of the source width with distinct destinations below `target_width`, including for zero/identity operators. It is not a fermionic mode permutation or a partial trace. Zero operators retain register widths. Sparse operations support more than 64 qubits; multiplication can still expand to the product of the input term counts. + +When a tensor product or embedding requires Pauli terms, widths beyond BitVec's representable length return `UnsupportedOperation` before its length assertion can panic. Zero operators require no bit vectors and can retain larger declared widths. + +```rust +use cqlib_core::qis::{Hamiltonian, QisError}; +use num_complex::Complex64; + +let x = Hamiltonian::from_pauli("X".parse()?); +let z = Hamiltonian::from_pauli("Z".parse()?); +let sum = x.try_add(&z)?; +assert!(sum.try_mul(&sum)?.equiv(&Hamiltonian::identity(1).scaled(2.0.into())?, 0.0)?); + +let xz = x.try_mul(&z)?; // XZ = -iY +let minus_i_y = Hamiltonian::from_terms(1, [("Y".parse()?, Complex64::new(0.0, -1.0))])?; +assert!(xz.equiv(&minus_i_y, 0.0)?); +assert!(!xz.is_hermitian(0.0)?); + +// [X,Z], and X ⊗ Z: the rightmost character is qubit zero. +let commutator = xz.try_sub(&z.try_mul(&x)?)?; +assert!(commutator.equiv(&minus_i_y.scaled(2.0.into())?, 0.0)?); +assert!(x.tensor(&z)?.equiv(&Hamiltonian::from_pauli("XZ".parse()?), 0.0)?); +assert!(x.remap_qubits(&[2], 3)?.equiv(&Hamiltonian::from_pauli("XII".parse()?), 0.0)?); +# Ok::<(), QisError>(()) +``` + +Legacy `Add`, `scale`, `simplify` and constructors retain their contracts; prefer the checked interfaces for new callers. `PartialEq` compares stored structure: operators with different term orders, duplicates or phase distributions may be algebraically equivalent but unequal. Use `equiv` for that comparison. + ### Queries and conversion - `fn all_terms_commute(&self) -> bool`: determine whether all Pauli terms commute pairwise. When true, $e^{-iHt}$ can be decomposed term by term into an exact evolution. -- `fn to_matrix(&self) -> Result, QisError>`: return a dense $2^N \times 2^N$ matrix. When the qubit count of any term is inconsistent with `num_qubits`, `QisError::QubitMismatch` is returned. An empty operator returns the zero matrix of the corresponding dimension. +- `fn to_matrix(&self) -> Result, QisError>`: validate and return a dense $2^N \times 2^N$ matrix. Empty operators produce the corresponding zero matrix; zero-qubit input produces a 1×1 matrix. Dimensions, element count and byte count are checked before allocation. Size overflow or reservation failure returns `UnsupportedOperation`; malformed terms and nonfinite coefficients (including matrix-entry accumulation overflow) return the corresponding errors. Requires $O(4^N)$ memory and is intended for small systems. - `fn to_trotter_circuit(&self, time: f64, steps: usize, mode: TrotterMode) -> Result`: see [Trotter Evolution](8_evolution.md). - `fn to_evolution_circuit(&self, time: f64, steps: usize, mode: TrotterMode) -> Result`: see [Trotter Evolution](8_evolution.md).