qft and qpe

example/grover.py

@@ -4,11 +4,11 @@
 def oracle(v, target):
     n = q.num_qubits(v)
     for i in range(n):
-        if (target >> i) & 1 == 0:
+        if (target >> (n - i - 1)) & 1 == 0:
             v = q.apply(v, q.op.x(), [i])
     v = q.apply_seq(v, q.op.controlled_seq(q.op.z(), list(range(n))))
     for i in range(n):
-        if (target >> i) & 1 == 0:
+        if (target >> (n - i - 1)) & 1 == 0:
             v = q.apply(v, q.op.x(), [i])
     return v
 
@@ -40,7 +40,7 @@ def grover(n, target, iter, use_mps):
 
 if __name__ == "__main__":
     n = 4
-    target = 0b0110
+    target = q.util.str2binary("0111")
     iter = 2 ** (n // 2)
     print("# number of qbits = {}".format(n))
     print("# target state = {}".format(q.util.binary2str(n, target)))

example/qpe.py

@@ -0,0 +1,36 @@
+from copy import deepcopy
+
+import numpy as np
+import qailo as q
+
+
+def qpe(n, u, ev):
+    m = q.num_qubits(u)
+    assert q.num_qubits(ev) == n + m
+    v = deepcopy(ev)
+
+    for p in range(n):
+        v = q.apply(v, q.op.h(), [p])
+
+    cp = q.op.controlled(u)
+    rep = 1
+    for p in range(n):
+        for _ in range(rep):
+            # print(f"apply cu on {p} and {list(range(n,n+m))}")
+            v = q.apply(v, cp, [p] + list(range(n, n + m)))
+        rep = rep * 2
+
+    v = q.apply_seq(v, q.alg.qft.inverse_qft_seq(n))
+    return v
+
+
+if __name__ == "__main__":
+    n = 3
+    phi = 2 * np.pi * (1 / 3)
+    u = q.op.p(phi)
+    ev = q.sv.state_vector(n + 1)
+    ev = q.apply(ev, q.op.x(), [n])
+    v = qpe(n, u, ev)
+    prob = np.diag(q.op.matrix(q.op.trace(q.sv.pure_state(v), [n])).real)
+    for i in range(len(prob)):
+        print("{} {}".format(q.util.binary2str(n, i), prob[i]))

example/test_grover.py

@@ -5,7 +5,7 @@
 
 def test_simple():
     n = 4
-    target = 0b0000
+    target = q.util.str2binary("0000")
     iter = 2 ** (n // 2)
     prob = q.probability(grover(n, target, iter, False))
     assert prob[0] == approx(0.581704139709473)

example/test_qpe.py

@@ -0,0 +1,16 @@
+import numpy as np
+import qailo as q
+import qpe
+from pytest import approx
+
+
+def test_qpe():
+    n = 3
+    phi = 2 * np.pi * (1 / 8)
+    u = q.op.p(phi)
+    ev = q.sv.state_vector(n + 1)
+    ev = q.apply(ev, q.op.x(), [n])
+    v = qpe.qpe(n, u, ev)
+    prob = np.diag(q.op.matrix(q.op.trace(q.sv.pure_state(v), [n])).real)
+    assert prob[4] == approx(1)
+    assert prob[0] == approx(0)

src/qailo/__init__.py

@@ -1,10 +1,11 @@
-from . import mps, operator, state_vector, util
+from . import alg, mps, operator, state_vector, util
 from . import operator as op
 from . import state_vector as sv
 from ._version import version
 from .dispatch import apply, apply_seq, num_qubits, probability, vector
 
 __all__ = [
+    alg,
     mps,
     operator,
     state_vector,

src/qailo/alg/__init__.py

@@ -0,0 +1,6 @@
+from .qft import inverse_qft_seq, qft_seq
+
+__all__ = [
+    inverse_qft_seq,
+    qft_seq,
+]

src/qailo/alg/qft.py

@@ -0,0 +1,36 @@
+# ref: https://learn.qiskit.org/course/ch-algorithms/quantum-fourier-transform
+
+import numpy as np
+
+import qailo as q
+
+
+def qft_rotations_seq(n):
+    seq = []
+    if n == 0:
+        return seq
+    n -= 1
+    # print(f"H on [{n}]")
+    seq.append([q.op.h(), [n]])
+    for p in range(n):
+        # print(f"CP(pi/{2**(n-p)} on [{p}, {n}]")
+        seq.append([q.op.cp(np.pi / 2 ** (n - p)), [p, n]])
+    seq += qft_rotations_seq(n)
+    return seq
+
+
+def swap_registers_seq(n):
+    seq = []
+    for p in range(n // 2):
+        # print(f"swap on [{p}, {n-p-1}]")
+        seq.append([q.op.swap(), [p, n - p - 1]])
+    return seq
+
+
+def qft_seq(n):
+    """QFT on the first n qubits in circuit"""
+    return qft_rotations_seq(n) + swap_registers_seq(n)
+
+
+def inverse_qft_seq(n):
+    return q.op.inverse_seq(qft_seq(n))

src/qailo/operator/__init__.py

@@ -10,6 +10,7 @@
 )
 from .hconj import hconj
 from .identity import identity
+from .inverse import inverse_seq
 from .matrix import matrix
 from .multiply import multiply
 from .one_qubit import h, p, rx, ry, rz, s, t, x, y, z
@@ -39,6 +40,7 @@
     z,
     hconj,
     identity,
+    inverse_seq,
     is_close,
     is_hermitian,
     is_identity,

src/qailo/operator/controlled.py

@@ -9,7 +9,7 @@ def controlled(u):
     assert is_operator(u)
     m = num_qubits(u)
     n = m + 1
-    op = np.identity(2**n).reshape([2, 2**m, 2, 2**m])
+    op = np.identity(2**n, dtype=u.dtype).reshape([2, 2**m, 2, 2**m])
     op[1, :, 1, :] = matrix(u)
     return op.reshape((2,) * (2 * n))
 

src/qailo/operator/inverse.py

@@ -0,0 +1,10 @@
+from .hconj import hconj
+from .type import is_operator
+
+
+def inverse_seq(seq):
+    res = []
+    for p, pos in seq[::-1]:
+        assert is_operator(p)
+        res.append([hconj(p), pos])
+    return res

src/qailo/operator/trace.py

@@ -16,4 +16,4 @@ def trace(q, pos=None):
     ss = list(range(2 * n))
     for i in pos:
         ss[n + i] = ss[i]
-    return np.einsum(q, ss)
+    return np.einsum(q.reshape((2,) * (2 * n)), ss)

src/qailo/util/strops.py

@@ -1,5 +1,5 @@
 def binary2str(n, i):
-    return bin(i)[2:].zfill(n)
+    return bin(i)[2:].zfill(n)[::-1]
 
 
 def str2binary(s):
@@ -8,5 +8,5 @@ def str2binary(s):
     for i in range(n):
         assert s[i] == "0" or s[i] == "1"
         if s[i] == "1":
-            c += 2 ** (n - i - 1)
+            c += 2**i
     return c

test/alg/test_qft.py

@@ -0,0 +1,32 @@
+import numpy as np
+import qailo as q
+from pytest import approx
+
+
+def test_qft():
+    n = 3
+    target = 5
+
+    # generate result of qft of target
+    v = q.sv.state_vector(n)
+    for p in range(n):
+        v = q.apply(v, q.op.h(), [p])
+    v = q.apply(v, q.op.p(target * np.pi / 4), [0])
+    v = q.apply(v, q.op.p(target * np.pi / 2), [1])
+    v = q.apply(v, q.op.p(target * np.pi), [2])
+
+    # apply inverse qft
+    v = q.apply_seq(v, q.alg.inverse_qft_seq(n))
+
+    v = q.vector(v)
+    for i in range(2**n):
+        if i == target:
+            print("* {} {}".format(q.util.binary2str(n, i), v[i]))
+            assert v[i] == approx(1)
+        else:
+            print("  {} {}".format(q.util.binary2str(n, i), v[i]))
+            assert v[i] == approx(0)
+
+
+if __name__ == "__main__":
+    test_qft()

test/util/test_binary2str.py

@@ -2,10 +2,14 @@
 
 
 def test_binary2str():
-    assert binary2str(4, 0b1011) == "1011"
-    assert binary2str(6, 0b1110) == "001110"
+    assert binary2str(4, 0b1011) == "1011"[::-1]
+    assert binary2str(6, 0b1110) == "001110"[::-1]
 
     assert str2binary("0000") == 0
-    assert str2binary("0010") == 2
+    assert str2binary("0010") == 4
     assert str2binary("0110") == 6
-    assert str2binary("1000") == 8
+    assert str2binary("1000") == 1
+
+    n = 6
+    for i in range(2**n):
+        assert str2binary(binary2str(n, i)) == i