fix: correct circuit and simulation semantics for beta6

Cargo.toml

@@ -8,7 +8,7 @@ members = [
 resolver = "2"
 
 [workspace.package]
-version = "0.1.0-beta.5"
+version = "0.1.0-beta.6"
 edition = "2024"
 rust-version = "1.89"
 license = "Apache-2.0"

crates/binding-c/Cargo.toml

@@ -9,7 +9,7 @@ license.workspace = true
 crate-type = ["rlib", "staticlib", "cdylib"]
 
 [dependencies]
-cqlib-core = { path = "../cqlib-core", version = "0.1.0-beta.5" }
+cqlib-core = { path = "../cqlib-core", version = "0.1.0-beta.6" }
 libc = "0.2"
 bitvec.workspace = true
 ndarray.workspace = true

crates/binding-c/src/circuit/parameter.rs

@@ -445,7 +445,8 @@ pub extern "C" fn param_replace(
     box_param(unsafe { (*a).inner.replace(symbol, replacement) })
 }
 
-/// Substitute several symbols in one call and simplify the result.
+/// Simultaneously substitute symbols and simplify the result. Replacement
+/// expressions are not substituted again; the original parameter is unchanged.
 /// `names` and `values` are parallel arrays of length `len`. Returns a new
 /// parameter, or NULL on NULL input, non-UTF-8 names or NULL elements.
 #[unsafe(no_mangle)]

crates/binding-c/src/qis/density_matrix_noise.rs

@@ -682,6 +682,10 @@ pub extern "C" fn density_matrix_noise_probabilities_with_readout_len(
 /// Copy the probability distribution including readout errors configured in
 /// the noise model for `qubits` into `buffer`. Returns 0 on success.
 ///
+/// `qubits` contains storage positions. For a circuit-bound state, noise is
+/// looked up using the circuit ID at each position; unbound states use identity
+/// IDs. The full output distribution retains storage order.
+///
 /// `qubits` may be NULL when `qubit_len` is 0, in which case no readout
 /// error is applied and the result matches `density_matrix_noise_probabilities`.
 #[unsafe(no_mangle)]

crates/binding-c/tests/parameter_test.rs

@@ -22,6 +22,63 @@ fn eval(ptr: *const CParameter, bindings: &str) -> f64 {
     param_evaluate(ptr, cstr(bindings).as_ptr())
 }
 
+#[test]
+fn param_substitution_swaps_symbols_simultaneously() {
+    let x = param_symbol(cstr("x").as_ptr());
+    let y = param_symbol(cstr("y").as_ptr());
+    let expression = param_parse(cstr("x + 2*y").as_ptr());
+    let name_x = cstr("x");
+    let name_y = cstr("y");
+    for (names, values) in [
+        (
+            [name_x.as_ptr(), name_y.as_ptr()],
+            [y.cast_const(), x.cast_const()],
+        ),
+        (
+            [name_y.as_ptr(), name_x.as_ptr()],
+            [x.cast_const(), y.cast_const()],
+        ),
+    ] {
+        for _ in 0..100 {
+            let result = param_substitute_many(expression, names.as_ptr(), values.as_ptr(), 2);
+            assert!(!result.is_null());
+            let value = eval(result, "x:1,y:3");
+            param_free(result);
+            assert_eq!(value, 5.0);
+        }
+    }
+    assert_eq!(eval(expression, "x:1,y:3"), 7.0);
+    for ptr in [x, y, expression] {
+        param_free(ptr);
+    }
+}
+
+#[test]
+fn param_substitution_reaches_polynomial_symbols() {
+    let x = param_symbol(cstr("x").as_ptr());
+    let y = param_symbol(cstr("y").as_ptr());
+    let square = param_mul(x, x);
+    let sum = param_add(square, x);
+    let one = param_from_double(1.0);
+    let expression = param_add(sum, one);
+    let name = cstr("x");
+    let replaced = param_replace(expression, name.as_ptr(), y);
+    let substituted = param_substitute_many(
+        expression,
+        [name.as_ptr()].as_ptr(),
+        [y.cast_const()].as_ptr(),
+        1,
+    );
+    for result in [replaced, substituted] {
+        assert!(!result.is_null());
+        assert_eq!(eval(result, "y:2"), 7.0);
+        param_free(result);
+    }
+    for ptr in [x, y, square, sum, one, expression] {
+        param_free(ptr);
+    }
+}
+
 #[test]
 fn param_construction_and_constants() {
     let theta = param_symbol(cstr("theta").as_ptr());

crates/binding-c/tests/qis_noise_test.rs

@@ -2,7 +2,9 @@
 
 #![allow(dead_code)]
 
-use binding_c::circuit::{circuit_cx, circuit_free, circuit_h, circuit_new, circuit_x};
+use binding_c::circuit::{
+    circuit_cx, circuit_free, circuit_from_qubits, circuit_h, circuit_new, circuit_x,
+};
 use binding_c::cqlib_string_free;
 use binding_c::device::{
     NOISE_BIT_FLIP, noise_model_add_readout, noise_model_add_single_qubit, noise_model_free,
@@ -216,6 +218,42 @@ fn test_readout_error_difference() {
     noise_model_free(nm);
 }
 
+#[test]
+fn test_readout_uses_circuit_ids_at_storage_positions() {
+    for ids in [[0, 1], [1, 0], [2, 0], [9, 2]] {
+        let circuit = circuit_from_qubits(ids.as_ptr(), ids.len());
+        assert!(!circuit.is_null());
+        let model = noise_model_new();
+        assert_eq!(noise_model_add_readout(model, ids[0], 0.0, 1.0), 0);
+        let sim = density_matrix_noise_from_circuit(circuit, model);
+        assert!(!sim.is_null());
+        assert_eq!(density_matrix_noise_probabilities_with_readout_len(sim), 4);
+        for (positions, expected) in [
+            (&[0][..], [0.0, 1.0, 0.0, 0.0]),
+            (&[1][..], [1.0, 0.0, 0.0, 0.0]),
+            (&[1, 0][..], [0.0, 1.0, 0.0, 0.0]),
+            (&[][..], [1.0, 0.0, 0.0, 0.0]),
+        ] {
+            let mut actual = [0.0; 4];
+            assert_eq!(
+                density_matrix_noise_probabilities_with_readout(
+                    sim,
+                    actual.as_mut_ptr(),
+                    actual.len(),
+                    positions.as_ptr(),
+                    positions.len()
+                ),
+                0
+            );
+            assert_eq!(actual, expected);
+        }
+        assert_eq!(read_probs(sim, 4), vec![1.0, 0.0, 0.0, 0.0]);
+        density_matrix_noise_free(sim);
+        noise_model_free(model);
+        circuit_free(circuit);
+    }
+}
+
 #[test]
 fn test_multi_qubit_readout() {
     // Mirrors core `test_multi_qubit_readout_error`.

crates/binding-c/tests/qis_test.rs

@@ -165,6 +165,61 @@ fn test_density_matrix_partial_trace() {
     density_matrix_free(dm);
 }
 
+#[test]
+fn density_matrix_duplicate_gate_targets_preserve_state() {
+    let dm = density_matrix_new(2);
+    assert!(!dm.is_null());
+    assert_eq!(density_matrix_apply_h(dm, 0), 0);
+    let mut before = vec![num_complex::Complex64::default(); 16];
+    assert_eq!(density_matrix_data(dm, before.as_mut_ptr(), 16), 0);
+    for q in [0, 1] {
+        assert_eq!(density_matrix_apply_cx(dm, q, q), -8);
+        assert_eq!(density_matrix_apply_cz(dm, q, q), -8);
+        let mut after = before.clone();
+        assert_eq!(density_matrix_data(dm, after.as_mut_ptr(), 16), 0);
+        assert_eq!(after, before);
+    }
+    assert_eq!(density_matrix_apply_cx(dm, 0, 1), 0);
+    density_matrix_free(dm);
+}
+
+#[test]
+fn density_matrix_duplicate_kraus_targets_preserve_state() {
+    let state: Vec<_> = (1..=4)
+        .map(|i| num_complex::Complex64::new(i as f64 / 30.0_f64.sqrt(), 0.0))
+        .collect();
+    let dm = density_matrix_from_state(2, state.as_ptr(), state.len());
+    assert!(!dm.is_null());
+    let vector = [4.0, -1.0, 2.0, 3.0];
+    let operator: Vec<_> = (0..4)
+        .flat_map(|r| {
+            (0..4).map(move |c| {
+                num_complex::Complex64::new(
+                    if r == c { 1.0 } else { 0.0 } - 2.0 * vector[r] * vector[c] / 30.0,
+                    0.0,
+                )
+            })
+        })
+        .collect();
+    let mut before = vec![num_complex::Complex64::default(); 16];
+    assert_eq!(density_matrix_data(dm, before.as_mut_ptr(), 16), 0);
+    for q in [0, 1] {
+        let targets = [q, q];
+        assert_eq!(
+            density_matrix_apply_kraus(dm, operator.as_ptr(), 1, 4, targets.as_ptr(), 2),
+            -8
+        );
+        let mut after = before.clone();
+        assert_eq!(density_matrix_data(dm, after.as_mut_ptr(), 16), 0);
+        assert_eq!(after, before);
+    }
+    assert_eq!(
+        density_matrix_apply_kraus(dm, operator.as_ptr(), 1, 4, [0, 1].as_ptr(), 2),
+        0
+    );
+    density_matrix_free(dm);
+}
+
 #[test]
 fn test_stabilizer_new_and_clifford() {
     let ss = stabilizer_new(2);

crates/binding-c/tests/test_circuit.c

@@ -165,11 +165,98 @@ static void test_extended_gates(void) {
     circuit_free(circuit);
 }
 
+static void test_circuit_gate_operation_matrix(void) {
+    const uint32_t orders[][2] = {{0, 1}, {1, 0}, {7, 2}};
+    const uint32_t targets[] = {0, 1};
+    for (size_t order = 0; order < 3; ++order) {
+        for (int parameterized = 0; parameterized < 2; ++parameterized) {
+            CCircuit* body = circuit_new(0);
+            assert(body != NULL);
+            assert(circuit_add_qubits(body, orders[order], 2) == 0);
+            /* Build on dense ids, then map into the explicit storage order. */
+            CCircuit* definition = circuit_new(2);
+            assert(definition != NULL);
+            CParameter* theta = param_parse("theta");
+            assert(theta != NULL);
+            if (parameterized) {
+                assert(circuit_rx_param(definition, 0, theta) == 0);
+            } else {
+                assert(circuit_x(definition, 0) == 0);
+            }
+            assert(circuit_compose(body, definition, orders[order], 2) == 0);
+            circuit_free(definition);
+            assert(circuit_set_global_phase(body, 0.13) == 0);
+            CCircuitGate* gate = circuit_to_gate(body, "first_gate");
+            assert(gate != NULL);
+            for (size_t application = 0; application < 2; ++application) {
+                double angle = application == 0 ? 0.37 : -0.61;
+                CParameter* argument = param_from_double(angle);
+                const CParameter* arguments[] = {argument};
+                assert(argument != NULL);
+                CCircuit* circuit = circuit_new(2);
+                assert(circuit != NULL);
+                assert(circuit_circuit_gate(circuit, gate, targets, 2,
+                                            parameterized ? arguments : NULL,
+                                            parameterized ? 1 : 0) == 0);
+                /* CParameter arguments are stored as parameter-table indices.
+                 * Resolve them in the owning circuit before extracting the raw
+                 * DAG Operation, which has no access to that parameter table.
+                 * The gate definition still contains theta; matrix() must bind it. */
+                CCircuit* resolved = circuit_assign_params(circuit, "");
+                assert(resolved != NULL && circuit_num_parameters(resolved) == 0);
+                circuit_free(circuit);
+                circuit = resolved;
+                CCircuitDag* dag = circuit_dag_from_circuit(circuit);
+                assert(dag != NULL && circuit_dag_op_nodes_len(dag) == 1);
+                uint32_t node;
+                assert(circuit_dag_op_nodes(dag, &node, 1) == 0);
+                COperation* operation = circuit_dag_operation(dag, node);
+                assert(operation != NULL);
+                size_t rows = 0, cols = 0;
+                assert(operation_matrix_len(operation) == 16);
+                assert(operation_matrix_dims(operation, &rows, &cols) == 0);
+                assert(rows == 4 && cols == 4);
+                double local[32], full[32];
+                assert(operation_matrix(operation, local, 32) == 16);
+                assert(circuit_to_matrix(circuit, NULL, 0, full, 32) == 16);
+                for (size_t row = 0; row < 4; ++row) {
+                    for (size_t col = 0; col < 4; ++col) {
+                        for (int local_basis = 0; local_basis < 2; ++local_basis) {
+                            size_t flip = local_basis ? 2 : 1;
+                            double re = 0, im = 0;
+                            if (parameterized) {
+                                if (row == col) re = cos(angle / 2);
+                                if (row == (col ^ flip)) im = -sin(angle / 2);
+                            } else if (row == (col ^ flip)) {
+                                re = 1;
+                            }
+                            double expected_re = re * cos(0.13) - im * sin(0.13);
+                            double expected_im = re * sin(0.13) + im * cos(0.13);
+                            const double* matrix = local_basis ? local : full;
+                            size_t index = 2 * (row * 4 + col);
+                            assert(fabs(matrix[index] - expected_re) < 1e-12);
+                            assert(fabs(matrix[index + 1] - expected_im) < 1e-12);
+                        }
+                    }
+                }
+                operation_free(operation);
+                circuit_dag_free(dag);
+                circuit_free(circuit);
+                param_free(argument);
+            }
+            circuit_gate_free(gate);
+            param_free(theta);
+            circuit_free(body);
+        }
+    }
+}
+
 int main(void) {
     test_basic_circuit();
     test_errors();
     test_symbolic_parameters();
     test_extended_gates();
+    test_circuit_gate_operation_matrix();
     printf("binding-c circuit tests passed\n");
     return 0;
 }