Doublecheck your inverse QFT circuit. I get a distribution similar to yours if I put QFT instead of its inverse in the second part of the phase estimation algorithm. If the correct inverse QFT is used, the distribution is very concentrated around 5.
Here is Qiskit code:
from qiskit import QuantumRegister, ClassicalRegister, QuantumCircuit
from numpy import pi
qreg_phi = QuantumRegister(1, 'phi')
qreg_q = QuantumRegister(5, 'q')
creg_c = ClassicalRegister(5, 'c')
circuit = QuantumCircuit(qreg_phi, qreg_q, creg_c)
theta = 5.039153255477287 * pi / 16
# Prepare |1> in qreg_phi
# |1> is the eigenvector of p(theta)
circuit.x(qreg_phi[0])
# Start phase estimation for U = p(theta)
# Hadamard
circuit.h(qreg_q[0])
circuit.h(qreg_q[1])
circuit.h(qreg_q[2])
circuit.h(qreg_q[3])
circuit.h(qreg_q[4])
# Apply controlled U's
circuit.cp(theta, qreg_q[0], qreg_phi[0])
circuit.cp(2 * theta, qreg_q[1], qreg_phi[0])
circuit.cp(4 * theta, qreg_q[2], qreg_phi[0])
circuit.cp(8 * theta, qreg_q[3], qreg_phi[0])
circuit.cp(16 * theta, qreg_q[4], qreg_phi[0])
# Uncompute QFT
circuit.h(qreg_q[4])
circuit.cp(-pi / 2, qreg_q[4], qreg_q[3])
circuit.h(qreg_q[3])
circuit.cp(-pi / 4, qreg_q[4], qreg_q[2])
circuit.cp(-pi / 2, qreg_q[3], qreg_q[2])
circuit.h(qreg_q[2])
circuit.cp(-pi / 8, qreg_q[4], qreg_q[1])
circuit.cp(-pi / 4, qreg_q[3], qreg_q[1])
circuit.cp(-pi / 2, qreg_q[2], qreg_q[1])
circuit.h(qreg_q[1])
circuit.cp(-pi / 16, qreg_q[4], qreg_q[0])
circuit.cp(-pi / 8, qreg_q[3], qreg_q[0])
circuit.cp(-pi / 4, qreg_q[2], qreg_q[0])
circuit.cp(-pi / 2, qreg_q[1], qreg_q[0])
circuit.h(qreg_q[0])
# Swap qubits
circuit.swap(qreg_q[0], qreg_q[4])
circuit.swap(qreg_q[1], qreg_q[3])
circuit.reset(qreg_phi)
circuit.measure(qreg_q, creg_c)
UPD. Here is an example with arbitrary number of qubits
def QPE(U, phi, n):
p = U.num_qubits
qreg_phi = QuantumRegister(p, 'phi')
qreg_q = QuantumRegister(n, 'q')
circuit = QuantumCircuit(qreg_phi, qreg_q)
circuit.initialize(phi, qreg_phi)
circuit.h(qreg_q)
for i in range(n):
circuit.append(U.power(2**i).control(1), [qreg_q[i], qreg_phi])
for i in range(n):
for j in range(i):
circuit.cp(-pi / 2**(i - j), qreg_q[n - 1 - j], qreg_q[n - 1 - i])
circuit.h(qreg_q[n - i - 1])
circuit.swap(qreg_q[:n//2], qreg_q[-1:-(n//2)-1:-1])
return circuit
theta = 5.039153255477287 * pi / 16
U = PhaseGate(theta)
n = 5
circuit = QPE(U, [0.0, 1.0], n)
qreg_q = circuit.qregs[1]
creg_x = ClassicalRegister(n, 'x')
circuit.add_register(creg_x)
circuit.measure(qreg_q, creg_x)