#  Copyright 2023 Google LLC
#
#  Licensed under the Apache License, Version 2.0 (the "License");
#  you may not use this file except in compliance with the License.
#  You may obtain a copy of the License at
#
#      https://www.apache.org/licenses/LICENSE-2.0
#
#  Unless required by applicable law or agreed to in writing, software
#  distributed under the License is distributed on an "AS IS" BASIS,
#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#  See the License for the specific language governing permissions and
#  limitations under the License.

Unary Iteration

Given an array of potential operations, for example:

ops = [X(i) for i in range(5)]

we would like to select an operation to apply:

n = 4 --> apply ops[4]

If \(n\) is a quantum integer, we need to apply the transformation

\[ |n \rangle |\psi\rangle \rightarrow |n\rangle \, \mathrm{ops}_n \cdot |\psi\rangle \]

The simplest conceptual way to do this is to use a “total control” quantum circuit where you introduce a multi-controlled operation for each of the len(ops) possible n values.

import cirq
from cirq.contrib.svg import SVGCircuit
import numpy as np
from typing import *

import operator
import cirq._compat
import itertools

Total Control

Here, we’ll use Sympy’s boolean logic to show how total control works. We perform an And( ... ) for each possible bit pattern. We use an Xnor on each selection bit to toggle whether it’s a positive or negative control (filled or open circle in quantum circuit diagrams).

In this example, we indeed consider \(X_n\) as our potential operations and toggle bits in the target register according to the total control.

import sympy as S
import sympy.logic.boolalg as slb

def total_control(selection, target):
    """Toggle bits in `target` depending on `selection`."""
    print(f"Selection is {selection}")
    
    for n, trial in enumerate(itertools.product((0, 1), repeat=len(selection))):
        print(f"Step {n}, apply total control: {trial}")
        target[n] ^= slb.And(*[slb.Xnor(s, t) for s, t in zip(selection, trial)])
          
        if target[n] == S.true:
            print(f"  -> At this stage, {n}= and our output bit is set")

        
selection = [0, 0, 0]
target = [False]*8
total_control(selection, target)    
print()
print("Target:")
print(target)

Selection is [0, 0, 0]
Step 0, apply total control: (0, 0, 0)
  -> At this stage, 0= and our output bit is set
Step 1, apply total control: (0, 0, 1)
Step 2, apply total control: (0, 1, 0)
Step 3, apply total control: (0, 1, 1)
Step 4, apply total control: (1, 0, 0)
Step 5, apply total control: (1, 0, 1)
Step 6, apply total control: (1, 1, 0)
Step 7, apply total control: (1, 1, 1)

Target:
[True, False, False, False, False, False, False, False]

Note that our target register shows we have indeed applied \(X_\mathrm{0b010}\). Try changing selection to other bit patterns and notice how it changes.

Of course, we don’t know what state the selection register will be in. We can use sympy’s support for symbolic boolean logic to verify our gadget for all possible selection inputs.

selection = [S.Symbol(f's{i}') for i in range(3)]
target = [S.false for i in range(2**len(selection)) ]
total_control(selection, target)

print()
print("Target:")
for n, t in enumerate(target):
    print(f'{n}= {t}')
    
tc_target = target.copy()

Selection is [s0, s1, s2]
Step 0, apply total control: (0, 0, 0)
Step 1, apply total control: (0, 0, 1)
Step 2, apply total control: (0, 1, 0)
Step 3, apply total control: (0, 1, 1)
Step 4, apply total control: (1, 0, 0)
Step 5, apply total control: (1, 0, 1)
Step 6, apply total control: (1, 1, 0)
Step 7, apply total control: (1, 1, 1)

Target:
0= ~s0 & ~s1 & ~s2
1= s2 & ~s0 & ~s1
2= s1 & ~s0 & ~s2
3= s1 & s2 & ~s0
4= s0 & ~s1 & ~s2
5= s0 & s2 & ~s1
6= s0 & s1 & ~s2
7= s0 & s1 & s2

As expected, the “not pattern” (where ~ is boolean not) matches the binary representations of n.

Unary Iteration with segment trees

A segment tree is a data structure that allows logrithmic-time querying of intervals. We use a segment tree where each interval is length 1 and comprises all the n integers we may select.

It is defined recursively by dividing the input interval into two half-size intervals until the left limit meets the right limit.

def segtree(ctrl, selection, target, depth, left, right):
    """Toggle bits in `target` depending on `selection` using a recursive segment tree."""
    print(f'depth={depth} left={left} right={right}', end=' ')
    
    if left == (right - 1):
        # Leaf of the recusion.
        print(f'n={n} ctrl={ctrl}')
        target[left] ^= ctrl
        return 
    print()
    
    assert depth < len(selection)
    mid = (left + right) >> 1
    
    # Recurse left interval
    new_ctrl = slb.And(ctrl, slb.Not(selection[depth]))
    segtree(ctrl=new_ctrl, selection=selection, target=target, depth=depth+1, left=left, right=mid)
    
    # Recurse right interval
    new_ctrl = slb.And(ctrl, selection[depth])
    segtree(ctrl=new_ctrl, selection=selection, target=target, depth=depth+1, left=mid, right=right)
    
    # Quantum note:
    # instead of throwing away the first value of `new_ctrl` and re-anding
    # with selection, we can just invert the first one (but only if `ctrl` is active)
    # new_ctrl ^= ctrl

selection = [S.Symbol(f's{i}') for i in range(3)]
target = [S.false for i in range(2**len(selection)) ]
segtree(S.true, selection, target, 0, 0, 2**len(selection))

print()
print("Target:")
for n, t in enumerate(target):
    print(f'n={n} {slb.simplify_logic(t)}')

depth=0 left=0 right=8 
depth=1 left=0 right=4 
depth=2 left=0 right=2 
depth=3 left=0 right=1 n=7 ctrl=~s0 & ~s1 & ~s2
depth=3 left=1 right=2 n=7 ctrl=s2 & ~s0 & ~s1
depth=2 left=2 right=4 
depth=3 left=2 right=3 n=7 ctrl=s1 & ~s0 & ~s2
depth=3 left=3 right=4 n=7 ctrl=s1 & s2 & ~s0
depth=1 left=4 right=8 
depth=2 left=4 right=6 
depth=3 left=4 right=5 n=7 ctrl=s0 & ~s1 & ~s2
depth=3 left=5 right=6 n=7 ctrl=s0 & s2 & ~s1
depth=2 left=6 right=8 
depth=3 left=6 right=7 n=7 ctrl=s0 & s1 & ~s2
depth=3 left=7 right=8 n=7 ctrl=s0 & s1 & s2

Target:
n=0 ~s0 & ~s1 & ~s2
n=1 s2 & ~s0 & ~s1
n=2 s1 & ~s0 & ~s2
n=3 s1 & s2 & ~s0
n=4 s0 & ~s1 & ~s2
n=5 s0 & s2 & ~s1
n=6 s0 & s1 & ~s2
n=7 s0 & s1 & s2

print(f"{'n':3s} | {'segtree':18s} | {'total control':18s} | same?")
for n, (t1, t2) in enumerate(zip(target, tc_target)):
    t1 = slb.simplify_logic(t1)
    print(f'{n:3d} | {str(t1):18s} | {str(t2):18s} | {str(t1==t2)}')

n   | segtree            | total control      | same?
  0 | ~s0 & ~s1 & ~s2    | ~s0 & ~s1 & ~s2    | True
  1 | s2 & ~s0 & ~s1     | s2 & ~s0 & ~s1     | True
  2 | s1 & ~s0 & ~s2     | s1 & ~s0 & ~s2     | True
  3 | s1 & s2 & ~s0      | s1 & s2 & ~s0      | True
  4 | s0 & ~s1 & ~s2     | s0 & ~s1 & ~s2     | True
  5 | s0 & s2 & ~s1      | s0 & s2 & ~s1      | True
  6 | s0 & s1 & ~s2      | s0 & s1 & ~s2      | True
  7 | s0 & s1 & s2       | s0 & s1 & s2       | True

Quantum Circuit

We can translate the boolean logic to reversible, quantum logic. It is instructive to start from the suboptimal total control quantum circuit for comparison purposes. We can build this as in the sympy boolean-logic case by adding controlled X operations to the target signature, with the controls on the selection signature toggled on or off according to the binary representation of the selection index.

Let us first build a GateWithRegisters object to implement the circuit

import cirq
from functools import cached_property
from qualtran import Signature, GateWithRegisters, QUInt

class TotallyControlledNot(GateWithRegisters):
    
    def __init__(self, selection_bitsize: int, target_bitsize: int, control_bitsize: int = 1):
        self._selection_bitsize = selection_bitsize
        self._target_bitsize = target_bitsize
        self._control_bitsize = control_bitsize

    @cached_property
    def signature(self) -> Signature:
        return Signature(
            [
                *Signature.build(control=self._control_bitsize),
                *Signature.build(selection=self._selection_bitsize),
                *Signature.build(target=self._target_bitsize)
            ]
        )

    def decompose_from_registers(self, **qubit_regs: Sequence[cirq.Qid]) -> cirq.OP_TREE:
        num_controls = self._control_bitsize + self._selection_bitsize
        for target_bit in range(self._target_bitsize):
            bit_pattern = QUInt(self._selection_bitsize).to_bits(target_bit)
            control_values = [1]*self._control_bitsize + list(bit_pattern)
            yield cirq.X.controlled(
                num_controls=num_controls,
                control_values=control_values
                ).on(
                    *qubit_regs["control"], 
                    *qubit_regs["selection"],
                    qubit_regs["target"][-(target_bit+1)])
        

import qualtran.cirq_interop.testing as cq_testing
tc_not = TotallyControlledNot(3, 5)
tc = cq_testing.GateHelper(tc_not)
cirq.Circuit((cirq.decompose_once(tc.operation)))
SVGCircuit(cirq.Circuit(cirq.decompose_once(tc.operation)))

Tests for Correctness

We can use a full statevector simulation to compare the desired statevector to the one generated by the unary iteration circuit for each basis state.

selection_bitsize = 3
target_bitsize = 5
for n in range(target_bitsize):
    # Initial qubit values
    qubit_vals = {q: 0 for q in tc.all_qubits}
    # All controls 'on' to activate circuit
    qubit_vals.update({c: 1 for c in tc.quregs['control']})
    # Set selection according to `n`
    qubit_vals.update(zip(tc.quregs['selection'], QUInt(selection_bitsize).to_bits(n)))

    initial_state = [qubit_vals[x] for x in tc.all_qubits]
    final_state = [qubit_vals[x] for x in tc.all_qubits]
    final_state[-(n+1)] = 1
    cq_testing.assert_circuit_inp_out_cirqsim(
        tc.circuit, tc.all_qubits, initial_state, final_state
    )
    print(f'n={n} checked!')

n=0 checked!
n=1 checked!
n=2 checked!
n=3 checked!
n=4 checked!

Towards a segment tree

Next let’s see how we can reduce the circuit to the observe the tree structure. First let’s recall what we are trying to do with the controlled not. Given a selection integer (say 3 = 011), we want to toggle the bit in the target register to on if the qubit 1 and 2 are set to on in the selection register.

# The selection bits [1-3] are set according to binary representation of the number 3 (011)
initial_state = [1, 0, 1, 1, 0, 0, 0, 0, 0]
final_state = [1, 0, 1, 1, 0, 1, 0, 0, 0]
actual, should_be = cq_testing.get_circuit_inp_out_cirqsim(
        tc.circuit, tc.all_qubits, initial_state, final_state
    )
print("simulated: ", actual)
print("expected : ", should_be)

simulated:  101101000
expected :  101101000

Now what is important to note is that we can remove many repeated controlled operations by using ancilla qubits to flag what part of the circuit we need to apply, this works because we know the bit pattern of nearby integers is very similar.

A circuit demonstrating this for our example is given below.

from qualtran.bloqs.mcmt import And

selection_bitsize = 2
target_bitsize = 4
qubits = cirq.LineQubit(0).range(1 + selection_bitsize * 2 + target_bitsize)
circuit = cirq.Circuit()
circuit.append(
    [
        And(1, 0).on(qubits[0], qubits[1], qubits[2]),
        And(1, 0).on(qubits[2], qubits[3], qubits[4]),
        cirq.CX(qubits[4], qubits[8]),
        cirq.CNOT(qubits[2], qubits[4]),
        cirq.CX(qubits[4], qubits[7]),
        And().adjoint().on(qubits[2], qubits[3], qubits[4]),
        cirq.CNOT(qubits[0], qubits[2]),
        And(1, 0).on(qubits[2], qubits[3], qubits[4]),
        cirq.CX(qubits[4], qubits[6]),
        cirq.CNOT(qubits[2], qubits[4]),
        cirq.CX(qubits[4], qubits[5]),
        And().adjoint().on(qubits[2], qubits[3], qubits[4]),
        And().adjoint().on(qubits[0], qubits[1], qubits[2]),
    ]
)

SVGCircuit(circuit)

Reading from left to right we first check the control is set to on and the selection qubit is off, if both these conditions are met then the ancilla qubit is now set to 1. The next control checks if the previous condition was met and likewise the second selection index is also off. At this point if both these conditions are met we must be indexing 0 as the first two qubits are set to off (00), otherwise we know that we want to apply X to qubit 1 so we perform a CNOT operation to flip the bit value in the second ancilla qubit, before returning back up the circuit. Now if the left half of the circuit was not applied (i.e. the first selection register was set to 1) then the CNOT between the control qubit and the first ancilla qubit causes the ancilla qubit to toggle on. This triggers the right side of the circuit, which now performs the previously described operations to figure out if the lowest bit is set. Combining these two then yields the expected controlled X operation.

Below we check the circuit is giving the expected behaviour.

initial_state = [1, 0, 0, 0, 0, 0, 0, 0, 0]
target_indx = 3
sel_bits = QUInt(selection_bitsize).to_bits(target_indx)
sel_indices = [i for i in range(1, 2*selection_bitsize+1, 2)]
initial_state[sel_indices[0]] = sel_bits[0]
initial_state[sel_indices[1]] = sel_bits[1]
result = cirq.Simulator(dtype=np.complex128).simulate(
    circuit, initial_state=initial_state
)
actual = result.dirac_notation(decimals=2)[1:-1]
print("simulated: {}, index set in string {}".format(actual, len(qubits)-1-target_indx))

simulated: 110101000, index set in string 5

Extending the above idea to larger ranges of integers is relatively straightforward. For example consider the next simplest case of \(L=8 = 2^3\). The circuit above takes care of the last two bits and can be duplicated. For the extra bit we just need to add a additional AND operations, and a CNOT to switch between the original range [0,3] or the new range [4,7] depending on whether the new selection register is off or on respectively. This procedure can be repeated and we can begin to notice the recursive tree-like structure.

This structure is just the segtree described previously for boolean logic and this gives is the basic idea of unary iteration, which uses L-1 AND operations. Below the ApplyXToLthQubit builds the controlled Not operation using the UnaryIterationGate as a base class which defines the decompose_from_registers method appropriately to recursively construct the unary iteration circuit.

Note below a different ordering of ancilla and selection qubits is taken to what was used in the simpler L=4 example.

from qualtran import QAny, Register, Register, BQUInt
from qualtran.bloqs.multiplexers.unary_iteration_bloq import UnaryIterationGate
from functools import cached_property



class ApplyXToLthQubit(UnaryIterationGate):
    def __init__(self, selection_bitsize: int, target_bitsize: int, control_bitsize: int = 1):
        self._selection_bitsize = selection_bitsize
        self._target_bitsize = target_bitsize
        self._control_bitsize = control_bitsize

    @cached_property
    def control_registers(self) -> Tuple[Register, ...]:
        return (Register('control', QAny(self._control_bitsize)),)

    @cached_property
    def selection_registers(self) -> Tuple[Register, ...]:
        return (Register('selection', BQUInt(self._selection_bitsize, self._target_bitsize)),)

    @cached_property
    def target_registers(self) -> Tuple[Register, ...]:
        return (Register('target', QAny(self._target_bitsize)),)

    def nth_operation(
        self, context, control: cirq.Qid, selection: int, target: Sequence[cirq.Qid]
    ) -> cirq.OP_TREE:
        return cirq.CNOT(control, target[-(selection + 1)])

import qualtran.cirq_interop.testing as cq_testing
selection_bitsize = 3
target_bitsize = 5

g = cq_testing.GateHelper(
    ApplyXToLthQubit(selection_bitsize, target_bitsize))
SVGCircuit(cirq.Circuit(cirq.decompose_once(g.operation)))

Tests for Correctness

We can use a full statevector simulation to check again that the optimized circuit produces the expected result.

from qualtran import QUInt

for n in range(target_bitsize):
    # Initial qubit values
    qubit_vals = {q: 0 for q in g.all_qubits}
    # All controls 'on' to activate circuit
    qubit_vals.update({c: 1 for c in g.quregs['control']})
    # Set selection according to `n`
    qubit_vals.update(zip(g.quregs['selection'], QUInt(selection_bitsize).to_bits(n)))

    initial_state = [qubit_vals[x] for x in g.all_qubits]
    qubit_vals[g.quregs['target'][-(n + 1)]] = 1
    final_state = [qubit_vals[x] for x in g.all_qubits]
    cq_testing.assert_circuit_inp_out_cirqsim(
        g.decomposed_circuit, g.all_qubits, initial_state, final_state
    )
    print(f'n={n} checked!')

/usr/local/google/home/mpharrigan/qualtran/conda-311/lib/python3.11/site-packages/cotengra/hyperoptimizers/hyper.py:57: UserWarning: Couldn't find `optuna`, `cmaes`, `baytune (btb)`, `chocolate`, or `nevergrad` so will use completely random sampling in place of hyper-optimization.
  warnings.warn(
/usr/local/google/home/mpharrigan/qualtran/conda-311/lib/python3.11/site-packages/cotengra/hyperoptimizers/hyper.py:76: UserWarning: Couldn't find `optuna`, `cmaes`, `baytune (btb)`, `chocolate`, or `nevergrad` so will use completely random sampling in place of hyper-optimization.
  warnings.warn(

n=0 checked!
n=1 checked!

n=2 checked!
n=3 checked!

n=4 checked!