Level 04 · Quantum Computing

Oracles and Search

Phase oracles, Deutsch–Jozsa, Bernstein–Vazirani and Grover search.

Level 4 · Quantum Computing: Oracles and Search

Your simulator can now hold a register of qubits, run gates on any of them and measure. This level uses it for the first algorithms that beat every classical program at something. The setting is a black box: a function f on bit strings that you may ask but never read, and the cost of an algorithm is the number of questions, or queries, it asks. One quantum query tells a constant function from a balanced one, or reads a hidden string of n bits, where a classical program needs many. Grover's search finds a marked item among N in about \frac{\pi}{4}\sqrt N queries instead of about N/2.

The first half builds the tools: the two ways to turn f into a gate (the phase oracle and the bit oracle), phase kickback, which turns one into the other, the Hadamard gate on every qubit, and the sealed box, an oracle you may use once. Then come Deutsch–Jozsa and Bernstein–Vazirani, with a glance at Simon. The second half is Grover's search in small steps: inversion about the mean, one round worked by hand, the success probability round by round, the picture in a plane that explains it, how many rounds to run, and what to do when a run misses.

The other paths meet these ideas. The Algorithms path counts comparisons in the worst case; here queries are counted the same way, and that path's bits, XOR and linear search are the raw material. The Hadamard gate on every qubit is computed with the same sum-and-difference butterflies as the fast Fourier transform of the Signals & Control path. Expected numbers of attempts and seeded random draws come from the Statistics path. numpy is available, and the examples use it for small full matrices.

A function you may only ask

A Boolean function f takes n bits and answers 0 or 1. The inputs where it answers 1 are its marked inputs. Written out for every input, it is a truth table truth of length 2^n: truth[x] is f(x), with the input read as a binary number and qubit 0 as the leftmost, most significant bit, as in Level 3. In 3 bits, |011\rangle is index 3. Here f marks the 3-bit strings with exactly two 1s:

n = 3
def f(bits):                      # bits is a string like "011", qubit 0 first
    return 1 if bits.count("1") == 2 else 0

truth = [f(format(x, f"0{n}b")) for x in range(2 ** n)]
print(truth)
print([format(x, f"0{n}b") for x in range(2 ** n) if truth[x]])
x = 6
print(format(x, "03b"), [(x >> (n - 1 - q)) & 1 for q in range(n)])
[0, 0, 0, 1, 0, 1, 1, 0]
['011', '101', '110']
110 [1, 1, 0]

format(x, "03b") keeps the leading zeros that bin(x) drops, and bit q of x is (x >> (n - 1 - q)) & 1 because qubit 0 sits highest. A query is one evaluation of f. A problem hands you a truth table when you are meant to simulate a circuit, and a sealed box when you are meant to ask it questions.

Common mistakes: (x >> q) & 1, which puts qubit 0 on the right. bin(x)[2:], which loses the leading zeros. Confusing the input x (a number) with its bit string.

The phase oracle

The simplest way to put f into a circuit leaves every basis state where it is and only changes its sign:

O_f\,|x\rangle = (-1)^{f(x)}\,|x\rangle

So O_f is a diagonal matrix, -1 for each marked input and +1 elsewhere. Applied to the even spread of 3 qubits, with the function above:

import numpy as np
truth = [0, 0, 0, 1, 0, 1, 1, 0]
O = np.diag([(-1) ** t for t in truth])
print(np.diag(O).tolist())
state = np.full(8, 1 / np.sqrt(8))
after = O @ state
print(np.round(after * np.sqrt(8), 4).tolist())
print(np.allclose(np.abs(after) ** 2, np.abs(state) ** 2), np.allclose(O @ O, np.eye(8)))
[1, 1, 1, -1, 1, -1, -1, 1]
[1.0, 1.0, 1.0, -1.0, 1.0, -1.0, -1.0, 1.0]
True True

The marked inputs 3, 5 and 6 carry a minus sign, and no probability changed: |-a|^2 = |a|^2. On a basis state, O_f gives back the same state times \pm1, a global phase (Level 2). On a superposition it writes relative signs, which a measurement right after cannot see. Every algorithm of this level is a way to turn those signs into something measurable, by interference. O_f is unitary and its own inverse: a second application removes every sign.

Common mistakes: putting 0 on the unmarked inputs instead of +1, which erases their amplitudes. Flipping the bits of the marked inputs (that is X, not a sign). Expecting a measurement straight after the oracle to show the marks.

Practice: Minus Signs on the Marked Rows

The bit oracle

A gate must be reversible, but f usually is not: the OR of two bits answers 1 for 01, 10 and 11, so the answer alone does not say which input it came from. The standard fix keeps the input and XORs the answer into one extra qubit, the target y:

U_f\,|x, y\rangle = |x,\ y\oplus f(x)\rangle \qquad (\oplus = \text{XOR})

The input register is qubits 0 to n-1 and the target is the last qubit, qubit n. So |x, y\rangle is index 2 * x + y: the target is the lowest bit. A target that starts at 0 ends as f(x), one that starts at 1 ends as 1 - f(x), and applying U_f twice gives back |x, y\rangle, since y \oplus f \oplus f = y. Here is the matrix for the OR of two bits, built column by column: column 2 * x + y has its single 1 in row 2 * x + (y ^ f(x)).

import numpy as np
def bit_oracle_matrix(truth):
    size = 2 * len(truth)
    U = np.zeros((size, size))
    for x, fx in enumerate(truth):
        for y in (0, 1):
            U[2 * x + (y ^ fx), 2 * x + y] = 1
    return U

U = bit_oracle_matrix([0, 1, 1, 1])          # f(x) = x0 OR x1
print([format(i, "03b") + "->" + format(int(np.argmax(U[:, i])), "03b") for i in range(8)])
print(np.allclose(U.T @ U, np.eye(8)), np.allclose(U @ U, np.eye(8)))
['000->000', '001->001', '010->011', '011->010', '100->101', '101->100', '110->111', '111->110']
True True

Read each string as register, then target: 100->101 is |10, 0\rangle \mapsto |10, 1\rangle, because f(10) = 1. Inputs with f(x) = 1 swap their two target states; inputs with f(x) = 0 stay put. Every row and every column holds a single 1: a permutation matrix, which is always unitary. On a superposition, U_f computes nothing new; it only moves each amplitude to its new index.

Aside (optional): where do oracles come from? On real hardware an oracle is not a magic part. It is compiled from an ordinary circuit for f, rebuilt from reversible gates such as CNOT and the Toffoli gate (which flips a target when two controls are both 1), usually with spare scratch qubits that must be returned to 0 afterwards. For a real search, such as testing a candidate key, that circuit is long, and every gate adds a little error. Counting queries hides that cost, which is one reason the speed-ups of this level are easier to state than to collect.

Common mistakes: putting the target first (index y * 2 ** n + x). Writing f(x) into the target instead of XOR-ing it in. Reading the register's bits from the right.

Practice: The Wire That Remembers an Answer, Is This Wiring an Oracle?

Phase kickback

Put the target in |-\rangle = (|0\rangle - |1\rangle)/\sqrt2 and follow one basis input |x\rangle by hand. If f(x) = 0, nothing happens. If f(x) = 1, the two parts of the target swap:

|x\rangle\,\tfrac{1}{\sqrt2}\big(|1\rangle - |0\rangle\big) = -\,|x\rangle\,\tfrac{1}{\sqrt2}\big(|0\rangle - |1\rangle\big)

Both cases in one line:

U_f\,|x\rangle|-\rangle = (-1)^{f(x)}\,|x\rangle|-\rangle

The target comes out unchanged and the answer lands on the register, as a sign. That is phase kickback. In matrix terms, the oracle applies X to the target exactly when f(x) = 1, and |-\rangle is an eigenvector of X with eigenvalue -1. With |+\rangle, an eigenvector with eigenvalue +1, there is nothing to kick back. The simulator agrees, for the OR function:

import numpy as np
def bit_oracle_matrix(truth):
    size = 2 * len(truth)
    U = np.zeros((size, size))
    for x, fx in enumerate(truth):
        for y in (0, 1):
            U[2 * x + (y ^ fx), 2 * x + y] = 1
    return U

r = 1 / np.sqrt(2)
targets = {"-": np.array([r, -r]), "+": np.array([r, r])}
truth = [0, 1, 1, 1]
U = bit_oracle_matrix(truth)
for x in range(4):
    e = np.zeros(4)
    e[x] = 1
    signs = []
    for name, t in targets.items():
        before = np.kron(e, t)
        after = U @ before
        signs.append(name + ("  same" if np.allclose(after, before) else " minus" if np.allclose(after, -before) else " other"))
    print(format(x, "02b"), "f =", truth[x], signs)
00 f = 0 ['-  same', '+  same']
01 f = 1 ['- minus', '+  same']
10 f = 1 ['- minus', '+  same']
11 f = 1 ['- minus', '+  same']

Because U_f is linear, a whole register \sum_x \psi_x|x\rangle next to a target |-\rangle becomes \sum_x (-1)^{f(x)}\psi_x|x\rangle, still next to an unchanged |-\rangle: the bit oracle has become the phase oracle on the register. The product state |\psi\rangle\otimes|-\rangle (Level 3's Kronecker product) has amplitude \psi_x/\sqrt2 at index 2 * x and -\psi_x/\sqrt2 at 2 * x + 1. A constant f = 1 multiplies everything by -1, a global phase that nothing can detect.

Common mistakes: a target in |1\rangle instead of |-\rangle: on a superposition the register then gets entangled with the target instead of picking up clean signs. Expecting the target to change. Forgetting the minus sign on the odd indices when building |\psi\rangle\otimes|-\rangle.

Practice: The Answer Lands on the Question

A Hadamard on every wire

Almost every algorithm of this level starts and ends with H on every qubit, written H^{\otimes n}. On |0\dots0\rangle it gives the even spread |s\rangle, every amplitude 1/\sqrt{2^n}. On any basis state |x\rangle it gives an even spread with signs:

H^{\otimes n}|x\rangle = \frac{1}{\sqrt{2^n}}\sum_{z=0}^{2^n-1}(-1)^{x\cdot z}\,|z\rangle, \qquad x\cdot z = \text{the number of positions where } x \text{ and } z \text{ both have a } 1

In Python x\cdot z is bin(x & z).count("1"). The reason: for one bit, H|b\rangle = (|0\rangle + (-1)^b|1\rangle)/\sqrt2, and the Kronecker product multiplies one such factor per qubit, so a minus sign appears once for every position where x and z both have a 1. Checked for |101\rangle with the full matrix:

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
H3 = np.kron(np.kron(H, H), H)
x = 5                                                  # |101>
formula = [(-1) ** bin(x & z).count("1") / np.sqrt(8) for z in range(8)]
print(np.round(H3[:, x] * np.sqrt(8), 4).tolist())
print(np.allclose(H3[:, x], formula), np.allclose(H3 @ H3, np.eye(8)))
print(np.round(H3[:, 0] * np.sqrt(8), 4).tolist())
[1.0, -1.0, 1.0, -1.0, -1.0, 1.0, -1.0, 1.0]
True True
[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

Column x of H^{\otimes n} is the signed spread of |x\rangle, and column 0 is all plus. H^{\otimes n} is its own inverse, so it also turns the signed spread of x back into |x\rangle. Every qubit gets the same gate, so the bit order does not matter here.

The full matrix has 4^n entries, 268 million at 14 qubits. Level 3 showed the way out: a gate on one qubit only mixes pairs of indices that differ in that qubit's bit. For H on qubit q of n, with step = 2 ** (n - 1 - q) (Level 3's bit mask), each pair i, i + step (where i has that bit 0) holding a, b becomes (a+b)/\sqrt2, (a-b)/\sqrt2: the sum-and-difference butterfly of the FFT. Once per qubit, that costs n\cdot2^n operations instead of 4^n. The pairs are visible inside the full matrix of H on qubit 0 of 3, where step is 4:

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
I = np.eye(2)
H0 = np.kron(np.kron(H, I), I)                         # H on qubit 0 of 3
for row in (1, 5):
    nonzero = [j for j in range(8) if abs(H0[row, j]) > 1e-12]
    print("row", row, "uses columns", nonzero, "with weights", np.round(H0[row, nonzero] * np.sqrt(2), 4).tolist())
row 1 uses columns [1, 5] with weights [1.0, 1.0]
row 5 uses columns [1, 5] with weights [1.0, -1.0]

New amplitude 1 is (a_1 + a_5)/\sqrt2 and new amplitude 5 is (a_1 - a_5)/\sqrt2; nothing else enters.

Common mistakes: (-1) ** (x * z), the product of the two numbers instead of the shared 1 bits. Forgetting the factor 1/\sqrt{2^n}, or applying 1/\sqrt2 at every butterfly and 1/\sqrt{2^n} again at the end. Overwriting v[i] before it is used for v[i + step].

Practice: A Hadamard on Every Wire

Sealed boxes: one question, one reading

Some problems give you a sealed box instead of a truth table: hardware you may use once. The pattern is always the same.

  1. You prepare a state yourself, a normalised list of amplitudes. Gates you want before the oracle you apply yourself, with your simulator, and hand over the result.
  2. The box applies its oracle once, then a short list of gates you choose, such as ("h", q), ("x", q) or ("z", q). These are the only gates after the oracle.
  3. The box measures every qubit and returns only the outcome, a bit string with qubit 0 first. You never see an amplitude.
  4. A second use raises TooManyQueries.

The whole skill is choosing the state and the gates so that the characters you need are certain. Here is a toy box holding the bit oracle of the OR function, built with full matrices and Level 2's seeded draw. First a classical question, then a second use, then a register in |01\rangle with the target in |-\rangle, for six different seeds:

import numpy as np, random

class TooManyQueries(Exception):
    pass

def make_box(truth, seed):
    n = len(truth).bit_length() - 1                    # input bits; the box has n + 1 qubits
    size = 2 * len(truth)
    U = np.zeros((size, size))
    for x, fx in enumerate(truth):
        for y in (0, 1):
            U[2 * x + (y ^ fx), 2 * x + y] = 1
    gate = {"h": np.array([[1, 1], [1, -1]]) / np.sqrt(2), "x": np.array([[0, 1], [1, 0]]), "z": np.diag([1, -1])}
    rng, used = random.Random(seed), [False]
    def box(state, gates=()):
        if used[0]:
            raise TooManyQueries("this box answers only once")
        used[0] = True
        v = U @ np.array(state, dtype=complex)
        for name, q in gates:
            v = np.kron(np.kron(np.eye(2 ** q), gate[name]), np.eye(2 ** (n - q))) @ v
        p = np.abs(v) ** 2
        r, total = rng.random(), 0.0
        for i in range(size):
            total += p[i]
            if total > r:
                return format(i, f"0{n + 1}b")
        return format(max(i for i in range(size) if p[i] > 1e-12), f"0{n + 1}b")
    return box

box = make_box([0, 1, 1, 1], seed=1)
state = [0] * 8
state[4] = 1                                           # |10, 0>
print(box(state))
try:
    box(state)
except TooManyQueries as e:
    print("second use:", e)

r = 1 / np.sqrt(2)
readings = []
for seed in range(6):
    state = [0] * 8
    state[2], state[3] = r, -r                         # |01> next to |->
    readings.append(make_box([0, 1, 1, 1], seed)(state))
print(readings, [s[:2] for s in readings])
101
second use: this box answers only once
['011', '010', '011', '010', '010', '011'] ['01', '01', '01', '01', '01', '01']

|10, 0\rangle came back as "101": register 10, target f(10) = 1. That is one classical query and one bit of information. With the target in |-\rangle the register characters are 01 for every seed, while the last character is random, because |-\rangle has weight \tfrac12 on each of 0 and 1. So with a target in |-\rangle, read only the register's characters. A box that holds a phase oracle has no target: it acts on n qubits and every character is a register bit.

Common mistakes: comparing the whole string, target included. Calling the box a second time to double-check. Handing over a list that is not normalised. Reversing the string: qubit 0 is the first character.

Deutsch–Jozsa: constant or balanced

The promise: f on n bits is either constant (the same answer on every input) or balanced (1 on exactly half of the 2^n inputs, 0 on the other half). Which is it?

Classically, to be certain, you may need 2^{n-1} + 1 queries. Picture an adversary who answers 0 as long as possible: after 2^{n-1} zeros, f could still be constant 0, or balanced with all its ones among the inputs you have not asked. One more answer settles it, because a balanced f has only 2^{n-1} zeros. For n = 10 that is 513 queries; for n = 1 it is both f(0) and f(1) (that case is Deutsch's problem).

Quantumly, one query is enough:

  1. the register in the even spread, the target in |-\rangle;
  2. the oracle, which by kickback writes (-1)^{f(x)} on each |x\rangle;
  3. H on every register qubit;
  4. measure the register.

Why it works: after step 2 the register is \frac{1}{\sqrt{2^n}}\sum_x(-1)^{f(x)}|x\rangle. Row 0 of H^{\otimes n} is all +\frac{1}{\sqrt{2^n}}, so after step 3 the amplitude of |0\dots0\rangle is

\frac{1}{2^n}\sum_x(-1)^{f(x)},

the average of the signs. For a constant f the average is +1 or -1, and the register reads all zeros for sure. For a balanced f the plus and minus signs cancel exactly, and it never does. Here is the register on its own (kickback already applied) for three functions on 3 bits:

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
H3 = np.kron(np.kron(H, H), H)
spread = np.full(8, 1 / np.sqrt(8))
for name, truth in [("constant 1", [1] * 8), ("parity", [0, 1, 1, 0, 1, 0, 0, 1]), ("balanced", [1, 0, 0, 1, 1, 0, 1, 0])]:
    signs = np.array([(-1) ** t for t in truth])
    final = H3 @ (signs * spread)
    a0 = round(float(final[0]), 4) + 0.0
    print(f"{name:10}  average sign {float(signs.mean()):+.2f}  amplitude of 000 {a0:+.4f}  P(000) {a0 ** 2:.4f}")
constant 1  average sign -1.00  amplitude of 000 -1.0000  P(000) 1.0000
parity      average sign +0.00  amplitude of 000 +0.0000  P(000) 0.0000
balanced    average sign +0.00  amplitude of 000 +0.0000  P(000) 0.0000

The balanced functions send the register to other outcomes, never to 000. In the full circuit there are n + 1 qubits and the target stays in |-\rangle, so its character is random; only the first n characters carry the answer.

Common mistakes: reading the target's character too. Skipping the final Hadamards: the signs alone leave every outcome with probability 1/2^n. Preparing the target in |0\rangle or |1\rangle. Testing for "some zero" instead of "all zeros".

Practice: One Look at a Sealed Switch, Constant or Balanced, One Question Only, How Often Does the Register Read All Zeros?

Bernstein–Vazirani, and Simon briefly

Now the promise is that f is a parity: for a hidden string s,

f(x) = s\cdot x \bmod 2 = (s_0x_0 + s_1x_1 + \dots + s_{n-1}x_{n-1}) \bmod 2,

which is bin(s & x).count("1") % 2. Classically, a query reveals at most one bit (asking x = 100\dots0 reveals s_0), so finding s takes n queries. Quantumly, send the even spread through the phase oracle: the amplitude of |x\rangle becomes (-1)^{s\cdot x}/\sqrt{2^n}. Compare with the formula for H^{\otimes n}: that is exactly the signed spread H^{\otimes n}|s\rangle, column s of the matrix. For s = 110:

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
H3 = np.kron(np.kron(H, H), H)
s = int("110", 2)
signed = np.array([(-1) ** (bin(s & x).count("1") % 2) for x in range(8)]) / np.sqrt(8)
print(np.round(signed * np.sqrt(8), 4).tolist())
print(np.allclose(signed, H3[:, s]))
[1.0, 1.0, -1.0, -1.0, -1.0, -1.0, 1.0, 1.0]
True

And H^{\otimes n} undoes itself. The bits of s are numbered like qubits, s_0 first.

Simon's problem (optional) goes one step further. A function hides s so that f(x) = f(x') exactly when x' = x or x' = x \oplus s. Each run of Simon's circuit ends with a random string y with y\cdot s = 0 \bmod 2, and the quantum part is then done: what remains is linear algebra mod 2. Each y is one equation in the unknown bits of s. Store it as an integer; adding two equations mod 2 is XOR, and the sum is again an equation that every solution satisfies:

ys = ["1010", "0110"]
rows = [int(y, 2) for y in ys]
def solves(z):
    return all(bin(row & z).count("1") % 2 == 0 for row in rows)
print([format(z, "04b") for z in range(16) if solves(z)])
combo = rows[0] ^ rows[1]
print(format(combo, "04b"), all(bin(combo & z).count("1") % 2 == 0 for z in range(16) if solves(z)))
['0000', '0001', '1110', '1111']
1100 True

Two equations in four unknown bits leave four candidates. Trying all 2^n candidates is fine for 4 bits and hopeless for 40; Gaussian elimination does better, with XOR to subtract rows and nothing to divide by.

Common mistakes: computing s\cdot x as the product of two numbers. Forgetting the "mod 2". Returning s reversed, or as a number without its leading zeros.

Practice: The Mask Behind the Parity Lock, Twin Inputs in the Lab Notebook

Search: signs are not enough

Grover's problem: among N = 2^n items, M are marked by f, and nothing else is known. Classically you check items one by one; with one marked item that takes about N/2 queries on average.

Grover's search starts from the even spread |s\rangle, every amplitude 1/\sqrt N (made by H on every qubit), where a measurement finds a marked item with probability M/N. The phase oracle flips the signs of the marked amplitudes, which changes no probability. The running example of this half: N = 64, one marked item, number 42.

import numpy as np
N, marked = 64, [42]
a = np.full(N, 1 / np.sqrt(N))
b = np.array([-x if i in marked else x for i, x in enumerate(a)])
print(float(a[42]), float(b[42]), bool(a[42] ** 2 == b[42] ** 2))
print(float(a.mean()), round(float(b.mean()), 8))
0.125 -0.125 True
0.125 0.12109375

What did change is the mean: one minus sign pulls it down from 0.125 to 0.12109375. The marked item now sits far below the mean, and the next step turns that distance into size.

Inversion about the mean

Reflect every entry through the mean m, so that an entry d below the mean ends up d above it:

a_i' = 2m - a_i, \qquad m = \frac1N\sum_{i=0}^{N-1}a_i

For the list 2, -1, 4, 3 the mean is 2: the 2 stays, -1 (3 below the mean) becomes 5, 4 becomes 0 and 3 becomes 1. The mean is unchanged, and so is the sum of squares, 4 + 1 + 16 + 9 = 30 = 4 + 25 + 0 + 1. That is no accident: the operation is the unitary matrix

D = 2\,|s\rangle\langle s| - I, \qquad D_{ij} = \frac{2}{N} - \delta_{ij}

where \delta_{ij} is 1 when i = j and 0 otherwise. Why: \langle s|a\rangle = \frac{1}{\sqrt N}\sum_i a_i, so every entry of 2|s\rangle\langle s|a\rangle is 2\cdot\frac1N\sum_i a_i = 2m, and -I subtracts a.

import numpy as np
amps = np.array([2.0, -1.0, 4.0, 3.0])
N = len(amps)
s = np.full(N, 1 / np.sqrt(N))
D = 2 * np.outer(s, s) - np.eye(N)
print(np.round(D, 4).tolist())
out = D @ amps
print(np.round(out, 4).tolist(), float(amps.mean()), round(float(out.mean()), 4))
print(float((amps ** 2).sum()), round(float((out ** 2).sum()), 4), np.allclose(D @ D, np.eye(N)))
[[-0.5, 0.5, 0.5, 0.5], [0.5, -0.5, 0.5, 0.5], [0.5, 0.5, -0.5, 0.5], [0.5, 0.5, 0.5, -0.5]]
[2.0, 5.0, 0.0, 1.0] 2.0 2.0
30.0 30.0 True

D leaves |s\rangle alone (D|s\rangle = 2|s\rangle - |s\rangle = |s\rangle) and flips every vector orthogonal to |s\rangle, that is, every list whose mean is 0. So D is a reflection, and applying it twice changes nothing. It is called the diffusion step, or inversion about the mean.

Common mistakes: a - 2 * m, which gets every sign wrong. Recomputing the mean after some entries have already changed. Using the mean of the absolute values or of the squares.

Practice: Flip Every Number Around the Average

One round, worked by hand

One Grover round is the oracle, then the inversion, with the mean taken after the flip. Back to N = 64 with item 42 marked:

  • start: every amplitude 1/8 = 0.125;
  • oracle: item 42 becomes -0.125, and the mean becomes \frac{63\cdot0.125 - 0.125}{64} = 0.12109375;
  • inversion: every unmarked 0.125 becomes 2m - 0.125 = 0.1171875, and the marked -0.125 becomes 2m + 0.125 = 0.3671875.

The marked amplitude almost tripled, and its probability rose from 1/64 \approx 0.0156 to 0.3671875^2 \approx 0.1348. The 63 others shrank a little each, 63\cdot0.1171875^2 \approx 0.8652, and the total is still 1. The same round with full matrices, and the two steps in the wrong order for contrast:

import numpy as np
N, marked = 64, [42]
s = np.full(N, 1 / np.sqrt(N))
O = np.diag([-1.0 if x in marked else 1.0 for x in range(N)])
D = 2 * np.outer(s, s) - np.eye(N)
after = D @ (O @ s)                          # oracle first, then the inversion
print(round(float(after[42]), 4), round(float(after[0]), 4), round(float(after[42] ** 2), 4))
wrong = O @ (D @ s)                          # inversion first, then the oracle
print(round(float(wrong[42]), 4), round(float(wrong[42] ** 2), 4))
0.3672 0.1172 0.1348
-0.125 0.0156

Inverting first does nothing to the even spread, since D|s\rangle = |s\rangle, so the wrong order is just the oracle again: the probability stays 1/64. Note the order of the matrices: the round applies O first, so its matrix is D\,O.

Common mistakes: inverting before marking. Taking the mean before the flip. Marking with -abs(a) instead of multiplying by -1: an amplitude that is already negative must turn positive.

Practice: Mark, Then Flip Around the Average

Round after round

Repeat the round and watch the probability of a marked item, the success probability p_k after k rounds: the sum of |a_x|^2 over the marked x. With N = 64 the matrices are small, and one round is the fixed matrix G = D\,O. Press Run to see p_k for k = 0 to 25:

import numpy as np
N, marked = 64, [42]
s = np.full(N, 1 / np.sqrt(N))
O = np.diag([-1.0 if x in marked else 1.0 for x in range(N)])
D = 2 * np.outer(s, s) - np.eye(N)
G = D @ O                                    # one round as one matrix
state, probs = s.copy(), []
for k in range(26):
    probs.append(float(sum(state[x] ** 2 for x in marked)))
    state = G @ state
print([round(p, 4) for p in probs[:9]])
print(probs.index(max(probs)), round(max(probs), 4), [round(p, 4) for p in probs[12:14]])
plot(range(26), probs, kind="stem", title="N = 64, one marked item", xlabel="rounds k", ylabel="success probability")
[0.0156, 0.1348, 0.3439, 0.5914, 0.8164, 0.9635, 0.9966, 0.9074, 0.718]
6 0.9966 [0.0001, 0.0576]

The probability climbs to 0.9966 at k = 6, falls to almost nothing at k = 12 (0.0001), and climbs again. More rounds are not safer: past the peak every round makes things worse. The amplitudes stay real throughout, since every matrix involved is real. The curve is a squared sine, and the next section shows why.

Common mistakes: forgetting p_0, the probability before any round. Recording after the round only, which shifts the whole curve by one. Adding up the probability of only the first marked item when several are marked.

Practice: Watching the Odds Climb, Round by Round

The picture in a plane

Both steps treat all marked items alike, and all unmarked ones alike. So the state is always a combination of two fixed unit vectors:

  • |w\rangle, the even spread over the M marked items (each 1/\sqrt M, the others 0);
  • |r\rangle, the even spread over the N - M unmarked items.

They are orthogonal, and the start is |s\rangle = \sin\theta\,|w\rangle + \cos\theta\,|r\rangle with

\sin\theta = \sqrt{M/N}, \qquad \theta = \arcsin\sqrt{M/N}.

Draw |r\rangle as the horizontal axis and |w\rangle as the vertical one: the start is the unit arrow at angle \theta above the horizontal. Now the two steps:

  • the oracle flips the |w\rangle part: a reflection in the horizontal axis;
  • D keeps the |s\rangle part and flips the part orthogonal to it: a reflection in the line through |s\rangle.

Two reflections in mirrors that meet at angle \theta make a rotation by 2\theta. So every round turns the arrow by 2\theta, after k rounds its angle is (2k+1)\theta, and the success probability is its squared height:

p_k = \sin^2\big((2k+1)\,\theta\big)

For N = 64, M = 1: \theta = \arcsin(1/8) \approx 0.1253 rad, about 7.2 degrees. The formula gives the curve you just simulated, and the arrow walks round the circle 0.25 rad at a time:

import math
N, M = 64, 1
theta = math.asin(math.sqrt(M / N))
print(round(theta, 4), round(math.degrees(theta), 2))
print([round(math.sin((2 * k + 1) * theta) ** 2, 4) for k in range(9)])
angles = [(2 * k + 1) * theta for k in range(13)]
plot([math.cos(t) for t in angles], [math.sin(t) for t in angles], kind="scatter", title="the state after k = 0..12 rounds", xlabel="unmarked part", ylabel="marked part")
0.1253 7.18
[0.0156, 0.1348, 0.3439, 0.5914, 0.8164, 0.9635, 0.9966, 0.9074, 0.718]

The best moment is when the angle reaches \pi/2: (2k+1)\theta = \pi/2, or k = \frac{\pi}{4\theta} - \frac12. When few items are marked, \theta \approx \sqrt{M/N}, which gives the well-known rule of thumb

k \approx \left\lfloor \frac{\pi}{4}\sqrt{N/M} \right\rfloor

For N = 64 it says 6, where the peak is. The number of queries grows like \sqrt N: a million items with one marked need about 785 rounds, not half a million checks.

Common mistakes: math.sin((2 * k + 1) * theta ** 2) or forgetting to square the sine. Using M / N without the square root. Degrees where math.sin expects radians.

Practice: A Trillion Items Without a Simulator, A Billion Amplitudes in Two Numbers

How many rounds?

The rule of thumb rounds twice, once in \theta \approx \sqrt{M/N} and once in the floor, and it can miss by one. The angle moves 2\theta per round and lands on either side of \pi/2; which side is closer decides. Take N = 512 and M = 35:

import math
def curve(N, M, K):
    theta = math.asin(math.sqrt(M / N))
    return [round(math.sin((2 * k + 1) * theta) ** 2, 4) for k in range(K + 1)]

print(curve(512, 35, 4), math.floor(math.pi / 4 * math.sqrt(512 / 35)))
print(curve(8, 6, 2), math.floor(math.pi / 4 * math.sqrt(8 / 6)))
print(curve(64, 38, 2), math.floor(math.pi / 4 * math.sqrt(64 / 38)))
print(curve(16, 8, 3))
[0.0684, 0.5082, 0.9397, 0.9231, 0.4754] 3
[0.75, 0.0, 0.75] 0
[0.5937, 0.2319, 0.9047] 1
[0.5, 0.5, 0.5, 0.5]

For N = 512, M = 35 the rule says 3 rounds, but the peak is at 2. Many marked items are a second trap. When M > N/2 the start angle is already above \pi/4, so the first round overshoots: with 6 of 8 items marked, \theta = \pi/3, p_0 = 0.75 and one round takes the success to 0. The best plan there is to measure at once, k = 0. Between M/N = \tfrac12 and about 0.62 the rule still says one round, which makes things worse: for N = 64, M = 38, p_0 = 38/64 \approx 0.594 but p_1 \approx 0.232. With exactly half marked, \theta = \pi/4 and every p_k is \tfrac12: rounds change nothing, and the smallest k is as good as any.

Common mistakes: trusting the rule when many items are marked. Forgetting that k = 0, measuring at once, is a plan too. round instead of math.floor in the rule.

Practice: How Many Rounds Before You Look?, Stirred Too Long

When a run misses: retries and one generator

Even at the peak, p_k is rarely exactly 1, and the measurement can miss. A miss costs only time: check the measured item with one more, classical query, and if it is not marked, run the whole search again from a fresh start state. If one attempt succeeds with probability p, independently of the others, the number of attempts until the first success has mean 1/p (the geometric distribution of the Statistics path). A seeded simulation agrees:

import random
rng = random.Random(5)
p = 0.25
counts = []
for trial in range(20000):
    attempts = 1
    while rng.random() >= p:
        attempts += 1
    counts.append(attempts)
print(round(sum(counts) / len(counts), 2), 1 / p)
4.02 4.0

With p = 0.25 that is 4 attempts on average. Each attempt costs its k rounds plus the check.

Simulated measurements must be reproducible, so they use Level 2's seeded draw: r = rng.random(), then the first outcome whose running total of probabilities is greater than r. When a search is repeated, create the generator once for the whole run and keep drawing from it. A fresh random.Random(seed) in every attempt restarts the same sequence, so every attempt draws the same r, measures the same outcome, and misses forever:

import random
print([round(random.Random(3).random(), 4) for attempt in range(4)])
rng = random.Random(3)
print([round(rng.random(), 4) for attempt in range(4)])
[0.238, 0.238, 0.238, 0.238]
[0.238, 0.5442, 0.37, 0.6039]

Common mistakes: creating the generator inside the retry loop. Continuing the next attempt from the old state instead of a fresh even spread. Reading the answer off the amplitudes, which a real chip never shows you.

Practice: Try, Check, Try Again, The Vault That Answers in Signs

The inversion as a circuit

A chip cannot apply the N\times N matrix D in one go; it runs gates. Two facts from the Hadamard section give a short circuit: H^{\otimes n}|0\dots0\rangle = |s\rangle, and H^{\otimes n} is its own inverse. So

D = 2|s\rangle\langle s| - I = H^{\otimes n}\,\big(2|0\dots0\rangle\langle0\dots0| - I\big)\,H^{\otimes n}.

The middle matrix is diagonal: +1 for |0\dots0\rangle and -1 for every other basis state. Checked for 3 qubits:

import numpy as np
H = np.array([[1, 1], [1, -1]]) / np.sqrt(2)
H3 = np.kron(np.kron(H, H), H)
s = np.full(8, 1 / np.sqrt(8))
D = 2 * np.outer(s, s) - np.eye(8)
middle = np.diag([1.0] + [-1.0] * 7)
print(np.allclose(H3 @ middle @ H3, D), np.round(D[0], 4).tolist())
True [-0.75, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25]

The middle matrix is -1 times "flip the sign of |0\dots0\rangle only", and a factor -1 on a whole circuit is a global phase (Level 2), so either one will do. Flipping the sign of a single basis state is the job of a multi-controlled Z, which flips only |1\dots1\rangle; an X on every qubit swaps |0\dots0\rangle with |1\dots1\rangle. Put together, about 4n gates replace a matrix with 4^n entries. On a chip, then, Grover's search is: Hadamards on every qubit, rounds of (oracle circuit, diffusion circuit), and one measurement.

Common mistakes: Hadamards on one side only. Rejecting a circuit whose matrix is -D: it acts on every state like D, up to a global phase. Flipping |1\dots1\rangle without the X gates around it, which reflects about the wrong state.

Practice: The Average Flip, Built From Gates

Practice

These problems are required on this path at this level; the path adds one problem from each of the level's topics, picked for you and kept.

Checklist

You can now:

  • Write a Boolean function as a truth table, with qubit 0 as the leftmost bit of the input.
  • Build the phase oracle O_f|x\rangle = (-1)^{f(x)}|x\rangle and the bit oracle U_f|x, y\rangle = |x, y\oplus f(x)\rangle, with the target as the last qubit at index 2 * x + y.
  • Explain phase kickback: a target in |-\rangle turns the bit oracle into the phase oracle on the register, and a target in |+\rangle changes nothing.
  • Use H^{\otimes n}|x\rangle = 2^{-n/2}\sum_z(-1)^{x\cdot z}|z\rangle, and apply H qubit by qubit to pairs of amplitudes instead of building the full matrix.
  • Plan a one-query experiment with a sealed box: the state you prepare, the gates after the oracle, and which characters of the outcome to read.
  • Explain why Deutsch–Jozsa's all-zeros amplitude is the average of the signs, why Bernstein–Vazirani reads s in one query, and why the classical tester needs 2^{n-1} + 1 queries.
  • Run Grover rounds, oracle first and then the inversion about the mean taken after the flip, and record the success probability round by round.
  • Explain each round as a turn by 2\theta in a plane, use p_k = \sin^2((2k+1)\theta) and the rule \lfloor\frac{\pi}{4}\sqrt{N/M}\rfloor, and spot when the rule is wrong.
  • Repeat a failed search with one seeded generator for the whole run, expecting 1/p attempts on average.