Project Locksport II, Shor's algorithm

Project Locksport II

Phase-estimation and factoring, Shor's algorithm 

ImageFX

https://github.com/barleyjohn/IBM-Quantum-Learning-Tutorials/blob/main/Phase-estimation%20and%20factoring.ipynb

test/study AI on Quantum Simulator https://blogbarley.blogspot.com/2024/08/ai-on-quantum-simulator-project.html

The phase estimation problem This section explains the phase estimation problem. We'll begin with a short discussion of the spectral theorem from linear algebra, and then move on to a statement of the phase estimation problem itself.

Spectral theorem The spectral theorem is an important fact from linear algebra that states that matrices of a certain type, called normal matrices, can be expressed in a simple and useful way. We'll only need this theorem for unitary matrices in this lesson, but in later lessons we'll apply it to Hermitian matrices as well.

https://learning.quantum.ibm.com/course/fundamentals-of-quantum-algorithms/phase-estimation-and-factoring

In [ ]:

from qiskit import __version__ print(__version__)

In [ ]:

from qiskit.quantum_info import Statevector from qiskit.quantum_info.operators import Operator from math import pi, cos, sin

In [ ]:

psi1 = Statevector([cos(pi / 8), sin(pi / 8)]) psi2 = Statevector([cos(5 * pi / 8), sin(5 * pi / 8)]) # When given a Statevector input, the Operator function returns the outer # product of that state vector with itself — or, in other words, the # product of the vector times its conjugate transpose. H = Operator(psi1) - Operator(psi2) display(H.draw("latex"))

In [ ]:

from qiskit import QuantumCircuit from qiskit_aer import AerSimulator from qiskit.visualization import plot_histogram from math import pi, cos, sin

In [ ]:

theta = 0.7 # Can be changed to any value between 0 and 1 qc = QuantumCircuit(2, 1) # Prepare the eigenvector, which is the |1> state qc.x(1) qc.barrier() # Implement the estimation procedure qc.h(0) qc.cp(2 * pi * theta, 0, 1) qc.h(0) qc.barrier() # Perform the final measurement qc.measure(0, 0) # Draw the circuit display(qc.draw(output="mpl"))

In [ ]:

result = AerSimulator().run(qc, shots = 1000).result() statistics = result.get_counts() display(plot_histogram(statistics))

In [ ]:

print(f"cos(pi * {theta})**2 = {cos(pi * theta) ** 2:.4f}") print(f"sin(pi * {theta})**2 = {sin(pi * theta) ** 2:.4f}")

theta = 0.7 qc = QuantumCircuit(3, 2)

Prepare the eigenvector

qc.x(2) qc.barrier()

The initial Hadamard gates

qc.h(0) qc.h(1) qc.barrier()

The controlled unitary gates

qc.cp(2 * pi * theta, 0, 2) qc.cp(2 * pi * (2 * theta), 1, 2) qc.barrier()

An implementation of the inverse of the two-qubit QFT

qc.swap(0, 1) qc.h(0) qc.cp(-pi / 2, 0, 1) qc.h(1) qc.barrier()

And finally the measurements

qc.measure([0, 1], [0, 1]) display(qc.draw(output="mpl"))

In [ ]:

result = AerSimulator().run(qc).result() statistics = result.get_counts() display(plot_histogram(statistics))

In [ ]:

from qiskit.circuit.library import QFT display(QFT(4).decompose().draw(output="mpl"))

In [ ]:

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister from qiskit_aer import AerSimulator from qiskit.circuit.library import QFT

In [ ]:

theta = 0.7 m = 3 # Number of control qubits control = QuantumRegister(m, name="X") target = QuantumRegister(1, name="Y") output = ClassicalRegister(m, name="Z") qc = QuantumCircuit(control, target, output) # Prepare the eigenvector qc.x(target) qc.barrier() # Perform phase estimation for index, qubit in enumerate(control): qc.h(qubit) for _ in range(2**index): qc.cp(2 * pi * theta, qubit, target) qc.barrier() # Perform the inverse quantum Fourier transform qc.compose( QFT(m, inverse=True).decompose(), inplace=True ) qc.barrier() # Measure everything qc.measure(range(m), range(m)) display(qc.draw(output="mpl"))

In [ ]:

result = AerSimulator().run(qc).result() statistics = result.get_counts() display(plot_histogram(statistics))

In [ ]:

most_probable = int(max(statistics, key=statistics.get),2) print(f"Most probable output: {most_probable}") print(f"Estimated theta: {most_probable/2**m}")

Shor's algorithm Now we'll turn our attention to the integer factorization problem, and see how it can be solved efficiently on a quantum computer using phase estimation. The algorithm we'll obtain is Shor's algorithm for integer factorization. Shor didn't describe his algorithm specifically in terms of phase estimation, but it is a natural and intuitive way to explain how it works. In the two subsections that follow we'll describe the two main parts of Shor's algorithm.

In the first subsection we'll define an intermediate problem known as the order-finding problem and see how phase estimation provides a solution to this problem.

In the second subsection we'll explain how an efficient solution to the order-finding problem also gives us an efficient solution to the integer factorization problem. When a solution to one problem provides a solution to another problem like this, we say that the second problem reduces to the first — so in this case we're reducing integer factorization to order finding. This part of Shor's algorithm doesn't make use of quantum computing at all, it's completely classical.

In [ ]:

from qiskit import QuantumCircuit, QuantumRegister, ClassicalRegister, transpile from qiskit.quantum_info.operators import Operator from qiskit.circuit.library import UnitaryGate from qiskit.circuit.library import QFT from qiskit_aer import AerSimulator import numpy as np from math import gcd, floor, log from fractions import Fraction import random

In [ ]:

def mod_mult_gate(b,N): if gcd(b,N)>1: print(f"Error: gcd({b},{N}) > 1") else: n = floor(log(N-1,2)) + 1 U = np.full((2**n,2**n),0) for x in range(N): U[b*x % N][x] = 1 for x in range(N,2**n): U[x][x] = 1 G = UnitaryGate(U) G.name = f"M_{b}" return G

In [ ]:

def order_finding_circuit(a,N): if gcd(a,N)>1: print(f"Error: gcd({a},{N}) > 1") else: n = floor(log(N-1,2)) + 1 m = 2*n control = QuantumRegister(m, name = "X") target = QuantumRegister(n, name = "Y") output = ClassicalRegister(m, name = "Z") circuit = QuantumCircuit(control, target, output) # Initialize the target register to the state |1> circuit.x(m) # Add the Hadamard gates and controlled versions of the # multiplication gates for k, qubit in enumerate(control): circuit.h(k) b = pow(a,2**k,N) circuit.compose( mod_mult_gate(b,N).control(), qubits = [qubit] + list(target), inplace=True) # Apply the inverse QFT to the control register circuit.compose( QFT(m, inverse=True), qubits=control, inplace=True) # Measure the control register circuit.measure(control, output) return circuit

In [ ]:

order_finding_circuit(2,11).draw(output = "mpl")

In [ ]:

def find_order(a,N): if gcd(a,N)>1: print(f"Error: gcd({a},{N}) > 1") else: n = floor(log(N-1,2)) + 1 m = 2*n circuit = order_finding_circuit(a,N) transpiled_circuit = transpile(circuit,AerSimulator()) while True: result = AerSimulator().run( transpiled_circuit, shots=1, memory=True).result() y = int(result.get_memory()[0],2) r = Fraction(y/2**m).limit_denominator(N).denominator if pow(a,r,N)==1: break return r

In [ ]:

N = 11 a = 2 print(f"The order of {a} modulo {N} is {find_order(a,N)}.")

In [ ]:

N = 39 FACTOR_FOUND = False # First we'll check to see if N is even or a nontrivial power. # Order finding won't help for factoring a *prime* power, but # we can easily find a nontrivial factor of *any* nontrivial # power, whether prime or not. if N % 2 == 0: print("Even number") d = 2 FACTOR_FOUND = True else: for k in range(2,round(log(N,2))+1): d = int(round(N ** (1/k))) if d**k == N: FACTOR_FOUND = True print("Number is a power") break # Now we'll iterate until a nontrivial factor of N is found. while not FACTOR_FOUND: a = random.randint(2,N-1) d = gcd(a,N) if d>1: FACTOR_FOUND = True print(f"Lucky guess of {a} modulo {N}") else: r = find_order(a,N) print(f"The order of {a} modulo {N} is {r}") if r % 2 == 0: x = pow(a,r//2,N) - 1 d = gcd(x,N) if d>1: FACTOR_FOUND = True print(f"Factor found: {d}")

ImageFX