Arithmetic (adders, comparators, modular)

A register is a group of qubits read as one binary number. Quantum arithmetic adds, subtracts and compares these numbers while keeping every superposition intact: an adder applied to a register holding |0> + |1> returns |c> + |c+1>, both sums at once. These operations are the building blocks of Shor's algorithm and of many oracles.

Every arithmetic operation is reversible, so it simply moves each basis state to another basis state. Dense-Evolution applies it exactly that way, as a permutation of the statevector amplitudes.

Step 1. Load two numbers into registers

import numpy as np
import dense_evolution as de

qasm = """OPENQASM 2.0;
include "qelib1.inc";
qreg q[5];
x q[0];
x q[1];
x q[3];
"""
circ = de.QASMParser().parse(qasm)
sim = de.DenseSVSimulator(5)
sim.run_circuit_jit(circ)
sv0 = sim.get_statevector()
print(format(int(np.argmax(np.abs(sv0))), '05b'))

11010

Register a is qubits [1, 0] and register b is qubits [4, 3, 2], each listed least significant bit first. The x gates set a = 3 (both qubits on) and b = 2 (only its middle bit, qubit 3, on). The printed bitstring is q0 q1 q2 q3 q4.

Step 2. Add a into b

sv1 = de.add_registers(sv0, 5, [1, 0], [4, 3, 2])
print(format(int(np.argmax(np.abs(sv1))), '05b'))

11101

add_registers maps |a, b> to |a, a + b>. Register a keeps 3; register b now holds 5 (qubits 2 and 4 on). b has one qubit more than a, so the sum can never overflow.

Step 3. Subtract it back

sv2 = de.subtract_registers(sv1, 5, [1, 0], [4, 3, 2])
print(format(int(np.argmax(np.abs(sv2))), '05b'))

11010

subtract_registers maps |a, b> to |a, b - a> and is the exact inverse of add_registers: the state returns to Step 1. If b < a the result wraps around and the most significant qubit of b is 1, which is how a subtraction tells which number is larger.

Step 4. Add a fixed number to a superposition

qasm = """OPENQASM 2.0;
include "qelib1.inc";
qreg q[3];
h q[2];
"""
circ = de.QASMParser().parse(qasm)
sim = de.DenseSVSimulator(3)
sim.run_circuit_jit(circ)
sv0 = sim.get_statevector()
sv1 = de.add_constant(sv0, 3, [2, 1, 0], 3)
print(np.round(np.asarray(sv1), 4))

[0.    +0.j 0.    +0.j 0.    +0.j 0.7071+0.j 0.7071+0.j 0.    +0.j
 0.    +0.j 0.    +0.j]

h q[2] puts the register [2, 1, 0] in (|0> + |1>)/√2. Adding the classical number 3 gives (|3> + |4>)/√2: both values are shifted, and nothing is measured. c = 1 is an increment and c = -1 a decrement.

Step 5. Compare a register with a number

qasm = """OPENQASM 2.0;
include "qelib1.inc";
qreg q[4];
h q[2];
x q[1];
"""
circ = de.QASMParser().parse(qasm)
sim = de.DenseSVSimulator(4)
sim.run_circuit_jit(circ)
sv0 = sim.get_statevector()
sv1 = de.compare_constant(sv0, 4, [2, 1], 3, 3, '<')
print(np.round(np.asarray(sv1), 4))

[0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.7071+0.j
 0.7071+0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j
 0.    +0.j 0.    +0.j 0.    +0.j 0.    +0.j]

Register b is qubits [2, 1] and holds (|2> + |3>)/√2. compare_constant flips the output qubit 3 where b < 3 is true: only the |2> branch, so the answer is now entangled with the register (amplitudes at 0101 and 0110). The registers are never changed, only the output qubit. op can be '<', '>', '<=', '>=', '==' or '!='; compare_registers does the same between two registers.

Step 6. Modular exponentiation, the core of Shor's algorithm

qasm = """OPENQASM 2.0;
include "qelib1.inc";
qreg q[6];
h q[0];
h q[1];
x q[5];
"""
circ = de.QASMParser().parse(qasm)
sim = de.DenseSVSimulator(6)
sim.run_circuit_jit(circ)
sv0 = sim.get_statevector()
sv1 = de.power_mod(sv0, 6, [1, 0], [5, 4, 3, 2], 7, 15)
p = np.abs(np.asarray(sv1)) ** 2
print([format(int(k), '06b') for k in np.flatnonzero(p > 1e-9)])

['000001', '010111', '100100', '111101']

Register x is qubits [1, 0], put in an equal superposition of 0, 1, 2, 3 by the two h gates. Register y is qubits [5, 4, 3, 2], set to 1 by x q[5]. power_mod maps |x, y> to |x, y * 7^x mod 15>, so each branch now pairs x with 7^x mod 15: 1, 7, 4, 13. Reading the four bitstrings, x is the first two digits and y the last four, least significant last. This is the state Shor's algorithm measures to find the period of 7^x mod 15.

The same module has the steps that build it: add_registers_mod and add_constant_mod ((a + b) mod N), multiply_add_mod (b + a*x mod N) and multiply_mod (x -> a*x mod N in place, for a coprime with N; with N = 2^n and odd a it is plain multiplication). The modular functions accept controls, a list of qubits that must all be 1 for the operation to act.

Step 7. The same additions as gate circuits

import numpy as np
import dense_evolution as de
from dense_evolution.circuits.arithmetic_qasm import cuccaro_adder_qasm

circ = de.QASMParser().parse(cuccaro_adder_qasm(2))
sim = de.DenseSVSimulator(circ.n_qubits)
sim.set_initial_state(np.eye(2 ** circ.n_qubits)[0b110100])
sim.run_circuit_jit(circ.to_tuples())
print(format(int(np.argmax(np.abs(sim.get_statevector()))), '06b'))

111010

The functions above move amplitudes directly, which is exact but is not a circuit you could run on a device or put through a noise model. cuccaro_adder_qasm(n) returns the ripple-carry adder of Cuccaro et al. as OPENQASM 2.0: registers a[2], b[3] and one ancilla x[1], each least significant bit first. The input 110100 is a = 3, b = 2; the output 111010 is a = 3, b = 5, ancilla back to 0. The circuit uses 4 Toffoli and 9 CNOT gates (2n and 4n + 1).

The same module has draper_adder_qasm(n) (addition in Fourier space, no ancilla), constant_adder_qasm(c, n) (adds a fixed number) and modular_constant_adder_qasm(c, N, n) (adds a fixed number modulo N). Each one gives the same state as the matching function above: add_registers, add_constant, add_constant_mod. cmult_mod_qasm(a, N, n) and controlled_ua_qasm(a, N, n) go one level up, to the controlled modular multiplication of Shor's algorithm, and match multiply_add_mod and multiply_mod.

Step 8. Shor's algorithm on 11 qubits

from fractions import Fraction
from math import gcd
from dense_evolution.circuits.shor import shor_order_finding

y, m = shor_order_finding(7, 15, rng=2)
r = Fraction(y, 2 ** m).limit_denominator(15).denominator
y, r, gcd(7 ** (r // 2) - 1, 15), gcd(7 ** (r // 2) + 1, 15)

(192, 4, 3, 5)

Shor's algorithm factors N by finding the order r of a number a: the smallest r with a^r mod N = 1. shor_order_finding runs the quantum part with the gate circuits of Step 7 and returns one measured integer y out of 2^m (m = 8 for N = 15). y / 2^m = 192 / 256 = 3/4, so the order of 7 modulo 15 is 4, and gcd(7^2 - 1, 15) and gcd(7^2 + 1, 15) give the factors 3 and 5. Each run is one measurement, so another seed can return 0 or 1/2, which only gives a divisor of r; in practice the run is repeated. shor_order_finding(7, 15, noise_model="depolarizing", p=0.01) runs the same algorithm with noise after each controlled multiplication (Details).

---

Details

Why a permutation. A reversible function f on basis states has the unitary U|x> = |f(x)> (Vedral, Barenco, Ekert, Quantum networks for elementary arithmetic operations, quant-ph/9511018, Sect. II). Applying it as an index permutation is exact and touches each amplitude once.

Checked against the gate-level circuits. The tests build the circuits from the papers out of engine gates and compare them with these functions on every basis state:

• add_registers against the ripple-carry adder with one ancilla (Cuccaro, Draper, Kutin, Moulton, quant-ph/0410184, MAJ/UMA construction), for 1, 2 and 3-bit numbers;
• add_constant against the transform adder in Fourier space with a classical addend (Draper, quant-ph/0008033, Sect. 5; Beauregard, quant-ph/0205095, Sect. 2.1);
• compare_registers(..., '<') against the comparator: complement a, compute only the high bit of the sum, undo (Cuccaro et al., Sect. 4.3);
• compare_constant(..., '<') against the most significant qubit after the reverse φADD(c) (Beauregard, Sect. 2.1, Fig. 4), for c < 2^n as in the paper.
• add_constant_mod against the modular adder built from adders, a subtraction of N, the overflow qubit and its reset (Beauregard, Sect. 2.2, Fig. 5);
• multiply_add_mod against n controlled modular additions of 2^i*a mod N (Vedral et al., Sect. III.C; Beauregard, Sect. 2.3, Eq. 2);
• multiply_mod against the controlled-U_a circuit: multiply-add, controlled swap, inverse multiply-add by a^-1 mod N (Beauregard, Sect. 2.3, Fig. 7);
• power_mod against the chain of controlled multiplications by a^(2^i) mod N (Vedral et al., Sect. III.D).

Modular domain. The modular functions act on values below N, as in the papers. Basis states holding a value >= N are left unchanged, so every function stays a permutation and therefore unitary.

Register sizes. With len(b) == len(a) + 1 the sum is exact (Vedral et al., Sect. III.A); with len(b) == len(a) it is addition modulo 2^n (Cuccaro et al., Sect. 4.1).

Gate circuits in arithmetic_qasm. The OPENQASM 2.0 generators follow the figures of the papers: MAJ and UMA as user-defined gates (Cuccaro et al., Fig. 1-2, the two-CNOT UMA), the adder of Fig. 3-4; the QFT adder without swaps (Draper); φADD(c) and φADD(c)MOD(N) (Beauregard, Fig. 3 and 5), here without the two controls. CMULT(a)MOD(N) is n doubly controlled φADD(2^i a mod N)MOD(N) (Fig. 6), the controlled-U_a is CMULT(a), a controlled swap and the inverse of CMULT(a^-1) (Fig. 7): 2n + 3 qubits in total. shor_order_finding uses one control qubit, measured and reset after each controlled-U_{a^(2^k)}, with the inverse QFT done semiclassically (Fig. 8); the measurement is sampled in Python between circuit runs, since QASMParser skips measure, reset and if, so one QASM program cannot feed a measurement back into later gates. Shor under noise. With noise_model and p, one stochastic trajectory of NoiseModel.apply_to_sv acts on all 11 qubits after each of the 8 controlled multiplications. For a = 7, N = 15, 40 runs per row (seeds 0-39), fraction of outcomes on the four ideal peaks {0, 64, 128, 192} (a uniform guess would give 4/256):

Depolarizing p	On a peak
0	40/40
0.01	37/40
0.05	32/40
0.2	8/40

tests/unit/test_arithmetic_qasm.py compares each one with the permutation functions on random states over every basis state of its domain. The modular adder needs 0 <= b < N and 0 <= c < N, as in the paper; on b >= N it is outside its contract.

arithmetic

Reversible arithmetic on qubit registers, applied to a statevector.

Every operation here is a bijection f on computational basis states, so its unitary is the permutation U|x> = |f(x)> (Vedral, Barenco, Ekert, quant-ph/9511018, Sect. II, Eqs. 2-3). It is applied directly as an index permutation of the statevector, which is exact and costs one gather over the 2^n amplitudes. The gate-level circuits from the papers (Cuccaro et al. quant-ph/0410184 ripple-carry adder, Draper quant-ph/0008033 transform adder) are the references the tests compare against.

A register is a list of qubit indices, least significant bit first. Qubit 0 is the most significant bit of the statevector index, as everywhere in the engine.

add_registers

add_registers(sv, n_qubits, a, b)

|a, b> -> |a, (a + b) mod 2^len(b)>.

With len(b) = len(a) + 1 the sum never overflows (Vedral et al., Sect. III.A, Eq. 9); with len(b) = len(a) it is addition modulo 2^n (Cuccaro et al., Sect. 4.1).

Source code in dense_evolution/circuits/arithmetic.py

def add_registers(sv, n_qubits, a, b):
    """
    |a, b> -> |a, (a + b) mod 2^len(b)>.

    With len(b) = len(a) + 1 the sum never overflows (Vedral et al.,
    Sect. III.A, Eq. 9); with len(b) = len(a) it is addition modulo 2^n
    (Cuccaro et al., Sect. 4.1).
    """
    _check(n_qubits, a, b)
    idx, va = _register_values(n_qubits, a)
    _, vb = _register_values(n_qubits, b)
    target = _with_register(n_qubits, idx, b, (va + vb) % (1 << len(b)))
    return _permute(sv, target)

subtract_registers

subtract_registers(sv, n_qubits, a, b)

|a, b> -> |a, (b - a) mod 2^len(b)>, the inverse of add_registers.

If b < a the result is 2^len(b) - (a - b) and, with len(b) = len(a) + 1, its most significant qubit is 1 (Vedral et al., Sect. III.A).

Source code in dense_evolution/circuits/arithmetic.py

def subtract_registers(sv, n_qubits, a, b):
    """
    |a, b> -> |a, (b - a) mod 2^len(b)>, the inverse of add_registers.

    If b < a the result is 2^len(b) - (a - b) and, with len(b) = len(a) + 1,
    its most significant qubit is 1 (Vedral et al., Sect. III.A).
    """
    _check(n_qubits, a, b)
    idx, va = _register_values(n_qubits, a)
    _, vb = _register_values(n_qubits, b)
    target = _with_register(n_qubits, idx, b, (vb - va) % (1 << len(b)))
    return _permute(sv, target)

add_constant

add_constant(sv, n_qubits, b, c)

|b> -> |(b + c) mod 2^len(b)> for a classical integer c.

c = 1 is increment, c = -1 decrement. Same operation as Draper's transform adder with classical addend (Sect. 5) and Beauregard's phiADD(c) (quant-ph/0205095, Sect. 2.1).

Source code in dense_evolution/circuits/arithmetic.py

def add_constant(sv, n_qubits, b, c):
    """
    |b> -> |(b + c) mod 2^len(b)> for a classical integer c.

    c = 1 is increment, c = -1 decrement. Same operation as Draper's
    transform adder with classical addend (Sect. 5) and Beauregard's
    phiADD(c) (quant-ph/0205095, Sect. 2.1).
    """
    _check(n_qubits, b)
    idx, vb = _register_values(n_qubits, b)
    target = _with_register(n_qubits, idx, b, (vb + int(c)) % (1 << len(b)))
    return _permute(sv, target)

compare_registers

compare_registers(sv, n_qubits, a, b, out, op='<')

|a, b, z> -> |a, b, z XOR [a op b]>, inputs unchanged.

op = '<' is the comparator of Cuccaro et al. (Sect. 4.3): the high bit of a - b, which is 1 if and only if a < b. The other comparisons follow by swapping the operands and negating the output qubit.

Source code in dense_evolution/circuits/arithmetic.py

def compare_registers(sv, n_qubits, a, b, out, op='<'):
    """
    |a, b, z> -> |a, b, z XOR [a op b]>, inputs unchanged.

    op = '<' is the comparator of Cuccaro et al. (Sect. 4.3): the high bit of
    a - b, which is 1 if and only if a < b. The other comparisons follow by
    swapping the operands and negating the output qubit.
    """
    _check(n_qubits, a, b, [out])
    idx, va = _register_values(n_qubits, a)
    _, vb = _register_values(n_qubits, b)
    return _permute(sv, _flip_if(n_qubits, idx, out, _comparison(op)(va, vb)))

compare_constant

compare_constant(sv, n_qubits, b, c, out, op='<')

|b, z> -> |b, z XOR [b op c]> for a classical integer c.

op = '<' is the most significant qubit after the reverse phiADD(c) of Beauregard (quant-ph/0205095, Sect. 2.1, Fig. 4).

Source code in dense_evolution/circuits/arithmetic.py

def compare_constant(sv, n_qubits, b, c, out, op='<'):
    """
    |b, z> -> |b, z XOR [b op c]> for a classical integer c.

    op = '<' is the most significant qubit after the reverse phiADD(c) of
    Beauregard (quant-ph/0205095, Sect. 2.1, Fig. 4).
    """
    _check(n_qubits, b, [out])
    idx, vb = _register_values(n_qubits, b)
    return _permute(sv, _flip_if(n_qubits, idx, out, _comparison(op)(vb, int(c))))

add_registers_mod

add_registers_mod(sv, n_qubits, a, b, N)

|a, b> -> |a, (a + b) mod N> for 0 <= a, b < N (Vedral et al., Sect. III.B, Eq. 10). Basis states with a >= N or b >= N are left unchanged, which keeps the map a permutation.

Source code in dense_evolution/circuits/arithmetic.py

def add_registers_mod(sv, n_qubits, a, b, N):
    """
    |a, b> -> |a, (a + b) mod N> for 0 <= a, b < N (Vedral et al., Sect. III.B,
    Eq. 10). Basis states with a >= N or b >= N are left unchanged, which
    keeps the map a permutation.
    """
    _check(n_qubits, a, b)
    _check_modulus(N, a, b)
    idx, va = _register_values(n_qubits, a)
    _, vb = _register_values(n_qubits, b)
    ok = (va < N) & (vb < N)
    return _permute(sv, _with_register(n_qubits, idx, b, np.where(ok, (va + vb) % N, vb)))

add_constant_mod

add_constant_mod(sv, n_qubits, b, c, N, controls=())

|b> -> |(b + c) mod N> for 0 <= b < N, applied only where every control qubit is 1: the doubly controlled phiADD(c)MOD(N) of Beauregard (quant-ph/0205095, Sect. 2.2, Fig. 5). Basis states with b >= N are left unchanged.

Source code in dense_evolution/circuits/arithmetic.py

def add_constant_mod(sv, n_qubits, b, c, N, controls=()):
    """
    |b> -> |(b + c) mod N> for 0 <= b < N, applied only where every control
    qubit is 1: the doubly controlled phiADD(c)MOD(N) of Beauregard
    (quant-ph/0205095, Sect. 2.2, Fig. 5). Basis states with b >= N are left
    unchanged.
    """
    _check(n_qubits, b, *([list(controls)] if controls else []))
    _check_modulus(N, b)
    idx, vb = _register_values(n_qubits, b)
    ok = _active(n_qubits, idx, controls) & (vb < N)
    return _permute(sv, _with_register(n_qubits, idx, b, np.where(ok, (vb + int(c)) % N, vb)))

multiply_add_mod

multiply_add_mod(sv, n_qubits, x, b, a, N, controls=())

|x, b> -> |x, (b + ax) mod N> for 0 <= b < N, applied only where every control qubit is 1: CMULT(a)MOD(N) of Beauregard (Sect. 2.3, Fig. 6), built in the paper from n controlled modular additions of 2^ia mod N (Vedral et al., Sect. III.C). Basis states with b >= N are left unchanged.

Source code in dense_evolution/circuits/arithmetic.py

def multiply_add_mod(sv, n_qubits, x, b, a, N, controls=()):
    """
    |x, b> -> |x, (b + a*x) mod N> for 0 <= b < N, applied only where every
    control qubit is 1: CMULT(a)MOD(N) of Beauregard (Sect. 2.3, Fig. 6),
    built in the paper from n controlled modular additions of 2^i*a mod N
    (Vedral et al., Sect. III.C). Basis states with b >= N are left unchanged.
    """
    _check(n_qubits, x, b, *([list(controls)] if controls else []))
    _check_modulus(N, b)
    idx, vx = _register_values(n_qubits, x)
    _, vb = _register_values(n_qubits, b)
    ok = _active(n_qubits, idx, controls) & (vb < N)
    new = (vb + (int(a) % N) * (vx % N)) % N
    return _permute(sv, _with_register(n_qubits, idx, b, np.where(ok, new, vb)))

multiply_mod

multiply_mod(sv, n_qubits, x, a, N, controls=())

|x> -> |a*x mod N> in place for 0 <= x < N and gcd(a, N) = 1, applied only where every control qubit is 1 (Vedral et al., Sect. II, Eqs. 4-6; Beauregard, Sect. 2.3, Fig. 7). With N = 2^len(x) and odd a this is plain multiplication modulo 2^n. Basis states with x >= N are left unchanged.

Source code in dense_evolution/circuits/arithmetic.py

def multiply_mod(sv, n_qubits, x, a, N, controls=()):
    """
    |x> -> |a*x mod N> in place for 0 <= x < N and gcd(a, N) = 1, applied only
    where every control qubit is 1 (Vedral et al., Sect. II, Eqs. 4-6;
    Beauregard, Sect. 2.3, Fig. 7). With N = 2^len(x) and odd a this is plain
    multiplication modulo 2^n. Basis states with x >= N are left unchanged.
    """
    _check(n_qubits, x, *([list(controls)] if controls else []))
    _check_modulus(N, x)
    if np.gcd(int(a), N) != 1:
        raise ValueError(f"a={a} and N={N} must be coprime for an in-place multiplication")
    idx, vx = _register_values(n_qubits, x)
    ok = _active(n_qubits, idx, controls) & (vx < N)
    return _permute(sv, _with_register(n_qubits, idx, x, np.where(ok, (int(a) * vx) % N, vx)))

power_mod

power_mod(sv, n_qubits, x, y, a, N)

|x, y> -> |x, y * a^x mod N> for 0 <= y < N and gcd(a, N) = 1. With y = 1 this is the modular exponentiation |x, 1> -> |x, a^x mod N> of Shor's algorithm (Vedral et al., Eq. 1 and Sect. III.D). Basis states with y >= N are left unchanged.

Source code in dense_evolution/circuits/arithmetic.py

def power_mod(sv, n_qubits, x, y, a, N):
    """
    |x, y> -> |x, y * a^x mod N> for 0 <= y < N and gcd(a, N) = 1. With y = 1
    this is the modular exponentiation |x, 1> -> |x, a^x mod N> of Shor's
    algorithm (Vedral et al., Eq. 1 and Sect. III.D). Basis states with
    y >= N are left unchanged.
    """
    _check(n_qubits, x, y)
    _check_modulus(N, y)
    if np.gcd(int(a), N) != 1:
        raise ValueError(f"a={a} and N={N} must be coprime")
    idx, vx = _register_values(n_qubits, x)
    _, vy = _register_values(n_qubits, y)
    powers = np.array([pow(int(a), int(e), N) for e in range(1 << len(x))])
    ok = vy < N
    return _permute(sv, _with_register(n_qubits, idx, y, np.where(ok, (vy * powers[vx]) % N, vy)))

arithmetic_qasm

Reversible arithmetic as gate circuits, written as OPENQASM 2.0.

arithmetic.py applies each operation as a basis permutation |x> -> |f(x)>: exact, but not a circuit. The functions here return the gate-level circuits of the original papers as OPENQASM 2.0 text, to be read with QASMParser and run, noised, counted or exported like any other circuit. Each one is checked exhaustively against the matching permutation in arithmetic.py.

Registers are separate qreg declarations, least significant bit first (b[0] is the lowest bit), so the parser's flattened qubit list matches the register lists arithmetic.py takes.

References

Cuccaro, Draper, Kutin, Moulton, quant-ph/0410184 (ripple-carry adder). Draper, quant-ph/0008033 (addition on a quantum computer, QFT adder). Beauregard, quant-ph/0205095 (circuit for Shor's algorithm using 2n+3 qubits).

cuccaro_adder_qasm

cuccaro_adder_qasm(n)

Ripple-carry adder |a, b, 0> -> |a, a + b, 0> of Cuccaro et al. (quant-ph/0410184, Fig. 1-4): MAJ gates carry the sum up, UMA gates undo them and write each sum bit.

Registers: a[n], b[n+1] (b[n] receives the carry out, so the sum never overflows) and one ancilla x[1], initially 0 and returned to 0. 2n Toffoli gates and 4n + 1 CNOTs. Matches arithmetic.add_registers(sv, 2n + 2, a, b).

Source code in dense_evolution/circuits/arithmetic_qasm.py

def cuccaro_adder_qasm(n):
    """
    Ripple-carry adder |a, b, 0> -> |a, a + b, 0> of Cuccaro et al.
    (quant-ph/0410184, Fig. 1-4): MAJ gates carry the sum up, UMA gates
    undo them and write each sum bit.

    Registers: `a[n]`, `b[n+1]` (b[n] receives the carry out, so the sum
    never overflows) and one ancilla `x[1]`, initially 0 and returned to 0.
    2n Toffoli gates and 4n + 1 CNOTs. Matches
    `arithmetic.add_registers(sv, 2n + 2, a, b)`.
    """
    _check_n(n)
    lines = [f"qreg a[{n}];", f"qreg b[{n + 1}];", "qreg x[1];",
             "maj x[0], b[0], a[0];"]
    lines += [f"maj a[{i - 1}], b[{i}], a[{i}];" for i in range(1, n)]
    lines.append(f"cx a[{n - 1}], b[{n}];")
    lines += [f"uma a[{i - 1}], b[{i}], a[{i}];" for i in range(n - 1, 0, -1)]
    lines.append("uma x[0], b[0], a[0];")
    return _HEADER + _MAJ + _UMA + "\n".join(lines) + "\n"

draper_adder_qasm

draper_adder_qasm(n)

QFT adder |a, b> -> |a, a + b> of Draper (quant-ph/0008033): b is taken to Fourier space, each bit of a adds its phase with controlled rotations, and the inverse QFT brings the sum back. No ancilla, no carry gates.

Registers: a[n], b[n+1] (b[n] holds the carry, so the sum never overflows). Matches arithmetic.add_registers(sv, 2n + 1, a, b).

Source code in dense_evolution/circuits/arithmetic_qasm.py

def draper_adder_qasm(n):
    """
    QFT adder |a, b> -> |a, a + b> of Draper (quant-ph/0008033): b is
    taken to Fourier space, each bit of a adds its phase with controlled
    rotations, and the inverse QFT brings the sum back. No ancilla, no
    carry gates.

    Registers: `a[n]`, `b[n+1]` (b[n] holds the carry, so the sum never
    overflows). Matches `arithmetic.add_registers(sv, 2n + 1, a, b)`.
    """
    _check_n(n)
    m = n + 1
    lines = [f"qreg a[{n}];", f"qreg b[{m}];"] + _qft("b", m)
    for i in range(n):
        lines += _phi_add("b", m, 2 ** i, (f"a[{i}]",))
    lines += _iqft("b", m)
    return _HEADER + "\n".join(lines) + "\n"

constant_adder_qasm

constant_adder_qasm(c, n)

Adds a classical integer, |b> -> |(b + c) mod 2^n>: Beauregard's phiADD(c) (quant-ph/0205095, Fig. 3), one phase gate per qubit of b in Fourier space, between a QFT and an inverse QFT.

Register: b[n]. Matches arithmetic.add_constant(sv, n, b, c).

Source code in dense_evolution/circuits/arithmetic_qasm.py

def constant_adder_qasm(c, n):
    """
    Adds a classical integer, |b> -> |(b + c) mod 2^n>: Beauregard's
    phiADD(c) (quant-ph/0205095, Fig. 3), one phase gate per qubit of b in
    Fourier space, between a QFT and an inverse QFT.

    Register: `b[n]`. Matches `arithmetic.add_constant(sv, n, b, c)`.
    """
    _check_n(n)
    lines = [f"qreg b[{n}];"] + _qft("b", n) + _phi_add("b", n, int(c)) + _iqft("b", n)
    return _HEADER + "\n".join(lines) + "\n"

modular_constant_adder_qasm

modular_constant_adder_qasm(c, N, n)

Modular addition of a classical integer, |b, 0> -> |(b + c) mod N, 0> for 0 <= b < N and 0 <= c < N: Beauregard's phiADD(c)MOD(N) (quant-ph/0205095, Fig. 5, here without the two controls), between a QFT and an inverse QFT.

Registers: b[n+1] (one bit more than N needs, used as the sign of b + c - N) and one ancilla t[1], initially 0 and returned to 0. Requires N < 2^n. On inputs with b < N it matches arithmetic.add_constant_mod(sv, n + 2, b[:n], c, N); inputs with b >= N are outside the circuit's contract.

Source code in dense_evolution/circuits/arithmetic_qasm.py

def modular_constant_adder_qasm(c, N, n):
    """
    Modular addition of a classical integer, |b, 0> -> |(b + c) mod N, 0>
    for 0 <= b < N and 0 <= c < N: Beauregard's phiADD(c)MOD(N)
    (quant-ph/0205095, Fig. 5, here without the two controls), between a
    QFT and an inverse QFT.

    Registers: `b[n+1]` (one bit more than N needs, used as the sign of
    b + c - N) and one ancilla `t[1]`, initially 0 and returned to 0.
    Requires N < 2^n. On inputs with b < N it matches
    `arithmetic.add_constant_mod(sv, n + 2, b[:n], c, N)`; inputs with
    b >= N are outside the circuit's contract.
    """
    c, N = int(c), int(N)
    _check_modular(n, N)
    if not 0 <= c < N:
        raise ValueError(f"need 0 <= c < N, got c={c}, N={N}.")
    m = n + 1
    lines = [f"qreg b[{m}];", "qreg t[1];"] + _qft("b", m)
    lines += _phi_add_mod("b", m, c, N, "t[0]") + _iqft("b", m)
    return _HEADER + "\n".join(lines) + "\n"

cmult_mod_qasm

cmult_mod_qasm(a, N, n)

Controlled modular multiply-add, |c, x, b, 0> -> |c, x, (b + a*x) mod N, 0> when c = 1 and unchanged when c = 0, for 0 <= b < N: Beauregard's CMULT(a)MOD(N) (quant-ph/0205095, Fig. 6), n doubly controlled phiADD(2^i a mod N)MOD(N) gates (Fig. 5) controlled by c and x[i], between a QFT and an inverse QFT on b.

Registers: c[1], x[n], b[n+1] (b[n] is the sign bit, 0 on input and output) and one ancilla t[1], initially 0 and returned to 0. Requires N < 2^n. On inputs with b < N it matches arithmetic.multiply_add_mod(sv, 2n + 3, x, b[:n], a, N, controls=(c,)).

Source code in dense_evolution/circuits/arithmetic_qasm.py

def cmult_mod_qasm(a, N, n):
    """
    Controlled modular multiply-add, |c, x, b, 0> -> |c, x, (b + a*x) mod N, 0>
    when c = 1 and unchanged when c = 0, for 0 <= b < N: Beauregard's
    CMULT(a)MOD(N) (quant-ph/0205095, Fig. 6), n doubly controlled
    phiADD(2^i a mod N)MOD(N) gates (Fig. 5) controlled by c and x[i],
    between a QFT and an inverse QFT on b.

    Registers: `c[1]`, `x[n]`, `b[n+1]` (b[n] is the sign bit, 0 on input
    and output) and one ancilla `t[1]`, initially 0 and returned to 0.
    Requires N < 2^n. On inputs with b < N it matches
    `arithmetic.multiply_add_mod(sv, 2n + 3, x, b[:n], a, N, controls=(c,))`.
    """
    a, N = int(a), int(N)
    _check_modular(n, N)
    lines = ["qreg c[1];", f"qreg x[{n}];", f"qreg b[{n + 1}];", "qreg t[1];"]
    return _HEADER + "\n".join(lines + _cmult(a % N, N, n, "c[0]")) + "\n"

controlled_ua_qasm

controlled_ua_qasm(a, N, n)

Controlled modular multiplication in place, |c, x, 0, 0> -> |c, (ax) mod N, 0, 0> when c = 1, unchanged when c = 0, for 0 <= x < N and gcd(a, N) = 1: Beauregard's controlled-U_a (quant-ph/0205095, Fig. 7). CMULT(a)MOD(N) writes ax into b, a controlled swap exchanges x and b, and the inverse of CMULT(a^-1 mod N)MOD(N) clears b back to 0.

Registers as in cmult_mod_qasm: c[1], x[n], b[n+1], t[1], with b and t starting and ending at 0: 2n + 3 qubits, the count of the paper. Matches arithmetic.multiply_mod(sv, 2n + 3, x, a, N, controls=(c,)) on inputs with x < N and b = t = 0.

Source code in dense_evolution/circuits/arithmetic_qasm.py

def controlled_ua_qasm(a, N, n):
    """
    Controlled modular multiplication in place, |c, x, 0, 0> ->
    |c, (a*x) mod N, 0, 0> when c = 1, unchanged when c = 0, for 0 <= x < N
    and gcd(a, N) = 1: Beauregard's controlled-U_a (quant-ph/0205095, Fig. 7).
    CMULT(a)MOD(N) writes a*x into b, a controlled swap exchanges x and b,
    and the inverse of CMULT(a^-1 mod N)MOD(N) clears b back to 0.

    Registers as in `cmult_mod_qasm`: `c[1]`, `x[n]`, `b[n+1]`, `t[1]`, with
    b and t starting and ending at 0: 2n + 3 qubits, the count of the paper.
    Matches `arithmetic.multiply_mod(sv, 2n + 3, x, a, N, controls=(c,))` on
    inputs with x < N and b = t = 0.
    """
    a, N = int(a), int(N)
    _check_modular(n, N)
    if math.gcd(a, N) != 1:
        raise ValueError(f"a={a} has no inverse modulo N={N}.")
    swap = []
    for i in range(n):
        swap += [f"cx b[{i}], x[{i}];", f"ccx c[0], x[{i}], b[{i}];", f"cx b[{i}], x[{i}];"]
    lines = ["qreg c[1];", f"qreg x[{n}];", f"qreg b[{n + 1}];", "qreg t[1];"]
    lines += _cmult(a % N, N, n, "c[0]") + swap + _inverse(_cmult(pow(a, -1, N), N, n, "c[0]"))
    return _HEADER + "\n".join(lines) + "\n"

shor

Order finding for Shor's algorithm with one control qubit (Beauregard, quant-ph/0205095, Sect. 2.5, Fig. 8): 2n + 3 qubits in total.

The 2n control qubits of the textbook circuit are replaced by a single one, measured and reset after each controlled multiplication; the inverse QFT is done semiclassically, a phase gate on that qubit set by the outcomes already measured, then a Hadamard and the measurement. Each controlled-U_{a^(2^k)} is the gate circuit of controlled_ua_qasm.

shor_order_finding

shor_order_finding(a, N, rng=None, noise_model=None, p=0.0)

One run of quantum order finding for a modulo N, returns (y, m): a measured integer 0 <= y < 2^m with m = 2n, n = N.bit_length(). y / 2^m is close to s / r for the order r of a and a random s, so the continued-fraction expansion of y / 2^m recovers r with good probability (fractions.Fraction(y, 2**m).limit_denominator(N)).

Iteration k = m-1, ..., 0 prepares the control in |+>, applies controlled-U_{a^(2^k) mod N} to the work register (initially |1>), applies the phase -2 pi sum_l r_l / 2^(l-k+1) over the bits r_l already measured, a Hadamard, then measures and resets the control.

With noise_model (any model of NoiseModel.apply_to_sv, e.g. "depolarizing") and probability p, one stochastic trajectory of that channel acts on every qubit after each controlled multiplication, before the control is measured.

Source code in dense_evolution/circuits/shor.py

def shor_order_finding(a, N, rng=None, noise_model=None, p=0.0):
    """
    One run of quantum order finding for a modulo N, returns
    (y, m): a measured integer 0 <= y < 2^m with m = 2n, n = N.bit_length().
    y / 2^m is close to s / r for the order r of a and a random s, so the
    continued-fraction expansion of y / 2^m recovers r with good probability
    (`fractions.Fraction(y, 2**m).limit_denominator(N)`).

    Iteration k = m-1, ..., 0 prepares the control in |+>, applies
    controlled-U_{a^(2^k) mod N} to the work register (initially |1>),
    applies the phase -2 pi sum_l r_l / 2^(l-k+1) over the bits r_l already
    measured, a Hadamard, then measures and resets the control.

    With `noise_model` (any model of `NoiseModel.apply_to_sv`, e.g.
    "depolarizing") and probability `p`, one stochastic trajectory of that
    channel acts on every qubit after each controlled multiplication, before
    the control is measured.
    """
    from ..backends.statevector import DenseSVSimulator
    from ..noise import NoiseModel

    a, N = int(a), int(N)
    n = N.bit_length()
    m, nq = 2 * n, 2 * n + 3
    rng = np.random.default_rng(rng)
    sv = np.zeros((2, 2 ** (nq - 1)), dtype=np.complex128)
    sv[0, 1 << (nq - 2)] = 1.0
    bits = {}
    for k in range(m - 1, -1, -1):
        sv = _hadamard_on_first(sv)
        ak = pow(a, 2 ** k, N)
        if ak != 1:
            circ = QASMParser().parse(controlled_ua_qasm(ak, N, n))
            sim = DenseSVSimulator(nq)
            sim.set_initial_state(sv.reshape(-1))
            sim.run_circuit_jit(circ.to_tuples())
            sv = np.asarray(sim.get_statevector()).reshape(2, -1)
        if noise_model is not None and p > 0:
            sv = np.asarray(NoiseModel.apply_to_sv(sv.reshape(-1).copy(), nq, noise_model, p, rng=rng)).reshape(2, -1)
        phase = -2 * np.pi * sum(bits[l] / 2 ** (l - k + 1) for l in bits)
        sv[1] *= np.exp(1j * phase)
        sv = _hadamard_on_first(sv)
        p1 = float(np.sum(np.abs(sv[1]) ** 2))
        bits[k] = int(rng.random() < p1)
        kept = sv[bits[k]] / np.sqrt(p1 if bits[k] else 1.0 - p1)
        sv = np.stack([kept, np.zeros_like(kept)])
    y = sum(r << (m - 1 - k) for k, r in bits.items())
    return y, m