Quantum Simulators and AQS

This is the last post in the series on quantum computing. You can check the previous ones here:

In the previous post, we discussed the applications of quantum algorithms that are in development and we touched upon the idea of testing these algorithm ideas. In the final post of this series, we will discuss about the ways that quantum programs can be analyzed before running in actual quantum computers.

Status Quo of Quantum Computing

Real-world Quantum Computers have been undergoing a huge boom during the last few years with hundreds of businesses and tens of countries investing in this technology. Due to the expensive R&D cost for Quantum Computing, the main business model for Quantum Computing has been through the provision of Cloud Computing services for access to this architecture. This model provides programming languages, tools, and documentation for implementing quantum computing applications in various ways.

Quantum Computer Development by Qubit Count – Yoke Development

Some of the big players in this space are Google with its Sycamore 53-qubit transmon Quantum Computer: they claimed Quantum Supremacy in 2019; IBM with its latest Osprey Chip with 433 transmon qubits: for which they claimed 0.999 fidelity; and QWave with its Advantage 5760 annealing qubit Quantum Computer used mainly for optimization problems. Other cloud providers are Microsoft through Azure, Amazon with AWS, and Xanadu with Xanadu Cloud.

The fast development of Quantum Computing has also started a shift in the perspective of its impact on society and technology. For instance, last year the U.S. Department of Commerce’s National Institute of Standards and Technology (NIST) started a competition and selected some post-quantum cryptographic algorithms that are resistant to quantum methods of attack in response to the concern that quantum computers may be able to break current cryptographic algorithms and pose a risk to the government and individuals. Similarly, many businesses have started integrating quantum computing to speed up certain computations or improve their efficiency. It is expected that the use of quantum computing will continue to grow in the near future with integration into many areas of research and development, but it is still very unlikely to be a replacement for current classical computing.

Quantum Simulators

Quantum computer simulators are software tools that are designed to simulate the behavior of quantum computers. They allow researchers to explore and understand the behavior of quantum systems without the need for a physical quantum computer. Quantum computer simulators work by using classical computers to simulate the behavior of quantum systems. These simulators use complex algorithms to simulate the behavior of quantum bits (qubits) and the interactions between them. Because quantum systems are highly sensitive to external factors, quantum computer simulators must take into account various sources of noise and error that can impact the accuracy of the simulation. Despite these challenges, quantum computer simulators are a valuable tool for researchers who want to explore the potential of quantum computing and develop new quantum algorithms.

The ArrayFire Quantum Simulator AQS

ArrayFire provides a tensor-based Quantum Computer Simulator that can be accelerated through GPU computations called the ArrayFire Quantum Simulator or AQS for short. AQS is a C++ library that provides the functionality to create, manipulate, visualize, and simulate quantum circuits with quick and accurate results. It contains implementations of common quantum algorithms to speed up the implementation of more complex programs that utilize them and ways to display such quantum circuits. It supports up to 30 qubits which is plenty for testing many applications.

To see how you can test a Quantum Algorithm with the ArrayFire Quantum Simulator, we can go through the example of simulating Shor’s Algorithm.

Simulating Shor’s Algorithm

Currently, the majority of internet communications are protected by RSA encryption. RSA is an asymmetric encryption method that takes advantage of the fact that the process of multiplying two large prime numbers is a much easier computation than factoring the product to recover the two prime numbers, i.e., it is an asymmetric computation.

The encryption consists of the two prime numbers $p, q$ , then the only visible information to other is their product $N = p q$ . If we want to break the encryption, our goal is to factor $N$ and find $p$ and $q$ . The first idea that might come to mind is to try multiple primes and check if they work. While this brute force approach works well for small prime numbers, for prime numbers on the order of $10^{25}$ it would take thousands of years for a supercomputer to compute it. Smarter methods still take a significant time that is not feasible with Classical Computers; however, Shor’s algorithm finds a way to do it much more efficiently on a Quantum Computer.

Shor’s Algorithm starts from the observation that the remainder of powers of a number $a^x \bmod N$ is a periodic function. For example, for $a=3$ and $N=7$

Here we can see that the pattern for the remainder of powers is $1,3,2,6,4,5$ which then repeats. So every 6 powers of $r$ , the remainder of the powers repeat, which means the period $r$ is 6. Therefore, any power that can be written as $a^{(rt + s)}$ for some values $t, s$ . The interesting property of this operation is that $a^r \equiv 1 \bmod N$ can be rearranged to $a^r - 1 \equiv 0 \bmod N$ . If r is an even number, we can factor out this expression:

$(a^{\frac{r}{2}} + 1)(a^{\frac{r}{2}} - 1) \equiv 0 \bmod N$

From this equation we can see that if $a$ is chosen such that $r$ is even, you can get two possible factors that are reasonable guesses for the primes $p$ and $q$ . Usually, these factors are not prime numbers, they are more likely multiples of $p$ or $q$ , so this process can be repeated by choosing a guess $a$ , finding the periodicity $r$ , and computing better guesses for the prime factors. The idea of converting the random guess $a$ into a better guess is the crux of Shor’s Algorithm.

The Quantum Computing part comes into computing the periodicity $r$ . Computing $r$ is an $O(N)$ operation on a classical computer, which is more computationally expensive than some brute force methods for finding the prime factors. However, with Quantum computers we can take advantage of the Quantum Fourier Transform to find the periodicity.

Implementation

We can implement this algorithm in the ArrayFire Quantum Simulator to test it. You can check the full example in the AQS Repository. You can learn more about how to use AQS in the documentation in the AQS Repo.

To start we first include the libraries we will use: the Arrayfire Quantum Simulator libraries, the vector and string data structures, the sort algorithm, and some console stream output. We first start by initializing the AQS library:

#include "quantum.h"
#include "quantum_visuals.h"
#include "quantum_algo.h"

#include <algorithm>
#include <iostream>
#include <string>
#include <vector>

int main(int argc, char** argv)
{
    aqs::initialize(argc, argv);
    //...
}

Let’s take the problem where our prime numbers are $p = 3, q=5$ so $N=15$ ; however we only know $N$ and we want to find out the numbers $p, q$ . For our initial guess, we choose $a=13$ :

int a = 13;
int N = 15;

The first thing we have to do is implement a way to compute the modulus/remainder of numbers through the quantum computer. We know that the binary states in a quantum computer can be represented as follows $|0000\rangle= |0\rangle, |1000\rangle = |1\rangle, |0100\rangle = |2\rangle \dots$ . Then our unitary gate $U$ must have the behavior $U|x\rangle = |ax \text{mod} N\rangle$ . This can be expressed by the following mapping:

Note that for $N=15$ , it is sufficient to analyze the numbers from $0$ to $15$ to find its factors, so 4 qubits suffice for this gate. It turns out, this gate can be implemented by a set of SWAP and X gates:

aqs::QCircuit U(4);

U << aqs::Swap(0, 1);
U << aqs::Swap(1, 2);
U << aqs::Swap(2, 3);

U << aqs::X(0);
U << aqs::X(1);
U << aqs::X(2);
U << aqs::X(3);

There are other ways to implement this gate $U$ for which you can vary $a$ or $N$ but it gets much more complicated. For now, we will use this implementation for our problem values. (If you want to test for other $a$ values, you can check them in the example code in the lambda function $\texttt{a2jmod15}$ ).

We need to apply the $U$ operation to the complete space of possible values, while also entangling the resulting operation by calculating $r$ . To achieve this, we can explore the whole set of possibilities by entangling the states of the modulus operation with some states we can measure. For this simulation, we will take the first 8 qubits for deducing some information of $r$ , while the other 4 qubits will be used as ancilla qubits for the $U$ gate to be applied.

The entanglement is done by first applying Hadamard gates on all the first 8 qubits which put those qubits into superposition. Then, Control Gates of variations of the gate $U$ are added which forces the Quantum Computer to explore all the possible combinations of applying $U$ gates. It turns out that the amount of information obtained from the circuit can be maximized by making the variations of the $U$ gate be powers of 2 of the $U$ gate, i.e. $U^{2^j}$ . All this process is encapsulated in the following code:

aqs::QCircuit qc(8 + 4);

for (uint32_t i = 0; i < 8; ++i)
    qc << aqs::H{i};

for (uint32_t i = 0; i < 8; ++i)
{
    aqs::QCircuit U_powerj(4);
    U_powerj << aqs::Gate(U, 0);

    for (uint32_t j = 0; j < i; ++j)
        U_powerj << aqs::Gate(U_powerj, 0);
 
    qc << aqs::ControlGate{U_powerj,  8 - 1 - i, 8,
                            "13^(2^" + std::to_string(i) + ") mod 15"};
}

The final part of the circuit is obtaining the value of $r$ which can be done by using the Quantum Phase Estimation Algorithm. In simple terms, what it does is find the $\theta$ from the operation:

$P|\psi\rangle = e^{2\pi i \theta} |\psi\rangle$

This phase rotation $\theta$ can be found by applying the inverse Quantum Fourier Transform $QFT^\dag$ :

qc << aqs::Gate(aqs::inverse_fourier_transform(8), 0, "QFT†");

The way we obtain $r$ from $\theta$ comes from

Simulating

To test if the algorithm works, we need to start a simulation of the circuit and measure the output of the qubits multiple times to get a distribution of the results. Before simulating the algorithm, we need to set the last qubit to the $|1\rangle$ and make sure it propagates to the whole quantum computer state. The reason for setting this qubit reason is that we want to make sure that the $U$ operator acts on a state that has an interesting repeating remainder pattern, that is, a state that does not produce 0 as the remainder of 13. Setting the last qubit to $|1\rangle$ assigns our starting number in the sequence to be 1. All of this can be done as follows:

aqs::QSimulator qs(qc.qubit_count());

qs.qubit(qs.qubit_count() - 1) = aqs::QState::one();

qs.generate_statevector();

qs.simulate(qc)

auto profile = qs.profile_measure(1000);

What this code does is create a Quantum Simulator, then set the initial conditions of this simulator for this simulator to have all $|0\rangle$ for the qubits except for the last one, and then update the internal state of the quantum computer (the state vector). The simulation of the circuit $\texttt{qc}$ is then done with the simulator $\texttt{qs}$ . The repetition of multiple measurements is done by executing the simulation with $\texttt{profile\_measure}$ . For this example, we effectively repeated the experiment 1000 times, from which we can get a rough estimate of the probable value for $r$ .

Of course, we are interested in the most probable result. To get a rank of the most probable values of $r$ , we need to process the data of the results. The value of $r$ is stored in the measurement of the qubit registers coming from the inverse QFT through the Quantum Estimation Algorithm. For this circuit, the result of the inverse QFT is encoded in the first 8 qubits and their states have the most significant qubit last, so we need to reverse the binary digit results. Finally, we want to sort those results in terms of the number of times they have been measured. This is done by the following code:

std::vector<std::pair<int, int>> results;
results.resize(256, {0, 0});

for (uint32_t i = 0; i < profile.size(); ++i)
{
    uint32_t value = 0;
    for (uint32_t j = 0; j < 8; ++j)
        value = value | (((i / 16) >> j) & 1) << (7 - j);

    results[value].first = value;
    results[value].second += profile[i];
}

std::sort(results.begin(), results.end(),
[](const std::pair<int, int>& a, const std::pair<int, int>& b) {
    return b.second < a.second;
});

The first part takes the 8 bits of the binary string, reverses it, and counts how many times the values were measured. The second part sorts those results in descending order. From there, we can take, for example, the top 4 guesses to determine $r$ :

for (int i = 0; i < 4; ++i)
{
    auto value = results[i].first;
    double angle = value / 256.0;
    auto fraction = approximate_fraction(angle, 15);
    auto r = fraction.second;
        
   std::cout << "------ Guess #" << i + 1 << "------\n"
             << "state = 0b" << binary_string(value, 8) << " = " << value << '\n'
             << "angle = " << value << " / 256" << " = " << angle << " ~ "
             << fraction.first << " / " << fraction.second << '\n'
             << "r = " << r << '\n';

    if (r % 2)
    {
        auto guess_factor1 = gcd(pow(a, r/2) + 1, N);
        auto guess_factor2 = gcd(pow(a, r/2) - 1, N);

        std::cout << "Factors: " << guess_factor1
                        << " , " << guess_factor2 << '\n';

        if (guess_factor1 != 1 && guess_factor2 != N)
            std::cout << "Found prime factor: " << guess_factor1 << '\n';
        if (guess_factor2 != 1 && guess_factor2 != N)
            std::cout << "Found prime factor: " << guess_factor2 << '\n';
    }    
    std::cout << std::endl;
}

To determine $r$ , what we did was find to approximate the phase angle $\theta$ , then with this angle, we can find a fraction of the phase rotation to a complete phase rotation of $2\pi$ from which the denominator will be our value $r$ .

After obtaining $r$ , we compute the guess factors from $a, r, N$ . The reason we add the Greatest Common Divisor $\texttt{gcd}$ function is to remove any other factors that are not factors of $N$ but may be attached to the guess factor. Thus, the gcd gives us our prime factors, and with one, you can obtain the other.

You can look at what the quantum circuit diagram looks by calling the function:

aqs::print_circuit_text_image(qc, qs);

Which gives the following Unicode image:

     ┌───┐                                                                                                                                                                            ┌──────┐
|0⟩──┤ H ├───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────█────────────┤ QFT† ├──
     └───┘                                                                                                                                                               │            │      │
     ┌───┐                                                                                                                                                               │            │      │
|0⟩──┤ H ├──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────█────────────────────┼────────────┤      ├──
     └───┘                                                                                                                                          │                    │            │      │
     ┌───┐                                                                                                                                          │                    │            │      │
|0⟩──┤ H ├─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────█────────────────────┼────────────────────┼────────────┤      ├──
     └───┘                                                                                                                     │                    │                    │            │      │
     ┌───┐                                                                                                                     │                    │                    │            │      │
|0⟩──┤ H ├────────────────────────────────────────────────────────────────────────────────────────────────█────────────────────┼────────────────────┼────────────────────┼────────────┤      ├──
     └───┘                                                                                                │                    │                    │                    │            │      │
     ┌───┐                                                                                                │                    │                    │                    │            │      │
|0⟩──┤ H ├───────────────────────────────────────────────────────────────────────────█────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────┤      ├──
     └───┘                                                                           │                    │                    │                    │                    │            │      │
     ┌───┐                                                                           │                    │                    │                    │                    │            │      │
|0⟩──┤ H ├──────────────────────────────────────────────────────█────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────┤      ├──
     └───┘                                                      │                    │                    │                    │                    │                    │            │      │
     ┌───┐                                                      │                    │                    │                    │                    │                    │            │      │
|0⟩──┤ H ├─────────────────────────────────█────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────┤      ├──
     └───┘                                 │                    │                    │                    │                    │                    │                    │            │      │
     ┌───┐                                 │                    │                    │                    │                    │                    │                    │            │      │
|0⟩──┤ H ├────────────█────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────────────┼────────────┤      ├──
     └───┘            │                    │                    │                    │                    │                    │                    │                    │            └──────┘
              ┌───────┴───────┐    ┌───────┴───────┐    ┌───────┴───────┐    ┌───────┴───────┐    ┌───────┴───────┐    ┌───────┴───────┐    ┌───────┴───────┐    ┌───────┴───────┐
|0⟩───────────┤ 13^(2^0)mod15 ├────┤ 13^(2^1)mod15 ├────┤ 13^(2^2)mod15 ├────┤ 13^(2^3)mod15 ├────┤ 13^(2^4)mod15 ├────┤ 13^(2^5)mod15 ├────┤ 13^(2^6)mod15 ├────┤ 13^(2^7)mod15 ├──────────────
              │               │    │               │    │               │    │               │    │               │    │               │    │               │    │               │
              │               │    │               │    │               │    │               │    │               │    │               │    │               │    │               │
|0⟩───────────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├──────────────
              │               │    │               │    │               │    │               │    │               │    │               │    │               │    │               │
              │               │    │               │    │               │    │               │    │               │    │               │    │               │    │               │
|0⟩───────────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├──────────────
              │               │    │               │    │               │    │               │    │               │    │               │    │               │    │               │
              │               │    │               │    │               │    │               │    │               │    │               │    │               │    │               │
|1⟩───────────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├────┤               ├──────────────
              └───────────────┘    └───────────────┘    └───────────────┘    └───────────────┘    └───────────────┘    └───────────────┘    └───────────────┘    └───────────────┘

Testing

After running this, you should get an output like this:

------ Guess #1------
state = 0b11000000 = 192
angle = 192 / 256 = 0.75 ~ 3 / 4
r = 4
Factors: 5 , 3
Found prime factor: 5
Found prime factor: 3

------ Guess #2------
state = 0b00000000 = 0
angle = 0 / 256 = 0 ~ 0 / 1
r = 1

------ Guess #3------
state = 0b01000000 = 64
angle = 64 / 256 = 0.25 ~ 1 / 4
r = 4
Factors: 5 , 3
Found prime factor: 5
Found prime factor: 3

------ Guess #4------
state = 0b10000000 = 128
angle = 128 / 256 = 0.5 ~ 1 / 2
r = 2
Factors: 1 , 3
Found prime factor: 3

While there may be variations in the order of the guesses, they should be pretty much the same ones. From these results, we can see that the first guess gave us a value of $r=4$ which provides us with two prime factor guesses: 5 and 3 which are exactly the factors of $p, q$ of 15. Thus, we have found the prime factors of 15 and cracked the encryption.

Through this example, we have seen how you can simulate a quantum algorithm with a normal computer using the ArrayFire Quantum Simulator by using the process that would be followed for implementing this algorithm in an actual quantum computer. Therefore, we can see the utility of simulators like this for developing new algorithms and testing ideas for other quantum circuits.

Conclusion

Throughout this Quantum Computer series, we have seen how quantum computers represent a revolutionary approach to computation that promises to unlock unprecedented computing power. By going into some of the mathematics, we were able to describe how quantum computers operate on the principles of quantum mechanics, allowing them to perform calculations that would be infeasible with classical computers. Algorithms, such as Shor’s algorithm and Grover’s algorithm, are specifically designed to exploit this advantage. Although quantum computers are still in their infancy, significant progress has been made in recent years, and it is expected that quantum computing will continue to advance rapidly aiding in multiple areas of research including Finance, Biochemistry, and AI.

Finally, we saw that with the help of quantum computer simulators, researchers and developers can test and refine quantum algorithms before implementing them on physical quantum computers, accelerating the progress of quantum computing research and development. We saw how the ArrayFire Quantum Simulator (AQS) could be used to simulate and test these algorithms on your machines and better understand quantum computing. Quantum computing is an exciting and rapidly evolving field that has the potential to transform computing as we know it.

References

Quantum Computer Development by Qubit Count Image. Yoke Development. “Quantum Technologies 2021”.. https://s3.i-micronews.com/uploads/2021/06/YINTR21211-Quantum-Technologies-2021-Sample.pdf

Qiskit Textbook Shor’s Algorithm Implementation: https://qiskit.org/textbook/ch-algorithms/shor.html

Shor, P.W. (1994). Algorithms for quantum computation: discrete logarithms and factoring. Proceedings 35th Annual Symposium on Foundations of Computer Science, 124-134.

NIST. NIST Announces First Four Quantum-Resistant Cryptographic Algorithms. (2022). https://www.nist.gov/news-events/news/2022/07/nist-announces-first-four-quantum-resistant-cryptographic-algorithms

Google Sycamore Data Sheet, https://quantumai.google/hardware/datasheet/weber.pdf

IBM, Next Wave Quantum Centric SuperComputing, https://research.ibm.com/blog/next-wave-quantum-centric-supercomputing

DWave, Products & Solutions, https://www.dwavesys.com/solutions-and-products/systems/

The State Of Quantum Computing: Future, Present, Past (forbes.com)

Nielsen, M., & Chuang, I. (2010). Quantum Computation and Quantum Information: 10th Anniversary Edition. Cambridge: Cambridge University Press. doi:10.1017/CBO9780511976667

ArrayFire Quantum Simulator. https://github.com/arrayfire/afQuantumSim

$x$	$a^x$	$a^x\text{ mod} N$
0	1	1
1	3	3
2	9	2
3	27	6
4	81	4
5	243	5
6	729	1
7	2187	3
…	…	…

$\|x\rangle$	$U\|x\rangle = \|ax \text{mod} N\rangle$ for $a=13, N=15$
$\|0\rangle$	$\|15\rangle$
$\|1\rangle$	$\|13\rangle$
$\|2\rangle$	$\|11\rangle$
$\dots$	$\dots$
$\|15\rangle$	$\|0\rangle$