From ec1b9d6f0df0cd5934ede45de9fbdfacdf1c66ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?=E5=81=A5=E5=81=A5?= Date: Tue, 13 Oct 2026 11:01:55 +0800 Subject: [PATCH] feat(c): expose fermionic and checked Pauli algebra --- Cargo.lock | 1 + crates/binding-c/Cargo.toml | 1 + crates/binding-c/README.md | 13 +- crates/binding-c/src/qis/algebra.rs | 80 +++ crates/binding-c/src/qis/fermion.rs | 599 +++++++++++++++++ crates/binding-c/src/qis/hamiltonian.rs | 20 +- .../binding-c/src/qis/hamiltonian_algebra.rs | 302 +++++++++ crates/binding-c/src/qis/mod.rs | 5 + crates/binding-c/src/qis/pauli.rs | 30 +- crates/binding-c/tests/qis_algebra_test.rs | 630 ++++++++++++++++++ crates/binding-c/tests/test_qis_algebra.c | 249 +++++++ docs/api/cn/c/3_qis/0_overview.md | 2 + docs/api/cn/c/3_qis/10_fermion.md | 179 +++++ docs/api/cn/c/3_qis/6_pauli.md | 12 +- docs/api/cn/c/3_qis/7_hamiltonian.md | 92 ++- docs/api/en/c/3_qis/0_overview.md | 2 + docs/api/en/c/3_qis/10_fermion.md | 179 +++++ docs/api/en/c/3_qis/6_pauli.md | 12 +- docs/api/en/c/3_qis/7_hamiltonian.md | 92 ++- tools/c-sdk/CMakeLists.txt | 55 +- tools/c-sdk/README.md | 9 +- tools/ci.py | 3 +- 22 files changed, 2512 insertions(+), 55 deletions(-) create mode 100644 crates/binding-c/src/qis/algebra.rs create mode 100644 crates/binding-c/src/qis/fermion.rs create mode 100644 crates/binding-c/src/qis/hamiltonian_algebra.rs create mode 100644 crates/binding-c/tests/qis_algebra_test.rs create mode 100644 crates/binding-c/tests/test_qis_algebra.c create mode 100644 docs/api/cn/c/3_qis/10_fermion.md create mode 100644 docs/api/en/c/3_qis/10_fermion.md diff --git a/Cargo.lock b/Cargo.lock index 0dfc9dcb..3bf21508 100644 --- a/Cargo.lock +++ b/Cargo.lock @@ -167,6 +167,7 @@ checksum = "72b3254f16251a8381aa12e40e3c4d2f0199f8c6508fbecb9d91f575e0fbb8c6" name = "binding-c" version = "0.1.0-beta.4" dependencies = [ + "bitvec", "cbindgen", "cqlib-core", "libc", diff --git a/crates/binding-c/Cargo.toml b/crates/binding-c/Cargo.toml index 35b237cc..312365b9 100644 --- a/crates/binding-c/Cargo.toml +++ b/crates/binding-c/Cargo.toml @@ -11,6 +11,7 @@ crate-type = ["rlib", "staticlib", "cdylib"] [dependencies] cqlib-core = { path = "../cqlib-core", version = "0.1.0-beta.4" } libc = "0.2" +bitvec.workspace = true ndarray.workspace = true num-complex.workspace = true rustworkx-core.workspace = true diff --git a/crates/binding-c/README.md b/crates/binding-c/README.md index d5904499..a39a669e 100644 --- a/crates/binding-c/README.md +++ b/crates/binding-c/README.md @@ -121,4 +121,15 @@ gcc -I crates/binding-c/include \ -o test_circuit ./test_circuit # prints "binding-c circuit tests passed" -``` \ No newline at end of file +``` + +## Checked sparse algebra + +The new `fermionic_op_*`, `jordan_wigner_*`, and checked `hamiltonian_*` functions +expose fermionic products, excitation generators and sparse Pauli arithmetic. +They return an error code and deliver independent objects through output pointers; +errors leave outputs unchanged. Inputs are borrowed and coefficients are never +approximately pruned unless requested. See [Fermionic algebra](../../docs/api/en/c/3_qis/10_fermion.md) +and [Hamiltonian](../../docs/api/en/c/3_qis/7_hamiltonian.md). + +The staged SDK consumer tests include `test_qis_algebra.c` for both static and shared linking. diff --git a/crates/binding-c/src/qis/algebra.rs b/crates/binding-c/src/qis/algebra.rs new file mode 100644 index 00000000..cf20f8ec --- /dev/null +++ b/crates/binding-c/src/qis/algebra.rs @@ -0,0 +1,80 @@ +// 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. + +//! Internal helpers for the checked algebra ABI. Legacy QIS error mapping is unchanged. + +use cqlib_core::qis::QisError; + +pub(crate) fn algebra_err_code(err: &QisError) -> i32 { + match err { + QisError::ModeOutOfBounds { .. } | QisError::IndexOutOfBounds { .. } => -2, + QisError::ModeMismatch { .. } + | QisError::QubitMismatch { .. } + | QisError::DimensionMismatch { .. } + | QisError::NonFiniteCoefficient { .. } + | QisError::InvalidOccupation { .. } + | QisError::InvalidParameterValue(_) + | QisError::NotHermitian => -8, + _ => super::qis_err_code(err), + } +} + +/// C callers must provide valid, aligned storage, live handles and exclusive +/// mutable outputs. NULL is allowed for empty arrays. Check the Rust slice size +/// bound before constructing a slice, including for malicious length values. +pub(crate) unsafe fn input_slice<'a, T>(ptr: *const T, len: usize) -> Result<&'a [T], i32> { + if len == 0 { + return Ok(&[]); + } + if ptr.is_null() { + return Err(-1); + } + if len > (isize::MAX as usize) / std::mem::size_of::() { + return Err(-8); + } + Ok(unsafe { std::slice::from_raw_parts(ptr, len) }) +} + +/// Only called after validating the output pointer. Write once on success. +pub(crate) unsafe fn write_result(out: *mut T, result: Result) -> i32 { + match result { + Ok(value) => { + unsafe { *out = value }; + 0 + } + Err(err) => algebra_err_code(&err), + } +} + +/// Check the same BitVec materialization limit used by core's checked algebra. +pub(crate) fn check_pauli_width(width: usize) -> Result<(), i32> { + if width > bitvec::slice::BitSlice::::MAX_BITS { + Err(-7) + } else { + Ok(()) + } +} + +/// Return 2^width only when the side length is representable. Do not truncate +/// usize widths to u32: zero sums may legally declare widths as large as usize::MAX. +pub(crate) fn dense_matrix_dim(width: usize) -> Option { + 1usize.checked_shl(u32::try_from(width).ok()?) +} + +/// Check element count and signed addressable byte size before allocation or +/// writing C buffers, matching Hamiltonian::to_matrix's allocation bounds. +pub(crate) fn dense_matrix_len(width: usize) -> Option { + let dim = dense_matrix_dim(width)?; + let len = dim.checked_mul(dim)?; + let bytes = len.checked_mul(std::mem::size_of::())?; + (bytes <= isize::MAX as usize).then_some(len) +} diff --git a/crates/binding-c/src/qis/fermion.rs b/crates/binding-c/src/qis/fermion.rs new file mode 100644 index 00000000..b5784c2b --- /dev/null +++ b/crates/binding-c/src/qis/fermion.rs @@ -0,0 +1,599 @@ +// 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. + +//! Checked fermionic algebra. All outputs and mutable inputs are unchanged on +//! failure. Object outputs are owned copies freed with `fermionic_op_free`; +//! inputs are borrowed. Multiplication applies the right operand first. + +use super::CHamiltonian; +use super::algebra::{algebra_err_code, check_pauli_width, input_slice, write_result}; +use cqlib_core::qis::fermion::mapping::{ + FermionToQubitMapper, JordanWignerMapper, OccupationStateMapper, +}; +use cqlib_core::qis::{FermionicLadder, FermionicOp}; +use libc::size_t; +use num_complex::Complex64; + +/// Opaque owned fermionic operator. No internal storage is exposed to C. +/// All checked functions return 0 on success and leave outputs unchanged on error. +/// Pass empty object output slots; results must be freed with the matching free +/// function. Caller storage must be live, aligned and large enough; output slots +/// must not overlap one another or handle internals. Inputs are borrowed. +/// Do not read a handle concurrently with mutation. Empty arrays permit NULL. +pub struct CFermionicOp { + pub inner: FermionicOp, +} + +/// Creation ladder tag. +pub const FERMIONIC_CREATE: u8 = 0; +/// Annihilation ladder tag. +pub const FERMIONIC_ANNIHILATE: u8 = 1; + +/// One ladder on a zero-based mode. Unknown kind tags return -8. +#[repr(C)] +#[derive(Clone, Copy)] +pub struct CFermionicLadder { + pub mode: size_t, + pub kind: u8, +} + +/// Borrowed ordered product and coefficient; the rightmost ladder acts first. +/// A NULL ladders pointer with len=0 denotes the identity product. +#[repr(C)] +pub struct CFermionicTerm { + pub ladders: *const CFermionicLadder, + pub len: size_t, + pub coefficient: Complex64, +} + +unsafe fn ladders(ptr: *const CFermionicLadder, len: size_t) -> Result, i32> { + unsafe { input_slice(ptr, len) }? + .iter() + .map(|ladder| match ladder.kind { + FERMIONIC_CREATE => Ok(FermionicLadder::Create(ladder.mode)), + FERMIONIC_ANNIHILATE => Ok(FermionicLadder::Annihilate(ladder.mode)), + _ => Err(-8), + }) + .collect() +} + +fn owned(inner: FermionicOp) -> *mut CFermionicOp { + Box::into_raw(Box::new(CFermionicOp { inner })) +} + +/// Release an owned operator. NULL is allowed. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_free(ptr: *mut CFermionicOp) { + if !ptr.is_null() { + unsafe { + drop(Box::from_raw(ptr)); + } + } +} + +/// Construct zero, retaining the declared width including zero modes. +/// Returns -1 for NULL out; writes an owned handle only on success. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_new(num_modes: size_t, out: *mut *mut CFermionicOp) -> i32 { + if out.is_null() { + return -1; + } + unsafe { + *out = owned(FermionicOp::new(num_modes)); + } + 0 +} + +/// Construct identity, retaining the declared width including zero modes. +/// Returns -1 for NULL out; writes an owned handle only on success. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_identity(num_modes: size_t, out: *mut *mut CFermionicOp) -> i32 { + if out.is_null() { + return -1; + } + unsafe { + *out = owned(FermionicOp::identity(num_modes)); + } + 0 +} + +/// Construct total particle number, retaining the declared width including zero modes. +/// Returns -1 for NULL out; writes an owned handle only on success. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_total_number(num_modes: size_t, out: *mut *mut CFermionicOp) -> i32 { + if out.is_null() { + return -1; + } + unsafe { + *out = owned(FermionicOp::total_number(num_modes)); + } + 0 +} + +/// Construct a mode creation operator. Returns -1 for NULL out, -2 for invalid mode. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_creation( + num_modes: size_t, + mode: size_t, + out: *mut *mut CFermionicOp, +) -> i32 { + if out.is_null() { + return -1; + } + unsafe { write_result(out, FermionicOp::creation(num_modes, mode).map(owned)) } +} + +/// Construct a mode annihilation operator. Returns -1 for NULL out, -2 for invalid mode. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_annihilation( + num_modes: size_t, + mode: size_t, + out: *mut *mut CFermionicOp, +) -> i32 { + if out.is_null() { + return -1; + } + unsafe { write_result(out, FermionicOp::annihilation(num_modes, mode).map(owned)) } +} + +/// Construct a mode number operator. Returns -1 for NULL out, -2 for invalid mode. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_number( + num_modes: size_t, + mode: size_t, + out: *mut *mut CFermionicOp, +) -> i32 { + if out.is_null() { + return -1; + } + unsafe { write_result(out, FermionicOp::number(num_modes, mode).map(owned)) } +} + +/// Construct a weighted ordered product. Empty product is identity; even +/// zero coefficients validate all modes. Returns -1 NULL, -2 mode, -8 tag/coefficient. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_from_term( + num_modes: size_t, + term: *const CFermionicLadder, + len: size_t, + coefficient: Complex64, + out: *mut *mut CFermionicOp, +) -> i32 { + if out.is_null() { + return -1; + } + let term = match unsafe { ladders(term, len) } { + Ok(v) => v, + Err(e) => return e, + }; + unsafe { + write_result( + out, + FermionicOp::from_term(num_modes, term, coefficient).map(owned), + ) + } +} + +/// Construct a sum of borrowed terms, merging identical products and removing +/// only exact zeros. NULL terms with len=0 creates zero. No input is consumed. +/// Returns -1 NULL, -2 invalid mode, -8 invalid tag/coefficient or length. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_from_terms( + num_modes: size_t, + terms: *const CFermionicTerm, + len: size_t, + out: *mut *mut CFermionicOp, +) -> i32 { + if out.is_null() { + return -1; + } + let terms = match unsafe { input_slice(terms, len) } { + Ok(v) => v, + Err(e) => return e, + }; + let mut result = FermionicOp::new(num_modes); + for term in terms { + let product = match unsafe { ladders(term.ladders, term.len) } { + Ok(v) => v, + Err(e) => return e, + }; + if let Err(err) = result.add_term(product, term.coefficient) { + return algebra_err_code(&err); + } + } + unsafe { + *out = owned(result); + } + 0 +} + +/// Add an ordered term atomically. Inputs are copied. Empty product is identity. +/// Returns -1 NULL, -2 mode, -8 tag/nonfinite coefficient (including merge overflow). +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_add_term( + ptr: *mut CFermionicOp, + term: *const CFermionicLadder, + len: size_t, + coefficient: Complex64, +) -> i32 { + if ptr.is_null() { + return -1; + } + let term = match unsafe { ladders(term, len) } { + Ok(v) => v, + Err(e) => return e, + }; + match unsafe { &mut (*ptr).inner }.add_term(term, coefficient) { + Ok(()) => 0, + Err(err) => algebra_err_code(&err), + } +} + +/// Construct G=i(T-T†). Target creators retain input order; source annihilators +/// use reverse input order. No sorting or factorial. Returns -1 NULL, -2 mode, +/// -8 empty/unequal lists, duplicates or overlap. Inputs are borrowed. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_excitation_generator( + num_modes: size_t, + sources: *const size_t, + sources_len: size_t, + targets: *const size_t, + targets_len: size_t, + out: *mut *mut CFermionicOp, +) -> i32 { + if out.is_null() { + return -1; + } + let sources = match unsafe { input_slice(sources, sources_len) } { + Ok(v) => v, + Err(e) => return e, + }; + let targets = match unsafe { input_slice(targets, targets_len) } { + Ok(v) => v, + Err(e) => return e, + }; + unsafe { + write_result( + out, + FermionicOp::excitation_generator(num_modes, sources, targets).map(owned), + ) + } +} + +/// Read declared mode count. Returns -1 for NULL input or output. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_num_modes(ptr: *const CFermionicOp, out: *mut size_t) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { + *out = inner.num_modes(); + } + 0 +} + +/// Read stored product count. Returns -1 for NULL input or output. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_num_terms(ptr: *const CFermionicOp, out: *mut size_t) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { + *out = inner.num_terms(); + } + 0 +} + +/// Read structural zero (0/1); use equiv for algebraic zero. Returns -1 for NULL input or output. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_is_zero(ptr: *const CFermionicOp, out: *mut u8) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { + *out = u8::from(inner.is_zero()); + } + 0 +} + +/// Return an owned independent copy. +/// Returns -1 for NULL input or output. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_copy(ptr: *const CFermionicOp, out: *mut *mut CFermionicOp) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { write_result(out, Ok(owned(inner.clone()))) } +} + +/// Return an owned conjugate transpose, reversing products and exchanging ladder kinds. +/// Returns -1 for NULL input or output. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_adjoint( + ptr: *const CFermionicOp, + out: *mut *mut CFermionicOp, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { write_result(out, Ok(owned(inner.adjoint()))) } +} + +/// Return an owned normal-ordered expansion with contraction terms and fermionic signs. +/// Returns -1 NULL; normal ordering may return -8 on coefficient overflow. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_normal_order( + ptr: *const CFermionicOp, + out: *mut *mut CFermionicOp, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { write_result(out, inner.normal_order().map(owned)) } +} + +/// Checked add, returning an independent operator without normal ordering. +/// Returns -1 NULL, -8 mode mismatch/nonfinite arithmetic. Right operand acts first in mul. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_add( + left: *const CFermionicOp, + right: *const CFermionicOp, + out: *mut *mut CFermionicOp, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*left).inner.try_add(&(*right).inner).map(owned)) } +} + +/// Checked sub, returning an independent operator without normal ordering. +/// Returns -1 NULL, -8 mode mismatch/nonfinite arithmetic. Right operand acts first in mul. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_sub( + left: *const CFermionicOp, + right: *const CFermionicOp, + out: *mut *mut CFermionicOp, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*left).inner.try_sub(&(*right).inner).map(owned)) } +} + +/// Checked mul, returning an independent operator without normal ordering. +/// Returns -1 NULL, -8 mode mismatch/nonfinite arithmetic. Right operand acts first in mul. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_mul( + left: *const CFermionicOp, + right: *const CFermionicOp, + out: *mut *mut CFermionicOp, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*left).inner.try_mul(&(*right).inner).map(owned)) } +} + +/// Return a scaled copy. Returns -1 NULL, -8 nonfinite factor or result. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_scaled( + ptr: *const CFermionicOp, + factor: Complex64, + out: *mut *mut CFermionicOp, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*ptr).inner.scaled(factor).map(owned)) } +} + +/// Remove coefficients with magnitude <= atol in place. Returns -1 NULL, +/// -8 for negative/nonfinite atol. Failure leaves the operator unchanged. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_chop(ptr: *mut CFermionicOp, atol: f64) -> i32 { + if ptr.is_null() { + return -1; + } + match unsafe { &mut (*ptr).inner }.chop(atol) { + Ok(()) => 0, + Err(e) => algebra_err_code(&e), + } +} + +/// Structural equality including mode count; writes 0 or 1. Returns -1 NULL. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_equals( + left: *const CFermionicOp, + right: *const CFermionicOp, + out: *mut u8, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { + *out = u8::from((*left).inner == (*right).inner); + } + 0 +} + +/// Algebraic equivalence after normal ordering; writes 0 or 1. Returns -1 NULL, +/// -8 mode mismatch, invalid atol or arithmetic overflow. atol=0 is exact. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_equiv( + left: *const CFermionicOp, + right: *const CFermionicOp, + atol: f64, + out: *mut u8, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { + write_result( + out, + (*left).inner.equiv(&(*right).inner, atol).map(u8::from), + ) + } +} + +/// Algebraic Hermiticity test; writes 0 or 1. Returns -1 NULL, -8 invalid atol/overflow. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_is_hermitian( + ptr: *const CFermionicOp, + atol: f64, + out: *mut u8, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*ptr).inner.is_hermitian(atol).map(u8::from)) } +} + +/// Read one stored product's length in deterministic structural order. +/// Returns -1 NULL, -2 invalid term index. No internal storage is borrowed. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_term_len( + ptr: *const CFermionicOp, + index: size_t, + out: *mut size_t, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let Some((term, _)) = (unsafe { &(*ptr).inner }).terms().nth(index) else { + return -2; + }; + unsafe { + *out = term.len(); + } + 0 +} + +/// Copy a term's ladders and coefficient in deterministic structural order. +/// Returns -1 NULL required pointer, -2 index, -8 insufficient capacity. +/// NULL buffer is allowed for an empty product. Writes no partial output on error. +#[unsafe(no_mangle)] +pub extern "C" fn fermionic_op_get_term( + ptr: *const CFermionicOp, + index: size_t, + buffer: *mut CFermionicLadder, + capacity: size_t, + out_coefficient: *mut Complex64, +) -> i32 { + if ptr.is_null() || out_coefficient.is_null() { + return -1; + } + let Some((term, coefficient)) = (unsafe { &(*ptr).inner }).terms().nth(index) else { + return -2; + }; + if capacity < term.len() { + return -8; + } + if !term.is_empty() && buffer.is_null() { + return -1; + } + for (i, ladder) in term.iter().enumerate() { + let kind = match ladder { + FermionicLadder::Create(_) => FERMIONIC_CREATE, + FermionicLadder::Annihilate(_) => FERMIONIC_ANNIHILATE, + }; + unsafe { + *buffer.add(i) = CFermionicLadder { + mode: ladder.mode(), + kind, + }; + } + } + unsafe { + *out_coefficient = coefficient; + } + 0 +} + +/// Map a borrowed operator to an owned CHamiltonian; free with hamiltonian_free. +/// Non-Hermitian inputs are valid. Mode p maps to qubit p, with Z on lower modes. +/// Returns -1 NULL, -7 unrepresentable width, -8 nonfinite arithmetic. +#[unsafe(no_mangle)] +pub extern "C" fn jordan_wigner_map_operator( + ptr: *const CFermionicOp, + out: *mut *mut CHamiltonian, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + // Zero operators need no bit vectors. All other inputs materialize a Pauli + // identity before mapping; guard that allocation like checked core algebra. + if !inner.is_zero() + && let Err(e) = check_pauli_width(inner.num_modes()) + { + return e; + } + unsafe { + write_result( + out, + JordanWignerMapper + .map_operator(inner) + .map(|inner| Box::into_raw(Box::new(CHamiltonian { inner }))), + ) + } +} + +/// Jordan-Wigner uses one qubit per mode. Returns -1 for NULL out. +#[unsafe(no_mangle)] +pub extern "C" fn jordan_wigner_num_qubits(num_modes: size_t, out: *mut size_t) -> i32 { + if out.is_null() { + return -1; + } + unsafe { + *out = JordanWignerMapper.num_qubits(num_modes); + } + 0 +} + +/// Encode occupations as 0/1 bytes: bits[p] addresses qubit p (LSB is qubit 0). +/// Empty input is vacuum. Returns -1 NULL required buffer, -8 invalid occupations, +/// length or insufficient capacity. Buffers may overlap: input is read before output. +#[unsafe(no_mangle)] +pub extern "C" fn jordan_wigner_map_occupation( + occupation: *const u8, + len: size_t, + out_bits: *mut u8, + capacity: size_t, +) -> i32 { + if len > 0 && (occupation.is_null() || out_bits.is_null()) { + return -1; + } + if capacity < len { + return -8; + } + let occupation = match unsafe { input_slice(occupation, len) } { + Ok(v) => v, + Err(e) => return e, + }; + let bits = match JordanWignerMapper.map_occupation(occupation) { + Ok(v) => v, + Err(e) => return algebra_err_code(&e), + }; + for (i, bit) in bits.into_iter().enumerate() { + unsafe { + *out_bits.add(i) = u8::from(bit); + } + } + 0 +} diff --git a/crates/binding-c/src/qis/hamiltonian.rs b/crates/binding-c/src/qis/hamiltonian.rs index 25a18be3..dec76e59 100644 --- a/crates/binding-c/src/qis/hamiltonian.rs +++ b/crates/binding-c/src/qis/hamiltonian.rs @@ -14,6 +14,7 @@ use crate::circuit::CCircuit; use crate::error::CqlibError; +use crate::qis::algebra::{dense_matrix_dim, dense_matrix_len}; use crate::qis::{CHamiltonian, CPauliString, CStatevector, qis_err_code}; use cqlib_core::qis::evolution::TrotterMode; use cqlib_core::qis::{Hamiltonian, Observable}; @@ -141,19 +142,22 @@ pub extern "C" fn hamiltonian_scale(ptr: *mut CHamiltonian, factor_re: f64, fact 0 } -/// Return the side length of the dense matrix (`2^N`), or 0 for NULL. +/// Return the side length of the dense matrix (`2^N`), or 0 for NULL or +/// an unrepresentable side length. A representable side does not guarantee +/// that the complete dense matrix can be allocated. #[unsafe(no_mangle)] pub extern "C" fn hamiltonian_to_matrix_len(ptr: *const CHamiltonian) -> usize { if ptr.is_null() { return 0; } - 1usize << unsafe { (*ptr).inner.num_qubits } + dense_matrix_dim(unsafe { (*ptr).inner.num_qubits }).unwrap_or(0) } /// Copy the dense `2^N x 2^N` matrix into `buffer` (row-major, /// `Complex64` values). `len` must equal /// `hamiltonian_to_matrix_len(ptr) * hamiltonian_to_matrix_len(ptr)`. -/// Returns 0 on success, or a negative error code. +/// Returns 0 on success, -7 for unrepresentable dense allocation size or +/// allocation failure, -8 for an incorrect buffer length. Errors leave the buffer unchanged. #[unsafe(no_mangle)] pub extern "C" fn hamiltonian_to_matrix( ptr: *const CHamiltonian, @@ -164,6 +168,12 @@ pub extern "C" fn hamiltonian_to_matrix( return CqlibError::NullPtr as i32; } let wrapper = unsafe { &*ptr }; + let Some(expected_len) = dense_matrix_len(wrapper.inner.num_qubits) else { + return CqlibError::SimulationError as i32; + }; + if len != expected_len { + return CqlibError::InvalidParam as i32; + } let matrix = match wrapper.inner.to_matrix() { Ok(matrix) => matrix, Err(err) => return qis_err_code(&err), @@ -171,9 +181,9 @@ pub extern "C" fn hamiltonian_to_matrix( if matrix.len() != len { return CqlibError::InvalidParam as i32; } - let flat: Vec = matrix.iter().cloned().collect(); + // Checked core conversion returns a contiguous row-major matrix. unsafe { - std::ptr::copy_nonoverlapping(flat.as_ptr(), buffer, len); + std::ptr::copy_nonoverlapping(matrix.as_ptr(), buffer, len); } 0 } diff --git a/crates/binding-c/src/qis/hamiltonian_algebra.rs b/crates/binding-c/src/qis/hamiltonian_algebra.rs new file mode 100644 index 00000000..9a1752f4 --- /dev/null +++ b/crates/binding-c/src/qis/hamiltonian_algebra.rs @@ -0,0 +1,302 @@ +// 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. + +//! Checked Pauli sums. Inputs are borrowed and outputs owned; errors do not +//! overwrite outputs. Only exact zeros are removed unless chop is requested. + +use super::algebra::{algebra_err_code, check_pauli_width, input_slice, write_result}; +use super::{CHamiltonian, CPauliString}; +use cqlib_core::qis::Hamiltonian; +use libc::size_t; +use num_complex::Complex64; + +/// Borrowed Pauli term for checked batch construction. Pauli handles are copied. +#[repr(C)] +pub struct CHamiltonianTerm { + pub pauli: *const CPauliString, + pub coefficient: Complex64, +} + +fn owned(inner: Hamiltonian) -> *mut CHamiltonian { + Box::into_raw(Box::new(CHamiltonian { inner })) +} + +/// Construct identity, including scalar one on zero qubits. Returns -1 NULL out, +/// -7 width beyond the BitVec representation limit. Free with hamiltonian_free. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_identity(num_qubits: size_t, out: *mut *mut CHamiltonian) -> i32 { + if out.is_null() { + return -1; + } + if let Err(e) = check_pauli_width(num_qubits) { + return e; + } + unsafe { + *out = owned(Hamiltonian::identity(num_qubits)); + } + 0 +} + +/// Construct a checked, phase-normalized sum, merging duplicate Paulis. +/// Empty input retains num_qubits, including zero. Inputs are copied. +/// Returns -1 NULL term/out/Pauli, -8 width mismatch, nonfinite coefficient or overflow. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_from_terms( + num_qubits: size_t, + terms: *const CHamiltonianTerm, + len: size_t, + out: *mut *mut CHamiltonian, +) -> i32 { + if out.is_null() { + return -1; + } + let terms = match unsafe { input_slice(terms, len) } { + Ok(v) => v, + Err(e) => return e, + }; + let mut copied = Vec::with_capacity(terms.len()); + for term in terms { + if term.pauli.is_null() { + return -1; + } + copied.push((unsafe { &(*term.pauli).inner }.clone(), term.coefficient)); + } + unsafe { write_result(out, Hamiltonian::from_terms(num_qubits, copied).map(owned)) } +} + +/// Return an owned independent structural copy. +/// Free with hamiltonian_free. Returns -1 NULL; copies storage without normalization. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_copy(ptr: *const CHamiltonian, out: *mut *mut CHamiltonian) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { write_result(out, Ok(owned(inner.clone()))) } +} + +/// Return an owned validated copy with phases absorbed and exact zeros removed. +/// Free with hamiltonian_free. Returns -1 NULL, -8 invalid input or merge overflow. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_normalized( + ptr: *const CHamiltonian, + out: *mut *mut CHamiltonian, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { write_result(out, inner.normalized().map(owned)) } +} + +/// Return an owned normalized Hermitian copy; imaginary residuals <=1e-10 are projected to zero. +/// Free with hamiltonian_free. Returns -1 NULL, -8 invalid input or non-Hermitian projection. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_normalized_hermitian( + ptr: *const CHamiltonian, + out: *mut *mut CHamiltonian, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { write_result(out, inner.normalized_hermitian().map(owned)) } +} + +/// Return an owned conjugate transpose after absorbing Pauli phases. +/// Free with hamiltonian_free. Returns -1 NULL, -8 invalid input or merge overflow. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_adjoint( + ptr: *const CHamiltonian, + out: *mut *mut CHamiltonian, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let inner = unsafe { &(*ptr).inner }; + unsafe { write_result(out, inner.adjoint().map(owned)) } +} + +/// Checked add, retaining every nonzero coefficient. +/// Returns an owned result. -1 NULL, -8 width mismatch/nonfinite arithmetic, +/// -7 unrepresentable tensor width. No input is modified. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_add( + left: *const CHamiltonian, + right: *const CHamiltonian, + out: *mut *mut CHamiltonian, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*left).inner.try_add(&(*right).inner).map(owned)) } +} + +/// Checked sub, retaining every nonzero coefficient. +/// Returns an owned result. -1 NULL, -8 width mismatch/nonfinite arithmetic, +/// -7 unrepresentable tensor width. No input is modified. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_sub( + left: *const CHamiltonian, + right: *const CHamiltonian, + out: *mut *mut CHamiltonian, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*left).inner.try_sub(&(*right).inner).map(owned)) } +} + +/// Checked mul, retaining every nonzero coefficient. Right operand acts first. +/// Returns an owned result. -1 NULL, -8 width mismatch/nonfinite arithmetic, +/// -7 unrepresentable tensor width. No input is modified. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_mul( + left: *const CHamiltonian, + right: *const CHamiltonian, + out: *mut *mut CHamiltonian, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*left).inner.try_mul(&(*right).inner).map(owned)) } +} + +/// Checked tensor, retaining every nonzero coefficient. Right operand occupies low qubits, left high qubits. +/// Returns an owned result. -1 NULL, -8 width mismatch/nonfinite arithmetic, +/// -7 unrepresentable tensor width. No input is modified. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_tensor( + left: *const CHamiltonian, + right: *const CHamiltonian, + out: *mut *mut CHamiltonian, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*left).inner.tensor(&(*right).inner).map(owned)) } +} + +/// Return a checked scaled copy. Returns -1 NULL, -8 nonfinite input/factor/result. +/// Unlike legacy hamiltonian_scale, all stored coefficients are validated. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_scaled( + ptr: *const CHamiltonian, + factor: Complex64, + out: *mut *mut CHamiltonian, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*ptr).inner.scaled(factor).map(owned)) } +} + +/// Embed mapping[source]=destination. Destinations must be distinct and in range; +/// unused target qubits are identity. Not a fermionic mode permutation. +/// Returns -1 NULL, -8 invalid mapping/input, -7 unrepresentable target width. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_remap_qubits( + ptr: *const CHamiltonian, + mapping: *const size_t, + len: size_t, + target_width: size_t, + out: *mut *mut CHamiltonian, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + let mapping = match unsafe { input_slice(mapping, len) } { + Ok(v) => v, + Err(e) => return e, + }; + unsafe { + write_result( + out, + (*ptr).inner.remap_qubits(mapping, target_width).map(owned), + ) + } +} + +/// Normalize then remove coefficients with magnitude <=atol in place. +/// Returns -1 NULL, -8 negative/nonfinite atol or invalid input. Failure is atomic. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_chop(ptr: *mut CHamiltonian, atol: f64) -> i32 { + if ptr.is_null() { + return -1; + } + match unsafe { &mut (*ptr).inner }.chop(atol) { + Ok(()) => 0, + Err(e) => algebra_err_code(&e), + } +} + +/// Compare normalized coefficient differences to atol; writes 0 or 1. +/// Returns -1 NULL, -8 invalid atol, width mismatch or nonfinite arithmetic. +/// atol is a coefficient tolerance, not an operator norm bound. atol=0 is exact. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_equiv( + left: *const CHamiltonian, + right: *const CHamiltonian, + atol: f64, + out: *mut u8, +) -> i32 { + if left.is_null() || right.is_null() || out.is_null() { + return -1; + } + unsafe { + write_result( + out, + (*left).inner.equiv(&(*right).inner, atol).map(u8::from), + ) + } +} + +/// Test equivalence with the adjoint; writes 0 or 1. For normalized coefficients +/// this bounds twice the imaginary magnitude. Does not project coefficients. +/// Returns -1 NULL, -8 invalid atol/input or arithmetic overflow. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_is_hermitian( + ptr: *const CHamiltonian, + atol: f64, + out: *mut u8, +) -> i32 { + if ptr.is_null() || out.is_null() { + return -1; + } + unsafe { write_result(out, (*ptr).inner.is_hermitian(atol).map(u8::from)) } +} + +/// Read one stored term by index, returning an independent Pauli handle and +/// coefficient. Free the Pauli with pauli_string_free; it may have zero width. +/// Both outputs are unchanged on error. Returns -1 NULL, -2 invalid index. +#[unsafe(no_mangle)] +pub extern "C" fn hamiltonian_get_term( + ptr: *const CHamiltonian, + index: size_t, + out_pauli: *mut *mut CPauliString, + out_coefficient: *mut Complex64, +) -> i32 { + if ptr.is_null() || out_pauli.is_null() || out_coefficient.is_null() { + return -1; + } + let Some((pauli, coefficient)) = (unsafe { &(*ptr).inner }).terms.get(index) else { + return -2; + }; + let copied = Box::into_raw(Box::new(CPauliString { + inner: pauli.clone(), + })); + unsafe { + *out_pauli = copied; + *out_coefficient = *coefficient; + } + 0 +} diff --git a/crates/binding-c/src/qis/mod.rs b/crates/binding-c/src/qis/mod.rs index 8628c2fc..ebd47e89 100644 --- a/crates/binding-c/src/qis/mod.rs +++ b/crates/binding-c/src/qis/mod.rs @@ -14,10 +14,13 @@ #![allow(clippy::not_unsafe_ptr_arg_deref)] +mod algebra; pub mod density_matrix; pub mod density_matrix_noise; pub mod entropy; +pub mod fermion; pub mod hamiltonian; +pub mod hamiltonian_algebra; pub mod measurement; pub mod metrics; pub mod pauli; @@ -103,7 +106,9 @@ pub(crate) fn qis_err_code(err: &cqlib_core::qis::QisError) -> i32 { pub use density_matrix::*; pub use density_matrix_noise::*; pub use entropy::*; +pub use fermion::*; pub use hamiltonian::*; +pub use hamiltonian_algebra::*; pub use measurement::*; pub use metrics::*; pub use pauli::*; diff --git a/crates/binding-c/src/qis/pauli.rs b/crates/binding-c/src/qis/pauli.rs index 7a88e9a7..305dde8e 100644 --- a/crates/binding-c/src/qis/pauli.rs +++ b/crates/binding-c/src/qis/pauli.rs @@ -13,8 +13,9 @@ //! C ABI for `PauliString` and single-qubit Pauli operators. use crate::error::CqlibError; +use crate::qis::algebra::{dense_matrix_dim, dense_matrix_len}; use crate::qis::{CPauliString, qis_err_code}; -use cqlib_core::qis::{Pauli, PauliString}; +use cqlib_core::qis::{Hamiltonian, Pauli, PauliString}; use std::collections::HashMap; use std::ffi::{CStr, CString}; use std::os::raw::c_char; @@ -184,22 +185,22 @@ pub extern "C" fn pauli_string_to_string(ptr: *const CPauliString) -> *mut c_cha } } -/// Flatten the dense matrix representation into `buffer`. -/// -/// The matrix has dimensions `2^N × 2^N` and is stored in row-major order -/// as `Complex64` values. Use `pauli_string_matrix_len` first to obtain -/// the side length (`2^N`). Returns 0 on success. +/// Return the side length (`2^N`) of the row-major dense matrix, or 0 for +/// NULL or an unrepresentable side length. A representable side does not +/// guarantee that the complete dense matrix can be allocated. #[unsafe(no_mangle)] pub extern "C" fn pauli_string_matrix_len(ptr: *const CPauliString) -> usize { if ptr.is_null() { return 0; } let wrapper = unsafe { &*ptr }; - 1usize << wrapper.inner.num_qubits + dense_matrix_dim(wrapper.inner.num_qubits).unwrap_or(0) } /// Copy the dense matrix into `buffer` (row-major, `Complex64` values). /// `len` must equal `pauli_string_matrix_len(ptr) * pauli_string_matrix_len(ptr)`. +/// Returns -7 for unrepresentable dense allocation size or allocation failure, +/// -8 for an incorrect buffer length. Errors leave the buffer unchanged. #[unsafe(no_mangle)] pub extern "C" fn pauli_string_matrix( ptr: *const CPauliString, @@ -210,14 +211,19 @@ pub extern "C" fn pauli_string_matrix( return CqlibError::NullPtr as i32; } let wrapper = unsafe { &*ptr }; - let matrix = wrapper.inner.to_matrix(); - let dim = 1usize << wrapper.inner.num_qubits; - if matrix.len() != len || len != dim * dim { + let Some(expected_len) = dense_matrix_len(wrapper.inner.num_qubits) else { + return CqlibError::SimulationError as i32; + }; + if len != expected_len { return CqlibError::InvalidParam as i32; } - let flat: Vec = matrix.iter().cloned().collect(); + let matrix = match Hamiltonian::from_pauli(wrapper.inner.clone()).to_matrix() { + Ok(matrix) => matrix, + Err(err) => return qis_err_code(&err), + }; + // Checked core conversion returns a contiguous row-major matrix. unsafe { - std::ptr::copy_nonoverlapping(flat.as_ptr(), buffer, len); + std::ptr::copy_nonoverlapping(matrix.as_ptr(), buffer, len); } 0 } diff --git a/crates/binding-c/tests/qis_algebra_test.rs b/crates/binding-c/tests/qis_algebra_test.rs new file mode 100644 index 00000000..3cd0b35f --- /dev/null +++ b/crates/binding-c/tests/qis_algebra_test.rs @@ -0,0 +1,630 @@ +// This code is part of Cqlib. +// +// (C) Copyright China Telecom Quantum Group 2026 +// +// This code is licensed under the Apache License, Version 2.0. You may +// obtain a copy of this license in the LICENSE.txt file in the root directory +// of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +// +// Any modifications or derivative works of this code must retain this +// copyright notice, and modified files need to carry a notice indicating +// that they have been altered from the originals. + +use binding_c::{circuit::*, qis::*}; +use num_complex::Complex64 as C; +use std::ptr::{null, null_mut}; + +struct F(*mut CFermionicOp); +impl F { + fn make(call: impl FnOnce(*mut *mut CFermionicOp) -> i32) -> Self { + let mut out = null_mut(); + assert_eq!(call(&mut out), 0); + assert!(!out.is_null()); + Self(out) + } + fn mapped(&self) -> H { + H::make(|out| jordan_wigner_map_operator(self.0, out)) + } + fn equiv(&self, other: &Self) -> bool { + let mut result = 99; + assert_eq!(fermionic_op_equiv(self.0, other.0, 0.0, &mut result), 0); + result == 1 + } +} +impl Drop for F { + fn drop(&mut self) { + fermionic_op_free(self.0); + } +} +struct H(*mut CHamiltonian); +impl H { + fn make(call: impl FnOnce(*mut *mut CHamiltonian) -> i32) -> Self { + let mut out = null_mut(); + assert_eq!(call(&mut out), 0); + assert!(!out.is_null()); + Self(out) + } + fn from_pauli(text: &std::ffi::CStr, coefficient: C) -> Self { + let p = pauli_string_parse(text.as_ptr()); + assert!(!p.is_null()); + let terms = [CHamiltonianTerm { + pauli: p, + coefficient, + }]; + let h = Self::make(|out| { + hamiltonian_from_terms(pauli_string_num_qubits(p), terms.as_ptr(), 1, out) + }); + pauli_string_free(p); + h + } + fn equiv(&self, other: &Self) -> bool { + let mut result = 99; + assert_eq!(hamiltonian_equiv(self.0, other.0, 0.0, &mut result), 0); + result == 1 + } + fn matrix(&self) -> Vec { + let dim = hamiltonian_to_matrix_len(self.0); + let mut result = vec![C::default(); dim * dim]; + assert_eq!( + hamiltonian_to_matrix(self.0, result.as_mut_ptr(), result.len()), + 0 + ); + result + } +} +impl Drop for H { + fn drop(&mut self) { + hamiltonian_free(self.0); + } +} + +#[test] +fn empty_spaces_and_owned_snapshots() { + let zero = F::make(|out| fermionic_op_from_terms(0, null(), 0, out)); + let identity = F::make(|out| fermionic_op_from_term(0, null(), 0, C::new(1.0, 0.0), out)); + let total = F::make(|out| fermionic_op_total_number(0, out)); + assert!(zero.equiv(&total)); + assert!(!zero.equiv(&identity)); + let mut n = 99; + assert_eq!(fermionic_op_num_modes(zero.0, &mut n), 0); + assert_eq!(n, 0); + assert_eq!(fermionic_op_num_terms(identity.0, &mut n), 0); + assert_eq!(n, 1); + assert_eq!(fermionic_op_term_len(identity.0, 0, &mut n), 0); + assert_eq!(n, 0); + let mut coefficient = C::default(); + assert_eq!( + fermionic_op_get_term(identity.0, 0, null_mut(), 0, &mut coefficient), + 0 + ); + assert_eq!(coefficient, C::new(1.0, 0.0)); + let h = identity.mapped(); + let mut pauli = null_mut(); + assert_eq!( + hamiltonian_get_term(h.0, 0, &mut pauli, &mut coefficient), + 0 + ); + drop(h); + assert_eq!(pauli_string_num_qubits(pauli), 0); + pauli_string_free(pauli); + let scalar = H::make(|out| hamiltonian_identity(0, out)); + assert_eq!(scalar.matrix(), vec![C::new(1.0, 0.0)]); + let huge_zero = F::make(|out| fermionic_op_new(usize::MAX, out)); + let mapped = huge_zero.mapped(); + assert_eq!(hamiltonian_num_qubits(mapped.0), usize::MAX); + let huge_identity = F::make(|out| fermionic_op_identity(usize::MAX, out)); + let mut out = scalar.0; + assert_eq!(jordan_wigner_map_operator(huge_identity.0, &mut out), -7); + assert_eq!(out, scalar.0); + assert_eq!(hamiltonian_identity(usize::MAX, &mut out), -7); + assert_eq!(out, scalar.0); +} + +#[test] +fn anticommutation_normal_order_and_structural_equality() { + let a = F::make(|out| fermionic_op_annihilation(2, 0, out)); + let creator = F::make(|out| fermionic_op_creation(2, 0, out)); + let ac = F::make(|out| fermionic_op_mul(a.0, creator.0, out)); + let ca = F::make(|out| fermionic_op_mul(creator.0, a.0, out)); + let sum = F::make(|out| fermionic_op_add(ac.0, ca.0, out)); + let identity = F::make(|out| fermionic_op_identity(2, out)); + assert!(sum.equiv(&identity)); + let mut same = 99; + assert_eq!(fermionic_op_equals(sum.0, identity.0, &mut same), 0); + assert_eq!(same, 0); + let ordered = F::make(|out| fermionic_op_normal_order(sum.0, out)); + assert_eq!(fermionic_op_equals(ordered.0, identity.0, &mut same), 0); + assert_eq!(same, 1); + let square = F::make(|out| fermionic_op_mul(a.0, a.0, out)); + assert_eq!(fermionic_op_is_zero(square.0, &mut same), 0); + assert_eq!(same, 0); + let ordered = F::make(|out| fermionic_op_normal_order(square.0, out)); + assert_eq!(fermionic_op_is_zero(ordered.0, &mut same), 0); + assert_eq!(same, 1); + let other = F::make(|out| fermionic_op_creation(2, 1, out)); + let ao = F::make(|out| fermionic_op_mul(a.0, other.0, out)); + let oa = F::make(|out| fermionic_op_mul(other.0, a.0, out)); + let cross = F::make(|out| fermionic_op_add(ao.0, oa.0, out)); + let zero = F::make(|out| fermionic_op_new(2, out)); + assert!(cross.equiv(&zero)); + let adjoint = F::make(|out| fermionic_op_adjoint(a.0, out)); + assert_eq!(fermionic_op_equals(adjoint.0, creator.0, &mut same), 0); + assert_eq!(same, 1); + let copy = F::make(|out| fermionic_op_copy(creator.0, out)); + assert_eq!(fermionic_op_chop(copy.0, 1.0), 0); + assert!(!copy.equiv(&creator)); +} + +#[test] +fn excitation_order_signs_and_number_conservation() { + let source = [0, 1]; + let target = [2, 3]; + let generator = F::make(|out| { + fermionic_op_excitation_generator(4, source.as_ptr(), 2, target.as_ptr(), 2, out) + }); + let reversed = [3, 2]; + let swapped = F::make(|out| { + fermionic_op_excitation_generator(4, source.as_ptr(), 2, reversed.as_ptr(), 2, out) + }); + let negative = F::make(|out| fermionic_op_scaled(generator.0, C::new(-1.0, 0.0), out)); + assert!(swapped.equiv(&negative)); + let mut yes = 99; + assert_eq!(fermionic_op_is_hermitian(generator.0, 0.0, &mut yes), 0); + assert_eq!(yes, 1); + let number = F::make(|out| fermionic_op_total_number(4, out)); + let ng = F::make(|out| fermionic_op_mul(number.0, generator.0, out)); + let gn = F::make(|out| fermionic_op_mul(generator.0, number.0, out)); + let commutator = F::make(|out| fermionic_op_sub(ng.0, gn.0, out)); + let zero = F::make(|out| fermionic_op_new(4, out)); + assert!(commutator.equiv(&zero)); + let mut out = generator.0; + for (sources, targets, expected) in [ + (&[][..], &[][..], -8), + (&source[..], &target[..1], -8), + (&[0, 0][..], &target[..], -8), + (&source[..], &[1, 3][..], -8), + (&source[..], &[2, 4][..], -2), + ] { + assert_eq!( + fermionic_op_excitation_generator( + 4, + sources.as_ptr(), + sources.len(), + targets.as_ptr(), + targets.len(), + &mut out + ), + expected + ); + assert_eq!(out, generator.0); + } +} + +#[test] +fn jordan_wigner_signs_occupations_and_nonhermitian_results() { + let creator = F::make(|out| fermionic_op_creation(1, 0, out)); + let annihilator = F::make(|out| fermionic_op_annihilation(1, 0, out)); + let creation_matrix = creator.mapped().matrix(); + assert_eq!( + creation_matrix, + vec![C::default(), C::default(), C::new(1.0, 0.0), C::default()] + ); + assert_eq!( + annihilator.mapped().matrix(), + vec![C::default(), C::new(1.0, 0.0), C::default(), C::default()] + ); + let number = F::make(|out| fermionic_op_number(1, 0, out)); + assert_eq!( + number.mapped().matrix(), + vec![C::default(), C::default(), C::default(), C::new(1.0, 0.0)] + ); + let high = F::make(|out| fermionic_op_creation(2, 1, out)); + let matrix = high.mapped().matrix(); + assert_eq!(matrix[2 * 4], C::new(1.0, 0.0)); + assert_eq!(matrix[3 * 4 + 1], C::new(-1.0, 0.0)); + let mut bits = [99; 4]; + let occupation = [1, 0, 1]; + assert_eq!( + jordan_wigner_map_occupation(occupation.as_ptr(), 3, bits.as_mut_ptr(), 4), + 0 + ); + assert_eq!(bits, [1, 0, 1, 99]); + let before = bits; + assert_eq!( + jordan_wigner_map_occupation([1, 2, 0].as_ptr(), 3, bits.as_mut_ptr(), 4), + -8 + ); + assert_eq!(bits, before); + assert_eq!( + jordan_wigner_map_occupation(occupation.as_ptr(), 3, bits.as_mut_ptr(), 2), + -8 + ); + assert_eq!(bits, before); + assert_eq!(jordan_wigner_map_occupation(null(), 0, null_mut(), 0), 0); + assert_eq!( + jordan_wigner_map_occupation(bits.as_ptr(), 3, bits.as_mut_ptr(), 4), + 0 + ); + let mut n = 99; + assert_eq!(jordan_wigner_num_qubits(3, &mut n), 0); + assert_eq!(n, 3); +} + +#[test] +fn checked_pauli_phases_tensor_and_remap() { + let x = H::from_pauli(c"X", C::new(1.0, 0.0)); + let z = H::from_pauli(c"Z", C::new(1.0, 0.0)); + let minus_i_y = H::from_pauli(c"Y", C::new(0.0, -1.0)); + let xz = H::make(|out| hamiltonian_mul(x.0, z.0, out)); + assert!(xz.equiv(&minus_i_y)); + let dagger = H::make(|out| hamiltonian_adjoint(xz.0, out)); + let negative = H::make(|out| hamiltonian_scaled(xz.0, C::new(-1.0, 0.0), out)); + assert!(dagger.equiv(&negative)); + let sum = H::make(|out| hamiltonian_add(x.0, z.0, out)); + let difference = H::make(|out| hamiltonian_sub(sum.0, z.0, out)); + assert!(difference.equiv(&x)); + let tensor = H::make(|out| hamiltonian_tensor(x.0, z.0, out)); + let expected = H::from_pauli(c"XZ", C::new(1.0, 0.0)); + assert!(tensor.equiv(&expected)); + let embed = H::make(|out| hamiltonian_remap_qubits(x.0, [2].as_ptr(), 1, 3, out)); + let expected = H::from_pauli(c"XII", C::new(1.0, 0.0)); + assert!(embed.equiv(&expected)); + let mut out = embed.0; + assert_eq!( + hamiltonian_remap_qubits(tensor.0, [0, 0].as_ptr(), 2, 3, &mut out), + -8 + ); + assert_eq!( + hamiltonian_remap_qubits(x.0, [3].as_ptr(), 1, 3, &mut out), + -8 + ); + assert_eq!(hamiltonian_remap_qubits(x.0, null(), 0, 3, &mut out), -8); + assert_eq!(out, embed.0); + let zero = H::make(|out| hamiltonian_from_terms(0, null(), 0, out)); + let scalar_embed = H::make(|out| hamiltonian_remap_qubits(zero.0, null(), 0, 3, out)); + assert_eq!(hamiltonian_num_qubits(scalar_embed.0), 3); + let phased = H::from_pauli(c"iX", C::new(0.0, -1.0)); + assert!(phased.equiv(&x)); +} + +#[test] +fn nonzero_coefficients_projection_and_legacy_semantics() { + let tiny = 1e-20; + let creator = F::make(|out| fermionic_op_creation(1, 0, out)); + let scaled = F::make(|out| fermionic_op_scaled(creator.0, C::new(tiny, 0.0), out)); + let mut n = 99; + assert_eq!(fermionic_op_num_terms(scaled.0, &mut n), 0); + assert_eq!(n, 1); + let mapped = scaled.mapped(); + assert_eq!(hamiltonian_num_terms(mapped.0), 2); + let h = H::from_pauli(c"X", C::new(tiny, 0.0)); + let normalized = H::make(|out| hamiltonian_normalized(h.0, out)); + assert_eq!(hamiltonian_num_terms(normalized.0), 1); + let copy = H::make(|out| hamiltonian_copy(h.0, out)); + assert_eq!(hamiltonian_simplify(copy.0), 0); + assert_eq!(hamiltonian_num_terms(copy.0), 0); + assert_eq!(hamiltonian_num_terms(h.0), 1); + assert_eq!(hamiltonian_chop(normalized.0, tiny), 0); + assert_eq!(hamiltonian_num_terms(normalized.0), 0); + let residual = H::from_pauli(c"Z", C::new(2.0, 4e-11)); + let mut yes = 99; + assert_eq!(hamiltonian_is_hermitian(residual.0, 0.0, &mut yes), 0); + assert_eq!(yes, 0); + assert_eq!(hamiltonian_is_hermitian(residual.0, 8e-11, &mut yes), 0); + assert_eq!(yes, 1); + let projected = H::make(|out| hamiltonian_normalized_hermitian(residual.0, out)); + let real = H::from_pauli(c"Z", C::new(2.0, 0.0)); + assert!(projected.equiv(&real)); + let imaginary = H::from_pauli(c"Z", C::new(2.0, 2e-10)); + let mut out = real.0; + assert_eq!(hamiltonian_normalized_hermitian(imaginary.0, &mut out), -8); + assert_eq!(out, real.0); +} + +#[test] +fn invalid_coefficients_atomicity_and_array_contracts() { + let creator = F::make(|out| fermionic_op_creation(1, 0, out)); + let snapshot = F::make(|out| fermionic_op_copy(creator.0, out)); + let ladder = CFermionicLadder { + mode: 0, + kind: FERMIONIC_CREATE, + }; + let mut out = creator.0; + for bad in [f64::NAN, f64::INFINITY, f64::NEG_INFINITY] { + assert_eq!( + fermionic_op_add_term(creator.0, &ladder, 1, C::new(bad, 0.0)), + -8 + ); + assert!(creator.equiv(&snapshot)); + assert_eq!( + fermionic_op_scaled(creator.0, C::new(0.0, bad), &mut out), + -8 + ); + assert_eq!(out, creator.0); + } + for bad in [-1.0, f64::NAN, f64::INFINITY] { + assert_eq!(fermionic_op_chop(creator.0, bad), -8); + assert!(creator.equiv(&snapshot)); + let mut result = 99; + assert_eq!( + fermionic_op_equiv(creator.0, snapshot.0, bad, &mut result), + -8 + ); + assert_eq!(result, 99); + } + let invalid = CFermionicLadder { mode: 0, kind: 255 }; + assert_eq!( + fermionic_op_from_term(1, &invalid, 1, C::default(), &mut out), + -8 + ); + let invalid = CFermionicLadder { + mode: 1, + kind: FERMIONIC_CREATE, + }; + assert_eq!( + fermionic_op_from_term(1, &invalid, 1, C::default(), &mut out), + -2 + ); + assert_eq!( + fermionic_op_from_term(1, &ladder, usize::MAX, C::default(), &mut out), + -8 + ); + assert_eq!(fermionic_op_from_terms(1, null(), 1, &mut out), -1); + let mut buffer = [CFermionicLadder { mode: 99, kind: 99 }; 2]; + let mut coefficient = C::new(99.0, 99.0); + assert_eq!( + fermionic_op_get_term(creator.0, 0, buffer.as_mut_ptr(), 0, &mut coefficient), + -8 + ); + assert_eq!(buffer[0].mode, 99); + assert_eq!(coefficient, C::new(99.0, 99.0)); + assert_eq!( + fermionic_op_get_term(creator.0, 9, buffer.as_mut_ptr(), 2, &mut coefficient), + -2 + ); + assert_eq!( + fermionic_op_get_term(creator.0, 0, null_mut(), 2, &mut coefficient), + -1 + ); + assert_eq!( + fermionic_op_get_term(creator.0, 0, buffer.as_mut_ptr(), 2, &mut coefficient), + 0 + ); + assert_eq!((buffer[0].mode, buffer[0].kind), (0, FERMIONIC_CREATE)); + assert_eq!(buffer[1].mode, 99); + buffer[0].kind = FERMIONIC_ANNIHILATE; + assert_eq!(buffer[0].kind, FERMIONIC_ANNIHILATE); + assert!(creator.equiv(&snapshot)); + let huge = F::make(|out| fermionic_op_from_term(1, &ladder, 1, C::new(f64::MAX, 0.0), out)); + let mut out = creator.0; + assert_eq!(fermionic_op_add(huge.0, huge.0, &mut out), -8); + assert_eq!( + fermionic_op_add_term(huge.0, &ladder, 1, C::new(f64::MAX, 0.0)), + -8 + ); + assert_eq!(fermionic_op_scaled(huge.0, C::new(2.0, 0.0), &mut out), -8); + assert_eq!(out, creator.0); + assert_eq!(fermionic_op_creation(1, 1, &mut out), -2); + let wrong_width = F::make(|out| fermionic_op_new(2, out)); + assert_eq!(fermionic_op_mul(creator.0, wrong_width.0, &mut out), -8); + assert_eq!(out, creator.0); +} + +#[test] +fn checked_hamiltonian_validates_legacy_storage_without_overwriting_outputs() { + let p = pauli_string_parse(c"X".as_ptr()); + let h = H(hamiltonian_new(1)); + assert_eq!(hamiltonian_add_term(h.0, p, f64::NAN, 0.0), 0); // legacy behavior + pauli_string_free(p); + let snapshot = unsafe { (*h.0).inner.clone() }; + let valid = H::from_pauli(c"Z", C::new(1.0, 0.0)); + let mut out = valid.0; + assert_eq!(hamiltonian_normalized(h.0, &mut out), -8); + assert_eq!(hamiltonian_scaled(h.0, C::new(1.0, 0.0), &mut out), -8); + assert_eq!(hamiltonian_add(h.0, valid.0, &mut out), -8); + assert_eq!(hamiltonian_chop(h.0, 0.0), -8); + assert_eq!(out, valid.0); + assert!(unsafe { (&(*h.0).inner.terms)[0].1.re.is_nan() }); + assert_eq!(unsafe { (*h.0).inner.terms.len() }, snapshot.terms.len()); + let large = H::from_pauli(c"X", C::new(f64::MAX, 0.0)); + assert_eq!(hamiltonian_add(large.0, large.0, &mut out), -8); + assert_eq!(hamiltonian_scaled(large.0, C::new(2.0, 0.0), &mut out), -8); + let p = pauli_string_parse(c"X".as_ptr()); + let terms = [ + CHamiltonianTerm { + pauli: p, + coefficient: C::new(f64::MAX, 0.0), + }, + CHamiltonianTerm { + pauli: p, + coefficient: C::new(f64::MAX, 0.0), + }, + ]; + assert_eq!(hamiltonian_from_terms(1, terms.as_ptr(), 2, &mut out), -8); + assert_eq!(hamiltonian_from_terms(2, terms.as_ptr(), 1, &mut out), -8); + assert_eq!( + hamiltonian_from_terms(1, terms.as_ptr(), usize::MAX, &mut out), + -8 + ); + pauli_string_free(p); + for bad in [-1.0, f64::NAN, f64::INFINITY] { + let mut yes = 99; + assert_eq!(hamiltonian_chop(valid.0, bad), -8); + assert_eq!(hamiltonian_is_hermitian(valid.0, bad, &mut yes), -8); + assert_eq!(yes, 99); + } + let mut copied = null_mut(); + let mut coefficient = C::new(99.0, 0.0); + assert_eq!( + hamiltonian_get_term(valid.0, 9, &mut copied, &mut coefficient), + -2 + ); + assert!(copied.is_null()); + assert_eq!(coefficient.re, 99.0); + assert_eq!( + hamiltonian_get_term(valid.0, 0, &mut copied, &mut coefficient), + 0 + ); + assert_eq!(pauli_string_set_pauli(copied, 0, PAULI_X), 0); + pauli_string_free(copied); + assert_eq!( + valid.matrix(), + vec![ + C::new(1.0, 0.0), + C::default(), + C::default(), + C::new(-1.0, 0.0) + ] + ); +} + +#[test] +fn excitation_to_evolution_and_numeric_phase_regressions() { + let generator = + F::make(|out| fermionic_op_excitation_generator(2, [0].as_ptr(), 1, [1].as_ptr(), 1, out)); + let mapped = generator.mapped(); + let theta: f64 = 0.37; + let evolution = hamiltonian_to_evolution_circuit(mapped.0, theta, 1, TROTTER_MODE_SECOND_ORDER); + assert!(!evolution.is_null()); + let state = [C::default(), C::new(1.0, 0.0), C::default(), C::default()]; + let sv = statevector_from_state(2, state.as_ptr(), 4); + assert_eq!(statevector_apply_circuit(sv, evolution), 0); + let mut data = [C::default(); 4]; + assert_eq!(statevector_data(sv, data.as_mut_ptr(), 4), 0); + assert!((data[1] - C::new(theta.cos(), 0.0)).norm() < 1e-12); + assert!((data[2] - C::new(theta.sin(), 0.0)).norm() < 1e-12); + assert!(data[0].norm() < 1e-12 && data[3].norm() < 1e-12); + statevector_free(sv); + circuit_free(evolution); + let phase = circuit_new(1); + assert_eq!(circuit_set_global_phase(phase, 0.37), 0); + let sv = statevector_from_circuit(phase); + let mut data = [C::default(); 2]; + assert_eq!(statevector_data(sv, data.as_mut_ptr(), 2), 0); + assert!((data[0] - C::from_polar(1.0, 0.37)).norm() < 1e-12); + statevector_free(sv); + circuit_free(phase); + let a = param_from_double(1e-20); + let b = param_from_double(2.0); + let product = param_mul(a, b); + assert_eq!(param_evaluate(product, null()), 2e-20); + param_free(product); + param_free(a); + param_free(b); +} + +#[test] +fn all_checked_entry_points_reject_required_null_pointers() { + assert_eq!(fermionic_op_new(1, null_mut()), -1); + assert_eq!(fermionic_op_identity(1, null_mut()), -1); + assert_eq!(fermionic_op_total_number(1, null_mut()), -1); + assert_eq!(fermionic_op_creation(1, 1, null_mut()), -1); + assert_eq!(fermionic_op_annihilation(1, 1, null_mut()), -1); + assert_eq!(fermionic_op_number(1, 1, null_mut()), -1); + assert_eq!( + fermionic_op_from_term(1, null(), 1, C::default(), null_mut()), + -1 + ); + assert_eq!(fermionic_op_from_terms(1, null(), 1, null_mut()), -1); + assert_eq!( + fermionic_op_add_term(null_mut(), null(), 1, C::default()), + -1 + ); + assert_eq!( + fermionic_op_excitation_generator(1, null(), 1, null(), 1, null_mut()), + -1 + ); + assert_eq!(fermionic_op_num_modes(null(), null_mut()), -1); + assert_eq!(fermionic_op_num_terms(null(), null_mut()), -1); + assert_eq!(fermionic_op_is_zero(null(), null_mut()), -1); + assert_eq!(fermionic_op_copy(null(), null_mut()), -1); + assert_eq!(fermionic_op_adjoint(null(), null_mut()), -1); + assert_eq!(fermionic_op_normal_order(null(), null_mut()), -1); + assert_eq!(fermionic_op_add(null(), null(), null_mut()), -1); + assert_eq!(fermionic_op_sub(null(), null(), null_mut()), -1); + assert_eq!(fermionic_op_mul(null(), null(), null_mut()), -1); + assert_eq!(fermionic_op_scaled(null(), C::default(), null_mut()), -1); + assert_eq!(fermionic_op_chop(null_mut(), 0.0), -1); + assert_eq!(fermionic_op_equals(null(), null(), null_mut()), -1); + assert_eq!(fermionic_op_equiv(null(), null(), 0.0, null_mut()), -1); + assert_eq!(fermionic_op_is_hermitian(null(), 0.0, null_mut()), -1); + assert_eq!(fermionic_op_term_len(null(), 1, null_mut()), -1); + assert_eq!( + fermionic_op_get_term(null(), 1, null_mut(), 1, null_mut()), + -1 + ); + assert_eq!(jordan_wigner_map_operator(null(), null_mut()), -1); + assert_eq!(jordan_wigner_num_qubits(1, null_mut()), -1); + assert_eq!(jordan_wigner_map_occupation(null(), 1, null_mut(), 1), -1); + assert_eq!(hamiltonian_identity(1, null_mut()), -1); + assert_eq!(hamiltonian_from_terms(1, null(), 1, null_mut()), -1); + assert_eq!(hamiltonian_copy(null(), null_mut()), -1); + assert_eq!(hamiltonian_normalized(null(), null_mut()), -1); + assert_eq!(hamiltonian_normalized_hermitian(null(), null_mut()), -1); + assert_eq!(hamiltonian_adjoint(null(), null_mut()), -1); + assert_eq!(hamiltonian_add(null(), null(), null_mut()), -1); + assert_eq!(hamiltonian_sub(null(), null(), null_mut()), -1); + assert_eq!(hamiltonian_mul(null(), null(), null_mut()), -1); + assert_eq!(hamiltonian_tensor(null(), null(), null_mut()), -1); + assert_eq!(hamiltonian_scaled(null(), C::default(), null_mut()), -1); + assert_eq!( + hamiltonian_remap_qubits(null(), null(), 1, 1, null_mut()), + -1 + ); + assert_eq!(hamiltonian_chop(null_mut(), 0.0), -1); + assert_eq!(hamiltonian_equiv(null(), null(), 0.0, null_mut()), -1); + assert_eq!(hamiltonian_is_hermitian(null(), 0.0, null_mut()), -1); + assert_eq!(hamiltonian_get_term(null(), 1, null_mut(), null_mut()), -1); +} + +#[test] +fn wide_sparse_registers_reject_dense_sizes_without_panicking_or_writing() { + let bits = usize::BITS as usize; + // Exercise side-length, element-count and byte-count overflow separately. + for width in [ + bits, + bits + 1, + bits - 1, + bits / 2, + bits / 2 - 1, + bits / 2 - 2, + ] { + let identity = F::make(|out| fermionic_op_identity(width, out)); + let mapped = identity.mapped(); + let mut pauli = null_mut(); + let mut coefficient = C::default(); + assert_eq!( + hamiltonian_get_term(mapped.0, 0, &mut pauli, &mut coefficient), + 0 + ); + let expected_dim = if width < bits { 1usize << width } else { 0 }; + assert_eq!(hamiltonian_to_matrix_len(mapped.0), expected_dim); + assert_eq!(pauli_string_matrix_len(pauli), expected_dim); + let sentinel = C::new(99.0, -99.0); + let mut buffer = sentinel; + assert_eq!(hamiltonian_to_matrix(mapped.0, &mut buffer, 1), -7); + assert_eq!(buffer, sentinel); + assert_eq!(pauli_string_matrix(pauli, &mut buffer, 1), -7); + assert_eq!(buffer, sentinel); + pauli_string_free(pauli); + } + let zero = F::make(|out| fermionic_op_new(usize::MAX, out)).mapped(); + assert_eq!(hamiltonian_to_matrix_len(zero.0), 0); + let mut buffer = C::new(99.0, -99.0); + assert_eq!(hamiltonian_to_matrix(zero.0, &mut buffer, 1), -7); + assert_eq!(buffer, C::new(99.0, -99.0)); + // A valid width with a tiny buffer must be rejected before dense allocation. + let identity = H::make(|out| hamiltonian_identity(20, out)); + let mut pauli = null_mut(); + let mut coefficient = C::default(); + assert_eq!( + hamiltonian_get_term(identity.0, 0, &mut pauli, &mut coefficient), + 0 + ); + assert_eq!(hamiltonian_to_matrix(identity.0, &mut buffer, 1), -8); + assert_eq!(pauli_string_matrix(pauli, &mut buffer, 1), -8); + assert_eq!(buffer, C::new(99.0, -99.0)); + pauli_string_free(pauli); + assert_eq!(hamiltonian_to_matrix_len(null()), 0); + assert_eq!(pauli_string_matrix_len(null()), 0); +} diff --git a/crates/binding-c/tests/test_qis_algebra.c b/crates/binding-c/tests/test_qis_algebra.c new file mode 100644 index 00000000..474f48c6 --- /dev/null +++ b/crates/binding-c/tests/test_qis_algebra.c @@ -0,0 +1,249 @@ +// 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. + +#include +#include +#include +#include +#include + +#include "cqlib_c.h" + +_Static_assert(sizeof(Complex64) == 2 * sizeof(double), "complex layout"); +_Static_assert(offsetof(Complex64, im) == sizeof(double), "complex imaginary offset"); +_Static_assert(offsetof(CFermionicLadder, kind) == sizeof(size_t), "ladder tag offset"); +_Static_assert(offsetof(CFermionicTerm, len) == sizeof(void*), "term length offset"); +_Static_assert(offsetof(CHamiltonianTerm, coefficient) == sizeof(void*), + "Pauli coefficient offset"); + +static void test_fermions(void) { + CFermionicOp *a = NULL, *c = NULL, *ac = NULL, *ca = NULL; + CFermionicOp *sum = NULL, *identity = NULL, *ordered = NULL, *copy = NULL; + CFermionicOp *number = NULL, *total = NULL, *zero = NULL, *scaled = NULL; + CFermionicOp *difference = NULL, *adjoint = NULL; + uint8_t yes = 99; + size_t n = 99; + assert(fermionic_op_annihilation(1, 0, &a) == 0); + assert(fermionic_op_creation(1, 0, &c) == 0); + assert(fermionic_op_mul(a, c, &ac) == 0); + assert(fermionic_op_mul(c, a, &ca) == 0); + assert(fermionic_op_add(ac, ca, &sum) == 0); + assert(fermionic_op_identity(1, &identity) == 0); + assert(fermionic_op_equiv(sum, identity, 0, &yes) == 0 && yes == 1); + assert(fermionic_op_equals(sum, identity, &yes) == 0 && yes == 0); + assert(fermionic_op_normal_order(sum, &ordered) == 0); + assert(fermionic_op_equals(ordered, identity, &yes) == 0 && yes == 1); + assert(fermionic_op_num_modes(c, &n) == 0 && n == 1); + assert(fermionic_op_num_terms(c, &n) == 0 && n == 1); + assert(fermionic_op_term_len(c, 0, &n) == 0 && n == 1); + CFermionicLadder buffer[2] = {{99, 99}, {99, 99}}; + Complex64 coefficient = {99, 99}; + assert(fermionic_op_get_term(c, 0, buffer, 0, &coefficient) == -8); + assert(buffer[0].mode == 99 && coefficient.re == 99); + assert(fermionic_op_get_term(c, 0, buffer, 2, &coefficient) == 0); + assert(buffer[0].mode == 0 && buffer[0].kind == FERMIONIC_CREATE); + assert(buffer[1].mode == 99); + assert(fermionic_op_copy(c, ©) == 0); + assert(fermionic_op_chop(copy, 1) == 0); + assert(fermionic_op_is_zero(copy, &yes) == 0 && yes == 1); + assert(fermionic_op_number(1, 0, &number) == 0); + assert(fermionic_op_total_number(1, &total) == 0); + assert(fermionic_op_equiv(number, total, 0, &yes) == 0 && yes == 1); + assert(fermionic_op_new(1, &zero) == 0); + assert(fermionic_op_scaled(c, (Complex64){1e-20, 0}, &scaled) == 0); + assert(fermionic_op_num_terms(scaled, &n) == 0 && n == 1); + assert(fermionic_op_sub(c, c, &difference) == 0); + assert(fermionic_op_equiv(difference, zero, 0, &yes) == 0 && yes == 1); + assert(fermionic_op_adjoint(a, &adjoint) == 0); + assert(fermionic_op_equals(adjoint, c, &yes) == 0 && yes == 1); + assert(fermionic_op_is_hermitian(number, 0, &yes) == 0 && yes == 1); + CFermionicLadder creator = {0, FERMIONIC_CREATE}; + CFermionicTerm terms[] = {{&creator, 1, {1, 0}}, {&creator, 1, {2, 0}}}; + CFermionicOp *batch = NULL, *single = NULL; + assert(fermionic_op_from_terms(1, terms, 2, &batch) == 0); + assert(fermionic_op_from_term(1, &creator, 1, (Complex64){3, 0}, &single) == 0); + assert(fermionic_op_equiv(batch, single, 0, &yes) == 0 && yes == 1); + assert(fermionic_op_add_term(batch, &creator, 1, (Complex64){-3, 0}) == 0); + assert(fermionic_op_is_zero(batch, &yes) == 0 && yes == 1); + CFermionicOp* out = c; + assert(fermionic_op_creation(1, 1, &out) == -2 && out == c); + creator.kind = 255; + assert(fermionic_op_from_term(1, &creator, 1, (Complex64){0, 0}, &out) == -8 && out == c); + assert(fermionic_op_scaled(c, (Complex64){NAN, 0}, &out) == -8 && out == c); + assert(fermionic_op_chop(c, -1) == -8); + assert(fermionic_op_equals(c, adjoint, &yes) == 0 && yes == 1); + fermionic_op_free(batch); + fermionic_op_free(single); + fermionic_op_free(a); + fermionic_op_free(c); + fermionic_op_free(ac); + fermionic_op_free(ca); + fermionic_op_free(sum); + fermionic_op_free(identity); + fermionic_op_free(ordered); + fermionic_op_free(copy); + fermionic_op_free(number); + fermionic_op_free(total); + fermionic_op_free(zero); + fermionic_op_free(scaled); + fermionic_op_free(difference); + fermionic_op_free(adjoint); + fermionic_op_free(NULL); +} + +static void test_pauli_algebra(void) { + CPauliString *px = pauli_string_parse("X"), *pz = pauli_string_parse("Z"); + CPauliString *py = pauli_string_parse("Y"), *pxz = pauli_string_parse("XZ"); + assert(px && pz && py && pxz); + CHamiltonianTerm tx = {px, {1, 0}}, tz = {pz, {1, 0}}; + CHamiltonianTerm ty = {py, {0, -1}}, txz = {pxz, {1, 0}}; + CHamiltonian *x = NULL, *z = NULL, *y = NULL, *expected = NULL; + CHamiltonian *product = NULL, *tensor = NULL, *sum = NULL, *difference = NULL; + CHamiltonian *adjoint = NULL, *scaled = NULL, *copy = NULL, *normalized = NULL; + CHamiltonian *hermitian = NULL, *remapped = NULL, *identity = NULL; + uint8_t yes = 99; + assert(hamiltonian_from_terms(1, &tx, 1, &x) == 0); + assert(hamiltonian_from_terms(1, &tz, 1, &z) == 0); + assert(hamiltonian_from_terms(1, &ty, 1, &y) == 0); + assert(hamiltonian_from_terms(2, &txz, 1, &expected) == 0); + pauli_string_free(px); + pauli_string_free(pz); + pauli_string_free(py); + pauli_string_free(pxz); + assert(hamiltonian_mul(x, z, &product) == 0); + assert(hamiltonian_equiv(product, y, 0, &yes) == 0 && yes == 1); + assert(hamiltonian_tensor(x, z, &tensor) == 0); + assert(hamiltonian_equiv(tensor, expected, 0, &yes) == 0 && yes == 1); + assert(hamiltonian_add(x, z, &sum) == 0); + assert(hamiltonian_sub(sum, z, &difference) == 0); + assert(hamiltonian_equiv(difference, x, 0, &yes) == 0 && yes == 1); + assert(hamiltonian_adjoint(product, &adjoint) == 0); + assert(hamiltonian_scaled(product, (Complex64){-1, 0}, &scaled) == 0); + assert(hamiltonian_equiv(adjoint, scaled, 0, &yes) == 0 && yes == 1); + assert(hamiltonian_copy(x, ©) == 0); + assert(hamiltonian_chop(copy, 1) == 0 && hamiltonian_num_terms(copy) == 0); + assert(hamiltonian_num_terms(x) == 1); + assert(hamiltonian_normalized(x, &normalized) == 0); + assert(hamiltonian_normalized_hermitian(x, &hermitian) == 0); + assert(hamiltonian_is_hermitian(product, 0, &yes) == 0 && yes == 0); + assert(hamiltonian_is_hermitian(hermitian, 0, &yes) == 0 && yes == 1); + const size_t mapping[] = {1}; + assert(hamiltonian_remap_qubits(x, mapping, 1, 2, &remapped) == 0); + assert(hamiltonian_identity(0, &identity) == 0); + CPauliString* snapshot = NULL; + Complex64 coefficient = {99, 99}; + assert(hamiltonian_get_term(identity, 0, &snapshot, &coefficient) == 0); + hamiltonian_free(identity); + assert(pauli_string_num_qubits(snapshot) == 0 && coefficient.re == 1); + pauli_string_free(snapshot); + CHamiltonian* out = x; + assert(hamiltonian_add(x, tensor, &out) == -8 && out == x); + assert(hamiltonian_normalized_hermitian(product, &out) == -8 && out == x); + assert(hamiltonian_identity(SIZE_MAX, &out) == -7 && out == x); + assert(hamiltonian_equiv(x, z, NAN, &yes) == -8 && yes == 1); + hamiltonian_free(x); + hamiltonian_free(z); + hamiltonian_free(y); + hamiltonian_free(expected); + hamiltonian_free(product); + hamiltonian_free(tensor); + hamiltonian_free(sum); + hamiltonian_free(difference); + hamiltonian_free(adjoint); + hamiltonian_free(scaled); + hamiltonian_free(copy); + hamiltonian_free(normalized); + hamiltonian_free(hermitian); + hamiltonian_free(remapped); +} + +static void test_mapping_and_evolution(void) { + size_t sources[] = {0}, targets[] = {1}, width = 99; + CFermionicOp* generator = NULL; + CHamiltonian* mapped = NULL; + assert(fermionic_op_excitation_generator(2, sources, 1, targets, 1, &generator) == 0); + assert(jordan_wigner_num_qubits(2, &width) == 0 && width == 2); + assert(jordan_wigner_map_operator(generator, &mapped) == 0); + uint8_t occupation[] = {1, 0}, bits[] = {99, 99}; + assert(jordan_wigner_map_occupation(occupation, 2, bits, 2) == 0 && bits[0] == 1 && + bits[1] == 0); + occupation[1] = 2; + assert(jordan_wigner_map_occupation(occupation, 2, bits, 2) == -8 && bits[0] == 1 && + bits[1] == 0); + assert(jordan_wigner_map_occupation(NULL, 0, NULL, 0) == 0); + const double theta = 0.37; + CCircuit* evolution = + hamiltonian_to_evolution_circuit(mapped, theta, 1, TROTTER_MODE_SECOND_ORDER); + assert(evolution); + Complex64 state[] = {{0, 0}, {1, 0}, {0, 0}, {0, 0}}, data[4]; + CStatevector* sv = statevector_from_state(2, state, 4); + assert(sv); + assert(statevector_apply_circuit(sv, evolution) == 0); + assert(statevector_data(sv, data, 4) == 0); + assert(fabs(data[1].re - cos(theta)) < 1e-12 && fabs(data[1].im) < 1e-12); + assert(fabs(data[2].re - sin(theta)) < 1e-12 && fabs(data[2].im) < 1e-12); + assert(hypot(data[0].re, data[0].im) < 1e-12 && hypot(data[3].re, data[3].im) < 1e-12); + statevector_free(sv); + circuit_free(evolution); + hamiltonian_free(mapped); + fermionic_op_free(generator); + CCircuit* phase = circuit_new(1); + assert(circuit_set_global_phase(phase, theta) == 0); + sv = statevector_from_circuit(phase); + assert(sv && statevector_data(sv, data, 2) == 0); + assert(fabs(data[0].re - cos(theta)) < 1e-12 && fabs(data[0].im - sin(theta)) < 1e-12); + statevector_free(sv); + circuit_free(phase); +} + +static void test_wide_sparse_registers(void) { + const size_t bits = sizeof(size_t) * CHAR_BIT; + const size_t widths[] = {bits, bits + 1, bits - 1, bits / 2, bits / 2 - 1, bits / 2 - 2}; + for (size_t i = 0; i < sizeof(widths) / sizeof(widths[0]); ++i) { + CFermionicOp* identity = NULL; + CHamiltonian* mapped = NULL; + CPauliString* pauli = NULL; + Complex64 coefficient = {0, 0}, buffer = {99, -99}; + assert(fermionic_op_identity(widths[i], &identity) == 0); + assert(jordan_wigner_map_operator(identity, &mapped) == 0); + assert(hamiltonian_get_term(mapped, 0, &pauli, &coefficient) == 0); + const size_t dim = widths[i] < bits ? (size_t)1 << widths[i] : 0; + assert(hamiltonian_to_matrix_len(mapped) == dim); + assert(pauli_string_matrix_len(pauli) == dim); + assert(hamiltonian_to_matrix(mapped, &buffer, 1) == -7); + assert(buffer.re == 99 && buffer.im == -99); + assert(pauli_string_matrix(pauli, &buffer, 1) == -7); + assert(buffer.re == 99 && buffer.im == -99); + pauli_string_free(pauli); + hamiltonian_free(mapped); + fermionic_op_free(identity); + } + CFermionicOp* zero = NULL; + CHamiltonian* mapped = NULL; + assert(fermionic_op_new(SIZE_MAX, &zero) == 0); + assert(jordan_wigner_map_operator(zero, &mapped) == 0); + assert(hamiltonian_to_matrix_len(mapped) == 0); + Complex64 buffer = {99, -99}; + assert(hamiltonian_to_matrix(mapped, &buffer, 1) == -7); + assert(buffer.re == 99 && buffer.im == -99); + hamiltonian_free(mapped); + fermionic_op_free(zero); +} + +int main(void) { + test_fermions(); + test_pauli_algebra(); + test_mapping_and_evolution(); + test_wide_sparse_registers(); + puts("binding-c checked algebra tests passed"); + return 0; +} diff --git a/docs/api/cn/c/3_qis/0_overview.md b/docs/api/cn/c/3_qis/0_overview.md index 18e0c9fe..3807ce50 100644 --- a/docs/api/cn/c/3_qis/0_overview.md +++ b/docs/api/cn/c/3_qis/0_overview.md @@ -84,3 +84,5 @@ int main(void) { ``` 噪声模型句柄 `CNoiseModel` 的构造与噪声登记见 [NoiseModel](../2_device/4_noise.md)。 + +新增的[费米子代数与 Jordan–Wigner](10_fermion.md)和 [checked Hamiltonian 代数](7_hamiltonian.md#带校验的稀疏-pauli-代数)使用“错误码 + 输出指针”接口；构造结果通过输出指针返回，判断通过字节输出写入 0/1。旧 QIS 接口保持原有约定。 diff --git a/docs/api/cn/c/3_qis/10_fermion.md b/docs/api/cn/c/3_qis/10_fermion.md new file mode 100644 index 00000000..888a3d8d --- /dev/null +++ b/docs/api/cn/c/3_qis/10_fermion.md @@ -0,0 +1,179 @@ +# 费米子代数与 Jordan–Wigner（C） + +`CFermionicOp` 是有序升降算符乘积的稀疏复数线性组合。模式从零编号；乘积 `A B` 表示 B 先作用。相同有序乘积自动合并，仅删除精确零；不会隐式正规序或截断小系数。零算符和恒等算符都保留声明的模式数，支持零模式。 + +## ABI 与数据 + +所有返回 `int32_t` 的新接口以 `0` 表示成功。对象通过 `**out` 返回，标量通过输出指针返回；失败不改变输出或原对象。对象结果独立拥有，使用 `fermionic_op_free` 释放；输入借用并复制。调用方传入空输出槽，并提供有效、对齐、存活的指针；输出存储不得相互重叠或与句柄内部存储重叠。同一句柄被修改时，不可并发读取。数组长度为零允许 NULL；非空数组必须指向足够存储。 + +错误码：`-1` 必需指针为空，`-2` 模式或项索引越界，`-8` 非法标签、模式数不匹配、非有限系数或算术溢出、非法容差、容量不足或无法表示的数组长度，`-7` 不支持的 Pauli 表示规模。判断通过 `uint8_t *out` 写入 0/1，成功返回值始终为 0。 + +```c +typedef struct CFermionicOp CFermionicOp; +#define FERMIONIC_CREATE 0 +#define FERMIONIC_ANNIHILATE 1 +typedef struct CFermionicLadder { + size_t mode; + uint8_t kind; +} CFermionicLadder; +typedef struct CFermionicTerm { + const CFermionicLadder *ladders; + size_t len; + Complex64 coefficient; +} CFermionicTerm; +``` + +## 构造 + +```c +int32_t fermionic_op_new(size_t num_modes, CFermionicOp **out); + +int32_t fermionic_op_identity(size_t num_modes, CFermionicOp **out); + +int32_t fermionic_op_from_term(size_t num_modes, + const CFermionicLadder *term, + size_t len, + Complex64 coefficient, + CFermionicOp **out); + +int32_t fermionic_op_from_terms(size_t num_modes, + const CFermionicTerm *terms, + size_t len, + CFermionicOp **out); + +int32_t fermionic_op_creation(size_t num_modes, size_t mode, CFermionicOp **out); + +int32_t fermionic_op_annihilation(size_t num_modes, size_t mode, CFermionicOp **out); + +int32_t fermionic_op_number(size_t num_modes, size_t mode, CFermionicOp **out); + +int32_t fermionic_op_total_number(size_t num_modes, CFermionicOp **out); + +int32_t fermionic_op_excitation_generator(size_t num_modes, + const size_t *sources, + size_t sources_len, + const size_t *targets, + size_t targets_len, + CFermionicOp **out); +``` + +- `new` 创建零，`identity` 创建一；`from_term` 的空乘积为 identity，`from_terms` 的空求和为 zero。批量输入及嵌套升降算符数组只在调用期间借用。即使系数为零，也校验标签和模式。 +- `creation/annihilation/number` 分别构造 a†、a、a†a；`total_number` 为所有模式 number 之和。 +- 激发生成元 `G=i(T−T†)`，其中 T 的 target 创建算符按输入顺序排列，source 湮灭算符按逆序排列。两侧必须非空、等长、各自无重复且互不重叠。奇排列改变符号，无排序、阶乘或自旋轨道约定。 + +## 生命周期与读取 + +```c +void fermionic_op_free(CFermionicOp *ptr); + +int32_t fermionic_op_copy(const CFermionicOp *ptr, CFermionicOp **out); + +int32_t fermionic_op_num_modes(const CFermionicOp *ptr, size_t *out); + +int32_t fermionic_op_num_terms(const CFermionicOp *ptr, size_t *out); + +int32_t fermionic_op_is_zero(const CFermionicOp *ptr, uint8_t *out); + +int32_t fermionic_op_equals(const CFermionicOp *left, + const CFermionicOp *right, + uint8_t *out); + +int32_t fermionic_op_term_len(const CFermionicOp *ptr, size_t index, size_t *out); + +int32_t fermionic_op_get_term(const CFermionicOp *ptr, + size_t index, + CFermionicLadder *buffer, + size_t capacity, + Complex64 *out_coefficient); +``` + +`free(NULL)` 合法；`copy` 返回独立副本。`is_zero` 检查存储是否为空，`equals` 比较存储结构和模式数（模式数不同返回 false）。代数比较使用 `equiv`。 + +先调用 `term_len` 获取项长度，再分配 buffer 并调用 `get_term`；capacity 可以大于所需长度，额外元素保持原值。空乘积允许 NULL buffer。按 core 的确定性结构顺序读取；修改算符后应重新查询长度和索引。输出数组是副本，没有需要释放的内部指针。 + +## 代数 + +```c +int32_t fermionic_op_add_term(CFermionicOp *ptr, + const CFermionicLadder *term, + size_t len, + Complex64 coefficient); + +int32_t fermionic_op_add(const CFermionicOp *left, + const CFermionicOp *right, + CFermionicOp **out); + +int32_t fermionic_op_sub(const CFermionicOp *left, + const CFermionicOp *right, + CFermionicOp **out); + +int32_t fermionic_op_mul(const CFermionicOp *left, + const CFermionicOp *right, + CFermionicOp **out); + +int32_t fermionic_op_scaled(const CFermionicOp *ptr, + Complex64 factor, + CFermionicOp **out); + +int32_t fermionic_op_adjoint(const CFermionicOp *ptr, CFermionicOp **out); + +int32_t fermionic_op_normal_order(const CFermionicOp *ptr, CFermionicOp **out); + +int32_t fermionic_op_chop(CFermionicOp *ptr, double atol); + +int32_t fermionic_op_equiv(const CFermionicOp *left, + const CFermionicOp *right, + double atol, + uint8_t *out); + +int32_t fermionic_op_is_hermitian(const CFermionicOp *ptr, double atol, uint8_t *out); +``` + +`add_term/chop` 原地修改，其余对象运算返回副本。`mul` 不自动正规序；`adjoint` 共轭系数、逆转乘积并交换升降种类。`normal_order` 将创建算符放在湮灭算符前，每组按模式升序，包含反对易收缩项。 + +`chop` 删除模长 ≤ atol 的系数。`equiv` 对差正规序后比较系数，要求模式数相同；`is_hermitian` 与伴随比较。atol 必须有限且非负，精确比较显式传入 0。它是系数容差，不是矩阵算符范数误差界。 + +## Jordan–Wigner + +```c +int32_t jordan_wigner_map_operator(const CFermionicOp *ptr, CHamiltonian **out); + +int32_t jordan_wigner_num_qubits(size_t num_modes, size_t *out); + +int32_t jordan_wigner_map_occupation(const uint8_t *occupation, + size_t len, + uint8_t *out_bits, + size_t capacity); +``` + +无状态映射无需句柄。mode p 对应 qubit p，Z 链覆盖所有更低模式：creation = Z链(X−iY)/2，annihilation = Z链(X+iY)/2。qubit 0 为最低有效位，也是显示的 Pauli 字符串最右字符。 + +`map_operator` 返回独立 `CHamiltonian`，使用 `hamiltonian_free` 释放；非厄米输入合法，厄米性要求由期望值和演化等下游接口检查。映射不构造稠密矩阵，但 k 个升降算符的乘积在合并前可能展开为 2^k 项。 + +`map_occupation` 仅接受 0/1，占据长度决定模式数；输出字节 bits[p] 对应 qubit p。capacity 至少为 len，输入输出允许重叠；空输入为零模式真空。 + +## 示例 + +```c +#include "cqlib_c.h" + +int main(void) { + const size_t sources[] = {0}, targets[] = {1}; + CFermionicOp *g = NULL; + CHamiltonian *h = NULL; + if (fermionic_op_excitation_generator(2, sources, 1, targets, 1, &g) != 0) { + return 1; + } + if (jordan_wigner_map_operator(g, &h) != 0) { + fermionic_op_free(g); + return 1; + } + CCircuit *u = hamiltonian_to_evolution_circuit( + h, 0.37, 1, TROTTER_MODE_SECOND_ORDER); + const int result = u == NULL ? 1 : 0; + circuit_free(u); + hamiltonian_free(h); + fermionic_op_free(g); + return result; +} +``` diff --git a/docs/api/cn/c/3_qis/6_pauli.md b/docs/api/cn/c/3_qis/6_pauli.md index 5079672e..ceae851b 100644 --- a/docs/api/cn/c/3_qis/6_pauli.md +++ b/docs/api/cn/c/3_qis/6_pauli.md @@ -132,13 +132,15 @@ int32_t pauli_mul_with_phase(uint8_t left, uint8_t right, uint8_t *out_pauli, Co ### pauli_string_matrix_len(ptr) -返回稠密矩阵的边长 `2^N`；NULL 返回 0。 +返回稠密矩阵的边长 `2^N`；NULL 或边长无法表示时返回 0。 ### pauli_string_matrix(ptr, buffer, len) 把 `2^N × 2^N` 稠密矩阵按行主序拷入 `buffer`（`Complex64` 值）。`len` 必须等于 `pauli_string_matrix_len(ptr) * pauli_string_matrix_len(ptr)`。 -返回：0 成功；-1 空指针；-8 `len` 与矩阵元素数不符。 +返回：0 成功；-1 空指针；-7 元素数或字节数无法表示，或分配失败；-8 `len` 与矩阵元素数不符。失败不改写缓冲区。 + +可表示的边长不保证完整矩阵可以分配；调用方应在计算 `dim * dim` 前检查溢出。 --- @@ -169,6 +171,7 @@ int32_t pauli_mul_with_phase(uint8_t left, uint8_t right, uint8_t *out_pauli, Co ## 示例 ```c +#include #include #include #include "cqlib_c.h" @@ -192,6 +195,11 @@ int main(void) { struct CPauliString *zz = pauli_string_parse("ZZI"); int32_t commutes = pauli_string_commutes_with(ps, zz); /* 0（比特 2 上 X 与 Z 反对易） */ + printf("%s: qubits=%" PRIuPTR ", p2=%u, p0=%u\n", text, n, + (unsigned)p2, (unsigned)p0); + printf("x_mask=%" PRIuPTR ", z_mask=%" PRIuPTR ", commutes=%" PRId32 "\n", + x_mask, z_mask, commutes); + /* 4. 两步式读取 8x8 稠密矩阵 */ uintptr_t dim = pauli_string_matrix_len(ps); /* 8 */ Complex64 *matrix = malloc(dim * dim * sizeof(Complex64)); diff --git a/docs/api/cn/c/3_qis/7_hamiltonian.md b/docs/api/cn/c/3_qis/7_hamiltonian.md index ad066c5f..27e9e7fd 100644 --- a/docs/api/cn/c/3_qis/7_hamiltonian.md +++ b/docs/api/cn/c/3_qis/7_hamiltonian.md @@ -99,9 +99,11 @@ Trotter-Suzuki 分解模式标签（`hamiltonian_to_trotter_circuit` 与 `hamilt ### hamiltonian_to_matrix_len(ptr) / hamiltonian_to_matrix(ptr, buffer, len) -两步式读取稠密矩阵：`*_len` 返回边长 `2^N`（NULL 返回 0），填充接口把 `2^N × 2^N` 矩阵按行主序拷入 `buffer`（`Complex64` 值）。`len` 必须等于边长的平方。 +两步式读取稠密矩阵：`*_len` 返回边长 `2^N`（NULL 或边长无法表示时返回 0），填充接口把 `2^N × 2^N` 矩阵按行主序拷入 `buffer`（`Complex64` 值）。`len` 必须等于边长的平方。 -返回（填充接口）：0 成功；-1 空指针；-8 `len` 与矩阵元素数不符。 +返回（填充接口）：0 成功；-1 空指针；-7 元素数或字节数无法表示，或分配失败；-8 `len` 与矩阵元素数不符。失败不改写缓冲区。 + +可表示的边长不保证完整矩阵可以分配；调用方应在计算 `dim * dim` 前检查溢出。 --- @@ -212,6 +214,7 @@ int32_t hamiltonian_variance_statevector(const struct CHamiltonian *ptr, ## 示例 ```c +#include #include #include #include "cqlib_c.h" @@ -230,6 +233,8 @@ int main(void) { uintptr_t nq = hamiltonian_num_qubits(ham); /* 2 */ uintptr_t nt = hamiltonian_num_terms(ham); /* 2 */ int32_t commute = hamiltonian_all_terms_commute(ham); /* 1（两单项对易） */ + printf("qubits=%" PRIuPTR ", terms=%" PRIuPTR ", commute=%" PRId32 "\n", + nq, nt, commute); hamiltonian_simplify(ham); /* 3. 两步式读取 4x4 稠密矩阵 */ @@ -280,3 +285,86 @@ int main(void) { ``` 哈密顿量作为可观测量用于误差缓解流程（见 [零噪声外推](../6_error_mitigation/1_zne.md)、[虚拟蒸馏](../6_error_mitigation/2_virtual_distillation.md)），其期望值可直接由态矢量模拟器计算（见 [态矢量](1_statevector.md)）。 + + +## 带校验的稀疏 Pauli 代数 + +`CHamiltonian` 同时表示一般非厄米稀疏算符。以下新接口采用与[费米子绑定](10_fermion.md)相同的错误码和输出约定：成功为 0，通过 `**out` 返回独立句柄，失败不改写输出或输入。输入 Pauli 句柄借用并复制，不消耗。保留旧 `from_pauli` 的所有权转移、原地 `scale` 和近零截断 `simplify` 行为。 + +```c +typedef struct CHamiltonianTerm { + const CPauliString *pauli; + Complex64 coefficient; +} CHamiltonianTerm; + +int32_t hamiltonian_identity(size_t num_qubits, CHamiltonian **out); + +int32_t hamiltonian_from_terms(size_t num_qubits, + const CHamiltonianTerm *terms, + size_t len, + CHamiltonian **out); + +int32_t hamiltonian_copy(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_normalized(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_normalized_hermitian(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_add(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_sub(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_mul(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_scaled(const CHamiltonian *ptr, + Complex64 factor, + CHamiltonian **out); + +int32_t hamiltonian_adjoint(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_tensor(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_remap_qubits(const CHamiltonian *ptr, + const size_t *mapping, + size_t len, + size_t target_width, + CHamiltonian **out); + +int32_t hamiltonian_chop(CHamiltonian *ptr, double atol); + +int32_t hamiltonian_equiv(const CHamiltonian *left, + const CHamiltonian *right, + double atol, + uint8_t *out); + +int32_t hamiltonian_is_hermitian(const CHamiltonian *ptr, double atol, uint8_t *out); + +int32_t hamiltonian_get_term(const CHamiltonian *ptr, + size_t index, + CPauliString **out_pauli, + Complex64 *out_coefficient); +``` + +| 接口 | 行为 | +| --- | --- | +| `identity/from_terms` | identity 支持零量子位标量一；空 terms 保留声明宽度。批量构造校验每项宽度和系数、吸收 Pauli 相位、合并重复项，仅删除精确零。 | +| `copy/get_term` | 独立存储副本。get_term 返回项的 Pauli 副本及系数，不自动规范化；Pauli 用 pauli_string_free 释放，零宽度也合法。索引错误为 -2。 | +| `normalized` | 返回经过校验和精确合并的副本，无近零截断。 | +| `normalized_hermitian` | 规范化后将绝对值 ≤1e-10 的虚部投影到零，更大的虚部返回 -8；保留非零实系数。 | +| `add/sub/mul/scaled/adjoint` | 调用 core checked 方法；校验输入与结果系数，非有限值或溢出返回 -8。二元运算要求宽度相同，mul 右侧先作用。 | +| `tensor` | 左侧 ⊗ 右侧，右侧占低位 qubit；宽度相加，溢出返回 -8。 | +| `remap_qubits` | mapping[source]=destination，长度等于输入宽度，目标唯一且小于 target_width；空闲 qubit 填 identity。不是费米模式置换。 | +| `chop` | 原地规范化后删除模长 ≤atol 的系数；失败不修改对象。 | +| `equiv/is_hermitian` | 写入 0/1；atol 有限且非负，0 为精确比较。equiv 要求宽度相同；厄米性检查等价于与伴随比较，不投影系数。 | + +对象结果用 `hamiltonian_free` 释放，输出槽须为空且有效。必需指针为空返回 -1；参数、宽度、系数或容差错误返回 -8；identity、tensor、remap 无法表示的 Pauli 位向量宽度返回 -7。零求和可以保留无需物化的较大声明宽度。 + +atol 是系数容差；规范化后 `is_hermitian` 使用 `2*abs(imag) <= atol`，与 `normalized_hermitian` 的投影阈值是两个独立契约。稀疏代数不构造稠密矩阵，也不隐式调用旧 simplify。 diff --git a/docs/api/en/c/3_qis/0_overview.md b/docs/api/en/c/3_qis/0_overview.md index cc94ee7b..b6f69def 100644 --- a/docs/api/en/c/3_qis/0_overview.md +++ b/docs/api/en/c/3_qis/0_overview.md @@ -84,3 +84,5 @@ int main(void) { ``` For constructing and populating the `CNoiseModel` handle, see [NoiseModel](../2_device/4_noise.md). + +New [fermionic algebra and Jordan–Wigner](10_fermion.md) and [checked Hamiltonian algebra](7_hamiltonian.md#checked-sparse-pauli-algebra) use error codes plus output pointers: constructors deliver handles through outputs and predicates write 0/1 bytes. Existing QIS functions retain their original conventions. diff --git a/docs/api/en/c/3_qis/10_fermion.md b/docs/api/en/c/3_qis/10_fermion.md new file mode 100644 index 00000000..44822fab --- /dev/null +++ b/docs/api/en/c/3_qis/10_fermion.md @@ -0,0 +1,179 @@ +# Fermionic algebra and Jordan–Wigner (C) + +`CFermionicOp` is a sparse complex sum of ordered ladder products. Modes are zero-based; in `A B`, B acts first. Identical ordered products merge and only exact zeros disappear. No implicit normal ordering or numerical pruning occurs. Zero and identity retain the declared mode count, including zero modes. + +## ABI and data + +New `int32_t` functions return 0 on success. Objects use `**out` and scalar queries use output pointers; errors leave outputs and mutable inputs unchanged. Results are independently owned and freed with `fermionic_op_free`; inputs are borrowed and copied. Pass an empty output slot and valid, aligned, live storage. Output slots must not overlap one another or handle internals, and a handle must not be read concurrently with mutation. Empty arrays allow NULL; nonempty arrays require sufficient storage. + +Errors: -1 required NULL pointer; -2 mode/term index out of bounds; -8 invalid tag, mode-count mismatch, nonfinite coefficient/arithmetic overflow, invalid tolerance, insufficient capacity or unrepresentable array length; -7 unsupported Pauli representation width. Predicates write 0/1 to `uint8_t *out`; their successful return code is always 0. + +```c +typedef struct CFermionicOp CFermionicOp; +#define FERMIONIC_CREATE 0 +#define FERMIONIC_ANNIHILATE 1 +typedef struct CFermionicLadder { + size_t mode; + uint8_t kind; +} CFermionicLadder; +typedef struct CFermionicTerm { + const CFermionicLadder *ladders; + size_t len; + Complex64 coefficient; +} CFermionicTerm; +``` + +## Construction + +```c +int32_t fermionic_op_new(size_t num_modes, CFermionicOp **out); + +int32_t fermionic_op_identity(size_t num_modes, CFermionicOp **out); + +int32_t fermionic_op_from_term(size_t num_modes, + const CFermionicLadder *term, + size_t len, + Complex64 coefficient, + CFermionicOp **out); + +int32_t fermionic_op_from_terms(size_t num_modes, + const CFermionicTerm *terms, + size_t len, + CFermionicOp **out); + +int32_t fermionic_op_creation(size_t num_modes, size_t mode, CFermionicOp **out); + +int32_t fermionic_op_annihilation(size_t num_modes, size_t mode, CFermionicOp **out); + +int32_t fermionic_op_number(size_t num_modes, size_t mode, CFermionicOp **out); + +int32_t fermionic_op_total_number(size_t num_modes, CFermionicOp **out); + +int32_t fermionic_op_excitation_generator(size_t num_modes, + const size_t *sources, + size_t sources_len, + const size_t *targets, + size_t targets_len, + CFermionicOp **out); +``` + +- `new` creates zero and `identity` creates one. An empty `from_term` product is identity; an empty `from_terms` sum is zero. Batch descriptors and nested ladder arrays are borrowed only during the call. Even zero-weight products validate tags and modes. +- `creation/annihilation/number` construct a†, a, and a†a; `total_number` sums the number operator over all modes. +- Excitation generates G=i(T−T†), with target creators in input order and source annihilators in reverse order. Lists must be nonempty, equally long, individually unique and disjoint. An odd permutation changes the sign. No sorting, factorial or spin-orbital convention is imposed. + +## Lifetime and queries + +```c +void fermionic_op_free(CFermionicOp *ptr); + +int32_t fermionic_op_copy(const CFermionicOp *ptr, CFermionicOp **out); + +int32_t fermionic_op_num_modes(const CFermionicOp *ptr, size_t *out); + +int32_t fermionic_op_num_terms(const CFermionicOp *ptr, size_t *out); + +int32_t fermionic_op_is_zero(const CFermionicOp *ptr, uint8_t *out); + +int32_t fermionic_op_equals(const CFermionicOp *left, + const CFermionicOp *right, + uint8_t *out); + +int32_t fermionic_op_term_len(const CFermionicOp *ptr, size_t index, size_t *out); + +int32_t fermionic_op_get_term(const CFermionicOp *ptr, + size_t index, + CFermionicLadder *buffer, + size_t capacity, + Complex64 *out_coefficient); +``` + +`free(NULL)` is valid; `copy` is independent. `is_zero` tests empty storage and `equals` compares stored structure and mode count (different mode counts return false). Use `equiv` for algebraic equality. + +Call `term_len`, allocate a buffer, then call `get_term`. Capacity may exceed the required length; trailing elements stay untouched. Empty products allow NULL buffer. Terms follow core deterministic structural order; requery lengths and indices after mutation. Output arrays are copies and contain no library-owned pointers to release. + +## Algebra + +```c +int32_t fermionic_op_add_term(CFermionicOp *ptr, + const CFermionicLadder *term, + size_t len, + Complex64 coefficient); + +int32_t fermionic_op_add(const CFermionicOp *left, + const CFermionicOp *right, + CFermionicOp **out); + +int32_t fermionic_op_sub(const CFermionicOp *left, + const CFermionicOp *right, + CFermionicOp **out); + +int32_t fermionic_op_mul(const CFermionicOp *left, + const CFermionicOp *right, + CFermionicOp **out); + +int32_t fermionic_op_scaled(const CFermionicOp *ptr, + Complex64 factor, + CFermionicOp **out); + +int32_t fermionic_op_adjoint(const CFermionicOp *ptr, CFermionicOp **out); + +int32_t fermionic_op_normal_order(const CFermionicOp *ptr, CFermionicOp **out); + +int32_t fermionic_op_chop(CFermionicOp *ptr, double atol); + +int32_t fermionic_op_equiv(const CFermionicOp *left, + const CFermionicOp *right, + double atol, + uint8_t *out); + +int32_t fermionic_op_is_hermitian(const CFermionicOp *ptr, double atol, uint8_t *out); +``` + +`add_term/chop` mutate in place; other object operations return copies. `mul` does not normal order. `adjoint` conjugates coefficients, reverses products and exchanges ladder kinds. `normal_order` places creators before annihilators with ascending modes within each group, including anticommutation contractions. + +`chop` removes coefficients with magnitude <= atol. `equiv` normal-orders the difference and compares coefficients, requiring equal mode counts; `is_hermitian` compares with the adjoint. atol must be finite and nonnegative; explicitly pass 0 for exact comparison. This is a coefficient tolerance, not a matrix operator-norm bound. + +## Jordan–Wigner + +```c +int32_t jordan_wigner_map_operator(const CFermionicOp *ptr, CHamiltonian **out); + +int32_t jordan_wigner_num_qubits(size_t num_modes, size_t *out); + +int32_t jordan_wigner_map_occupation(const uint8_t *occupation, + size_t len, + uint8_t *out_bits, + size_t capacity); +``` + +Mapping is stateless and needs no mapper handle. Mode p maps to qubit p with Z on lower modes: creation = Z-chain (X−iY)/2 and annihilation = Z-chain (X+iY)/2. Qubit 0 is the least significant bit and rightmost displayed Pauli character. + +`map_operator` returns an independent `CHamiltonian`, freed with `hamiltonian_free`. Non-Hermitian inputs are valid; downstream expectation and evolution functions check Hermiticity. Mapping builds no dense matrix, but a product of k ladders may expand to 2^k terms before merging. + +`map_occupation` accepts only 0/1 bytes. Input length defines the mode count; output bits[p] addresses qubit p. Capacity must be at least len; input and output may overlap. Empty input is zero-mode vacuum. + +## Example + +```c +#include "cqlib_c.h" + +int main(void) { + const size_t sources[] = {0}, targets[] = {1}; + CFermionicOp *g = NULL; + CHamiltonian *h = NULL; + if (fermionic_op_excitation_generator(2, sources, 1, targets, 1, &g) != 0) { + return 1; + } + if (jordan_wigner_map_operator(g, &h) != 0) { + fermionic_op_free(g); + return 1; + } + CCircuit *u = hamiltonian_to_evolution_circuit( + h, 0.37, 1, TROTTER_MODE_SECOND_ORDER); + const int result = u == NULL ? 1 : 0; + circuit_free(u); + hamiltonian_free(h); + fermionic_op_free(g); + return result; +} +``` diff --git a/docs/api/en/c/3_qis/6_pauli.md b/docs/api/en/c/3_qis/6_pauli.md index d42ac6f3..e94853e9 100644 --- a/docs/api/en/c/3_qis/6_pauli.md +++ b/docs/api/en/c/3_qis/6_pauli.md @@ -132,13 +132,15 @@ Returns: 1 when they commute; 0 otherwise; -1 for NULL handles; -8 when the qubi ### pauli_string_matrix_len(ptr) -Returns the side length of the dense matrix (`2^N`); 0 for NULL. +Returns the side length of the dense matrix (`2^N`); 0 for NULL or an unrepresentable side length. ### pauli_string_matrix(ptr, buffer, len) Copies the `2^N × 2^N` dense matrix into `buffer` in row-major order (`Complex64` values). `len` must equal `pauli_string_matrix_len(ptr) * pauli_string_matrix_len(ptr)`. -Returns: 0 on success; -1 for NULL pointers; -8 when `len` does not match the matrix element count. +Returns: 0 on success; -1 for NULL pointers; -7 for unrepresentable element/byte counts or allocation failure; -8 when `len` does not match the matrix element count. Errors leave the buffer unchanged. + +A representable side does not guarantee the complete matrix can be allocated. Check for overflow before computing `dim * dim`. --- @@ -169,6 +171,7 @@ Returns: 0 on success; -1 when `ptr` or `out` is NULL, or when either array is N ## Example ```c +#include #include #include #include "cqlib_c.h" @@ -192,6 +195,11 @@ int main(void) { struct CPauliString *zz = pauli_string_parse("ZZI"); int32_t commutes = pauli_string_commutes_with(ps, zz); /* 0 (X and Z anticommute on qubit 2) */ + printf("%s: qubits=%" PRIuPTR ", p2=%u, p0=%u\n", text, n, + (unsigned)p2, (unsigned)p0); + printf("x_mask=%" PRIuPTR ", z_mask=%" PRIuPTR ", commutes=%" PRId32 "\n", + x_mask, z_mask, commutes); + /* 4. Two-step read of the 8x8 dense matrix */ uintptr_t dim = pauli_string_matrix_len(ps); /* 8 */ Complex64 *matrix = malloc(dim * dim * sizeof(Complex64)); diff --git a/docs/api/en/c/3_qis/7_hamiltonian.md b/docs/api/en/c/3_qis/7_hamiltonian.md index 08f1b73b..6688f7ec 100644 --- a/docs/api/en/c/3_qis/7_hamiltonian.md +++ b/docs/api/en/c/3_qis/7_hamiltonian.md @@ -99,9 +99,11 @@ Returns: 1 when all terms commute; 0 otherwise; -1 for NULL. ### hamiltonian_to_matrix_len(ptr) / hamiltonian_to_matrix(ptr, buffer, len) -Two-step read of the dense matrix: `*_len` returns the side length `2^N` (0 for NULL), and the fill function copies the `2^N × 2^N` matrix into `buffer` in row-major order (`Complex64` values). `len` must equal the square of the side length. +Two-step read of the dense matrix: `*_len` returns the side length `2^N` (0 for NULL or an unrepresentable side length), and the fill function copies the `2^N × 2^N` matrix into `buffer` in row-major order (`Complex64` values). `len` must equal the square of the side length. -Returns (fill function): 0 on success; -1 for NULL pointers; -8 when `len` does not match the matrix element count. +Returns (fill function): 0 on success; -1 for NULL pointers; -7 for unrepresentable element/byte counts or allocation failure; -8 when `len` does not match the matrix element count. Errors leave the buffer unchanged. + +A representable side does not guarantee the complete matrix can be allocated. Check for overflow before computing `dim * dim`. --- @@ -212,6 +214,7 @@ Returns: a new `CCircuit*` on success (free with `circuit_free`); NULL on NULL i ## Example ```c +#include #include #include #include "cqlib_c.h" @@ -230,6 +233,8 @@ int main(void) { uintptr_t nq = hamiltonian_num_qubits(ham); /* 2 */ uintptr_t nt = hamiltonian_num_terms(ham); /* 2 */ int32_t commute = hamiltonian_all_terms_commute(ham); /* 1 (single-qubit terms commute) */ + printf("qubits=%" PRIuPTR ", terms=%" PRIuPTR ", commute=%" PRId32 "\n", + nq, nt, commute); hamiltonian_simplify(ham); /* 3. Two-step read of the 4x4 dense matrix */ @@ -280,3 +285,86 @@ int main(void) { ``` Hamiltonians serve as observables in the error-mitigation workflows (see [Zero-Noise Extrapolation](../6_error_mitigation/1_zne.md) and [Virtual Distillation](../6_error_mitigation/2_virtual_distillation.md)); their expectation values can be computed directly with the statevector simulator (see [Statevector](1_statevector.md)). + + +## Checked sparse Pauli algebra + +`CHamiltonian` also represents general non-Hermitian sparse operators. These new functions follow the error/output conventions of the [fermionic binding](10_fermion.md): success is 0, object outputs are independent handles delivered through `**out`, and errors leave outputs and inputs unchanged. Input Pauli handles are borrowed and copied. Legacy consuming `from_pauli`, in-place `scale` and approximate-pruning `simplify` remain unchanged. + +```c +typedef struct CHamiltonianTerm { + const CPauliString *pauli; + Complex64 coefficient; +} CHamiltonianTerm; + +int32_t hamiltonian_identity(size_t num_qubits, CHamiltonian **out); + +int32_t hamiltonian_from_terms(size_t num_qubits, + const CHamiltonianTerm *terms, + size_t len, + CHamiltonian **out); + +int32_t hamiltonian_copy(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_normalized(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_normalized_hermitian(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_add(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_sub(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_mul(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_scaled(const CHamiltonian *ptr, + Complex64 factor, + CHamiltonian **out); + +int32_t hamiltonian_adjoint(const CHamiltonian *ptr, CHamiltonian **out); + +int32_t hamiltonian_tensor(const CHamiltonian *left, + const CHamiltonian *right, + CHamiltonian **out); + +int32_t hamiltonian_remap_qubits(const CHamiltonian *ptr, + const size_t *mapping, + size_t len, + size_t target_width, + CHamiltonian **out); + +int32_t hamiltonian_chop(CHamiltonian *ptr, double atol); + +int32_t hamiltonian_equiv(const CHamiltonian *left, + const CHamiltonian *right, + double atol, + uint8_t *out); + +int32_t hamiltonian_is_hermitian(const CHamiltonian *ptr, double atol, uint8_t *out); + +int32_t hamiltonian_get_term(const CHamiltonian *ptr, + size_t index, + CPauliString **out_pauli, + Complex64 *out_coefficient); +``` + +| Functions | Behavior | +| --- | --- | +| `identity/from_terms` | Identity supports scalar one on zero qubits; empty terms retain width. Batch construction validates widths/coefficients, absorbs Pauli phases, merges duplicates and removes only exact zeros. | +| `copy/get_term` | Independent storage copies. get_term returns a stored Pauli copy and coefficient without normalizing; free the Pauli with pauli_string_free. Zero-width Paulis are valid. Bad indices return -2. | +| `normalized` | Validated, exactly merged copy without near-zero pruning. | +| `normalized_hermitian` | After normalization, projects imaginary parts of magnitude <=1e-10 to zero; larger parts return -8. Retains every nonzero real coefficient. | +| `add/sub/mul/scaled/adjoint` | Call core checked methods; validate input/result coefficients. Nonfinite values and overflow return -8. Binary arithmetic requires equal widths; mul applies the right operand first. | +| `tensor` | Left tensor right; right occupies low qubits. Widths add; arithmetic width overflow returns -8. | +| `remap_qubits` | mapping[source]=destination; length equals input width; destinations are unique and less than target_width. Unused qubits are identity. This is not a fermionic permutation. | +| `chop` | Normalize in place then remove coefficients with magnitude <=atol; failure is atomic. | +| `equiv/is_hermitian` | Write 0/1. atol must be finite/nonnegative; 0 is exact. equiv requires equal widths; Hermiticity compares with the adjoint without projection. | + +Free object results with `hamiltonian_free`; supply an empty, valid output slot. Required NULL pointers return -1; invalid parameters, widths, coefficients or tolerances return -8. Unrepresentable materialized Pauli widths in identity/tensor/remap return -7. Zero sums may retain larger widths without materializing bit vectors. + +atol is a coefficient tolerance. On normalized sums, `is_hermitian` bounds `2*abs(imag) <= atol`; this is independent of the projection threshold in `normalized_hermitian`. Sparse algebra builds no dense matrix and never implicitly calls legacy simplify. diff --git a/tools/c-sdk/CMakeLists.txt b/tools/c-sdk/CMakeLists.txt index befa07af..4fc74dd4 100644 --- a/tools/c-sdk/CMakeLists.txt +++ b/tools/c-sdk/CMakeLists.txt @@ -28,30 +28,33 @@ else() endif() foreach(kind static shared) - add_executable(test_${kind} test_circuit.c) - target_include_directories(test_${kind} PRIVATE "${CQLIB_SDK_ROOT}/include") - target_link_libraries(test_${kind} PRIVATE cqlib_${kind}) - # The C tests put API calls inside assert(); retain them even in Release. - if(MSVC) - target_compile_options(test_${kind} PRIVATE /UNDEBUG) - else() - target_compile_options(test_${kind} PRIVATE -UNDEBUG) - endif() - if(WIN32) - target_link_libraries(test_${kind} PRIVATE - ntdll ws2_32 bcrypt userenv advapi32) - elseif(APPLE) - target_link_libraries(test_${kind} PRIVATE c++ "-framework CoreFoundation") - else() - target_link_libraries(test_${kind} PRIVATE pthread dl m) - endif() - set_target_properties(test_${kind} PROPERTIES - BUILD_RPATH "${CQLIB_SDK_ROOT}/lib") - add_test(NAME c_${kind} COMMAND test_${kind}) + foreach(suite circuit qis_algebra) + set(test_target test_${suite}_${kind}) + add_executable(${test_target} test_${suite}.c) + set_target_properties(${test_target} PROPERTIES C_STANDARD 11 C_STANDARD_REQUIRED YES) + target_include_directories(${test_target} PRIVATE "${CQLIB_SDK_ROOT}/include") + target_link_libraries(${test_target} PRIVATE cqlib_${kind}) + # The C tests put API calls inside assert(); retain them even in Release. + if(MSVC) + target_compile_options(${test_target} PRIVATE /UNDEBUG) + else() + target_compile_options(${test_target} PRIVATE -UNDEBUG) + endif() + if(WIN32) + target_link_libraries(${test_target} PRIVATE + ntdll ws2_32 bcrypt userenv advapi32) + elseif(APPLE) + target_link_libraries(${test_target} PRIVATE c++ "-framework CoreFoundation") + else() + target_link_libraries(${test_target} PRIVATE pthread dl m) + endif() + set_target_properties(${test_target} PROPERTIES + BUILD_RPATH "${CQLIB_SDK_ROOT}/lib") + add_test(NAME c_${suite}_${kind} COMMAND ${test_target}) + if(WIN32 AND kind STREQUAL "shared") + add_custom_command(TARGET ${test_target} POST_BUILD + COMMAND ${CMAKE_COMMAND} -E copy_if_different + "$" "$") + endif() + endforeach() endforeach() - -if(WIN32) - add_custom_command(TARGET test_shared POST_BUILD - COMMAND ${CMAKE_COMMAND} -E copy_if_different - "$" "$") -endif() diff --git a/tools/c-sdk/README.md b/tools/c-sdk/README.md index 0ea4b178..ae465515 100644 --- a/tools/c-sdk/README.md +++ b/tools/c-sdk/README.md @@ -54,5 +54,10 @@ RPATH such as `@executable_path/../lib`; the dylib uses `@rpath/libbinding_c.dylib` as its install name. No checkout or Cargo build directory is needed by SDK consumers. -The supported circuit construction and parameter APIs are documented in the -generated header and illustrated in `tests/test_circuit.c`. +Circuit construction, parameters, fermionic algebra, Jordan-Wigner mapping and +checked sparse Pauli algebra are documented in the generated header and +illustrated in `tests/test_circuit.c` and `tests/test_qis_algebra.c`. Both suites +run against static and shared libraries. New checked algebra APIs return error +codes and deliver owned results through output pointers; failed calls leave +outputs unchanged. Initialize object output slots to NULL and free results with +their corresponding functions. diff --git a/tools/ci.py b/tools/ci.py index 0eb347ce..ccb5c4c6 100644 --- a/tools/ci.py +++ b/tools/ci.py @@ -105,7 +105,8 @@ def build_and_test(platform=None): shutil.copy2(ROOT / "LICENSE.txt", sdk) shutil.copy2(ROOT / "tools/c-sdk/README.md", sdk / "README.md") (sdk / "tests").mkdir() - shutil.copy2(ROOT / "crates/binding-c/tests/test_circuit.c", sdk / "tests") + for test in ("test_circuit.c", "test_qis_algebra.c"): + shutil.copy2(ROOT / "crates/binding-c/tests" / test, sdk / "tests") shutil.copy2(ROOT / "tools/c-sdk/CMakeLists.txt", sdk / "tests") if sys.platform == "darwin": dylib = sdk / "lib/libbinding_c.dylib"