Quantum GCD circuit: On reversibility and clearing ancillae Originally posted on PHYS, however, obviously it has more to do with CS

I am currently trying to implement a circuit for computing the greatest common divisor in the Quantum Computing Language. In my work, I partially proceed from the paper Quantum Circuits for GCD Computation with O(nlogn) Depth and O(n) Ancillae published by I. Markov and M. Saeedi in 2013, which may be found at arXiv http://arxiv.org/abs/1304.7516.

The paper introduces a binary GCD algorithm analogue and works out a design of the quantum circuit which is clarified in least details, however, in the same time it is somewhat confusing due to several statements that actually do not correspond with the rest of the work. First of all, let me enclose an image denoting design of the circuit. 

(Credit goes to the paper)

As you may see, the circuit entirely corresponds to the classical binary GCD algorithm. It takes advantage of quantum circuits for comparison, bit shift, addition and multiplication. Nevertheless, in our case those functions are controlled by state of ancillary qubits and thus, they may not regarded reversible on their own. It means, we have to account for a separate register of auxiliary qubits which cannot be erased till the end of our computation as it would simply violate desired reversibility otherwise. 

But in contrast to that, in the paper we may find at least three statements which do actually deny those simple facts. Let me cite: (About preceding works) The authors did not clear all zero-initialized ancillae which limits the applicability of their techniques. (About the circuit) The zero-initialized ancillae may be modified in-side a given circuit, but should be returned to zero at the end of computation to be reused.  (About implementation) To count ancillae, note that all computational ancillae are cleared inside each block. After the final multiplication block for R·B, one can copy (in log depth) the final GCD result to another n-qubit zero-initialized register and apply the whole circuit (except for copying the result) in reverse order to recover A,B, and zero-initialized ancillae. Given that all components use O(n) ancillae, the total number of ancillae remains linear.

I ask myself: How can be all ancillaes cleared at the end of each computational block while the binary GCD seems to be reversible only with regard to an auxiliary register of control qubits? It is apparent this question reduces to the problem of binary GCD reversibility which is actually something I am doubtful about. Maybe I lack of insight and the usage of other control bits is so obvious so as to omit it, however, in every case I would like to ask you about your view on the problem.
It looks like they just directly translated the algorithm into a circuit, so it's not too hard to reproduce.


Write out the binary gcd algorithm as a loop:

def gcd(a, b):
    let shift = 0
    while a != 0 or b != 0:
        let pa = a &amp; 1
        let pb = b &amp; 1
        let pc = a &gt;= b

        if pa == 0 and pb == 0:
            a &gt;&gt;= 1
            b &gt;&gt;= 1
            shift++
        elif pa == 0:
            a &gt;&gt;= 1
        elif pb == 0:
            b &gt;&gt;= 1
        elif pc:
            a -= b
            a &gt;&gt;= 1
        else:
            b -= a
            b &gt;&gt;= 1

    return (a | b) &lt;&lt; shift

Keep track of any intermediate values so you can uncompute them instead of leaking them. Also switch to a fixed number of iterations to avoid leaking the size of the inputs.

def gcd(a, b, len_a, len_b):
    let shift = 0
    # Note: &gt;&gt; and &lt;&lt; are now bit rotations instead of bit shifts
    for i in range(len_a + len_b):
        let g = a &gt; b
        let pa = (a &amp; 1) == 0
        let pb = (b &amp; 1) == 0

        if pa and pb and a != 0 and b != 0: shift++
        if pa: a &gt;&gt;= 1
        if pb: b &gt;&gt;= 1
        if not pa and not pb:
            if g:
                a -= b
                a &gt;&gt;= 1
            else:
                b -= a
                b &gt;&gt;= 1

        # We can't clear these yet, so just set them aside until later
        push g, pa, pb

    let result = (a | b) &lt;&lt; shift

    # Uncompute intermediate values
    for i in range(n + m):
        pop pb, pa, g

        if not pa and not pb:
            if g:
                a &lt;&lt;= 1
                a += b
            else:
                b &lt;&lt;= 1
                b += a
        if pb: b &lt;&lt;= 1
        if pa: a &lt;&lt;= 1
        if pa and pb and a != 0 and b != 0: shift--

        unlet pb = (b &amp; 1) == 0
        unlet pa = (a &amp; 1) == 0
        unlet g = a &gt; b

    unlet shift = 0
    return result

Convert to a circuit.

(Edit: oops, I forgot the shift part. You'll have to add that part in yourself.)




It takes $O(n)$ iterations to guarantee the inputs have been reduced all the way down, and each iteration requires $\Omega(\lg n)$ depth and $O(n)$ gates so the overall bound is $O(n^2)$ gates and $O(n \lg n)$ depth.

Now, all that being said, you need to be wary of the fact that the result is not totally coherent. Before doing interference-y stuff with it, I'd want to be a lot more careful about ensuring the stashed ancilla are independent of the computed gcd. Also you'd want to do any computation in the middle so that the result doesn't have to co-exist with the inputs.

But I hope this gives the general idea.