shor's Algorithm

from qiskit import QuantumCircuit, Aer, execute
from qiskit.visualization import plot_histogram
from math import gcd
import numpy as np

def qft_dagger(qc, n):
    """Apply the inverse QFT to the first n qubits in the circuit."""
    for qubit in range(n//2):
        qc.swap(qubit, n-qubit-1)
    for j in range(n):
        for m in range(j):
            qc.cp(-np.pi/float(2**(j-m)), m, j)
        qc.h(j)

def qpe_amod15(a):
    """Quantum Phase Estimation for finding the order of a modulo 15."""
    n_count = 4
    qc = QuantumCircuit(n_count + 4, n_count)
    
    # Initialize counting qubits in state |+>
    for q in range(n_count):
        qc.h(q)
    
    # Initialize auxiliary register in state |1>
    qc.x(3+n_count)
    
    # Apply the controlled-U operations
    for q in range(n_count):
        qc.append(c_amod15(a, 2**q), 
                  [q] + [i+n_count for i in range(4)])
    
    # Apply inverse QFT
    qft_dagger(qc, n_count)
    
    # Measure
    qc.measure(range(n_count), range(n_count))
    
    # Simulate the circuit
    backend = Aer.get_backend('qasm_simulator')
    counts = execute(qc, backend, shots=1024).result().get_counts()
    return counts

def c_amod15(a, power):
    """Controlled multiplication by a modulo 15."""
    U = QuantumCircuit(4)
    for iteration in range(power):
        if a in [2, 13]:
            U.swap(0, 1)
            U.swap(1, 2)
            U.swap(2, 3)
        if a in [4, 11]:
            U.swap(0, 2)
            U.swap(1, 3)
        if a in [7, 8]:
            U.swap(0, 3)
            U.swap(1, 2)
            U.swap(2, 3)
        if a in [13, 11, 8]:
            U.x(range(4))
    U = U.to_gate()
    U.name = f"{a}^{power} mod 15"
    c_U = U.control()
    return c_U

def shor(N):
    # Choose a random integer a such that 1 < a < N and gcd(a, N) = 1
    a = np.random.randint(2, N)
    while gcd(a, N) != 1:
        a = np.random.randint(2, N)

    # Perform Quantum Phase Estimation
    counts = qpe_amod15(a)

    # Analyze the results
    measurement = max(counts, key=counts.get)
    phase = int(measurement, 2) / (2**len(measurement))

    # Find r
    r = 1 / phase
    r = round(r)

    # Check if r is odd or if a^(r/2) ≡ -1 (mod N)
    if r % 2 == 1 or pow(a, r // 2, N) == N - 1:
        return None

    # Compute the factors using the Euclidean algorithm
    factor1 = gcd(pow(a, r // 2) + 1, N)
    factor2 = gcd(pow(a, r // 2) - 1, N)

    return factor1, factor2

# Test the Shor's algorithm
N = 15
factors = shor(N)
if factors is not None:
    print(f"The factors of {N} are {factors[0]} and {factors[1]}")
else:
    print("Failed to find the factors")

You can also try this code with Online Python Compiler

Run Code

#include <iostream>
#include <cmath>
#include <cstdlib>
#include <ctime>

using namespace std;

int gcd(int a, int b) {
    if (b == 0)
        return a;
    return gcd(b, a % b);
}

int power(int a, int b, int mod) {
    int result = 1;
    while (b > 0) {
        if (b & 1)
            result = (result * a) % mod;
        a = (a * a) % mod;
        b >>= 1;
    }
    return result;
}

int findOrder(int a, int N) {
    int r = 1;
    int powerMod = power(a, r, N);
    while (powerMod != 1 && r < N) {
        r++;
        powerMod = power(a, r, N);
    }
    return (powerMod == 1) ? r : -1;
}

int shor(int N) {
    srand(time(NULL));

    // Choose a random integer a such that 1 < a < N and gcd(a, N) = 1
    int a = rand() % (N - 2) + 2;
    while (gcd(a, N) != 1)
        a = rand() % (N - 2) + 2;

    // Find the order r of a modulo N
    int r = findOrder(a, N);

    // Check if we failed to find a useful r
    if (r == -1 || r % 2 == 1 || power(a, r / 2, N) == N - 1)
        return -1;

    // Compute the factors using the Euclidean algorithm
    int factor1 = gcd(power(a, r / 2, N) + 1, N);
    int factor2 = gcd(power(a, r / 2, N) - 1, N);

    return (factor1 == 1 || factor2 == 1) ? -1 : factor1;
}

int main() {
    int N = 15;
    int factor = shor(N);

    if (factor == -1)
        cout << "Failed to find the factors" << endl;
    else
        cout << "The factors of " << N << " are " << factor << " and " << N / factor << endl;

    return 0;
}

You can also try this code with Online C++ Compiler

Run Code

import java.util.Random;

public class Shor {
    public static int gcd(int a, int b) {
        if (b == 0)
            return a;
        return gcd(b, a % b);
    }

    public static int power(int a, int b, int mod) {
        int result = 1;
        while (b > 0) {
            if ((b & 1) != 0)
                result = (result * a) % mod;
            a = (a * a) % mod;
            b >>= 1;
        }
        return result;
    }

    public static int shor(int N) {
        Random rand = new Random();

        // Choose a random integer a such that 1 < a < N and gcd(a, N) = 1
        int a = rand.nextInt(N - 2) + 2;
        while (gcd(a, N) != 1)
            a = rand.nextInt(N - 2) + 2;

        // Find the order r of a modulo N
        int r = 1;
        while (power(a, r, N) != 1 && r < N) {
            r++;
        }

        // Check if we failed to find a useful r
        if (r == N || r % 2 == 1 || power(a, r / 2, N) == N - 1)
            return -1;

        // Compute the factors using the Euclidean algorithm
        int factor1 = gcd(power(a, r / 2, N) + 1, N);
        int factor2 = gcd(power(a, r / 2, N) - 1, N);

        return (factor1 == 1 || factor2 == 1) ? -1 : factor1;
    }

    public static void main(String[] args) {
        int N = 15;
        int factor = shor(N);

        if (factor == -1)
            System.out.println("Failed to find the factors");
        else
            System.out.println("The factors of " + N + " are " + factor + " and " + (N / factor));
    }
}

You can also try this code with Online Java Compiler

Run Code