-- Find r such that 7^r mod 15 = 1 -- See qiskit.org/textbook/ch-algorithms/shor.html -- 7 mod 15 (from qiskit code) U(psi,k) = rotate(psi, C,k,W,10,11, -- conditional swap bits 10 and 11 C,k,W,9,10, -- conditional swap bits 9 and 10 C,k,W,8,9, -- conditional swap bits 8 and 9 C,k,X,8, -- conditional not bit 8 C,k,X,9, -- conditional not bit 9 C,k,X,10, -- conditional not bit 10 C,k,X,11) -- conditional not bit 11 -- 12 quantum bits have 2^12 = 4096 eigenstates psi = zero(4096) -- initial state is eigenstate zero psi[1] = 1 -- start computing psi = rotate(psi, H,0, -- hadamard bit 0 H,1, -- hadamard bit 1 H,2, -- hadamard bit 2 H,3, -- hadamard bit 3 H,4, -- hadamard bit 4 H,5, -- hadamard bit 5 H,6, -- hadamard bit 6 H,7, -- hadamard bit 7 X,8) -- not bit 8 -- U^(2^0) U^(2^1) U^(2^2) ... U^(2^7) for(k,1,1, psi = U(psi,0)) for(k,1,2, psi = U(psi,1)) for(k,1,4, psi = U(psi,2)) for(k,1,8, psi = U(psi,3)) for(k,1,16, psi = U(psi,4)) for(k,1,32, psi = U(psi,5)) for(k,1,64, psi = U(psi,6)) for(k,1,128, psi = U(psi,7)) -- inverse QFT on bits 0-7 psi = rotate(psi,V,7) -- result is a probability distribution P = psi conj(psi) -- qubits 8-11 are don't care for(k,1,2048, P[k] = P[k] + P[k + 2048]) for(k,1,1024, P[k] = P[k] + P[k + 1024]) for(k,1, 512, P[k] = P[k] + P[k + 512]) for(k,1, 256, P[k] = P[k] + P[k + 256]) "Row 1 is the probability of observing the phase in row 2" n = 0 for(k,1,256,test(P[k],n=n+1)) test(n<2,n=2) M = zero(2,n) -- 2 rows, n columns n = 1 for(k,1,256,test(P[k],do(M[1,n]=float(P[k]),M[2,n]=(k-1)/256,n=n+1))) M "Observing 1/4 or 3/4 turns gives the correct answer r = 4"