From ec5d93b89254c96027b029e54400f42da5e0b218 Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Tue, 25 Aug 2026 02:31:52 +0000 Subject: [PATCH 1/6] Add the reversible integer arithmetic device kernels to primitives CDKM/Cuccaro ripple-carry family (register and constant add/subtract, >=-constant comparator) and the ancilla-free Draper QFT family (qft/iqft, Fourier-basis constant add/subtract, the comparator pair with its hand-written adjoint), all little-endian in-place kernels with hand-written inverses. Tests are exhaustive over all inputs for widths up to 5 via the superposition harness, pin every inverse by op-then-inverse identity, and hold the documented gate prices against the compiler with cudaq.estimate_resources: exactly 2n Toffolis for every CDKM operation (0 at K = 0), and zero Toffolis for the QFT family with exact controlled-r1 / r1 / h budgets (n(n-1) + n for the adder, n(n+1) + n + 1 per comparator side). Review fixes: document the 0 <= K <= 2^n precondition on cmp_ge_constant_qft (the borrow wraps mod 2^(n+1) for larger K) and the exact constant_bits/complement_bits length preconditions on add_constant and cmp_ge_constant; note in the module docstring that operands of any one op must be pairwise disjoint; version-scope the cudaq.adjoint limitation (broken as of CUDA-Q 0.15, still unresolved on the 0.16 pre-release) instead of stating it timelessly. Signed-off-by: Scott Thornton --- .../cudaq_algorithms/primitives/__init__.py | 29 +- .../primitives/_arithmetic.py | 333 +++++++++++++ tests/python/test_primitives_arithmetic.py | 437 ++++++++++++++++++ 3 files changed, 798 insertions(+), 1 deletion(-) create mode 100644 python/cudaq_algorithms/primitives/_arithmetic.py create mode 100644 tests/python/test_primitives_arithmetic.py diff --git a/python/cudaq_algorithms/primitives/__init__.py b/python/cudaq_algorithms/primitives/__init__.py index 11820a5..2c0d2c6 100644 --- a/python/cudaq_algorithms/primitives/__init__.py +++ b/python/cudaq_algorithms/primitives/__init__.py @@ -1,6 +1,6 @@ # SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. # SPDX-License-Identifier: Apache-2.0 -"""Fault-tolerant circuit primitives: unary iteration and QROM. +"""Fault-tolerant circuit primitives: unary iteration, QROM, arithmetic. ``unary_iteration_kernels`` mints unary iteration (Babbush et al., `arXiv:1805.03662`, Fig. 7) in its strictly unitary form — no @@ -19,15 +19,42 @@ minted kernels actually have, and the resource tests pin them against the compiler. +The reversible integer arithmetic device kernels (:mod:`._arithmetic`) +are the in-place little-endian adders and comparators the lookup-based +constructions compose with: the CDKM/Cuccaro ripple-carry family +(``add_register`` / ``subtract_register``, ``add_constant`` / +``subtract_constant``, ``cmp_ge_constant``) and the ancilla-free Draper +QFT family (``qft`` / ``iqft``, ``add_constant_qft`` / +``subtract_constant_qft``, the ``cmp_ge_constant_qft`` / +``cmp_ge_constant_qft_adj`` pair). Every inverse is hand-written and the +gate prices are compiler-pinned by the resource tests. + Import the subpackage directly (``from cudaq_algorithms.primitives import QROM``); nothing here is re-exported from the package root. """ +from ._arithmetic import (add_constant, add_constant_qft, add_register, + cmp_ge_constant, cmp_ge_constant_qft, + cmp_ge_constant_qft_adj, iqft, phase_add_constant, + qft, subtract_constant, subtract_constant_qft, + subtract_register) from ._qrom import QROM from ._unary_iteration import UnaryIterationKernels, unary_iteration_kernels __all__ = [ "QROM", "UnaryIterationKernels", + "add_constant", + "add_constant_qft", + "add_register", + "cmp_ge_constant", + "cmp_ge_constant_qft", + "cmp_ge_constant_qft_adj", + "iqft", + "phase_add_constant", + "qft", + "subtract_constant", + "subtract_constant_qft", + "subtract_register", "unary_iteration_kernels", ] diff --git a/python/cudaq_algorithms/primitives/_arithmetic.py b/python/cudaq_algorithms/primitives/_arithmetic.py new file mode 100644 index 0000000..4bb0297 --- /dev/null +++ b/python/cudaq_algorithms/primitives/_arithmetic.py @@ -0,0 +1,333 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +"""Reversible in-place integer arithmetic device kernels. + +Two families, both little-endian (``register[0]`` is the least-significant +bit, matching ``docs/conventions.md``): + +- **CDKM/Cuccaro ripple-carry** (`arXiv:quant-ph/0410184`): + ``add_register`` / ``subtract_register`` (in-place ``b <- b +/- a``), + ``add_constant`` / ``subtract_constant`` (constant loaded into a caller + provided work register), and ``cmp_ge_constant`` (a ``x >= K`` comparator + writing into an out qubit, leaving ``x`` untouched: MAJ sweep, copy the + carry, reverse MAJ sweep). Toffoli prices (documented contracts, pinned + by the resource tests in ``tests/python/test_primitives_arithmetic.py``): + ``add_register`` / ``subtract_register`` cost exactly ``2 n`` Toffolis + on ``n``-bit registers (``n`` in the MAJ sweep, ``n`` in the UMA sweep), + ``add_constant`` / ``subtract_constant`` inherit the same ``2 n`` (the + constant load is X-only), and ``cmp_ge_constant`` costs ``2 n`` for + ``K >= 1`` (MAJ sweep plus its reversal) and ``0`` for ``K = 0``. +- **Draper QFT arithmetic** (`arXiv:quant-ph/0008033`): ``qft`` / ``iqft`` + and ``add_constant_qft``, plus the ``cmp_ge_constant_qft`` / + ``cmp_ge_constant_qft_adj`` comparator pair. The QFT family needs *no* + work qubits — on a statevector simulator every extra work qubit doubles + the memory, and the CDKM constant ops need ``n + 1`` of them — and *no* + Toffolis: its price is rotations, ``n (n - 1)`` controlled-``r1`` plus + ``n`` free ``r1`` for ``add_constant_qft`` on ``n`` bits, and + ``n (n + 1)`` controlled-``r1`` plus ``n + 1`` free ``r1`` per side of + the ``cmp_ge_constant_qft`` pair (``K >= 1``; zero rotations at + ``K = 0``), also pinned by the resource tests. + +Registers passed to any one op must be pairwise disjoint: the kernels +never alias-check, and overlapping views (e.g. ``add_register(a, a, +carry)``, or a ``work`` view sharing qubits with ``target``) are +undefined. + +Every inverse is hand-written (``cudaq.adjoint`` is off-limits as of +CUDA-Q 0.15: cuda-quantum#4897/#4898, still unresolved as of the 0.16 +pre-release) and pinned by op-then-inverse identity tests in +``tests/python/test_primitives_arithmetic.py``, alongside exhaustive +value checks against classical integer arithmetic for widths up to 5. + +Kernel-language notes (see CLAUDE.md): guards are positive ``if`` blocks +(kernel ``return`` is silently ignored, cuda-quantum#4845), and all +classical indices are recomputed from loop variables rather than carried +as mutating accumulators. +""" + +from __future__ import annotations + +import cudaq + +__all__ = [ + "add_register", + "subtract_register", + "add_constant", + "subtract_constant", + "cmp_ge_constant", + "qft", + "iqft", + "phase_add_constant", + "add_constant_qft", + "subtract_constant_qft", + "cmp_ge_constant_qft", + "cmp_ge_constant_qft_adj", +] + +# ============================================================================ +# CDKM / Cuccaro ripple-carry family +# ============================================================================ +# +# MAJ(x, y, z) = cx(z, y); cx(z, x); ccx(x, y, z) +# UMA(x, y, z) = ccx(x, y, z); cx(z, x); cx(x, y) +# +# The adder chains MAJ(carry, b0, a0), MAJ(a0, b1, a1), ...; after the MAJ +# sweep the ripple carry-out sits on a[n-1]; the UMA sweep (in reverse) +# restores ``a`` and the carry ancilla and completes ``b <- a + b mod 2^n``. + + +@cudaq.kernel +def add_register(a: cudaq.qview, b: cudaq.qview, carry: cudaq.qview): + """``b <- (a + b) mod 2^n`` (CDKM). ``a`` and ``carry`` are restored. + + ``a`` and ``b`` have equal size ``n``; ``carry`` is a one-qubit view + that must be |0> on entry (it is returned to |0>). + """ + n = a.size() + if n > 0: + # MAJ sweep. + cx(a[0], b[0]) + cx(a[0], carry[0]) + x.ctrl(carry[0], b[0], a[0]) + for i in range(1, n): + cx(a[i], b[i]) + cx(a[i], a[i - 1]) + x.ctrl(a[i - 1], b[i], a[i]) + # UMA sweep (reverse). + for k in range(1, n): + i = n - k + x.ctrl(a[i - 1], b[i], a[i]) + cx(a[i], a[i - 1]) + cx(a[i - 1], b[i]) + x.ctrl(carry[0], b[0], a[0]) + cx(a[0], carry[0]) + cx(carry[0], b[0]) + + +@cudaq.kernel +def subtract_register(a: cudaq.qview, b: cudaq.qview, carry: cudaq.qview): + """``b <- (b - a) mod 2^n``: the exact gate-reversal of ``add_register``. + + Hand-written inverse (no ``cudaq.adjoint``); every gate of the adder is + self-inverse, so the reversed sequence is the inverse circuit. + """ + n = a.size() + if n > 0: + cx(carry[0], b[0]) + cx(a[0], carry[0]) + x.ctrl(carry[0], b[0], a[0]) + for i in range(1, n): + cx(a[i - 1], b[i]) + cx(a[i], a[i - 1]) + x.ctrl(a[i - 1], b[i], a[i]) + for k in range(1, n): + i = n - k + x.ctrl(a[i - 1], b[i], a[i]) + cx(a[i], a[i - 1]) + cx(a[i], b[i]) + x.ctrl(carry[0], b[0], a[0]) + cx(a[0], carry[0]) + cx(a[0], b[0]) + + +@cudaq.kernel +def add_constant(target: cudaq.qview, constant_bits: list[int], + work: cudaq.qview, carry: cudaq.qview): + """``target <- (target + K) mod 2^n`` (CDKM). + + ``constant_bits[k]`` is bit ``k`` of ``K`` (little-endian). Precondition: + ``constant_bits`` has length exactly ``n = target.size()``, i.e. + ``0 <= K < 2^n`` — a shorter or longer list is not truncated or padded + and the kernel is undefined. ``work`` (``n`` qubits) and ``carry`` (1 + qubit) must be |0> on entry and are returned to |0>: the constant is + X-loaded into ``work``, ripple-added, and X-unloaded. + """ + n = target.size() + for k in range(n): + if constant_bits[k] == 1: + x(work[k]) + add_register(work, target, carry) + for k in range(n): + if constant_bits[k] == 1: + x(work[k]) + + +@cudaq.kernel +def subtract_constant(target: cudaq.qview, constant_bits: list[int], + work: cudaq.qview, carry: cudaq.qview): + """``target <- (target - K) mod 2^n``: inverse of ``add_constant``.""" + n = target.size() + for k in range(n): + if constant_bits[k] == 1: + x(work[k]) + subtract_register(work, target, carry) + for k in range(n): + if constant_bits[k] == 1: + x(work[k]) + + +@cudaq.kernel +def cmp_ge_constant(x_reg: cudaq.qview, complement_bits: list[int], + k_is_zero: int, work: cudaq.qview, carry: cudaq.qview, + out: cudaq.qview): + """``out[0] ^= (x >= K)`` (CDKM), leaving ``x_reg`` unchanged. + + ``complement_bits`` are the little-endian bits of ``2^n - K``. + Precondition: ``complement_bits`` has length exactly ``n = + x_reg.size()`` and ``0 <= K <= 2^n`` (so ``2^n - K`` fits in ``n`` + bits; pass ``k_is_zero = 1`` and all-zero bits for ``K = 0``, which is + always true, and all-zero bits with ``k_is_zero = 0`` for ``K = 2^n``, + which is always false). Uses the ripple identity ``x >= K <=> carry_out(x + (2^n + - K))`` for ``K >= 1``: MAJ sweep, copy the carry-out (which sits on + the top work qubit) into ``out``, then reverse the MAJ sweep so + ``x_reg``, ``work`` and ``carry`` are all restored. + """ + if k_is_zero == 1: + x(out[0]) + else: + n = x_reg.size() + for k in range(n): + if complement_bits[k] == 1: + x(work[k]) + # MAJ sweep (a = work, b = x_reg). + cx(work[0], x_reg[0]) + cx(work[0], carry[0]) + x.ctrl(carry[0], x_reg[0], work[0]) + for i in range(1, n): + cx(work[i], x_reg[i]) + cx(work[i], work[i - 1]) + x.ctrl(work[i - 1], x_reg[i], work[i]) + cx(work[n - 1], out[0]) + # Reverse MAJ sweep (restores x_reg, work, carry). + for k in range(1, n): + i = n - k + x.ctrl(work[i - 1], x_reg[i], work[i]) + cx(work[i], work[i - 1]) + cx(work[i], x_reg[i]) + x.ctrl(carry[0], x_reg[0], work[0]) + cx(work[0], carry[0]) + cx(work[0], x_reg[0]) + for k in range(n): + if complement_bits[k] == 1: + x(work[k]) + + +# ============================================================================ +# Draper QFT family (no work qubits) +# ============================================================================ +# +# ``qft`` is the textbook circuit without the final bit-reversal swaps; +# after it, qubit t holds (|0> + exp(2 pi i x / 2^(t+1)) |1>)/sqrt(2), so a +# constant K is added by the single-qubit phases r1(2 pi K / 2^(t+1)). + + +@cudaq.kernel +def qft(reg: cudaq.qview): + """Quantum Fourier transform (no bit-reversal swaps; see module doc).""" + n = reg.size() + for j in range(n): + t = n - 1 - j + h(reg[t]) + for c in range(t): + r1.ctrl(3.141592653589793 / (1 << (t - c)), reg[c], reg[t]) + + +@cudaq.kernel +def iqft(reg: cudaq.qview): + """Inverse QFT: hand-written reversal of ``qft``.""" + n = reg.size() + for t in range(n): + for k in range(t): + c = t - 1 - k + r1.ctrl(-3.141592653589793 / (1 << (t - c)), reg[c], reg[t]) + h(reg[t]) + + +@cudaq.kernel +def phase_add_constant(reg: cudaq.qview, constant: int): + """``reg <- reg + K mod 2^n`` in the Fourier basis (between qft/iqft).""" + n = reg.size() + for t in range(n): + r1(6.283185307179586 * constant / (1 << (t + 1)), reg[t]) + + +@cudaq.kernel +def add_constant_qft(reg: cudaq.qview, constant: int): + """``reg <- (reg + K) mod 2^n`` (Draper; no work qubits).""" + qft(reg) + phase_add_constant(reg, constant) + iqft(reg) + + +@cudaq.kernel +def subtract_constant_qft(reg: cudaq.qview, constant: int): + """``reg <- (reg - K) mod 2^n``: inverse of ``add_constant_qft``.""" + qft(reg) + phase_add_constant(reg, -constant) + iqft(reg) + + +@cudaq.kernel +def _qft_extended(x_reg: cudaq.qview, msb: cudaq.qview): + """QFT over the (n+1)-bit register ``[x_reg, msb]`` (msb = bit n).""" + n = x_reg.size() + h(msb[0]) + for c in range(n): + r1.ctrl(3.141592653589793 / (1 << (n - c)), x_reg[c], msb[0]) + qft(x_reg) + + +@cudaq.kernel +def _iqft_extended(x_reg: cudaq.qview, msb: cudaq.qview): + """Inverse of ``_qft_extended``.""" + n = x_reg.size() + iqft(x_reg) + for k in range(n): + c = n - 1 - k + r1.ctrl(-3.141592653589793 / (1 << (n - c)), x_reg[c], msb[0]) + h(msb[0]) + + +@cudaq.kernel +def cmp_ge_constant_qft(x_reg: cudaq.qview, out: cudaq.qview, constant: int, + invert: int): + """Compute ``out[0] = (x >= K)`` (or ``x < K`` with ``invert = 1``). + + Precondition: ``0 <= K <= 2^n`` (``n = x_reg.size()``). For ``K > + 2^n`` the result is silently wrong: the subtraction acts on the + (n+1)-bit extension, so the borrow wraps mod ``2^(n+1)`` and ``out`` + no longer encodes ``x < K``. + + Draper-style: subtract ``K`` on the (n+1)-bit register ``[x_reg, + out]``; the MSB (``out``) becomes the borrow, i.e. ``x < K``. ``out`` + must be |0> on entry. Between this kernel and + ``cmp_ge_constant_qft_adj`` the low bits hold ``(x - K) mod 2^n`` — + callers may read ``out`` but must not consume ``x_reg`` until the + adjoint restores it. ``K = 0`` (with ``invert = 0``) yields the + constant-true comparator. + """ + n = x_reg.size() + if constant > 0: + _qft_extended(x_reg, out) + for t in range(n): + r1(-6.283185307179586 * constant / (1 << (t + 1)), x_reg[t]) + r1(-6.283185307179586 * constant / (1 << (n + 1)), out[0]) + _iqft_extended(x_reg, out) + if invert == 0: + x(out[0]) + + +@cudaq.kernel +def cmp_ge_constant_qft_adj(x_reg: cudaq.qview, out: cudaq.qview, + constant: int, invert: int): + """Hand-written inverse of ``cmp_ge_constant_qft``.""" + n = x_reg.size() + if invert == 0: + x(out[0]) + if constant > 0: + _qft_extended(x_reg, out) + r1(6.283185307179586 * constant / (1 << (n + 1)), out[0]) + for k in range(n): + t = n - 1 - k + r1(6.283185307179586 * constant / (1 << (t + 1)), x_reg[t]) + _iqft_extended(x_reg, out) diff --git a/tests/python/test_primitives_arithmetic.py b/tests/python/test_primitives_arithmetic.py new file mode 100644 index 0000000..ca60da6 --- /dev/null +++ b/tests/python/test_primitives_arithmetic.py @@ -0,0 +1,437 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +"""Exhaustive tests for the primitives arithmetic device kernels. + +Every operation is checked against classical integer arithmetic over *all* +inputs for register widths up to 5, and every inverse is pinned by an +op-then-inverse == identity test. Exhaustiveness comes from a superposition +harness: the input registers are prepared in a uniform superposition, so a +single statevector comparison validates the operation's action on every +computational-basis input at once (the operations are classical +permutations, so distinct inputs cannot interfere). + +The resource-contract tests at the bottom hold the module's documented +gate prices against the compiler: ``cudaq.estimate_resources`` counts the +operations actually synthesized, so the closed forms (``2 n`` Toffolis for +every CDKM operation, zero Toffolis and the exact ``r1`` budgets for the +Draper QFT family) are compiler facts, not emitter claims. +""" + +import numpy as np +import pytest + +import cudaq + +from cudaq_algorithms.primitives import _arithmetic as arith + +WIDTHS = [1, 2, 3, 4, 5] + + +def _basis(index: int, num_qubits: int) -> np.ndarray: + ket = np.zeros(1 << num_qubits, dtype=np.complex128) + ket[index] = 1.0 + return ket + + +# ---------------------------------------------------------------------- +# Register-register add / subtract (CDKM) +# ---------------------------------------------------------------------- +# +# Layout: a at bits [0, n), b at [n, 2n), carry at bit 2n. + + +@cudaq.kernel +def _run_add_register(n: int, subtract: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + for k in range(n): + h(a[k]) + h(b[k]) + if subtract == 0: + arith.add_register(a, b, carry) + else: + arith.subtract_register(a, b, carry) + + +@cudaq.kernel +def _run_add_then_subtract(n: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + for k in range(n): + h(a[k]) + h(b[k]) + arith.add_register(a, b, carry) + arith.subtract_register(a, b, carry) + + +def _register_op_expected(n: int, op) -> np.ndarray: + """Uniform superposition over (a, b) mapped through (a, op(a, b)).""" + expected = np.zeros(1 << (2 * n + 1), dtype=np.complex128) + norm = 1.0 / (1 << n) + for a in range(1 << n): + for b in range(1 << n): + expected[a + (op(a, b) % (1 << n)) * (1 << n)] += norm + return expected + + +@pytest.mark.parametrize("n", WIDTHS) +def test_add_register_all_inputs(n): + state = np.array(cudaq.get_state(_run_add_register, n, 0)) + np.testing.assert_allclose(state, + _register_op_expected(n, lambda a, b: a + b), + atol=1e-12) + + +@pytest.mark.parametrize("n", WIDTHS) +def test_subtract_register_all_inputs(n): + state = np.array(cudaq.get_state(_run_add_register, n, 1)) + np.testing.assert_allclose(state, + _register_op_expected(n, lambda a, b: b - a), + atol=1e-12) + + +@pytest.mark.parametrize("n", WIDTHS) +def test_subtract_register_inverts_add_register(n): + state = np.array(cudaq.get_state(_run_add_then_subtract, n)) + np.testing.assert_allclose(state, + _register_op_expected(n, lambda a, b: b), + atol=1e-12) + + +# ---------------------------------------------------------------------- +# Constant add / subtract (CDKM and Draper QFT) +# ---------------------------------------------------------------------- +# +# CDKM layout: target at [0, n), work at [n, 2n), carry at 2n. +# QFT layout: target only. + + +@cudaq.kernel +def _run_add_constant(n: int, bits: list[int], subtract: int, roundtrip: int): + target = cudaq.qvector(n) + work = cudaq.qvector(n) + carry = cudaq.qvector(1) + for k in range(n): + h(target[k]) + if subtract == 0: + arith.add_constant(target, bits, work, carry) + else: + arith.subtract_constant(target, bits, work, carry) + if roundtrip == 1: + if subtract == 0: + arith.subtract_constant(target, bits, work, carry) + else: + arith.add_constant(target, bits, work, carry) + + +@cudaq.kernel +def _run_add_constant_qft(n: int, constant: int, subtract: int, + roundtrip: int): + target = cudaq.qvector(n) + for k in range(n): + h(target[k]) + if subtract == 0: + arith.add_constant_qft(target, constant) + else: + arith.subtract_constant_qft(target, constant) + if roundtrip == 1: + if subtract == 0: + arith.subtract_constant_qft(target, constant) + else: + arith.add_constant_qft(target, constant) + + +def _constant_op_expected(n: int, shift: int, extra_qubits: int) -> np.ndarray: + expected = np.zeros(1 << (n + extra_qubits), dtype=np.complex128) + norm = 1.0 / np.sqrt(1 << n) + for value in range(1 << n): + expected[(value + shift) % (1 << n)] += norm + return expected + + +def _bits(value: int, n: int) -> list[int]: + return [(value >> k) & 1 for k in range(n)] + + +@pytest.mark.parametrize("n", WIDTHS) +def test_add_and_subtract_constant_all_inputs(n): + for constant in range(1 << n): + bits = _bits(constant, n) + added = np.array(cudaq.get_state(_run_add_constant, n, bits, 0, 0)) + np.testing.assert_allclose(added, + _constant_op_expected(n, constant, n + 1), + atol=1e-12) + subtracted = np.array(cudaq.get_state(_run_add_constant, n, bits, 1, + 0)) + np.testing.assert_allclose(subtracted, + _constant_op_expected(n, -constant, n + 1), + atol=1e-12) + roundtrip = np.array(cudaq.get_state(_run_add_constant, n, bits, 0, 1)) + np.testing.assert_allclose(roundtrip, + _constant_op_expected(n, 0, n + 1), + atol=1e-12) + + +@pytest.mark.parametrize("n", WIDTHS) +def test_add_and_subtract_constant_qft_all_inputs(n): + for constant in range(1 << n): + added = np.array( + cudaq.get_state(_run_add_constant_qft, n, constant, 0, 0)) + np.testing.assert_allclose(added, + _constant_op_expected(n, constant, 0), + atol=1e-10) + subtracted = np.array( + cudaq.get_state(_run_add_constant_qft, n, constant, 1, 0)) + np.testing.assert_allclose(subtracted, + _constant_op_expected(n, -constant, 0), + atol=1e-10) + roundtrip = np.array( + cudaq.get_state(_run_add_constant_qft, n, constant, 0, 1)) + np.testing.assert_allclose(roundtrip, + _constant_op_expected(n, 0, 0), + atol=1e-10) + + +# ---------------------------------------------------------------------- +# >= comparators (CDKM and Draper QFT) +# ---------------------------------------------------------------------- + + +@cudaq.kernel +def _run_cmp_ge_constant(n: int, bits: list[int], k_is_zero: int): + x_reg = cudaq.qvector(n) + work = cudaq.qvector(n) + carry = cudaq.qvector(1) + out = cudaq.qvector(1) + for k in range(n): + h(x_reg[k]) + arith.cmp_ge_constant(x_reg, bits, k_is_zero, work, carry, out) + + +@cudaq.kernel +def _run_cmp_ge_constant_qft(n: int, constant: int, invert: int, + uncompute: int): + x_reg = cudaq.qvector(n) + out = cudaq.qvector(1) + for k in range(n): + h(x_reg[k]) + arith.cmp_ge_constant_qft(x_reg, out, constant, invert) + if uncompute == 1: + arith.cmp_ge_constant_qft_adj(x_reg, out, constant, invert) + + +def _cmp_expected(n: int, predicate, extra_before_out: int) -> np.ndarray: + """|x> (work/carry |0>) |out = predicate(x)> over superposed x.""" + total = n + extra_before_out + 1 + expected = np.zeros(1 << total, dtype=np.complex128) + norm = 1.0 / np.sqrt(1 << n) + for value in range(1 << n): + expected[value + + (1 << (n + extra_before_out)) * int(predicate(value))] += norm + return expected + + +@pytest.mark.parametrize("n", WIDTHS) +def test_cmp_ge_constant_all_inputs(n): + for constant in range(1 << n): + complement = _bits((1 << n) - constant, n) if constant else [0] * n + state = np.array( + cudaq.get_state(_run_cmp_ge_constant, n, complement, + int(constant == 0))) + np.testing.assert_allclose(state, + _cmp_expected(n, lambda v: v >= constant, + n + 1), + atol=1e-12) + + +@pytest.mark.parametrize("n", WIDTHS) +def test_cmp_ge_constant_qft_all_inputs(n): + # Between the compute/adjoint pair the register is shifted by -K (a + # documented contract), so the compute-only check reads out through the + # shifted basis; the roundtrip check pins full restoration. + for constant in range(1 << n): + for invert in (0, 1): + state = np.array( + cudaq.get_state(_run_cmp_ge_constant_qft, n, constant, invert, + 0)) + expected = np.zeros(1 << (n + 1), dtype=np.complex128) + norm = 1.0 / np.sqrt(1 << n) + for value in range(1 << n): + flag = (value >= constant) if invert == 0 else (value + < constant) + shifted = (value - constant) % (1 << n) + expected[shifted + (int(flag) << n)] += norm + np.testing.assert_allclose(state, expected, atol=1e-10) + roundtrip = np.array( + cudaq.get_state(_run_cmp_ge_constant_qft, n, constant, invert, + 1)) + np.testing.assert_allclose(roundtrip, + _constant_op_expected(n, 0, 1), + atol=1e-10) + + +# ---------------------------------------------------------------------- +# Direct basis-state spot checks (readable, non-superposed) +# ---------------------------------------------------------------------- + + +@cudaq.kernel +def _spot_add(n: int, aval: int, bval: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + for k in range(n): + if ((aval >> k) & 1) == 1: + x(a[k]) + if ((bval >> k) & 1) == 1: + x(b[k]) + arith.add_register(a, b, carry) + + +def test_add_register_basis_spot_checks(): + for n, aval, bval in [(2, 1, 2), (3, 5, 6), (5, 21, 27)]: + state = np.array(cudaq.get_state(_spot_add, n, aval, bval)) + index = aval + (((aval + bval) % (1 << n)) << n) + np.testing.assert_allclose(state, _basis(index, 2 * n + 1), atol=1e-12) + + +# ---------------------------------------------------------------------- +# Compiler-pinned resource contracts +# ---------------------------------------------------------------------- +# +# The harnesses below pass the operation parameters as *runtime kernel +# arguments* (never source literals): constant-folded literals let the +# compiler specialize the circuit, and the pinned counts would then +# depend on the folding rather than the construction. +# +# CDKM derivations (n-bit registers, one Toffoli per MAJ and per UMA): +# - add_register / subtract_register: the MAJ sweep is one Toffoli at +# the carry step plus n - 1 in the ripple = n; the UMA sweep mirrors +# it = n. Total exactly 2 n. +# - add_constant / subtract_constant: the constant load/unload is X-only +# (free), so the price is the inner register adder's 2 n. +# - cmp_ge_constant (K >= 1): a MAJ sweep (n) plus its literal reversal +# (n) with a free CNOT carry-copy in between = 2 n; K = 0 is a bare X. +# +# Draper QFT derivations: no Toffolis at all. add_constant_qft is +# qft + phases + iqft = 2 * (n(n-1)/2) controlled-r1, n free r1 and 2 n +# H; each side of the cmp_ge_constant_qft pair is one extended +# (n+1)-bit QFT sandwich = n(n+1) controlled-r1, n + 1 free r1 and +# 2 (n + 1) H (K >= 1; K = 0 emits no rotations). + +_RESOURCES = pytest.mark.skipif( + not hasattr(cudaq, "estimate_resources"), + reason="cudaq.estimate_resources is not available in this CUDA-Q") + + +@cudaq.kernel +def _res_add_register(n: int, subtract: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + if subtract == 0: + arith.add_register(a, b, carry) + else: + arith.subtract_register(a, b, carry) + + +@cudaq.kernel +def _res_add_constant(n: int, bits: list[int], subtract: int): + target = cudaq.qvector(n) + work = cudaq.qvector(n) + carry = cudaq.qvector(1) + if subtract == 0: + arith.add_constant(target, bits, work, carry) + else: + arith.subtract_constant(target, bits, work, carry) + + +@cudaq.kernel +def _res_cmp_ge_constant(n: int, bits: list[int], k_is_zero: int): + x_reg = cudaq.qvector(n) + work = cudaq.qvector(n) + carry = cudaq.qvector(1) + out = cudaq.qvector(1) + arith.cmp_ge_constant(x_reg, bits, k_is_zero, work, carry, out) + + +@cudaq.kernel +def _res_add_constant_qft(n: int, constant: int): + target = cudaq.qvector(n) + arith.add_constant_qft(target, constant) + + +@cudaq.kernel +def _res_cmp_ge_constant_qft(n: int, constant: int, invert: int): + x_reg = cudaq.qvector(n) + out = cudaq.qvector(1) + arith.cmp_ge_constant_qft(x_reg, out, constant, invert) + + +def _toffolis(kernel, *args) -> int: + # A Toffoli is an x with two controls; ``count_controls`` is the + # arity-aware accessor (``count("ccx")`` matches nothing — the + # display name is not the lookup key and returns 0). + return cudaq.estimate_resources(kernel, *args).count_controls("x", 2) + + +@_RESOURCES +@pytest.mark.parametrize("n", WIDTHS) +def test_cdkm_register_ops_cost_exactly_2n_toffolis(n): + assert _toffolis(_res_add_register, n, 0) == 2 * n + assert _toffolis(_res_add_register, n, 1) == 2 * n + + +@_RESOURCES +@pytest.mark.parametrize("n", WIDTHS) +def test_cdkm_constant_ops_cost_exactly_2n_toffolis(n): + bits = _bits((1 << n) - 1, n) # worst-case load: every bit set + assert _toffolis(_res_add_constant, n, bits, 0) == 2 * n + assert _toffolis(_res_add_constant, n, bits, 1) == 2 * n + + +@_RESOURCES +@pytest.mark.parametrize("n", WIDTHS) +def test_cdkm_comparator_costs_exactly_2n_toffolis(n): + complement = _bits((1 << n) - 1, n) # K = 1 + assert _toffolis(_res_cmp_ge_constant, n, complement, 0) == 2 * n + # K = 0 short-circuits to a single X: no Toffolis. + assert _toffolis(_res_cmp_ge_constant, n, [0] * n, 1) == 0 + + +@_RESOURCES +def test_cdkm_toffoli_cost_grows_linearly_per_doubling(): + # The cost is exactly linear (2 n): each width doubling doubles it. + compiled = {n: _toffolis(_res_add_register, n, 0) for n in (2, 4, 8, 16)} + for n in (2, 4, 8): + assert compiled[2 * n] == 2 * compiled[n], compiled + + +@_RESOURCES +@pytest.mark.parametrize("n", WIDTHS) +def test_qft_add_constant_costs_rotations_not_toffolis(n): + resources = cudaq.estimate_resources(_res_add_constant_qft, n, 1) + assert resources.count_controls("x", 2) == 0 + assert resources.count_controls("r1", 1) == n * (n - 1) + assert resources.count_controls("r1", 0) == n + assert resources.count("h") == 2 * n + + +@_RESOURCES +@pytest.mark.parametrize("n", WIDTHS) +def test_qft_comparator_costs_rotations_not_toffolis(n): + for invert in (0, 1): + resources = cudaq.estimate_resources(_res_cmp_ge_constant_qft, n, 1, + invert) + assert resources.count_controls("x", 2) == 0 + assert resources.count_controls("r1", 1) == n * (n + 1) + assert resources.count_controls("r1", 0) == n + 1 + assert resources.count("h") == 2 * (n + 1) + assert resources.count("x") == (1 if invert == 0 else 0) + # K = 0 emits no rotations at all: the constant-true comparator is a + # bare X on the out qubit. + trivial = cudaq.estimate_resources(_res_cmp_ge_constant_qft, n, 0, 0) + assert trivial.count("r1") == 0 + assert trivial.count("x") == 1 From a412e3b6399e4a45f5aa38d70c5abc055c3b292a Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Fri, 11 Sep 2026 21:54:09 +0000 Subject: [PATCH 2/6] Add the register-register >= / > comparators to the CDKM arithmetic family MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit cmp_ge_register(a, b, carry, out) XOR-loads (a >= b) (unsigned, equal-width little-endian registers) into the out qubit and restores a, b and the carry ancilla: complement b in place (2^n - b = ~b + 1, carry-in set to 1), MAJ sweep, copy the top ripple carry into out, reverse the sweep and undo the complement. Because b doubles as the constant-load register, no n-qubit work register is needed — only the one-qubit carry — and the price is the family contract of exactly 2n Toffolis, pinned against cudaq.estimate_resources with runtime-argument harnesses at widths 1-5 and 8. cmp_gt_register is the free strict variant (a > b <=> not (b >= a): one X plus the >= comparator with the roles swapped), same price. Both are self-inverse (compute-copy- uncompute operand action, XOR-accumulated flag), pinned by apply-twice-is-identity tests on random entangled states, alongside exhaustive truth tables (widths 1-4, flag initially 0 and 1) and superposed-register statevector checks against dense NumPy references. This is the comparator the alias-sampling conditional (keep-value vs alt-value test) and the first-quantized momentum comparisons need (colleague review feedback; GAP_first_quantized_lcu.md). Signed-off-by: Scott Thornton --- .../cudaq_algorithms/primitives/__init__.py | 10 +- .../primitives/_arithmetic.py | 83 ++++++++++- tests/python/test_primitives_arithmetic.py | 140 ++++++++++++++++++ 3 files changed, 225 insertions(+), 8 deletions(-) diff --git a/python/cudaq_algorithms/primitives/__init__.py b/python/cudaq_algorithms/primitives/__init__.py index 2c0d2c6..263e1c6 100644 --- a/python/cudaq_algorithms/primitives/__init__.py +++ b/python/cudaq_algorithms/primitives/__init__.py @@ -23,7 +23,8 @@ are the in-place little-endian adders and comparators the lookup-based constructions compose with: the CDKM/Cuccaro ripple-carry family (``add_register`` / ``subtract_register``, ``add_constant`` / -``subtract_constant``, ``cmp_ge_constant``) and the ancilla-free Draper +``subtract_constant``, ``cmp_ge_constant``, and the register-register +comparators ``cmp_ge_register`` / ``cmp_gt_register``) and the ancilla-free Draper QFT family (``qft`` / ``iqft``, ``add_constant_qft`` / ``subtract_constant_qft``, the ``cmp_ge_constant_qft`` / ``cmp_ge_constant_qft_adj`` pair). Every inverse is hand-written and the @@ -35,8 +36,9 @@ from ._arithmetic import (add_constant, add_constant_qft, add_register, cmp_ge_constant, cmp_ge_constant_qft, - cmp_ge_constant_qft_adj, iqft, phase_add_constant, - qft, subtract_constant, subtract_constant_qft, + cmp_ge_constant_qft_adj, cmp_ge_register, + cmp_gt_register, iqft, phase_add_constant, qft, + subtract_constant, subtract_constant_qft, subtract_register) from ._qrom import QROM from ._unary_iteration import UnaryIterationKernels, unary_iteration_kernels @@ -50,6 +52,8 @@ "cmp_ge_constant", "cmp_ge_constant_qft", "cmp_ge_constant_qft_adj", + "cmp_ge_register", + "cmp_gt_register", "iqft", "phase_add_constant", "qft", diff --git a/python/cudaq_algorithms/primitives/_arithmetic.py b/python/cudaq_algorithms/primitives/_arithmetic.py index 4bb0297..86fb267 100644 --- a/python/cudaq_algorithms/primitives/_arithmetic.py +++ b/python/cudaq_algorithms/primitives/_arithmetic.py @@ -8,15 +8,22 @@ - **CDKM/Cuccaro ripple-carry** (`arXiv:quant-ph/0410184`): ``add_register`` / ``subtract_register`` (in-place ``b <- b +/- a``), ``add_constant`` / ``subtract_constant`` (constant loaded into a caller - provided work register), and ``cmp_ge_constant`` (a ``x >= K`` comparator + provided work register), ``cmp_ge_constant`` (a ``x >= K`` comparator writing into an out qubit, leaving ``x`` untouched: MAJ sweep, copy the - carry, reverse MAJ sweep). Toffoli prices (documented contracts, pinned - by the resource tests in ``tests/python/test_primitives_arithmetic.py``): + carry, reverse MAJ sweep), and the register-register comparators + ``cmp_ge_register`` / ``cmp_gt_register`` (``out ^= (a >= b)`` / + ``out ^= (a > b)`` on equal-width unsigned registers, same + MAJ-copy-reverse structure with ``b`` complemented in place, so no + ``n``-qubit work register at all — only the one-qubit carry). Toffoli + prices (documented contracts, pinned by the resource tests in + ``tests/python/test_primitives_arithmetic.py``): ``add_register`` / ``subtract_register`` cost exactly ``2 n`` Toffolis on ``n``-bit registers (``n`` in the MAJ sweep, ``n`` in the UMA sweep), ``add_constant`` / ``subtract_constant`` inherit the same ``2 n`` (the - constant load is X-only), and ``cmp_ge_constant`` costs ``2 n`` for - ``K >= 1`` (MAJ sweep plus its reversal) and ``0`` for ``K = 0``. + constant load is X-only), ``cmp_ge_constant`` costs ``2 n`` for + ``K >= 1`` (MAJ sweep plus its reversal) and ``0`` for ``K = 0``, and + ``cmp_ge_register`` / ``cmp_gt_register`` cost exactly ``2 n`` (the + complement and flag copy are X/CNOT-only). - **Draper QFT arithmetic** (`arXiv:quant-ph/0008033`): ``qft`` / ``iqft`` and ``add_constant_qft``, plus the ``cmp_ge_constant_qft`` / ``cmp_ge_constant_qft_adj`` comparator pair. The QFT family needs *no* @@ -55,6 +62,8 @@ "add_constant", "subtract_constant", "cmp_ge_constant", + "cmp_ge_register", + "cmp_gt_register", "qft", "iqft", "phase_add_constant", @@ -212,6 +221,70 @@ def cmp_ge_constant(x_reg: cudaq.qview, complement_bits: list[int], x(work[k]) +@cudaq.kernel +def cmp_ge_register(a: cudaq.qview, b: cudaq.qview, carry: cudaq.qview, + out: cudaq.qview): + """``out[0] ^= (a >= b)`` (CDKM, unsigned), leaving ``a``, ``b`` unchanged. + + ``a`` and ``b`` have equal size ``n >= 1``; ``carry`` is a one-qubit + view that must be |0> on entry (it is returned to |0>); ``out`` is the + one-qubit flag, XOR-loaded (any input state is allowed). ``a``, ``b``, + ``carry`` and ``out`` must be pairwise disjoint (module precondition). + + Uses the ripple identity ``a >= b <=> carry_out(a + (2^n - b))`` with + ``2^n - b = ~b + 1``: complement ``b`` in place and set the carry-in + (so ``b`` doubles as the constant-load work register — no extra + ``n``-qubit work is needed, unlike ``cmp_ge_constant``), MAJ sweep, + copy the carry-out (on the top ``b`` qubit) into ``out``, reverse the + MAJ sweep and undo the complement so ``a``, ``b`` and ``carry`` are + all restored. Because the flag is XOR-accumulated and the operand + action is compute-copy-uncompute, applying the kernel twice with the + same arguments is the identity (self-inverse). + """ + n = a.size() + if n > 0: + # 2^n - b = ~b + 1: complement b in place, set the carry-in to 1. + for k in range(n): + x(b[k]) + x(carry[0]) + # MAJ sweep of a + ~b + 1 (b accumulates the ripple carries). + cx(b[0], a[0]) + cx(b[0], carry[0]) + x.ctrl(carry[0], a[0], b[0]) + for i in range(1, n): + cx(b[i], a[i]) + cx(b[i], b[i - 1]) + x.ctrl(b[i - 1], a[i], b[i]) + # a >= b iff the sum carries out of the top bit. + cx(b[n - 1], out[0]) + # Reverse MAJ sweep (restores a, ~b and the carry-in). + for k in range(1, n): + i = n - k + x.ctrl(b[i - 1], a[i], b[i]) + cx(b[i], b[i - 1]) + cx(b[i], a[i]) + x.ctrl(carry[0], a[0], b[0]) + cx(b[0], carry[0]) + cx(b[0], a[0]) + x(carry[0]) + for k in range(n): + x(b[k]) + + +@cudaq.kernel +def cmp_gt_register(a: cudaq.qview, b: cudaq.qview, carry: cudaq.qview, + out: cudaq.qview): + """``out[0] ^= (a > b)`` (CDKM, unsigned), leaving ``a``, ``b`` unchanged. + + Free strict variant of ``cmp_ge_register`` via ``a > b <=> not + (b >= a)``: an X on ``out`` plus the ``>=`` comparator with the roles + swapped. Same preconditions, same ``2 n`` Toffoli price, and likewise + self-inverse. + """ + x(out[0]) + cmp_ge_register(b, a, carry, out) + + # ============================================================================ # Draper QFT family (no work qubits) # ============================================================================ diff --git a/tests/python/test_primitives_arithmetic.py b/tests/python/test_primitives_arithmetic.py index ca60da6..bfa3184 100644 --- a/tests/python/test_primitives_arithmetic.py +++ b/tests/python/test_primitives_arithmetic.py @@ -272,6 +272,123 @@ def test_cmp_ge_constant_qft_all_inputs(n): atol=1e-10) +# ---------------------------------------------------------------------- +# Register-register >= / > comparators (CDKM) +# ---------------------------------------------------------------------- +# +# Layout: a at bits [0, n), b at [n, 2n), carry at 2n, out at 2n + 1. + + +@cudaq.kernel +def _run_cmp_register_basis(n: int, aval: int, bval: int, flag_init: int, + strict: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + out = cudaq.qvector(1) + for k in range(n): + if ((aval >> k) & 1) == 1: + x(a[k]) + if ((bval >> k) & 1) == 1: + x(b[k]) + if flag_init == 1: + x(out[0]) + if strict == 0: + arith.cmp_ge_register(a, b, carry, out) + else: + arith.cmp_gt_register(a, b, carry, out) + + +@cudaq.kernel +def _run_cmp_register_superposed(n: int, strict: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + out = cudaq.qvector(1) + for k in range(n): + h(a[k]) + h(b[k]) + if strict == 0: + arith.cmp_ge_register(a, b, carry, out) + else: + arith.cmp_gt_register(a, b, carry, out) + + +@cudaq.kernel +def _run_cmp_register_twice(n: int, angles: list[float], strict: int, + apply_ops: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + out = cudaq.qvector(1) + for k in range(n): + ry(angles[k], a[k]) + ry(angles[n + k], b[k]) + cx(a[k], b[k]) # entangle a with b + ry(angles[2 * n], out[0]) + if apply_ops == 1: + if strict == 0: + arith.cmp_ge_register(a, b, carry, out) + arith.cmp_ge_register(a, b, carry, out) + else: + arith.cmp_gt_register(a, b, carry, out) + arith.cmp_gt_register(a, b, carry, out) + + +_CMP_REGISTER_OPS = [ + ("ge", 0, lambda a, b: a >= b), + ("gt", 1, lambda a, b: a > b), +] + + +@pytest.mark.parametrize("n", [1, 2, 3, 4]) +@pytest.mark.parametrize("name,strict,predicate", _CMP_REGISTER_OPS) +def test_cmp_register_truth_table_exhaustive(n, name, strict, predicate): + # All (a, b) pairs, flag initially 0 and 1: XOR semantics, operands + # (and carry) verified unchanged by the full statevector equality. + for aval in range(1 << n): + for bval in range(1 << n): + for flag_init in (0, 1): + state = np.array( + cudaq.get_state(_run_cmp_register_basis, n, aval, bval, + flag_init, strict)) + flag = flag_init ^ int(predicate(aval, bval)) + index = aval + (bval << n) + (flag << (2 * n + 1)) + np.testing.assert_allclose(state, + _basis(index, 2 * n + 2), + atol=1e-12) + + +@pytest.mark.parametrize("n", WIDTHS) +@pytest.mark.parametrize("name,strict,predicate", _CMP_REGISTER_OPS) +def test_cmp_register_all_inputs_superposed(n, name, strict, predicate): + state = np.array(cudaq.get_state(_run_cmp_register_superposed, n, strict)) + expected = np.zeros(1 << (2 * n + 2), dtype=np.complex128) + norm = 1.0 / (1 << n) + for aval in range(1 << n): + for bval in range(1 << n): + index = aval + (bval << n) + (int(predicate(aval, bval)) << + (2 * n + 1)) + expected[index] += norm + np.testing.assert_allclose(state, expected, atol=1e-12) + + +@pytest.mark.parametrize("n", WIDTHS) +@pytest.mark.parametrize("name,strict,predicate", _CMP_REGISTER_OPS) +def test_cmp_register_twice_is_identity(n, name, strict, predicate): + # Self-inverse contract: the operand action is compute-copy-uncompute + # and the flag is XOR-accumulated, so two applications are the + # identity — checked on a random entangled state (superposed flag + # included: XOR-loading is a permutation, so linearity carries it). + rng = np.random.default_rng(20260911 + 8 * n + strict) + angles = rng.uniform(0.1, 3.0, size=2 * n + 1).tolist() + twice = np.array( + cudaq.get_state(_run_cmp_register_twice, n, angles, strict, 1)) + reference = np.array( + cudaq.get_state(_run_cmp_register_twice, n, angles, strict, 0)) + np.testing.assert_allclose(twice, reference, atol=1e-12) + + # ---------------------------------------------------------------------- # Direct basis-state spot checks (readable, non-superposed) # ---------------------------------------------------------------------- @@ -314,6 +431,8 @@ def test_add_register_basis_spot_checks(): # (free), so the price is the inner register adder's 2 n. # - cmp_ge_constant (K >= 1): a MAJ sweep (n) plus its literal reversal # (n) with a free CNOT carry-copy in between = 2 n; K = 0 is a bare X. +# - cmp_ge_register / cmp_gt_register: the same MAJ sweep + reversal +# = 2 n; the b-complement, carry-in set and flag copy are X/CNOT-only. # # Draper QFT derivations: no Toffolis at all. add_constant_qft is # qft + phases + iqft = 2 * (n(n-1)/2) controlled-r1, n free r1 and 2 n @@ -357,6 +476,18 @@ def _res_cmp_ge_constant(n: int, bits: list[int], k_is_zero: int): arith.cmp_ge_constant(x_reg, bits, k_is_zero, work, carry, out) +@cudaq.kernel +def _res_cmp_register(n: int, strict: int): + a = cudaq.qvector(n) + b = cudaq.qvector(n) + carry = cudaq.qvector(1) + out = cudaq.qvector(1) + if strict == 0: + arith.cmp_ge_register(a, b, carry, out) + else: + arith.cmp_gt_register(a, b, carry, out) + + @cudaq.kernel def _res_add_constant_qft(n: int, constant: int): target = cudaq.qvector(n) @@ -401,6 +532,15 @@ def test_cdkm_comparator_costs_exactly_2n_toffolis(n): assert _toffolis(_res_cmp_ge_constant, n, [0] * n, 1) == 0 +@_RESOURCES +@pytest.mark.parametrize("n", WIDTHS + [8]) +def test_cmp_register_costs_exactly_2n_toffolis(n): + # Widths past the truth-table range (5 and 8) included: the count is + # a function of the runtime width argument, never a folded constant. + assert _toffolis(_res_cmp_register, n, 0) == 2 * n + assert _toffolis(_res_cmp_register, n, 1) == 2 * n + + @_RESOURCES def test_cdkm_toffoli_cost_grows_linearly_per_doubling(): # The cost is exactly linear (2 n): each width doubling doubles it. From fe37bd61759e608cc9bf1e69de74e9760f1242b7 Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Tue, 25 Aug 2026 02:36:15 +0000 Subject: [PATCH 3/6] Add the coherent alias-sampling PREPARE to primitives AliasSamplingPrepare: the Babbush-style (arXiv:1805.03662, Sec. III.D) PREPARE-with-garbage whose index-register marginal realizes the integer-Vose table distribution exactly (per-bin discretization bound 1.5 / (K 2^mu), derived not assumed), with a hand-written adjoint. The (alias, keep) lookup is the variants QROM under variant="auto": the ladder slice inside the garbage register is qrom.num_ladder wide, so num_garbage = m + 2mu + 4 + qrom.num_ladder (the historical 2m + 2mu + 4 whenever the plain select walk wins the pricing). No variant is pinned - every variant restores its ladder and self-inverts on the clean-ladder sector, which is what the adjoint's second lookup uses. The new qrom property exposes the priced lookup. Tests pin the index marginal against an independent brute-force branch enumeration, the table distribution against the ideal within the derived bound, prepare-then-adjoint identity, degenerate weights and validation messages, plus compiler-pinned resource contracts: each PREPARE side costs exactly qrom.toffoli_count + 4 (mu + 1) Toffolis (the lookup at the QROM's own reported price plus the two CDKM comparator adders) and num_index Fredkins, adjoint == forward gate-for-gate. Review fixes: pass variant/block_size through to the QROM constructor as keyword-only parameters (they shape num_garbage via qrom.num_ladder; QROM validation errors propagate as-is) and pin the variant-agnostic claim with a forced select_swap test asserting the exact marginal, the prepare-then-adjoint identity and the Toffoli cost identity; reject bool mu explicitly while keeping integral numpy ints; document that the per-bin garbage vectors are not mutually orthogonal (only the index amplitude magnitudes are guaranteed - exactly the qubitization <0|PREP' SELECT PREP|0> contract); say "the available QROM variants" instead of naming constructions this base does not ship; version-scope the cudaq.control nested-call rejection (through CUDA-Q 0.15, lifted in 0.16); mark the h-count assertion as incidental and the toffoli_count cross-check as a deliberate bookkeeping-vs-compiler pin. Signed-off-by: Scott Thornton --- .../cudaq_algorithms/primitives/__init__.py | 9 + .../primitives/_alias_sampling.py | 385 ++++++++++++++++++ .../python/test_primitives_alias_sampling.py | 300 ++++++++++++++ 3 files changed, 694 insertions(+) create mode 100644 python/cudaq_algorithms/primitives/_alias_sampling.py create mode 100644 tests/python/test_primitives_alias_sampling.py diff --git a/python/cudaq_algorithms/primitives/__init__.py b/python/cudaq_algorithms/primitives/__init__.py index 263e1c6..ed1fdae 100644 --- a/python/cudaq_algorithms/primitives/__init__.py +++ b/python/cudaq_algorithms/primitives/__init__.py @@ -30,10 +30,18 @@ ``cmp_ge_constant_qft_adj`` pair). Every inverse is hand-written and the gate prices are compiler-pinned by the resource tests. +``AliasSamplingPrepare`` composes both: the coherent alias-sampling +PREPARE-with-garbage of Babbush et al. (`arXiv:1805.03662`, Sec. III.D) +— integer Vose preprocessing, a QROM lookup of the (alias, keep) table +priced by ``variant="auto"``, and the CDKM comparator — realizing a +weighted index-register marginal exactly, up to a derived ``mu``-bit +discretization bound. + Import the subpackage directly (``from cudaq_algorithms.primitives import QROM``); nothing here is re-exported from the package root. """ +from ._alias_sampling import AliasSamplingPrepare from ._arithmetic import (add_constant, add_constant_qft, add_register, cmp_ge_constant, cmp_ge_constant_qft, cmp_ge_constant_qft_adj, cmp_ge_register, @@ -44,6 +52,7 @@ from ._unary_iteration import UnaryIterationKernels, unary_iteration_kernels __all__ = [ + "AliasSamplingPrepare", "QROM", "UnaryIterationKernels", "add_constant", diff --git a/python/cudaq_algorithms/primitives/_alias_sampling.py b/python/cudaq_algorithms/primitives/_alias_sampling.py new file mode 100644 index 0000000..193512b --- /dev/null +++ b/python/cudaq_algorithms/primitives/_alias_sampling.py @@ -0,0 +1,385 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +"""Coherent alias sampling: PREPARE-with-garbage for weighted indices. + +``AliasSamplingPrepare(weights, mu)`` mints the Babbush-style +(`arXiv:1805.03662`, Sec. III.D) state preparation whose *index-register +marginal* is the normalized weight distribution ``w / lambda`` up to a +``mu``-bit discretization: classical Vose preprocessing builds integer +keep/alias tables, and the circuit is + +1. uniform superposition over the ``K = 2^num_index`` bins (``weights`` + is zero-padded to a power of two; padded bins carry exactly zero + probability), +2. a :class:`~cudaq_algorithms.primitives.QROM` lookup of + ``(alias_k, keep_k)`` at bin ``k`` (``variant="auto"``: the lookup is + priced across the available QROM variants and the cheapest + construction is minted — its clean-ancilla ladder lives inside the + garbage register, so ``num_garbage`` follows the chosen variant), +3. a comparator of ``keep_k`` against a uniform ``mu``-bit reference + register (built from the CDKM register adders of ``_arithmetic``: + subtract on a ``mu+1``-bit extension, copy the borrow, add back), and +4. a controlled swap of the bin index with the alias when the reference + is not below ``keep_k``. + +PREPARE with garbage — read this before consuming +------------------------------------------------- + +The output is **not** ``sum_k sqrt(w_k / lambda) |k> |0...0>``: the +alias/keep/reference/flag registers stay *entangled* with the index (only +the amplitude *magnitudes* marginalize to ``w / lambda``). That is +exactly the coherent-alias-sampling contract: qubitization-style +consumers tolerate the garbage only because every use is the symmetric +sandwich ``PREPARE ... PREPARE^dagger`` — the garbage registers are +uncomputed by the hand-written ``adjoint_kernel()`` applied to the same +registers, never discarded, never reflected over while dirty, and never +consumed by any other circuit element in between. Anything that needs a +clean ``sum sqrt(p_k) |k>`` (e.g. amplitude arithmetic on the index +alone) must use a different preparation. + +Discretization: the integer tables represent bin probabilities exactly as +``W_k / (K 2^mu)`` (integer Vose is exact, and full bins are stored +self-aliased so the ``keep = 2^mu`` edge loses nothing); the only error +is rounding ``w_k / lambda`` to the integers ``W_k``, bounded per bin by +``1.5 / (K 2^mu)`` (half-ulp rounding plus at most one unit of residual +redistribution) — exposed as ``discretization_bound`` and derived, never +assumed, in the tests. + +Kernel signatures (little-endian; ``docs/conventions.md``): both +``kernel()`` and ``adjoint_kernel()`` are ``(index: qview, garbage: +qview)`` with ``index`` of width ``num_index`` and ``garbage`` of width +``num_garbage`` laid out as ``[alias(num_index) | keep(mu) | keep_pad | +ref(mu) | ref_pad | flag | ladder(qrom.num_ladder) | carry]`` — i.e. +``num_garbage = num_index + 2 mu + 3 + qrom.num_ladder + 1``, which is +``2 num_index + 2 mu + 4`` whenever the priced lookup is the plain +``"select"`` walk (``num_ladder = num_index``; always the case for the +small tables where select wins ``"auto"``'s pricing). Both registers +must be |0...0> on entry to ``kernel()``; ``adjoint_kernel()`` is the +literal gate-reversal (no ``cudaq.adjoint``, cuda-quantum#4897/#4898) and +returns them to |0...0>. The QROM lookup uncomputes itself on the clean +ladder sector the adjoint provides, whatever the variant (the ladder is +restored by the compute pass and untouched in between, so re-applying +the lookup XORs the table back out). These kernels call the QROM and +CDKM sub-kernels and are therefore NOT flat: do not place them under +``cudaq.control`` (control-variant generation rejects nested kernel +calls through CUDA-Q 0.15; the rejection is lifted in 0.16, but the +contract stands) — consumers control SELECT, not PREPARE. + +Cost (pinned by ``tests/python/test_primitives_alias_sampling.py`` +against the compiler): each of ``kernel()`` / ``adjoint_kernel()`` costs +exactly ``qrom.toffoli_count + 4 (mu + 1)`` Toffolis — the lookup at the +QROM's own reported price plus two CDKM register adders on the +``mu+1``-bit extension at ``2 (mu + 1)`` each — plus ``num_index`` +controlled swaps (Fredkins) for the alias swap. +""" + +from __future__ import annotations + +import math +from collections.abc import Sequence + +import cudaq +import numpy as np + +from ._arithmetic import add_register, subtract_register +from ._qrom import QROM +from ._unary_iteration import _retain + +__all__ = ["AliasSamplingPrepare"] + + +def _vose_tables(probabilities: Sequence[float], + mu: int) -> tuple[list[int], list[int], list[float]]: + """Integer Vose preprocessing at ``mu`` bits. + + Returns ``(keep, alias, table_probabilities)`` over the padded bins: + ``keep[k] in [0, 2^mu - 1]``, ``alias`` a bin index (full bins are + self-aliased), and ``table_probabilities[k] = W_k / (K 2^mu)`` — the + *exact* probability the circuit realizes for bin ``k``. + """ + num_bins = len(probabilities) + unit = 1 << mu + total = num_bins * unit + + # Integer weights summing exactly to K * 2^mu: round, then push the + # rounding residual back one unit at a time onto the bins whose + # rounding moved the furthest in the residual's direction (weights + # stay non-negative: a downward push targets a bin with W >= 1). + ideal = [p * total for p in probabilities] + weights = [int(round(v)) for v in ideal] + residual = total - sum(weights) + step = 1 if residual > 0 else -1 + order = sorted(range(num_bins), + key=lambda k: step * (ideal[k] - weights[k]), + reverse=True) + cursor = 0 + for _ in range(abs(residual)): + while step < 0 and weights[order[cursor % num_bins]] == 0: + cursor += 1 + weights[order[cursor % num_bins]] += step + cursor += 1 + + # Integer Vose: pair small bins (W < 2^mu) with large ones (W > 2^mu). + keep = [unit] * num_bins + alias = list(range(num_bins)) + remaining = list(weights) + small = [k for k in range(num_bins) if remaining[k] < unit] + large = [k for k in range(num_bins) if remaining[k] > unit] + while small and large: + s = small.pop() + l = large[-1] + keep[s] = remaining[s] + alias[s] = l + remaining[l] -= unit - remaining[s] + if remaining[l] <= unit: + large.pop() + if remaining[l] < unit: + small.append(l) + # Bins left at exactly 2^mu are full: stored self-aliased with the + # keep clamped into mu bits — self-aliasing makes the keep value + # irrelevant (both comparator branches land on the same bin), so the + # clamp is exact, not an approximation. + keep = [min(v, unit - 1) for v in keep] + table = [w / total for w in weights] + return keep, alias, table + + +class AliasSamplingPrepare: + """Coherent alias-sampling PREPARE (see the module docstring). + + Parameters + ---------- + weights + Non-negative, finite weights with a positive sum; zero entries + are allowed (their bins get exactly zero probability). Padded + with zero-weight bins to the next power of two. + mu + Keep-threshold precision in bits (>= 1); the per-bin + discretization error is bounded by ``discretization_bound``. + variant + Keyword-only; forwarded verbatim to the :class:`QROM` + constructor (default ``"auto"``: price the variants and mint + the cheapest). The chosen variant controls the garbage-register + size through ``qrom.num_ladder`` (see ``num_garbage``). + block_size + Keyword-only; forwarded verbatim to the :class:`QROM` + constructor (only meaningful with ``variant="select_swap"``; + ``None`` lets the QROM pick). Also shapes ``qrom.num_ladder`` + and hence ``num_garbage``. Invalid combinations raise the + QROM's own ``ValueError``. + + Garbage caveat: the per-bin garbage vectors ``|g_k>`` are NOT + mutually orthogonal (measured off-diagonal overlaps reach ~0.19), so + the reduced state on the index register is not ``diag(p_k)``; only + the amplitude *magnitudes* on the index register are guaranteed. + That is exactly what qubitization's ``<0| PREP^dagger SELECT PREP + |0>`` contract needs — index orthogonality kills the cross terms — + but consumers must not assume orthonormal garbage. + """ + + def __init__(self, + weights: Sequence[float], + mu: int, + *, + variant: str = "auto", + block_size: int | None = None) -> None: + values = [float(w) for w in weights] + if len(values) == 0: + raise ValueError("weights must be non-empty") + if any(not math.isfinite(w) for w in values): + raise ValueError("weights must be finite") + if any(w < 0.0 for w in values): + raise ValueError("weights must be non-negative") + lam = sum(values) + if not lam > 0.0: + raise ValueError("weights sum to zero: nothing to prepare") + if isinstance(mu, bool) or int(mu) != mu or mu < 1: + raise ValueError("mu must be a positive integer (mu >= 1)") + mu = int(mu) + + num_index = max(1, (len(values) - 1).bit_length()) + num_bins = 1 << num_index + padded = values + [0.0] * (num_bins - len(values)) + probabilities = [w / lam for w in padded] + keep, alias, table = _vose_tables(probabilities, mu) + + self._lam = lam + self._mu = mu + self._num_index = num_index + self._num_bins = num_bins + self._num_weights = len(values) + self._keep = tuple(keep) + self._alias = tuple(alias) + self._probabilities = np.asarray(probabilities) + self._table_probabilities = np.asarray(table) + + # QROM data: alias in the low num_index bits, keep above it. + # variant/block_size pass through verbatim (default "auto" + # prices the constructions); the chosen variant's ladder width + # shapes the garbage layout below. QROM validation errors + # propagate as-is. + self._qrom = QROM( + [alias[k] | (keep[k] << num_index) for k in range(num_bins)], + address_bits=num_index, + output_bits=num_index + mu, + variant=variant, + block_size=block_size) + self._build_kernels() + + # ------------------------------------------------------------------ + # Kernel construction + # ------------------------------------------------------------------ + + def _build_kernels(self) -> None: + # Garbage layout offsets (see the module docstring); everything + # unpacked into scalar locals (no tuple/self capture in kernels). + m = self._num_index + mu = self._mu + num_ladder = self._qrom.num_ladder + k0 = m # keep (keep_pad at k0 + mu) + r0 = m + mu + 1 # ref (ref_pad at r0 + mu) + flag = m + 2 * mu + 2 + l0 = m + 2 * mu + 3 + c0 = l0 + num_ladder # carry + qrom_kernel = self._qrom.kernel + + @cudaq.kernel + def primitives_alias_prepare(index: cudaq.qview, garbage: cudaq.qview): + # Uniform superpositions over bins and the mu-bit reference. + for b in range(m): + h(index[b]) + for b in range(mu): + h(garbage[r0 + b]) + # (alias_k, keep_k) lookup: output is garbage[0 : m + mu]. + qrom_kernel(index, garbage[l0:l0 + num_ladder], garbage[0:m + mu]) + # flag <- NOT (ref < keep): subtract keep on the (mu+1)-bit + # extension of ref, copy the borrow (MSB), add keep back. + subtract_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], + garbage[c0:c0 + 1]) + cx(garbage[r0 + mu], garbage[flag]) + x(garbage[flag]) + add_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], + garbage[c0:c0 + 1]) + # Swap the bin index with its alias on the flag. + for b in range(m): + swap.ctrl(garbage[flag], index[b], garbage[b]) + + @cudaq.kernel + def primitives_alias_prepare_adj(index: cudaq.qview, + garbage: cudaq.qview): + """Hand-written gate-reversal of ``primitives_alias_prepare``.""" + for j in range(m): + b = m - 1 - j + swap.ctrl(garbage[flag], index[b], garbage[b]) + subtract_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], + garbage[c0:c0 + 1]) + x(garbage[flag]) + cx(garbage[r0 + mu], garbage[flag]) + add_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], + garbage[c0:c0 + 1]) + # The QROM lookup uncomputes itself on the clean-ladder + # sector: the ladder is |0> here (restored by the compute + # pass, untouched since), so re-applying the lookup XORs the + # table back out of garbage[0 : m + mu] — the documented + # self-inverse contract of every QROM variant. + qrom_kernel(index, garbage[l0:l0 + num_ladder], garbage[0:m + mu]) + for j in range(mu): + b = mu - 1 - j + h(garbage[r0 + b]) + for j in range(m): + b = m - 1 - j + h(index[b]) + + _retain(primitives_alias_prepare, primitives_alias_prepare_adj) + self._prepare = primitives_alias_prepare + self._prepare_adj = primitives_alias_prepare_adj + + # ------------------------------------------------------------------ + # Inspection + # ------------------------------------------------------------------ + + @property + def lam(self) -> float: + """The weight one-norm ``lambda = sum_k w_k``.""" + return self._lam + + @property + def mu(self) -> int: + return self._mu + + @property + def num_index(self) -> int: + """Index register width (``log2`` of the padded bin count).""" + return self._num_index + + @property + def num_bins(self) -> int: + """Padded bin count ``2^num_index`` (>= ``len(weights)``).""" + return self._num_bins + + @property + def num_garbage(self) -> int: + """Garbage register width (layout in the module docstring).""" + return self._num_index + 2 * self._mu + 4 + self._qrom.num_ladder + + @property + def qrom(self) -> QROM: + """The minted ``(alias, keep)`` lookup: its ``variant``, + ``num_ladder`` and ``toffoli_count`` are the priced lookup facts + the PREPARE inherits.""" + return self._qrom + + @property + def ladder_offset(self) -> int: + """Offset of the ``qrom.num_ladder``-wide QROM ladder inside + garbage. + + The ladder qubits are clean (|0>) between ``kernel()`` and + ``adjoint_kernel()`` — unlike the rest of the garbage — so a + SELECT sandwiched between them may reuse + ``garbage[ladder_offset : ladder_offset + qrom.num_ladder]`` as + its own ladder. + """ + return self._num_index + 2 * self._mu + 3 + + @property + def keep(self) -> tuple[int, ...]: + """Per-bin mu-bit keep thresholds (full bins clamped, self-aliased).""" + return self._keep + + @property + def alias(self) -> tuple[int, ...]: + return self._alias + + @property + def probabilities(self) -> np.ndarray: + """The ideal padded distribution ``w / lambda`` (zero-padded).""" + return self._probabilities.copy() + + @property + def table_probabilities(self) -> np.ndarray: + """The exactly realized index marginal (integer-table exact).""" + return self._table_probabilities.copy() + + @property + def discretization_bound(self) -> float: + """Per-bin bound on |table - ideal| probability (see module doc).""" + return 1.5 / (self._num_bins * (1 << self._mu)) + + def __repr__(self) -> str: + return (f"AliasSamplingPrepare(weights={self._num_weights} " + f"(padded to {self.num_bins} bins), mu={self.mu}, " + f"lambda={self.lam:.6g}, index_qubits={self.num_index}, " + f"garbage_qubits={self.num_garbage})") + + # ------------------------------------------------------------------ + # Kernels + # ------------------------------------------------------------------ + + def kernel(self): + """PREPARE ``(index, garbage)`` — both |0...0> on entry.""" + return self._prepare + + def adjoint_kernel(self): + """PREPARE^dagger ``(index, garbage)`` — the exact inverse.""" + return self._prepare_adj diff --git a/tests/python/test_primitives_alias_sampling.py b/tests/python/test_primitives_alias_sampling.py new file mode 100644 index 0000000..f98b686 --- /dev/null +++ b/tests/python/test_primitives_alias_sampling.py @@ -0,0 +1,300 @@ +# SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES. All rights reserved. +# SPDX-License-Identifier: Apache-2.0 +"""Tests for the coherent alias-sampling PREPARE. + +Two-level validation: the circuit's index-register marginal (amplitude +magnitudes squared, summed over every garbage configuration) is pinned at +simulator precision against an *independent* brute-force enumeration of +the classical keep/alias branches — for bin k and reference sigma the +branch lands on k if sigma < keep[k] else on alias[k], each branch +carrying probability 1 / (K 2^mu) — and that exact table distribution is +pinned against the ideal w / lambda within the derived discretization +bound (integer weights W_k = round(p_k K 2^mu) plus at most one unit of +residual redistribution: |W_k - ideal| <= 3/2, so +|p_table - p_ideal| <= 1.5 / (K 2^mu) per bin; no bare constants). + +PREPARE followed by the hand-written PREPARE-dagger is the identity on the +full (index + garbage) register. + +The resource tests at the bottom hold the documented cost against the +compiler: each PREPARE (and its adjoint) costs exactly +``qrom.toffoli_count + 4 (mu + 1)`` Toffolis plus ``num_index`` Fredkins +— the lookup at the QROM's *own reported* price (whatever variant +``"auto"`` minted) plus two CDKM register adders on the ``mu+1``-bit +comparator extension at ``2 (mu + 1)`` Toffolis each. +""" + +import numpy as np +import pytest + +import cudaq + +from cudaq_algorithms.primitives import AliasSamplingPrepare + + +def _brute_force_marginal(prep: AliasSamplingPrepare) -> np.ndarray: + """Enumerate the classical (bin, reference) branches independently.""" + num_bins = prep.num_bins + unit = 1 << prep.mu + marginal = np.zeros(num_bins) + for k in range(num_bins): + for sigma in range(unit): + j = k if sigma < prep.keep[k] else prep.alias[k] + marginal[j] += 1.0 / (num_bins * unit) + return marginal + + +def _circuit_marginal(prep: AliasSamplingPrepare) -> np.ndarray: + prepare = prep.kernel() + num_index = prep.num_index + num_garbage = prep.num_garbage + + @cudaq.kernel + def run(): + index = cudaq.qvector(num_index) + garbage = cudaq.qvector(num_garbage) + prepare(index, garbage) + + state = np.array(cudaq.get_state(run)) + # Index register at qubits [0, num_index): its value is the low bits + # of the basis index, so the marginal sums |amplitude|^2 over the + # garbage (high) bits. + return np.sum(np.abs(state.reshape(-1, 1 << num_index))**2, axis=0) + + +def _padded_ideal(weights, num_bins: int) -> np.ndarray: + values = np.zeros(num_bins) + values[:len(weights)] = np.asarray(weights, dtype=float) + return values / values.sum() + + +CASES = [ + ([0.7, 0.2, 1.4], 2), + ([0.7, 0.2, 1.4], 4), + ([1.0, 3.5, 0.25, 2.0, 0.8], 3), + ([2.5], 3), # single weight: deterministic |0> index + ([0.0, 0.0, 7.0, 0.0], 3), # single nonzero weight + ([1.0] * 8, 2), # uniform: every bin full, exactly representable + ([0.9, 0.0, 2.2, 0.0, 0.4, 1.3, 0.0, 3.1], 2), # zero-weight entries +] + + +@pytest.mark.parametrize("weights,mu", CASES) +def test_alias_sampling_index_marginal(weights, mu): + prep = AliasSamplingPrepare(weights, mu) + brute = _brute_force_marginal(prep) + # Level 1: circuit vs the exact integer-table distribution, at + # simulator precision. + np.testing.assert_allclose(_circuit_marginal(prep), brute, atol=1e-10) + # The class's own table bookkeeping agrees with the enumeration. + np.testing.assert_allclose(prep.table_probabilities, brute, atol=1e-15) + # Level 2: table vs ideal within the derived mu-bit bound. + assert prep.discretization_bound == \ + 1.5 / (prep.num_bins * (1 << prep.mu)) + ideal = _padded_ideal(weights, prep.num_bins) + np.testing.assert_allclose(prep.probabilities, ideal, atol=1e-15) + assert np.max(np.abs(brute - ideal)) <= prep.discretization_bound + 1e-15 + + +def test_alias_sampling_larger_random_table(): + rng = np.random.default_rng(23) + weights = rng.random(16) * 3.0 + prep = AliasSamplingPrepare(weights, mu=2) + assert prep.num_index == 4 and prep.num_bins == 16 + brute = _brute_force_marginal(prep) + np.testing.assert_allclose(_circuit_marginal(prep), brute, atol=1e-10) + ideal = _padded_ideal(weights, 16) + assert np.max(np.abs(brute - ideal)) <= prep.discretization_bound + 1e-15 + + +def test_alias_sampling_exact_dyadic_weights(): + # w / lambda = [1/4, 1/4, 1/2, 0]: representable exactly at any mu, + # so the table distribution must match the ideal to fp precision. + prep = AliasSamplingPrepare([1.0, 1.0, 2.0], mu=2) + brute = _brute_force_marginal(prep) + np.testing.assert_allclose(brute, [0.25, 0.25, 0.5, 0.0], atol=1e-15) + np.testing.assert_allclose(_circuit_marginal(prep), brute, atol=1e-10) + + +def test_alias_sampling_uniform_is_exact(): + prep = AliasSamplingPrepare([1.0] * 8, mu=2) + np.testing.assert_allclose(_brute_force_marginal(prep), + np.full(8, 1.0 / 8.0), + atol=1e-15) + # Full bins are self-aliased with the keep clamped into mu bits. + assert prep.alias == tuple(range(8)) + + +@pytest.mark.parametrize("weights,mu", [([0.7, 0.2, 1.4], 2), + ([1.0, 3.5, 0.25, 2.0, 0.8], 3), + ([2.5], 3)]) +def test_alias_sampling_prepare_then_adjoint_is_identity(weights, mu): + prep = AliasSamplingPrepare(weights, mu) + prepare = prep.kernel() + unprepare = prep.adjoint_kernel() + num_index = prep.num_index + num_garbage = prep.num_garbage + + @cudaq.kernel + def run(): + index = cudaq.qvector(num_index) + garbage = cudaq.qvector(num_garbage) + prepare(index, garbage) + unprepare(index, garbage) + + state = np.array(cudaq.get_state(run)) + expected = np.zeros(1 << (num_index + num_garbage), dtype=np.complex128) + expected[0] = 1.0 + np.testing.assert_allclose(state, expected, atol=1e-10) + + +def test_alias_sampling_register_accounting(): + prep = AliasSamplingPrepare([1.0, 2.0, 3.0], mu=4) + assert prep.num_index == 2 + assert prep.num_bins == 4 + # [alias(m) | keep(mu) | keep_pad | ref(mu) | ref_pad | flag | + # ladder(qrom.num_ladder) | carry] = m + 2mu + 4 + num_ladder. + assert prep.num_garbage == 2 + 2 * 4 + 4 + prep.qrom.num_ladder + assert prep.ladder_offset == 2 + 2 * 4 + 3 + # A 4-entry table prices out to the plain select walk, whose ladder + # is one line per address bit — the layout then closes at the + # historical 2m + 2mu + 4. + assert prep.qrom.variant == "select" + assert prep.qrom.num_ladder == prep.num_index + assert prep.num_garbage == 2 * 2 + 2 * 4 + 4 + assert prep.lam == pytest.approx(6.0) + assert len(prep.keep) == 4 and len(prep.alias) == 4 + assert all(0 <= v < (1 << 4) for v in prep.keep) + assert all(0 <= a < 4 for a in prep.alias) + assert "mu=4" in repr(prep) + + +def test_alias_sampling_validation_raises(): + with pytest.raises(ValueError, match="weights must be non-empty"): + AliasSamplingPrepare([], mu=2) + with pytest.raises(ValueError, match="weights must be non-negative"): + AliasSamplingPrepare([1.0, -0.5], mu=2) + with pytest.raises(ValueError, match="weights must be finite"): + AliasSamplingPrepare([1.0, float("inf")], mu=2) + with pytest.raises(ValueError, match="sum to zero"): + AliasSamplingPrepare([0.0, 0.0], mu=2) + with pytest.raises(ValueError, match="mu must be a positive integer"): + AliasSamplingPrepare([1.0, 2.0], mu=0) + with pytest.raises(ValueError, match="mu must be a positive integer"): + AliasSamplingPrepare([1.0, 2.0], mu=1.5) + # bools are ints in Python but are rejected explicitly; integral + # numpy ints stay accepted. + with pytest.raises(ValueError, match="mu must be a positive integer"): + AliasSamplingPrepare([1.0, 2.0], mu=True) + assert AliasSamplingPrepare([1.0, 2.0], mu=np.int64(2)).mu == 2 + + +# ---------------------------------------------------------------------- +# Compiler-pinned resource contracts +# ---------------------------------------------------------------------- + +_RESOURCES = pytest.mark.skipif( + not hasattr(cudaq, "estimate_resources"), + reason="cudaq.estimate_resources is not available in this CUDA-Q") + + +def _prepare_resources(prep: AliasSamplingPrepare, adjoint: bool): + kernel = prep.adjoint_kernel() if adjoint else prep.kernel() + num_index = prep.num_index + num_garbage = prep.num_garbage + + @cudaq.kernel + def harness(): + index = cudaq.qvector(num_index) + garbage = cudaq.qvector(num_garbage) + kernel(index, garbage) + + return cudaq.estimate_resources(harness) + + +@_RESOURCES +@pytest.mark.parametrize("weights,mu", [([0.7, 0.2, 1.4], 2), + ([1.0, 3.5, 0.25, 2.0, 0.8], 3), + ([2.5], 3), ([1.0] * 8, 2)]) +@pytest.mark.parametrize("adjoint", [False, True]) +def test_alias_sampling_cost_is_qrom_price_plus_comparator( + weights, mu, adjoint): + # The lookup cost is not re-derived here: it is asserted consistent + # with the QROM's own reported count (whatever construction "auto" + # priced in), on top of which sit exactly the two CDKM adders of the + # comparator — 2 (mu + 1) Toffolis each on the (mu+1)-bit extension. + # The alias swap is num_index Fredkins (counted as cswap, not ccx — + # ``count_controls`` is arity-aware; ``count("ccx")`` matches + # nothing). + prep = AliasSamplingPrepare(weights, mu) + resources = _prepare_resources(prep, adjoint) + # Not a tautology: qrom.toffoli_count is classical bookkeeping in + # QROM, deliberately cross-pinned here against the compiled circuit. + assert resources.count_controls("x", 2) == \ + prep.qrom.toffoli_count + 4 * (mu + 1) + assert resources.count_controls("swap", 1) == prep.num_index + # The mu reference Hadamards and the num_index bin Hadamards. (An + # incidental property of the current QROM gate choices — the lookup + # happens to add no h gates — not a contract.) + assert resources.count("h") == prep.num_index + prep.mu + + +@_RESOURCES +def test_alias_sampling_forced_select_swap_variant(): + # Pins the docstring's variant-agnostic claim: force the lookup onto + # the select_swap construction through the new passthrough at the + # smallest simulable size (m = 2, mu = 2, B = 2: the ladder is + # 1 + 2 * (m + mu) = 9 lines, 21 qubits total) and re-assert the + # full PREPARE contract. + weights = [0.7, 0.2, 1.4, 0.5] + prep = AliasSamplingPrepare(weights, + mu=2, + variant="select_swap", + block_size=2) + assert prep.qrom.variant == "select_swap" + assert prep.num_index + prep.num_garbage <= 24 + # (i) Exact marginal amplitudes, against the independent enumeration. + brute = _brute_force_marginal(prep) + np.testing.assert_allclose(_circuit_marginal(prep), brute, atol=1e-10) + np.testing.assert_allclose(prep.table_probabilities, brute, atol=1e-15) + # (ii) PREPARE then the hand-written adjoint is the identity. + prepare = prep.kernel() + unprepare = prep.adjoint_kernel() + num_index = prep.num_index + num_garbage = prep.num_garbage + + @cudaq.kernel + def roundtrip(): + index = cudaq.qvector(num_index) + garbage = cudaq.qvector(num_garbage) + prepare(index, garbage) + unprepare(index, garbage) + + state = np.array(cudaq.get_state(roundtrip)) + expected = np.zeros(1 << (num_index + num_garbage), dtype=np.complex128) + expected[0] = 1.0 + np.testing.assert_allclose(state, expected, atol=1e-10) + # (iii) The Toffoli cost identity holds for this variant too. + for adjoint in (False, True): + resources = _prepare_resources(prep, adjoint) + assert resources.count_controls("x", 2) == \ + prep.qrom.toffoli_count + 4 * (prep.mu + 1) + + +def test_alias_sampling_passthrough_validation_propagates(): + # QROM's own validation errors surface unchanged. + with pytest.raises(ValueError, match="block_size"): + AliasSamplingPrepare([1.0, 2.0, 3.0], mu=2, block_size=2) + + +@_RESOURCES +def test_alias_sampling_adjoint_costs_exactly_the_forward_price(): + # The hand-written adjoint is the literal gate-reversal: same gate + # multiset, gate for gate. + prep = AliasSamplingPrepare([0.9, 0.0, 2.2, 0.0, 0.4, 1.3, 0.0, 3.1], 2) + forward = _prepare_resources(prep, adjoint=False) + backward = _prepare_resources(prep, adjoint=True) + for name, controls in (("x", 2), ("x", 1), ("swap", 1), ("x", 0)): + assert forward.count_controls(name, controls) == \ + backward.count_controls(name, controls) + assert forward.count("h") == backward.count("h") From 9c733627a5702dff348e236340a22992d52b9c6d Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Wed, 9 Sep 2026 16:17:51 +0000 Subject: [PATCH 4/6] Quote ket and absolute-value notation for the Sphinx build Bare |0...0> kets and the |table - ideal| absolute value in the AliasSamplingPrepare docstrings parse as unterminated RST substitution references and fail the warnings-as-errors docs build now that the class is on the API reference page. Same fix as the unary-iteration module received on the QROM PR. Signed-off-by: Scott Thornton --- .../primitives/_alias_sampling.py | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/python/cudaq_algorithms/primitives/_alias_sampling.py b/python/cudaq_algorithms/primitives/_alias_sampling.py index 193512b..5c5b3ef 100644 --- a/python/cudaq_algorithms/primitives/_alias_sampling.py +++ b/python/cudaq_algorithms/primitives/_alias_sampling.py @@ -54,9 +54,9 @@ ``2 num_index + 2 mu + 4`` whenever the priced lookup is the plain ``"select"`` walk (``num_ladder = num_index``; always the case for the small tables where select wins ``"auto"``'s pricing). Both registers -must be |0...0> on entry to ``kernel()``; ``adjoint_kernel()`` is the +must be ``|0...0>`` on entry to ``kernel()``; ``adjoint_kernel()`` is the literal gate-reversal (no ``cudaq.adjoint``, cuda-quantum#4897/#4898) and -returns them to |0...0>. The QROM lookup uncomputes itself on the clean +returns them to ``|0...0>``. The QROM lookup uncomputes itself on the clean ladder sector the adjoint provides, whatever the variant (the ladder is restored by the compute pass and untouched in between, so re-applying the lookup XORs the table back out). These kernels call the QROM and @@ -172,9 +172,9 @@ class AliasSamplingPrepare: mutually orthogonal (measured off-diagonal overlaps reach ~0.19), so the reduced state on the index register is not ``diag(p_k)``; only the amplitude *magnitudes* on the index register are guaranteed. - That is exactly what qubitization's ``<0| PREP^dagger SELECT PREP - |0>`` contract needs — index orthogonality kills the cross terms — - but consumers must not assume orthonormal garbage. + That is exactly what qubitization's ``<0|PREP^dagger SELECT PREP|0>`` + contract needs — index orthogonality kills the cross terms — but + consumers must not assume orthonormal garbage. """ def __init__(self, @@ -334,7 +334,7 @@ def ladder_offset(self) -> int: """Offset of the ``qrom.num_ladder``-wide QROM ladder inside garbage. - The ladder qubits are clean (|0>) between ``kernel()`` and + The ladder qubits are clean (``|0>``) between ``kernel()`` and ``adjoint_kernel()`` — unlike the rest of the garbage — so a SELECT sandwiched between them may reuse ``garbage[ladder_offset : ladder_offset + qrom.num_ladder]`` as @@ -363,7 +363,7 @@ def table_probabilities(self) -> np.ndarray: @property def discretization_bound(self) -> float: - """Per-bin bound on |table - ideal| probability (see module doc).""" + """Per-bin bound on ``|table - ideal|`` probability (see module doc).""" return 1.5 / (self._num_bins * (1 << self._mu)) def __repr__(self) -> str: @@ -377,7 +377,7 @@ def __repr__(self) -> str: # ------------------------------------------------------------------ def kernel(self): - """PREPARE ``(index, garbage)`` — both |0...0> on entry.""" + """PREPARE ``(index, garbage)`` — both ``|0...0>`` on entry.""" return self._prepare def adjoint_kernel(self): From 334e1608e0cb525e6dc2d2e517a056700d662923 Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Fri, 11 Sep 2026 22:14:20 +0000 Subject: [PATCH 5/6] Consume the family register comparator in the alias-sampling PREPARE Implements the reviewer-suggested consolidation: the PREPARE's keep-vs-reference conditional now calls the arithmetic family's cmp_ge_register (flag <- (ref >= keep) on the bare mu-bit operands) instead of the bespoke two-adder weave (subtract keep on a (mu+1)-bit extension, copy the borrow, negate, add back). - Toffoli count per PREPARE side drops from 4(mu+1) to 2mu: the cost identity is now qrom.toffoli_count + 2 mu Toffolis + num_index Fredkins, re-derived in the docstrings and re-pinned by the resource tests against cudaq.estimate_resources. - The (mu+1)-bit extension's two pad qubits (keep_pad, ref_pad) are gone: cmp_ge_register complements its second operand in place and needs only the one-qubit carry the layout already had. num_garbage shrinks by 2 to num_index + 2 mu + 2 + qrom.num_ladder, and ladder_offset moves down accordingly. - The hand-written adjoint stays the literal gate-reversal: cmp_ge_register is palindromic around its self-inverse flag copy, so re-applying it is the reversed gate sequence. All marginal/bound/identity tests pass unchanged; only the layout and cost accounting assertions are re-derived. Signed-off-by: Scott Thornton --- .../primitives/_alias_sampling.py | 60 +++++++++---------- .../python/test_primitives_alias_sampling.py | 36 +++++------ 2 files changed, 46 insertions(+), 50 deletions(-) diff --git a/python/cudaq_algorithms/primitives/_alias_sampling.py b/python/cudaq_algorithms/primitives/_alias_sampling.py index 5c5b3ef..92ea716 100644 --- a/python/cudaq_algorithms/primitives/_alias_sampling.py +++ b/python/cudaq_algorithms/primitives/_alias_sampling.py @@ -16,9 +16,10 @@ priced across the available QROM variants and the cheapest construction is minted — its clean-ancilla ladder lives inside the garbage register, so ``num_garbage`` follows the chosen variant), -3. a comparator of ``keep_k`` against a uniform ``mu``-bit reference - register (built from the CDKM register adders of ``_arithmetic``: - subtract on a ``mu+1``-bit extension, copy the borrow, add back), and +3. a comparator of a uniform ``mu``-bit reference register against + ``keep_k`` (the CDKM register comparator + :func:`~cudaq_algorithms.primitives.cmp_ge_register` of + ``_arithmetic``: ``flag <- (ref >= keep_k)``), and 4. a controlled swap of the bin index with the alias when the reference is not below ``keep_k``. @@ -48,10 +49,10 @@ Kernel signatures (little-endian; ``docs/conventions.md``): both ``kernel()`` and ``adjoint_kernel()`` are ``(index: qview, garbage: qview)`` with ``index`` of width ``num_index`` and ``garbage`` of width -``num_garbage`` laid out as ``[alias(num_index) | keep(mu) | keep_pad | -ref(mu) | ref_pad | flag | ladder(qrom.num_ladder) | carry]`` — i.e. -``num_garbage = num_index + 2 mu + 3 + qrom.num_ladder + 1``, which is -``2 num_index + 2 mu + 4`` whenever the priced lookup is the plain +``num_garbage`` laid out as ``[alias(num_index) | keep(mu) | ref(mu) | +flag | ladder(qrom.num_ladder) | carry]`` — i.e. +``num_garbage = num_index + 2 mu + 1 + qrom.num_ladder + 1``, which is +``2 num_index + 2 mu + 2`` whenever the priced lookup is the plain ``"select"`` walk (``num_ladder = num_index``; always the case for the small tables where select wins ``"auto"``'s pricing). Both registers must be ``|0...0>`` on entry to ``kernel()``; ``adjoint_kernel()`` is the @@ -67,10 +68,10 @@ Cost (pinned by ``tests/python/test_primitives_alias_sampling.py`` against the compiler): each of ``kernel()`` / ``adjoint_kernel()`` costs -exactly ``qrom.toffoli_count + 4 (mu + 1)`` Toffolis — the lookup at the -QROM's own reported price plus two CDKM register adders on the -``mu+1``-bit extension at ``2 (mu + 1)`` each — plus ``num_index`` -controlled swaps (Fredkins) for the alias swap. +exactly ``qrom.toffoli_count + 2 mu`` Toffolis — the lookup at the +QROM's own reported price plus one CDKM register comparator on the +``mu``-bit operands at ``2 mu`` — plus ``num_index`` controlled swaps +(Fredkins) for the alias swap. """ from __future__ import annotations @@ -81,7 +82,7 @@ import cudaq import numpy as np -from ._arithmetic import add_register, subtract_register +from ._arithmetic import cmp_ge_register from ._qrom import QROM from ._unary_iteration import _retain @@ -236,10 +237,10 @@ def _build_kernels(self) -> None: m = self._num_index mu = self._mu num_ladder = self._qrom.num_ladder - k0 = m # keep (keep_pad at k0 + mu) - r0 = m + mu + 1 # ref (ref_pad at r0 + mu) - flag = m + 2 * mu + 2 - l0 = m + 2 * mu + 3 + k0 = m # keep + r0 = m + mu # ref + flag = m + 2 * mu + l0 = m + 2 * mu + 1 c0 = l0 + num_ladder # carry qrom_kernel = self._qrom.kernel @@ -252,14 +253,10 @@ def primitives_alias_prepare(index: cudaq.qview, garbage: cudaq.qview): h(garbage[r0 + b]) # (alias_k, keep_k) lookup: output is garbage[0 : m + mu]. qrom_kernel(index, garbage[l0:l0 + num_ladder], garbage[0:m + mu]) - # flag <- NOT (ref < keep): subtract keep on the (mu+1)-bit - # extension of ref, copy the borrow (MSB), add keep back. - subtract_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], - garbage[c0:c0 + 1]) - cx(garbage[r0 + mu], garbage[flag]) - x(garbage[flag]) - add_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], - garbage[c0:c0 + 1]) + # flag <- (ref >= keep): the family register comparator + # (keep and ref are left untouched, the carry returns to |0>). + cmp_ge_register(garbage[r0:r0 + mu], garbage[k0:k0 + mu], + garbage[c0:c0 + 1], garbage[flag:flag + 1]) # Swap the bin index with its alias on the flag. for b in range(m): swap.ctrl(garbage[flag], index[b], garbage[b]) @@ -271,12 +268,11 @@ def primitives_alias_prepare_adj(index: cudaq.qview, for j in range(m): b = m - 1 - j swap.ctrl(garbage[flag], index[b], garbage[b]) - subtract_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], - garbage[c0:c0 + 1]) - x(garbage[flag]) - cx(garbage[r0 + mu], garbage[flag]) - add_register(garbage[k0:k0 + mu + 1], garbage[r0:r0 + mu + 1], - garbage[c0:c0 + 1]) + # The comparator is its own gate-reversal (a palindrome + # around the self-inverse flag copy), so re-applying it here + # IS the literal reversed gate sequence. + cmp_ge_register(garbage[r0:r0 + mu], garbage[k0:k0 + mu], + garbage[c0:c0 + 1], garbage[flag:flag + 1]) # The QROM lookup uncomputes itself on the clean-ladder # sector: the ladder is |0> here (restored by the compute # pass, untouched since), so re-applying the lookup XORs the @@ -320,7 +316,7 @@ def num_bins(self) -> int: @property def num_garbage(self) -> int: """Garbage register width (layout in the module docstring).""" - return self._num_index + 2 * self._mu + 4 + self._qrom.num_ladder + return self._num_index + 2 * self._mu + 2 + self._qrom.num_ladder @property def qrom(self) -> QROM: @@ -340,7 +336,7 @@ def ladder_offset(self) -> int: ``garbage[ladder_offset : ladder_offset + qrom.num_ladder]`` as its own ladder. """ - return self._num_index + 2 * self._mu + 3 + return self._num_index + 2 * self._mu + 1 @property def keep(self) -> tuple[int, ...]: diff --git a/tests/python/test_primitives_alias_sampling.py b/tests/python/test_primitives_alias_sampling.py index f98b686..6c40f05 100644 --- a/tests/python/test_primitives_alias_sampling.py +++ b/tests/python/test_primitives_alias_sampling.py @@ -18,10 +18,10 @@ The resource tests at the bottom hold the documented cost against the compiler: each PREPARE (and its adjoint) costs exactly -``qrom.toffoli_count + 4 (mu + 1)`` Toffolis plus ``num_index`` Fredkins +``qrom.toffoli_count + 2 mu`` Toffolis plus ``num_index`` Fredkins — the lookup at the QROM's *own reported* price (whatever variant -``"auto"`` minted) plus two CDKM register adders on the ``mu+1``-bit -comparator extension at ``2 (mu + 1)`` Toffolis each. +``"auto"`` minted) plus one CDKM register comparator +(``cmp_ge_register``) on the ``mu``-bit operands at ``2 mu`` Toffolis. """ import numpy as np @@ -152,16 +152,16 @@ def test_alias_sampling_register_accounting(): prep = AliasSamplingPrepare([1.0, 2.0, 3.0], mu=4) assert prep.num_index == 2 assert prep.num_bins == 4 - # [alias(m) | keep(mu) | keep_pad | ref(mu) | ref_pad | flag | - # ladder(qrom.num_ladder) | carry] = m + 2mu + 4 + num_ladder. - assert prep.num_garbage == 2 + 2 * 4 + 4 + prep.qrom.num_ladder - assert prep.ladder_offset == 2 + 2 * 4 + 3 + # [alias(m) | keep(mu) | ref(mu) | flag | + # ladder(qrom.num_ladder) | carry] = m + 2mu + 2 + num_ladder. + assert prep.num_garbage == 2 + 2 * 4 + 2 + prep.qrom.num_ladder + assert prep.ladder_offset == 2 + 2 * 4 + 1 # A 4-entry table prices out to the plain select walk, whose ladder - # is one line per address bit — the layout then closes at the - # historical 2m + 2mu + 4. + # is one line per address bit — the layout then closes at + # 2m + 2mu + 2. assert prep.qrom.variant == "select" assert prep.qrom.num_ladder == prep.num_index - assert prep.num_garbage == 2 * 2 + 2 * 4 + 4 + assert prep.num_garbage == 2 * 2 + 2 * 4 + 2 assert prep.lam == pytest.approx(6.0) assert len(prep.keep) == 4 and len(prep.alias) == 4 assert all(0 <= v < (1 << 4) for v in prep.keep) @@ -221,17 +221,17 @@ def test_alias_sampling_cost_is_qrom_price_plus_comparator( weights, mu, adjoint): # The lookup cost is not re-derived here: it is asserted consistent # with the QROM's own reported count (whatever construction "auto" - # priced in), on top of which sit exactly the two CDKM adders of the - # comparator — 2 (mu + 1) Toffolis each on the (mu+1)-bit extension. - # The alias swap is num_index Fredkins (counted as cswap, not ccx — - # ``count_controls`` is arity-aware; ``count("ccx")`` matches - # nothing). + # priced in), on top of which sits exactly the one CDKM register + # comparator (cmp_ge_register) — 2 mu Toffolis on the mu-bit + # operands. The alias swap is num_index Fredkins (counted as cswap, + # not ccx — ``count_controls`` is arity-aware; ``count("ccx")`` + # matches nothing). prep = AliasSamplingPrepare(weights, mu) resources = _prepare_resources(prep, adjoint) # Not a tautology: qrom.toffoli_count is classical bookkeeping in # QROM, deliberately cross-pinned here against the compiled circuit. assert resources.count_controls("x", 2) == \ - prep.qrom.toffoli_count + 4 * (mu + 1) + prep.qrom.toffoli_count + 2 * mu assert resources.count_controls("swap", 1) == prep.num_index # The mu reference Hadamards and the num_index bin Hadamards. (An # incidental property of the current QROM gate choices — the lookup @@ -244,7 +244,7 @@ def test_alias_sampling_forced_select_swap_variant(): # Pins the docstring's variant-agnostic claim: force the lookup onto # the select_swap construction through the new passthrough at the # smallest simulable size (m = 2, mu = 2, B = 2: the ladder is - # 1 + 2 * (m + mu) = 9 lines, 21 qubits total) and re-assert the + # 1 + 2 * (m + mu) = 9 lines, 19 qubits total) and re-assert the # full PREPARE contract. weights = [0.7, 0.2, 1.4, 0.5] prep = AliasSamplingPrepare(weights, @@ -278,7 +278,7 @@ def roundtrip(): for adjoint in (False, True): resources = _prepare_resources(prep, adjoint) assert resources.count_controls("x", 2) == \ - prep.qrom.toffoli_count + 4 * (prep.mu + 1) + prep.qrom.toffoli_count + 2 * prep.mu def test_alias_sampling_passthrough_validation_propagates(): From 1cc210b93888f9096f78f820e54c3114418b04c4 Mon Sep 17 00:00:00 2001 From: Scott Thornton Date: Thu, 17 Sep 2026 04:38:19 +0000 Subject: [PATCH 6/6] Address review: docstring, packing comment, offset names, inverse-QROM TODO Per smcardleQ's review of PR #48: - drop the instance-specific '~0.19' overlap figure from the garbage caveat - spell out the QROM-word packing convention (alias low bits, keep above) - rename terse offsets k0/r0/l0/c0 -> keep_off/ref_off/ladder_off/carry_off - add a TODO to expose the inverse QROM as an explicit qrom_inverse call, since measurement-based uncomputation makes the adjoint PREPARE cheaper No behavior change; alias-sampling suite 26/26. Signed-off-by: Scott Thornton --- .../primitives/_alias_sampling.py | 51 ++++++++++++------- 1 file changed, 34 insertions(+), 17 deletions(-) diff --git a/python/cudaq_algorithms/primitives/_alias_sampling.py b/python/cudaq_algorithms/primitives/_alias_sampling.py index 92ea716..49580f1 100644 --- a/python/cudaq_algorithms/primitives/_alias_sampling.py +++ b/python/cudaq_algorithms/primitives/_alias_sampling.py @@ -170,9 +170,9 @@ class AliasSamplingPrepare: QROM's own ``ValueError``. Garbage caveat: the per-bin garbage vectors ``|g_k>`` are NOT - mutually orthogonal (measured off-diagonal overlaps reach ~0.19), so - the reduced state on the index register is not ``diag(p_k)``; only - the amplitude *magnitudes* on the index register are guaranteed. + mutually orthogonal, so the reduced state on the index register is + not ``diag(p_k)``; only the amplitude *magnitudes* on the index + register are guaranteed. That is exactly what qubitization's ``<0|PREP^dagger SELECT PREP|0>`` contract needs — index orthogonality kills the cross terms — but consumers must not assume orthonormal garbage. @@ -214,7 +214,12 @@ def __init__(self, self._probabilities = np.asarray(probabilities) self._table_probabilities = np.asarray(table) - # QROM data: alias in the low num_index bits, keep above it. + # QROM data: each address k stores one (num_index + mu)-bit word + # packing both lookup outputs into a single QROM entry -- + # bits [0, num_index) = alias[k] (the alias bin index) + # bits [num_index, num_index+mu) = keep[k] (the keep threshold) + # so one lookup emits alias into garbage's alias register and keep + # into its keep register (see the garbage-layout docstring). # variant/block_size pass through verbatim (default "auto" # prices the constructions); the chosen variant's ladder width # shapes the garbage layout below. QROM validation errors @@ -237,11 +242,13 @@ def _build_kernels(self) -> None: m = self._num_index mu = self._mu num_ladder = self._qrom.num_ladder - k0 = m # keep - r0 = m + mu # ref - flag = m + 2 * mu - l0 = m + 2 * mu + 1 - c0 = l0 + num_ladder # carry + # Base offsets into the garbage register, matching the layout + # [alias(m) | keep(mu) | ref(mu) | flag | ladder | carry]. + keep_off = m # keep register (mu bits) + ref_off = m + mu # reference register (mu bits) + flag = m + 2 * mu # comparator flag (1 qubit) + ladder_off = m + 2 * mu + 1 # QROM ladder ancillas + carry_off = ladder_off + num_ladder # comparator carry (1 qubit) qrom_kernel = self._qrom.kernel @cudaq.kernel @@ -250,13 +257,16 @@ def primitives_alias_prepare(index: cudaq.qview, garbage: cudaq.qview): for b in range(m): h(index[b]) for b in range(mu): - h(garbage[r0 + b]) + h(garbage[ref_off + b]) # (alias_k, keep_k) lookup: output is garbage[0 : m + mu]. - qrom_kernel(index, garbage[l0:l0 + num_ladder], garbage[0:m + mu]) + qrom_kernel(index, garbage[ladder_off:ladder_off + num_ladder], + garbage[0:m + mu]) # flag <- (ref >= keep): the family register comparator # (keep and ref are left untouched, the carry returns to |0>). - cmp_ge_register(garbage[r0:r0 + mu], garbage[k0:k0 + mu], - garbage[c0:c0 + 1], garbage[flag:flag + 1]) + cmp_ge_register(garbage[ref_off:ref_off + mu], + garbage[keep_off:keep_off + mu], + garbage[carry_off:carry_off + 1], + garbage[flag:flag + 1]) # Swap the bin index with its alias on the flag. for b in range(m): swap.ctrl(garbage[flag], index[b], garbage[b]) @@ -271,17 +281,24 @@ def primitives_alias_prepare_adj(index: cudaq.qview, # The comparator is its own gate-reversal (a palindrome # around the self-inverse flag copy), so re-applying it here # IS the literal reversed gate sequence. - cmp_ge_register(garbage[r0:r0 + mu], garbage[k0:k0 + mu], - garbage[c0:c0 + 1], garbage[flag:flag + 1]) + cmp_ge_register(garbage[ref_off:ref_off + mu], + garbage[keep_off:keep_off + mu], + garbage[carry_off:carry_off + 1], + garbage[flag:flag + 1]) # The QROM lookup uncomputes itself on the clean-ladder # sector: the ladder is |0> here (restored by the compute # pass, untouched since), so re-applying the lookup XORs the # table back out of garbage[0 : m + mu] — the documented # self-inverse contract of every QROM variant. - qrom_kernel(index, garbage[l0:l0 + num_ladder], garbage[0:m + mu]) + # TODO: expose this as an explicit ``qrom_inverse`` call. Once + # measurement-based uncomputation lands, uncomputing the QROM + # is much cheaper than the forward lookup, making the adjoint + # PREPARE significantly cheaper than the forward pass. + qrom_kernel(index, garbage[ladder_off:ladder_off + num_ladder], + garbage[0:m + mu]) for j in range(mu): b = mu - 1 - j - h(garbage[r0 + b]) + h(garbage[ref_off + b]) for j in range(m): b = m - 1 - j h(index[b])