# QCGPU

QCGPU is a high performance, hardware accelerated quantum computer simulator written with Rust and OpenCL.

QCGPU is free and open source. The source code is available on GitHub. Please report any places where the documentation is unclear, or where examples could be added or improved by creating an issue.

## Help Wanted

Please request any functionality you would like, need or appreciate by creating an issue. If you would like to add functionality or fix issues, please create a pull request and I will get back to you (usually within the day).

Any other feedback, suggestions or even tiny bits of criticism are welcome. When in doubt, file an issue!

## API Docs

Alongside this guide, you may also want the documentation.

## License

QCGPU is licensed under the MIT License.

# Getting Started

## Installing The Requirements

To use QCGPU, you will need Rust and OpenCL installed.

The setup process for OpenCL will be different for every device. All apple devices (MacOS / OSX) will have OpenCL already installed. For some hints on how to install on linux, look at the AWS EC2 Install Instructions. There is also a good chance that it is installed already. Check that `clinfo` for some other diagnostic command will run.

Rust is very easy to install. Check out rustup.

## Adding The Dependency

To use the library with rust, you must add the following to your `cargo.toml` file:

    [dependencies]
    qcgpu = "0.1"

You should now be able to use the library by adding

    extern crate qcgpu;

to `lib.rs` or `main.rs`, depending on if you are writing an executable or a library.

# User Guide

This chapter covers the usage of the QCGPU library. It will also contain some information about the mathematics that each of the functions represents.

All aspects of the library are covered by this chapter, whereas complete examples (in the form of algorithm implementations) are available in the Algorithms Chapter

# Quantum Registers

All of the simulation is done through quantum registers. QCGPU provides a struct as a register, but that contains fields to do with the OpenCL buffers and related items, so the creation of registers should be done through the provided methods.

The library is optimized for complex superpositions, so the registers are all dense. This means that the number of qubits you can initialize is directly related to the capacity / available memory of the device.

The register struct is called `State` and is available through `qcgpu::State`.

To create a register, the easiest way to do it is with the `State::new` method. It takes two parameters, the number of qubits and the device to use. The device is given as a usize, and corresponds to the OpenCL device which the register will be on.

The following example creates a register with 5 qubits on the 1st device.

    # extern crate qcgpu;

    use qcgpu::State;

    # fn main() {
    let mut register = State::new(5, 0);
    # }

Notice that the register is mutable. This allows the register to change. Also, the device is 0 indexed.

The implementation is equivilent to the description of a state vector \( \lvert \psi \rangle \) with

\[ \lvert \psi \rangle = \sum_{j = 0}^{2^n - 1} \alpha_j \lvert j \rangle \]

where \(n\) is the number of qubits, \(\alpha_j\) is the amplitude and the state is \(j\) runs over all \(2^n\) basis states.

There is one other way to initialize a state. Given a bitstring, a register can be initialized in that state using the `State::from_bit_string` method. For example, to initialize a register with the value `0100` the following is used:

    # extern crate qcgpu;

    use qcgpu::State;

    # fn main() {
    let mut register = State::from_bit_string("|0100>", 0);
    # }

The second argument is the same as before. The register that is outputed from this method is equivilent to the state

\[ \lvert 0100 \rangle\]

# Quantum Gates

Gates are used to manipulate quantum registers and to implement quantum algorithms.

## Built In Gates

There are a number of gates built in to QCGPU. They can all be applied the same way:

    # #![allow(unused_variables)]
    #fn main() {
    use qcgpu::State;

    let mut state = State::new(5, 0);

    state.h(0); // Applies the hadamard (`h`) gate to the 0th qubit
    #}

`h` can be replaced with any of the following:

- The hadmard gate: h - `state.h(0);`
- The S gate: s - `state.s(0);`
- The T gate: t - `state.t(0);`
- The Pauli-X / NOT gate: x - `state.x(0);`
- The Pauli-Y gate: y - `state.y(0);`
- The Pauli-Z gate: z - `state.z(0);`
- The CNOT gate: cx - `state.cx(0, 1); // CNOT with control = 0, target = 1`
- The SWAP gate: swap - `state.swap(0,1); // Swaps the 0th and 1st qubit`
- The Toffoli gate: toffoli - `state.toffoli(0, 1, 2); // Toffoli with control1 = 0, control1 = 1, target = 2`

These are all shorthand methods for the application of arbitrary gates. For example, the application of a hadamard gate above is shorthand for

    # #![allow(unused_variables)]
    #fn main() {
    use qcgpu::gates::{h};
    use qcgpu::State;

    let mut state = State::new(5, 0);
    state.apply_gate(h(), 0);
    #}

You can also use any of the gates as controlled gates. For example, the application of the CNOT gate above is shorthand for

    # #![allow(unused_variables)]
    #fn main() {
    use qcgpu::gates::{x};
    use qcgpu::State;

    let mut state = State::new(5, 0);
    state.apply_controlled_gate(x(), 0, 1);
    #}

## User Defined Gates

Gates in QCGPU are represented by the `Gate` struct, available through `qcgpu::Gate`.

It is defined as follows:

    # #![allow(unused_variables)]
    #fn main() {
    extern crate num_complex;
    use num_complex::Complex32;

    struct Gate {
        a: Complex32,
        b: Complex32,
        c: Complex32,
        d: Complex32,
    }
    #}

To create your own gate, you will need to add the `num_complex` crate to your dependencies.

A gate is created as follows:

    # #![allow(unused_variables)]
    #fn main() {
    let x = Gate {
        Gate {
            a: Complex32::new(0.0, 0.0),
            b: Complex32::new(1.0, 0.0),
            c: Complex32::new(1.0, 0.0),
            d: Complex32::new(0.0, 0.0),
        }
    }
    #}

This corresponds to the matrix

\[x = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}\]

This can be applied using the same long hand method as above:

    # #![allow(unused_variables)]
    #fn main() {
    let mut state = State::new(1, 0);
    state.apply_gate(x, 0);
    #}

# Quantum Operations

There are a number of operations you can preform on quantum registers with QCGPU.

## Measurement

You can measure the register in two ways. You can either do a single measurement and return an integer with the measured value or you can measure multiple times and return a `HashMap<String, i32>`, with the key being the bitstring, and the value being the number of times it was measured.

The measurements don't collapse the state.

They are used as follows:

    # #![allow(unused_variables)]
    #fn main() {
    use qcgpu::State;

    let mut state = State::new(5, 0);
    state.measure(); // Returns an integer
    state.measure_many(1000); // Measures 1000 times, returns a HashMap<String, i32>
    #}

There is also a convenience method to measure the first \(n\) qubits in the register. Again, the state is not collapsed

    # #![allow(unused_variables)]
    #fn main() {
    use qcgpu::State;

    let mut state = State::new(5, 0);
    state.measure_first(3, 1000); // Measures the first 3 qubits 1000 times, returns a HashMap<String, i32>
    #}

## Probability

QCGPU provides another method for getting the probability of each outcome.

The probability is calculated for a state \(\lvert \psi \rangle = \sum_{j = 0}^{2^n - 1} \alpha_j \lvert j \rangle\),

\[P(j) = |\alpha_j|^2\]

The method `get_probabilities` returns a `Vec<f32>` with each of the values corresponding to \(|\alpha_j|^2\) for each index \(j\).

    # #![allow(unused_variables)]
    #fn main() {
    use qcgpu::State;

    let mut state = State::new(1, 0);
    state.h(0);
    state.get_probabilities(); // [0.5, 0.5]
    #}

# Examples

Here is a complete example of using QCGPU to simulate the bell state.

    extern crate qcgpu;

    use qcgpu::State;

    fn main() {
        // Create a new quantum register with 2 qubits on OpenCL device #0
        let mut register = State::new(2, 0);

        // Apply a hadamard gate to the first qubit
        register.h(0);

        // Apply a CNOT, with the control = 0, and target = 1
        register.cx(0, 1);

        // Measure the register 1000 times and print the output
        println!("Measured: {:?}", register.measure_many(1000));
    }

The output should be something like

    Measured: {"00": 516, "11": 484}

For more, non trivial examples, see the Algorithms

# Decoherence

QCGPU provides a way to easily simulate the effects of decoherence on the quantum computer. The effects are simulated by a random gate, corresponding to a random rotation around the \(z\) axis of the bloch sphere.

The angle of the rotation is a normal distrobuted value, with the varience as the strength factor `d`.

To avoid performance costs when not being used, decoherence can be enabled via a feature.

Change the dependency in `cargo.toml` to

    [dependencies]
    qcgpu = { version = "0.1", features = ["decoherence"] }

Now you can add decoherence to your simulator.

To set the amount of decoherence, the `set_decoherence` method is used. Its arguments are the strength of the decoherence. This will affect all following gate applications.

    # #![allow(unused_variables)]
    #fn main() {
    use qcgpu::State;

    let mut register = State::new(5, 0);
    register.set_decoherence(0.4);
    #}

You can also manually decohere the register based on the previously set strength value. This method is automatically called by all gate applications.

    # #![allow(unused_variables)]
    #fn main() {
    register.decohere();
    #}

# Algorithms

The following chapter details some implementations of non trivial quantum algorithms. They are based on the source code from libquantum along with the lecture notes from the University of Waterloo.

Note that some implementations may not yet be complete.

# Bernstein-Vazirani Algorithm

This algorithm finds a hidden integer \(a \in { 0, 1}^n\)from an oracle \(f_a\)which returns a bit \(a \cdot x \equiv \sum_i a_i x_i \mod 2\) for an input \(x \in {0,1}^n\).

A classical oracle returns \(f_a(x) = a \dot x \mod 2\), while the quantum oracle must be queried with superpositions of input \(x\)'s.

To solve this problem classically, the hidden integer can be found by checking the oracle with the inputs \(x = 1,2,\dots,2^i,2^{n-1}\), where each query reveals the \(i\)th bit of \(a\) (\(a_i\)). This is the optimal classical solution, and is O(n). Using a quantum oracle and the Bernstein-Vazirani algorithm, \(a\) can be found with just one query to the oracle.

## The Algorithm

1. Initialize \(n\) qubits in the state \(\lvert 0, \dots, 0\rangle\).
2. Apply the Hadamard gate \(H\) to each qubit.
3. Apply the inner product oracle.
4. Apply the Hadamard gate \(H\) to each qubit.
5. Measure the register

From this procedure, we find that the registers measured value is equal to that of the original hidden integer.

    extern crate qcgpu;

    use qcgpu::State;
    use qcgpu::gates::h;

    fn main() {
        let num_qubits = 16; // Number of qubits to use
        let a = 101; // Hidden integer, bitstring is 1100101

        // You should also make sure that a is representable with n qubits,
        // by settings a as a mod 2^n.

        // Bernstein-Vazirani algorithm
        let mut state = State::new(num_qubits, 1); // New quantum register, using the GPU.

        // Apply a hadamard gate to each qubit
        state.apply_all(h());

        // Apply the inner products oracle
        for i in 0..num_qubits {
            if a & (1 << i) != 0 {
                state.z(i as i32);
            }
            // Otherwise should apply identity gate, but computationally this doens't change the state.
        }

        // Apply hadamard gates before measuring
        state.apply_all(h());

        println!("Measurement Results: {:?}", state.measure_many(1000));
        // Measurement Results: {"0000000001100101": 1000}
    }

# Deutsch-Jozsa Algorithm

This algorithm was the first to show that quantum computers could have a speedup over classical computers.

Consider a function \(f(x)\) which takes an input of an \(n\)-bit string \(x\) and returns 0 or 1.

Suppose that \(f(x)\) is either a constant function that has the same value \(c \in {0, 1}, \forall x\), or a balanced function, where the value is 0 for half of the inputs, and 1 for the other half.

The Deutsch-Jozsa problem is to find whether \(f\) is constant or balanced, in as few function evaluations as possible.

Using classical computing, in the worst case, this requires \(2^{n-1}+1\) function evaluations. Using quantum computing, this can be done with just one function evaluation.

The function \(f\), to be used in a quantum computer, must be specified specified by an oracle circuit \(U_{f}\) such that \(U_{f} \lvert x \rangle = (-1)^{f(x)}\lvert x \rangle\).

## The Algorithm

1. Initialize \(n\) qubits in the state \(\lvert 0, \dots, 0\rangle\).
2. Apply the Hadamard gate \(H\) to each qubit.
3. Apply the oracle circuit \(U_f\).
4. Apply the Hadamard gate \(H\) to each qubit.
5. Measure each qubit. Let \(y = (y_{1}, \dots, y_{n})\) be the list of measurement outcomes.

From this procedure, we find that \(f\) is constant if \(y\) is the all zero string.

    extern crate qcgpu;

    use qcgpu::State;
    use qcgpu::gates::h;

    fn main() {
        // 3 qubits, f(x) = x_0 NOT x_1 x_2
        // Balanced
        let mut balanced_state = State::new(3, 1);

        balanced_state.apply_all(h());

        // Oracle U_f
        balanced_state.h(2);
        balanced_state.z(0);
        balanced_state.cx(1, 2);
        balanced_state.h(2);

        balanced_state.apply_all(h());

        println!(
            "{}",
            if balanced_state.measure() == 0 {
                "constant"
            } else {
                "balanced"
            }
        );

        // 3 qubits, f(x) = 0
        // Constant
        let mut constant_state = State::new(3, 1);

        constant_state.apply_all(h());

        // Oracle is equivilant to the identity gate, thus has no effect on the
        // state.

        constant_state.apply_all(h());

        println!(
            "{}",
            if constant_state.measure() == 0 {
                "constant"
            } else {
                "balanced"
            }
        );
    }

# Grovers Algorithm

Given an unstructured set \(N = {a_1, a_2,\dots,a_n}\), find a given element \(a_i \in N\).

This implementation looks for a given number \(target\) in the set \({0,1,\dots, \text{regwidth}}\).

See https://cs.uwaterloo.ca/~watrous/LectureNotes.html

    extern crate qcgpu;
    extern crate rand;

    use qcgpu::State;
    use std::f32::consts::PI;

    fn main() {
        let mut state = State::new(8, 1);

        let target = 5;
        let reg_width = 3;

        let num_inversions = ((PI / 4.0) * ((1 << reg_width) as f32).sqrt()) as i32;

        state.x(reg_width);

        for i in 0..(reg_width + 1) {
            state.h(i);
        }

        for _ in 0..(num_inversions) {
            iteration(&mut state, target, reg_width);
        }

        state.h(reg_width);

        println!("Measured: {:?}", state.measure_first(reg_width, 1000));
    }

    fn oracle(state: &mut State, target: i32, reg_width: i32) {
        for i in 0..reg_width {
            if get_bit(target, i) == 0 {
                state.x(i);
            }
        }

        state.toffoli(0, 1, reg_width + 1);

        let mut i = 1;
        while i < reg_width {
            state.toffoli(i, reg_width + i, reg_width + i + 1);
            i += 1;
        }

        state.cx(reg_width + i, reg_width);

        i = reg_width - 1;
        while i > 0 {
            state.toffoli(i, reg_width + i, reg_width + i + 1);
            i -= 1;
        }

        state.toffoli(0, 1, reg_width + 1);

        for i in 0..reg_width {
            if get_bit(target, i) == 0 {
                state.x(i);
            }
        }
    }

    fn inversion(state: &mut State, reg_width: i32) {
        for i in 0..reg_width {
            state.x(i);
        }

        state.h(reg_width - 1);

        if reg_width == 3 {
            state.toffoli(0, 1, 2);
        } else {
            state.toffoli(0, 1, reg_width + 1);

            let mut i = 1;
            while i < reg_width {
                state.toffoli(i, reg_width + i, reg_width + i + 1);
                i += 1;
            }

            state.cx(reg_width + i, reg_width - 1);

            i = reg_width - 2;
            while i > 0 {
                state.toffoli(i, reg_width + i, reg_width + i + 1);
                i -= 1;
            }

            state.toffoli(0, 1, reg_width + 1);
        }

        state.h(reg_width - 1);

        for i in 0..reg_width {
            state.x(i);
        }
    }

    fn iteration(state: &mut State, target: i32, reg_width: i32) {
        oracle(state, target, reg_width);

        for i in 0..reg_width {
            state.h(i);
        }

        inversion(state, reg_width);

        for i in 0..reg_width {
            state.h(i);
        }
    }

    /// Get the value of a bit
    fn get_bit(number: i32, n: i32) -> i32 {
        if number & (1 << n) != 0 {
            return 1;
        }
        0
    }

# Shors Algorithm

This algorithm finds the prime factors (\)u\) and \(v\)) of an odd, composite integer \(n\), that is not a prime power.

## The Algorithm

(pseudo code)

    Repeat
        Randomly choose \\(a \in \{ 2, \dots, n - 1 \}\\)
        Compute \\(d = gcd(a, n)\\)
        If \\(d \geq 2\\) then
            Return \\(u = d\\) and \\(v = n/d\\)
        Else // We know \\(a \in \mathbb{Z}^*_N\\)
            Let \\(r\\) be the order of \\(a\\) in \\(\mathbb{Z}^*_N\\) // Order finding algorithm
            If \\(r\\) is even then
                Compute \\(x = a^{r/2} - 1 (\mod n)\\)
                Compute \\(d = \gcd(x, n)\\)
                If \\(d \geq 2\\) then
                    Return \\(u = d\\) and \\(v = n/d\\)
    Until a value is returned.

See https://cs.uwaterloo.ca/~watrous/LectureNotes.html

    extern crate qcgpu;
    extern crate rand;

    use rand::{thread_rng, Rng};
    use qcgpu::{gcd, State};
    use qcgpu::gates::h;

    fn main() {
        let n = 15; // Number to factor
        println!("Factoring {}.", n);

        // Here we should check if the number is even, or if it is a power factor

        let mut rng = thread_rng();
        loop {
            let mut a = rng.gen_range(2, n); // Randomly choose \\(a \in \{ 2, \dots, n - 1 \}\\)

            let mut d = gcd(a, n);
            if d >= 2 {
                // Found the factors; No Quantum needed
                println!(
                    "Factors are {} and {} (No quantum algorithm used)",
                    d,
                    n / d
                );
                break;
            } else {
                // We know \\(a \in \mathbb{Z}^*_N\\)
                let r = find_order(a, n);
                if r % 2 == 0 {
                    let x = (a.pow((r as f32 / 2 as f32) as u32) - 1) % n;
                    d = gcd(x, n);
                    if d >= 2 {
                        println!("Factors are {} and {}", d, n / d);
                        break;
                    }
                } else {
                    println!("Period is odd");
                }
            }
        }
    }

## Order Finding

Given a positive integer \(n \geq 2\) and an element \(a \in \mathbb{Z}_n^* \), Find the order of \(a\) in \(\mathbb{Z}_n^*\).

Order finding is the only quantum part in shors algorithm.

\(\mathbb{Z}_n\) is the set of integers from \({0,\dots, n - 1}\) or \(\mathbb{z} \mod n\). The set \(\mathbb{Z}_n^* = {a \in \mathbb{Z}_n : \gcd(a,n) = 1}\)

The set \( \mathbb{Z}_n^* \) forms a group wth the multiplication modulo \(n\) operation. Thus, for \( a \in \mathbb{Z}_n^* \) , \(\exists b \in \mathbb{Z}_n^*\) that uniquely satisfies

\[ ab \equiv 1 \mod n \]

The order of a given element \(a \in \mathbb{Z}_n^*\) is the smallest positive integer \(r\) such that

\[ a^r \equiv 1 \mod n \]

    # #![allow(unused_variables)]
    #fn main() {
    fn find_order(a: i32, n: i32) -> i32 {
        unimplemented!()
    }
    #}

# Super Dense Coding

if Alice and Bob share a pair of entangled qubits, then Alice can encode two classical bits into her one entangled qubit, send it to Bob, and Bob can decode it with the help of his entangled qubit.

    extern crate qcgpu;

    use qcgpu::State;

    fn superdense(input: &str) -> i32 {
        let mut state = State::new(2, 0);
        let input_str = String::from(input);

        // Prepare the bell state
        state.h(0);
        state.cx(0, 1);

        // Alice prepares her qubit
        let alice = 1;
        if input_str.get(0..1) == Some("1") {
            state.z(alice);
        }
        if input_str.get(1..2) == Some("1") {
            state.x(alice);
        }

        println!("\nState after Alice prepares her qubit: \n{}", state);

        // Alice sends her qubit to Bob
        let bob = 0;
        state.cx(alice, bob);
        state.h(alice);

        println!(
            "\nState after Bob receives Alice's qubit and 'decodes' it: \n{}",
            state
        );

        state.measure()
    }

    fn main() {
        use std::io;

        println!("Two bit string to send:");
        let mut input = String::new();
        match io::stdin().read_line(&mut input) {
            Ok(_n) => {
                let result = superdense(input.as_str());
                println!("\nDecoded string is: {}", result);
            }
            Err(error) => println!("error: {}", error),
        }
    }