from qsystem import QSystem, Gate
from random import randint
from math import log2, ceil, floor, gcd
from IPython.display import display, Latex, Math
n = 15
Latex('Factoring $n = {}$'.format(n))
a = 0
while gcd(n, a) != 1 or a == 1:
a = randint(2, n)
Math('a = {}'.format(a))
s = ceil(log2(n+1))
Math('s = {}'.format(s))
and a quantum oracle 'POWN' that $$\left|x\right>\left|0\right> \xrightarrow{\text{POWN}} \left|x\right>\left|a^x (\text{mod}\,n)\right>$$
def pown(x):
x = x >> s
fx = pow(a, x, n)
return (x << s) | fx
def it():
for x in range(2**s):
yield x << s
pown_gate = Gate.from_func(pown, 2*s, it())
seed = randint(0,10000)
q = QSystem(s, seed)
print('init first register at')
print(q)
init first register at +1.000000000 |0000>
q.evol(gate='H', qbit=0, count=s)
print(q)
+1/sqrt(16) |0101> +1/sqrt(16) |0100> +1/sqrt(16) |0111> +1/sqrt(16) |0110> +1/sqrt(16) |0001> +1/sqrt(16) |0000> +1/sqrt(16) |0011> +1/sqrt(16) |0010> +1/sqrt(16) |1101> +1/sqrt(16) |1100> +1/sqrt(16) |1111> +1/sqrt(16) |1110> +1/sqrt(16) |1001> +1/sqrt(16) |1000> +1/sqrt(16) |1011> +1/sqrt(16) |1010>
print('add the second register')
q.add_ancillas(s)
q.apply(gate=pown_gate, qbit=0)
print(q)
add the second register +1/sqrt(16) |0101>|1000> +1/sqrt(16) |0100>|0001> +1/sqrt(16) |0110>|0100> +1/sqrt(16) |0000>|0001> +1/sqrt(16) |0010>|0100> +1/sqrt(16) |0111>|0010> +1/sqrt(16) |1101>|1000> +1/sqrt(16) |1100>|0001> +1/sqrt(16) |1111>|0010> +1/sqrt(16) |0001>|1000> +1/sqrt(16) |1110>|0100> +1/sqrt(16) |0011>|0010> +1/sqrt(16) |1001>|1000> +1/sqrt(16) |1000>|0001> +1/sqrt(16) |1011>|0010> +1/sqrt(16) |1010>|0100>
help to understand the algorithm
$$ {1\over\sqrt{2^s}}\sum_{x=0}^{2^{s}-1} \left|x\right>\left|a^x(\text{mod}\, n)\right> \xrightarrow{\text{measure}[s:2s]} \sqrt{r\over{2^s}}\sum_{i=0}^{{2^{s}\over r}-1} \left|ir + x_0\right>\left|a^{x_0}(\text{mod}\, n)\right> $$print('measure and remove the second register')
q.measure(s, s)
q.rm_ancillas()
print(q)
measure and remove the second register +1/sqrt(4) |1001> +1/sqrt(4) |0001> +1/sqrt(4) |1101> +1/sqrt(4) |0101>
q.qft(qbit_begin=0, qbit_end=s)
print(q)
-i/sqrt(4) |1100> -1/sqrt(4) |1000> +i/sqrt(4) |0100> +1/sqrt(4) |0000>
q.measure_all()
c = q.bits()
c = sum([m*2**i for m, i in zip(c, reversed(range(len(c))))])
mea = [c]
for _ in range(s-1):
seed = randint(0,10000)
q = QSystem(s, seed)
q.evol('H', 0, s) # 1
q.add_ancillas(s)
q.apply(pown_gate, 0) # 2
q.measure(s, s)
q.rm_ancillas()
q.qft(0, s) # 4
q.measure_all() # 5
c = q.bits()
c = sum([m*2**i for m, i in zip(c, reversed(range(len(c))))])
mea.append(c)
print('measurements results =', mea)
measurements results = [12, 0, 4, 12]
c = mea[0]
for m in mea:
c = gcd(c, m)
if c == 0:
print('repite the period-finding algorithm')
else:
r = 2**s/c
display(Latex('possible period $r = {}$'.format(r)))
if c == 0:
print('repite the period-finding algorithm')
elif r % 2 == 1:
print('go to step 1')
else:
p = gcd(int(a**(r/2)+1), n)
q = gcd(int(a**(r/2)-1), n)
display(Math('{}*{}={}'.format(p,q,p*q)))