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.

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!

Alongside this guide, you may also want the documentation.

QCGPU is licensed under the MIT License.

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.

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.

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

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\]

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

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);
#}

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);
#}

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

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>
#}

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]
#}

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

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();
#}

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.

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.

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}
}

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 by an oracle circuit \(U_{f}\) such that \(U_{f} \lvert x \rangle = (-1)^{f(x)}\lvert x \rangle\).

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"
        }
    );
}

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
}

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

(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");
            }
        }
    }
}

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!()
}
#}

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),
    }
}