subset_sum_oracle.py

from qiskit import QuantumCircuit, QuantumRegister, AncillaRegister
import numpy as np


def qft(qc: QuantumCircuit, reg: QuantumRegister):
    """Apply Quantum Fourier Transform to register."""
    n = len(reg)
    for i in range(n - 1, -1, -1):
        qc.h(reg[i])
        for j in range(i - 1, -1, -1):
            angle = np.pi / (2 ** (i - j))
            qc.cp(angle, reg[j], reg[i])


def inverse_qft(qc: QuantumCircuit, reg: QuantumRegister):
    """Apply inverse Quantum Fourier Transform to register."""
    n = len(reg)
    for i in range(n):
        for j in range(i):
            angle = -np.pi / (2 ** (i - j))
            qc.cp(angle, reg[j], reg[i])
        qc.h(reg[i])


class SubsetSumOracle:
    """
    Translates the Subset Sum problem into a quantum oracle for Grover's algorithm.
    Uses a simple comparator-based approach.
    """

    def __init__(self, numbers: list[int], target: int):
        self.numbers = numbers
        self.target = target
        self.n_items = len(numbers)

        positive_max_sum = sum(n for n in numbers if n > 0)
        negative_min_sum = sum(n for n in numbers if n < 0)

        if target > 0 and target > positive_max_sum:
            raise ValueError(f"Target ({self.target}) must be less or equal to the sum of all positive numbers ({positive_max_sum})")
        elif target < 0 and target < negative_min_sum:
            raise ValueError(f"Target ({self.target}) must be greater or equal to the sum of all negative numbers ({negative_min_sum})")

        self.abs_max_sum = max(positive_max_sum, -negative_min_sum)
        self.sum_bits = max(1, int(np.ceil(np.log2(self.abs_max_sum + 1))))

        # add one (sign) bit to be able to represent negative sums
        if negative_min_sum < 0:
            self.sum_bits += 1

    def build_oracle(self, arithmetic_mode) -> QuantumCircuit:
        """
        Builds the oracle.
        For each possible subset, computes the sum and checks if it equals the target.

        Each possible subset comes from Hadamard gates on the item register before this oracle (superposition):

               ┌───┐     ┌──────┐
        |0> ≡≡≡┤ H ├≡≡≡≡≡┤      ├  item register (in superposition of all possible subsets)
               └───┘     │      │
        |0> ≡≡≡≡≡≡≡≡≡≡≡≡≡┤  U_f ├  sum register (all qubits are 0)
               ┌───┐     │      │
        |1> ───┤ H ├─────┤      ├  ancilla qubit
               └───┘     └──────┘

        :arg arithmetic_mode: 'qft' or 'ripple' for the arithmetic operations
        """
        if arithmetic_mode not in ['qft', 'ripple']: raise ValueError(f"Unknown arithmetic mode: {arithmetic_mode}")

        items = QuantumRegister(self.n_items, 'items')
        sum_reg = QuantumRegister(self.sum_bits, 'sum')
        ancilla = AncillaRegister(1, 'anc')

        qc = QuantumCircuit(items, sum_reg, ancilla, name='SubsetSumOracle')

        # compute the sum: for each item, add its value to sum_reg if selected
        for i, value in enumerate(self.numbers):
            if arithmetic_mode == 'qft': self._controlled_add_qft(qc, items[i], sum_reg, value)
            elif arithmetic_mode == 'ripple': self._controlled_add_ripple(qc, items[i], sum_reg, value)
        if arithmetic_mode == 'qft': inverse_qft(qc, sum_reg)

        # check if sum == target and flip phase
        self._phase_flip_if_equal(qc, sum_reg, ancilla[0], self.target)

        # uncompute sum (subtract in reverse order)
        if arithmetic_mode == 'qft': qft(qc, sum_reg)
        for i in range(self.n_items - 1, -1, -1):
            if arithmetic_mode == 'qft': self._controlled_subtract_qft(qc, items[i], sum_reg, self.numbers[i])
            elif arithmetic_mode == 'ripple': self._controlled_subtract_ripple(qc, items[i], sum_reg, self.numbers[i])

        return qc

    def _controlled_add_ripple(self, qc: QuantumCircuit, ctrl, sum_reg: QuantumRegister, value: int):
        """Add value to sum_reg, controlled by ctrl qubit. Uses simple ripple approach."""
        if value == 0:
            return

        n = len(sum_reg)

        # for each bit position, determine if we need to flip based on value and carries
        for bit_pos in range(n):
            bit_val = (value >> bit_pos) & 1
            if bit_val == 1:
                # need to add 1 at this position, which may cause carries
                # flip this bit and propagate carry
                # the carry propagation: if we flip bit i and it was 1, we need to flip bit i+1, etc.

                # start from MSB and work down to handle carries correctly
                # use controlled operations with multiple controls for carry chain
                for j in range(n - 1, bit_pos, -1):
                    # flip bit j if ctrl is 1 AND all bits from bit_pos to j-1 are 1
                    controls = [ctrl] + [sum_reg[k] for k in range(bit_pos, j)]
                    qc.mcx(controls, sum_reg[j])

                # finally, flip the bit at bit_pos
                qc.cx(ctrl, sum_reg[bit_pos])

    def _controlled_subtract_ripple(self, qc: QuantumCircuit, ctrl, sum_reg: QuantumRegister, value: int):
        """Subtract value from sum_reg, controlled by ctrl qubit. Inverse of add."""
        if value == 0:
            return

        n = len(sum_reg)

        for bit_pos in range(n):
            bit_val = (value >> bit_pos) & 1
            if bit_val == 1:
                # inverse of addition: first flip the bit, then handle carries in reverse order
                qc.cx(ctrl, sum_reg[bit_pos])

                for j in range(bit_pos + 1, n):
                    controls = [ctrl] + [sum_reg[k] for k in range(bit_pos, j)]
                    qc.mcx(controls, sum_reg[j])

    def _controlled_add_qft(self, qc: QuantumCircuit, ctrl, sum_reg: QuantumRegister, value: int):
        """Add value to sum_reg, controlled by ctrl qubit. Uses Quantum Fourier Transform for addition."""
        if value == 0:
            return

        n = len(sum_reg)

        # apply controlled phase rotations to add 'value'
        # in QFT basis, adding a classical value 'a' means applying phase e^(2πi * a * k / 2^n)
        # to each basis state |k⟩
        for i in range(n):
            # calculate the cumulative phase for qubit i
            # qubit i (from LSB) contributes to phases based on bits 0..i of value
            angle = 0
            for j in range(i + 1):
                if (value >> j) & 1:
                    angle += np.pi / (2 ** (i - j))

            if angle != 0:
                qc.cp(angle, ctrl, sum_reg[i])

    def _controlled_subtract_qft(self, qc: QuantumCircuit, ctrl, sum_reg: QuantumRegister, value: int):
        """Subtract value from sum_reg using QFT. Inverse of _controlled_add_qft."""
        if value == 0:
            return

        n = len(sum_reg)

        for i in range(n):
            angle = 0
            for j in range(i + 1):
                if (value >> j) & 1:
                    angle -= np.pi / (2 ** (i - j))  # negative for subtraction

            if angle != 0:
                qc.cp(angle, ctrl, sum_reg[i])

    def _phase_flip_if_equal(self, qc: QuantumCircuit, sum_reg: QuantumRegister,
                             ancilla, target: int):
        """Apply phase flip if sum_reg equals target."""
        n = len(sum_reg)

        # flip qubits of sum_reg where the target bit is 0 (so all-ones means match)
        for i in range(n):
            if target & (1 << i) == 0:
                qc.x(sum_reg[i])

        # multi-controlled Z on ancilla in |-> state
        # the ancilla should already be in |-> from the main circuit
        qc.mcx(list(sum_reg), ancilla)

        # unflip sum_reg qubits that were just flipped
        for i in range(n):
            if target & (1 << i) == 0:
                qc.x(sum_reg[i])