4. Simple examples

All of the examples of this section are copied verbatim from the directory examples and are fully compilable. For convenience, the location of the source file is displayed in the first line of each example as a C++ comment. The examples are simple and demonstrate the main features of Quantum++. They cover only a small part of library functions, but enough to get the interested user started. For an extensive reference of all library functions, including various overloads, the user should consult the official API documentation. Additional examples (not discussed in this document) are located in examples.

Gates and states

Let us introduce the main objects used by Quantum++: gates, states and basic operations. Consider the example examples/gates_states.cpp with content listed below.

// Gates and states
// Source: ./examples/gates_states.cpp

#include <iostream>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    ket psi = 0_ket; // |0> state
    cmat U = gt.X;
    ket result = U * psi;

    std::cout << ">> The result of applying the bit-flip gate X on |0> is:\n";
    std::cout << disp(result) << '\n';

    psi = 10_ket; // |10> state
    U = gt.CNOT;  // Controlled-NOT
    result = U * psi;

    std::cout << ">> The result of applying the gate CNOT on |10> is:\n";
    std::cout << disp(result) << '\n';

    U = randU();
    std::cout << ">> Generating a random one-qubit gate U:\n";
    std::cout << disp(U) << '\n';

    result = applyCTRL(psi, U, {0}, {1}); // Controlled-U
    std::cout << ">> The result of applying the CTRL-U gate on |10> is:\n";
    std::cout << disp(result) << '\n';
}

A possible output is

>> The result of applying the bit-flip gate X on |0> is:
0
1
>> The result of applying the gate CNOT on |10> is:
0
0
0
1
>> Generating a random one-qubit gate U:
-0.50206 + 0.0354176i  -0.422949 - 0.753522i
-0.206807 - 0.838995i  -0.390495 + 0.317541i
>> The result of applying the CTRL-U gate on |10> is:
                   0
                   0
-0.50206 + 0.0354176i
-0.206807 - 0.838995i

In the line

using namespace qpp;

we bring the namespace qpp into the global namespace. In the line

ket psi = st.z0; // |0> state

we use the qpp::States singleton st to declare psi as the zero eigenvector |0> of the Z Pauli operator. In the line

cmat U = gt.X;

we use the qpp::Gates singleton gt and assign to U the bit flip gate gt.X. In the next line

ket result = U * psi;

we compute the result of the operation X|0>, and display the result |1> in the lines

std::cout << ">> The result of applying the bit-flip gate X on |0> is:\n";
std::cout << disp(result) << '\n';

Above we used the display manipulator disp(), which is especially useful when displaying complex matrices, as it displays the entries of the latter in the form $a+bi$, in contrast to the form $(a,b)$ used by the C++ standard library. The manipulator also accepts additional parameters that allows, e.g., setting to zero numbers smaller than some given value (useful to chop small values), and it is in addition overloaded for standard containers, iterators and C-style arrays.

In the line

psi = 10_ket;  // |10> state

we reassign to psi the state |10> constructed via the user-defined literal qpp::ket qpp::operator"" _ket() which transforms a strings of zero and ones to the corresponding multi-qubit ket (in the computational basis). We could have also used the Eigen 3 insertion operator

ket psi(4); // must specify the dimension before insertion of elements via <<
psi << 0, 0, 1, 0;

or the Quantum++ qpp::mket() function. In the next line

U = gt.CNOT; // Controlled-NOT

we re-assign to the gate U the Controlled-NOT with control as the first subsystem, and target as the last, using the global singleton gt. In the next line

result = U * psi;

we declare the ket result as the result of applying the Controlled-NOT gate to the state |10>, i.e., |11>. In the next two lines

std::cout << ">> Generating a random one-qubit gate U:\n";
std::cout << disp(U) << '\n';

we then display the result of the computation. Next, in the line

U = randU();

we generate a random unitary gate via the function qpp::randU(), then in the line

result = applyCTRL(psi, U, {0}, {1}); // Controlled-U

we apply the Controlled-U, with control as the first qubit and target as the second qubit, to the state psi. Finally, we display the result in the last two lines.

Measurements

Let us now complicate things a bit and introduce measurements. Consider the example examples/measurements1.cpp displayed below.

// Measurements
// Source: ./examples/measurements1.cpp

#include <iostream>
#include <tuple>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    ket psi = 00_ket;
    cmat U = gt.CNOT * kron(gt.H, gt.Id2);
    ket result = U * psi; // we have the Bell state (|00> + |11>) / sqrt(2)

    std::cout << ">> We just produced the Bell state:\n";
    std::cout << disp(result) << '\n';

    // apply a bit flip on the second qubit
    result = apply(result, gt.X, {1}); // we produced (|01> + |10>) / sqrt(2)
    std::cout << ">> We produced the Bell state:\n";
    std::cout << disp(result) << '\n';

    // measure the first qubit in the X basis
    auto [m, probs, states] = measure(result, gt.H, {0});
    std::cout << ">> Measurement result: " << m << '\n';
    std::cout << ">> Probabilities: ";
    std::cout << disp(probs, IOManipContainerOpts{}.set_sep(", ")) << '\n';
    std::cout << ">> Resulting states:\n";
    for (auto&& elem : states) {
        std::cout << disp(elem) << "\n\n";
    }
}

A possible output is

>> We just produced the Bell state:
0.707107
       0
       0
0.707107
>> We produced the Bell state:
       0
0.707107
0.707107
       0
>> Measurement result: 0
>> Probabilities: [0.5, 0.5]
>> Resulting states:
0.707107
0.707107

-0.707107
 0.707107

In the line

cmat U = gt.CNOT * kron(gt.H, gt.Id2);

we use the function qpp::kron() to create the tensor product (Kronecker product) of the Hadamard gate on the first qubit and identity on the second qubit, then we left-multiply it by the Controlled-NOT gate. In the line

ket result = U * psi; // we have the Bell state (|00> + |11>) / sqrt(2)

we compute the result of the operation $CNOT_{ab}(H\otimes I)|00>$, which is the Bell state $(|00> + |11>)/\sqrt{2}$. We display the result in the next two lines.

Next, in the line

result = apply(result, gt.X, {1}); // we produced (|01> + |10>) / sqrt(2)

we use the function qpp::apply() to apply the gate X on the second qubit of the previously produced Bell state. The function qpp::apply() takes as its third parameter a list of subsystems, and in our case {1} denotes the second subsystem, not the first (again we follow the C++ indexing convention).

The function qpp::apply(), as well as many other functions that we will encounter, have a variety of useful overloads, see the official API documentation for a detailed library reference. In the next two lines we display the newly created Bell state.

In the line

auto [m, probs, states] = measure(result, gt.H, {0});

we use the function qpp::measure() to perform a measurement of the first qubit (subsystem {0}) in the X basis. You may be confused by the apparition of gt.H, however this overload of the function qpp::measure() takes as its second parameter the measurement basis, specified as the columns of a complex matrix. In our case, the eigenvectors of the X operator are just the columns of the Hadamard matrix. As mentioned before, as all other library functions, qpp::measure() returns by value, hence it does not modify its argument. The return of qpp::measure is a tuple consisting of the measurement result, the outcome probabilities, and the possible output states. Technically qpp::measure() returns a std::tuple<> with 3 elements

std::tuple<qpp::idx, std::vector<double>, std::vector<qpp::cmat>>

The first element represents the measurement result, the second the possible output probabilities and the third the output output states. We use the C++17 structured binding to bind the results to m, probs, and states, respectively.

Finally we display the result. Note that IOManipContainerOpts{}.set_sep(", ") is a display formatting option for STL-like containers, defined in qpp/options.hpp

Quantum operations

In examples/quantum_operations.cpp displayed below

// Quantum operations
// Source: ./examples/quantum_operations.cpp

#include <iostream>
#include <vector>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    cmat rho = st.pb00; // projector onto the Bell state (|00> + |11>) / sqrt(2)
    std::cout << ">> Initial state:\n";
    std::cout << disp(rho) << '\n';

    // partial transpose of first subsystem
    cmat rhoTA = ptranspose(rho, {0});
    std::cout << ">> Eigenvalues of the partial transpose "
              << "of Bell-0 state are:\n";
    std::cout << disp(transpose(hevals(rhoTA))) << '\n';

    std::cout << ">> Measurement channel with 2 Kraus operators:\n";
    std::vector<cmat> Ks{st.pz0, st.pz1}; // 2 Kraus operators
    std::cout << disp(Ks[0]) << "\nand\n" << disp(Ks[1]) << '\n';

    std::cout << ">> Superoperator matrix of the channel:\n";
    std::cout << disp(kraus2super(Ks)) << '\n';

    std::cout << ">> Choi matrix of the channel:\n";
    std::cout << disp(kraus2choi(Ks)) << '\n';

    // apply the channel onto the first subsystem
    cmat rhoOut = apply(rho, Ks, {0});
    std::cout << ">> After applying the measurement channel "
              << "on the first qubit:\n";
    std::cout << disp(rhoOut) << '\n';

    // take the partial trace over the second subsystem
    cmat rhoA = ptrace(rhoOut, {1});
    std::cout << ">> After partially tracing down the second subsystem:\n";
    std::cout << disp(rhoA) << '\n';

    // compute the von-Neumann entropy
    double ent = entropy(rhoA);
    std::cout << ">> Entropy: " << ent << '\n';
}

we introduce quantum operations: quantum channels, as well as the partial trace and partial transpose operations. The output of this program is

>> Initial state:
0.5   0   0   0.5
  0   0   0     0
  0   0   0     0
0.5   0   0   0.5
>> Eigenvalues of the partial transpose of Bell-0 state are:
-0.5   0.5   0.5   0.5
>> Measurement channel with 2 Kraus operators:
1   0
0   0
and
0   0
0   1
>> Superoperator matrix of the channel:
1   0   0   0
0   0   0   0
0   0   0   0
0   0   0   1
>> Choi matrix of the channel:
1   0   0   0
0   0   0   0
0   0   0   0
0   0   0   1
>> After applying the measurement channel on the first qubit:
0.5   0   0     0
  0   0   0     0
  0   0   0     0
  0   0   0   0.5
>> After partially tracing down the second subsystem:
0.5     0
  0   0.5
>> Entropy: 1

The example should by now be self-explanatory. In the line

cmat rho = st.pb00; // projector onto the Bell state (|00> + |11>) / sqrt(2)

we define the input state rho as the projector onto the Bell state $(|00> + |11>)/\sqrt{2}$, then display it in the next two lines.

In the lines

cmat rhoTA = ptranspose(rho, {0});
std::cout << ">> Eigenvalues of the partial transpose "
              << "of Bell-0 state are:\n";
std::cout << disp(transpose(hevals(rhoTA))) << '\n';

we partially transpose the first qubit, then display the eigenvalues of the resulting matrix rhoTA.

Next, in the line

std::vector<cmat> Ks{st.pz0, st.pz1}; // 2 Kraus operators

we define a quantum channel Ks consisting of two Kraus operators: $|0><0|$ and $|1><1|$, then display them. Note that Quantum++ uses the std::vector<cmat> container to store the Kraus operators and define a quantum channel.

In the following lines we display the superoperator matrix as well as the Choi matrix of the channel Ks.

Next, in the line

 cmat rhoOut = apply(rho, Ks, {0});

we apply the channel Ks to the first qubit of the input state rho. We then display the output state rhoOut.

In the line

cmat rhoA = ptrace(rhoOut, {1});

we take the partial trace of the output state rhoOut. In the next lines we display the resulting state rhoA.

Finally, in the line

cmat rhoA = ptrace(rhoOut, {1});

we compute the von-Neumann entropy of the resulting state, then display it in the last line.

Timing

To facilitate simple timing tasks, Quantum++ provides a qpp::Timer<> class that uses by default a std::steady_clock. The clock type is passed as a template argument to qpp::Timer<>, hence it can be changed at compile time, e.g., with a std::high_resolution_clock.

The program examples/timing1.cpp listed below demonstrate the usage of qpp::Timer<>.

// Timing
// Source: ./examples/timing1.cpp

#include <iomanip>
#include <iostream>
#include <vector>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    std::cout << std::setprecision(8); // increase the default output precision

    // get the first codeword from Shor's [[9,1,3]] code
    ket c0 = codes.codeword(Codes::Type::SHOR_NINE_QUBIT, 0);

    Timer<> t;                           // declare and start a timer
    std::vector<idx> perm = randperm(9); // declare a random permutation
    ket c0perm = syspermute(c0, perm);   // permute the system
    t.toc();                             // stops the timer
    std::cout << ">> Permuting subsystems according to " << disp(perm, ", ");
    std::cout << "\n>> It took " << t << " seconds to permute the subsytems.\n";

    t.tic(); // restart the timer
    std::cout << ">> Inverse permutation: ";
    std::cout << disp(invperm(perm), ", ") << '\n';
    ket c0invperm = syspermute(c0perm, invperm(perm)); // permute again
    std::cout << ">> It took " << t.toc();
    std::cout << " seconds to un-permute the subsystems.\n";

    std::cout << ">> Norm difference: " << norm(c0invperm - c0) << '\n';
}

A possible output of this program is

>> Permuting subsystems according to [1, 0, 2, 3, 6, 7, 5, 8, 4]
>> It took 0.00016 seconds to permute the subsytems.
>> Inverse permutation: [1, 0, 2, 3, 8, 6, 4, 5, 7]
>> It took 8.4e-05 seconds to un-permute the subsystems.
>> Norm difference: 0

In the line

std::cout << std::setprecision(8); // increase the default output precision

we change the default output precision from to 4 8 decimals after the delimiter.

In the line

ket c0 = codes.codeword(Codes::Type::NINE_QUBIT_SHOR, 0);

we use the qpp::Codes singleton codes to retrieve in c0 the first codeword of Shor’s 9,1,3 quantum error correcting code.

In the following line

Timer<> t;                           // declare and start a timer

we declare an instance timer of the class qpp::Timer<> (which automatically starts the timer). In the next line

 std::vector<idx> perm = randperm(9); // declare a random permutation

we declare a random permutation perm via the function qpp::randperm(). In the line

 ket c0perm = syspermute(c0, perm);   // permute the system

we permute the codeword according to the permutation perm using the function qpp::syspermute() and store the result . Next, in the line

t.toc();                             // stops the timer

we stop the timer. In the lines

std::cout << ">> Permuting subsystems according to " << disp(perm, ", ");

we display the permutation, using an overloaded form of the qpp::disp() manipulator for C++ standard library containers. The latter takes a std::string as its second parameter to specify the delimiter between the elements of the container. In the following line

std::cout << "\n>> It took " << t << " seconds to permute the subsytems.\n";

we display the elapsed time using the ostream operator<<() operator overload for qpp::Timer<> instances.

Next, in the line

t.tic(); // restart the timer

we reset the timer, then display the inverse permutation of perm in the following two lines.

In the lines

ket c0invperm = syspermute(c0perm, invperm(perm)); // permute again
std::cout << ">> It took " << t.toc();
std::cout << " seconds to un-permute the subsystems.\n";

we permute the already permuted state c0perm according to the inverse permutation of perm, and store the result in c0invperm, then we display the elapsed time. Note that we used directly t.toc() in the stream insertion operator, as, for convenience, the member function qpp::Timer<>::toc() returns a const qpp::Timer&.

Finally, in the line

 std::cout << ">> Norm difference: " << norm(c0invperm - c0) << '\n';

we verify that by permuting and permuting again using the inverse permutation we recover the initial codeword, i.e., the norm difference is zero, via qpp::norm() function.

Input/output

We now introduce the input/output functions of Quantum++, as well as the input/output interfacing with MATLAB. The following programs save a matrix in both Quantum++ text format (preserving double precision) as well as in MATLAB format, then loads it back and tests that the norm difference between the saved/loaded matrix is zero. The code displayed below is self-explanatory.

// Input/output
// Source: ./examples/input_output.cpp

#include <fstream>
#include <iostream>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    // Quantum++ input/output
    cmat rho = randrho(256); // an 8 qubit density operator
    {
        std::ofstream fout("rho.txt");
        save(rho, fout); // save it
    }
    std::ifstream fin("rho.txt");
    cmat loaded_rho = load<cmat>(fin); // load it back
    // display the difference in norm, should be 0
    std::cout << ">> Norm difference load/save: ";
    std::cout << norm(loaded_rho - rho) << '\n';
}

The output of this program is

>> Norm difference load/save: 0

// MATLAB input/output
// Source: ./examples/matlab_io.cpp

#include <iostream>

#include "qpp/qpp.h"

#include "qpp/MATLAB/matlab.hpp" // must be explicitly included

int main() {
    using namespace qpp;

    // interfacing with MATLAB
    cmat rho = randrho(256); // an 8 qubit density operator
    save_MATLAB(rho, "rho.mat", "rho", "w");
    cmat loaded_rho = load_MATLAB<cmat>("rho.mat", "rho");
    // display the difference in norm, should be 0
    std::cout << ">> Norm difference MATLAB load/save: ";
    std::cout << norm(loaded_rho - rho) << '\n';
}

The output of this program is

>> Norm difference MATLAB load/save: 0

Note that in order to use the MATLAB input/output interface support, you need to explicitly include the header file qpp/MATLAB/matlab.hpp, and you also need to have MATLAB or the MATLAB compiler installed, otherwise the program fails to compile. See the file README.md for extensive details about compiling with MATLAB support.

Exceptions

Most Quantum++ functions throw exceptions in the case of unrecoverable errors, such as out-of-range input parameters, input/output errors etc. The exceptions are handled via the base class qpp::exception::Exception (note the additional inner namespace qpp::exception) and its descendants.

The base class qpp::exception::Exception is derived from the standard library exception class std::exception. The exception types are hard-coded inside the scoped enumeration (enum class) qpp::exception::Exception::Type. If you want to add more exceptions, augment the enumeration qpp::exception::Exception::Type and also modify accordingly the member function qpp::exception::Exception::construct_exception_msg_(), which constructs the exception message displayed via the overridden virtual function qpp::exception::Exception::what(). Listing  examples/exceptions.cpp displayed below illustrates the basics of exception handling in Quantum++.

// Exceptions
// Source: ./examples/exceptions.cpp

#include <exception>
#include <iostream>

#include "qpq/qpp.h"

int main() {
    using namespace qpp;

    cmat rho = randrho(16); // 4 qubits (subsystems)
    try {
        // the line below throws qpp::exception::SubsysMismatchDims
        double mInfo = qmutualinfo(rho, {0}, {4});
        std::cout << ">> Mutual information between first and last subsystem: ";
        std::cout << mInfo << '\n';
    } catch (const std::exception& e) {
        std::cout << ">> Exception caught: " << e.what() << '\n';
    }
}

The output of this program is

>> Exception caught: qpp::qmutualinfo(): Subsystems mismatch dimensions! [subsysB]

In the line

cmat rho = randrho(16); // 4 qubits (subsystems)

we define a random density matrix on four qubits (dimension 16). In the line

double mInfo = qmutualinfo(rho, {0}, {4});

we compute the mutual information between the first and the 5-th subsystem. The code will throw an exception of type qpp::exception::SubsysMismatchDim, as there are only four systems. We next catch and display the exception in the catch block

catch (const std::exception& e) {
        std::cout << ">> Exception caught: " << e.what() << '\n';
}

via the std::exception standard exception base class. We could have also used the Quantum++ exception base class qpp::exception::Exception, however using the std::exception allows catching other exceptions, not just of the type qpp::exception::Exception and its descendants.

Classical reversible logic

The library offers support for classical reversible logic, via the classes qpp::Dynamic_bitset and qpp::Bit_circuit defined in qpp/classes/reversible.hpp. The program examples/reversible1.cpp displayed below is self-explanatory.

PS: note that for the sake of performance, the qpp::Dynamic_bitset and qpp::Bit_circuit classes do not perform any error/exception checks on their arguments, so please make sure you pass semantically correct arguments to their member functions.

// Classical reversible circuits
// Source: ./examples/reversible1.cpp

#include <iostream>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    std::cout << ">> Classical reversible circuits. ";
    std::cout << "Bits are labeled from right to left,\n   ";
    std::cout << "i.e., bit zero is the least significant bit (rightmost).\n";

    Dynamic_bitset bits{4};                             // 4 classical bits
    std::cout << ">> Initial bitset: " << bits << '\n'; // display them

    bits.rand(); // randomize the bits
    std::cout << ">> After randomization: " << bits << '\n'; // display them

    Bit_circuit bc{bits}; // construct a bit circuit out of a bit set
    std::cout << ">> Bit circuit (constructed from the above bitset):\n"
              << bc << '\n';

    std::cout << ">> Apply X_0, followed by CNOT_02, CNOT_13 and TOF_013\n";
    bc.X(0);                               // apply a NOT gate on the first bit
    bc.CNOT(0, 2).CNOT(1, 3).TOF(0, 1, 3); // sequence operations

    std::cout << ">> Final bit circuit:\n" << bc << '\n';
    std::cout << ">> 3rd bit: " << bc.get(2) << '\n';
    std::cout << ">> CNOT count: " << bc.get_gate_count("CNOT") << '\n';
    std::cout << ">> CNOT depth: " << bc.get_gate_depth("CNOT") << '\n';
}

A possible output of this program is

>> Classical reversible circuits. Bits are labeled from right to left,
   i.e., bit zero is the least significant bit (rightmost).
>> Initial bitset:
	0000
>> After randomization:
	0110
>> Apply X_0, followed by CNOT_02, CNOT_13 and TOF_013
>> Final bit circuit:
	0011
>> 3rd bit: 0
>> CNOT count: 2
>> Reseted circuit:
	0000
>> CNOT count: 0

Quantum noise

Quantum++ offers support for quantum noise models. The following program examples/noise.cpp is self-explanatory.

// Quantum noise
// Source: ./examples/noise.cpp

#include <iostream>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    std::cout << ">> Depolarizing qubit noise acting on a random state\n";
    double p = 0.75;     // depolarizing probability (fully depolarizing)
    idx N = 10000;       // number of trials
    ket psi = randket(); // initial state
    cmat rho = cmat::Zero(2, 2);        // final density matrix
    QubitDepolarizingNoise noise{0.75}; // constructs a noise instance
    // compute the noise action; see qpp::NoiseBase::operator() for other
    // overloads for multi-partite states
    for (idx i = 0; i < N; ++i) {
        // apply the noise
        rho += prj(noise(psi));
    }
    rho /= static_cast<double>(N); // normalize the resulting density matrix
    std::cout << ">> Resulting state after " << N
              << " actions of depolarizing noise with p = " << p << ":\n";
    std::cout << disp(rho) << '\n';
}

A possible output of this program is

>> Depolarizing qubit noise acting on a random state
>> Resulting state after 10000 actions of depolarizing noise with p = 0.75:
                 0.498064   -6.97305e-05 - 0.00122461i
-6.97305e-05 + 0.00122461i                   0.501936

Quantum circuits

Quantum++ offers support for quantum circuits. The following program examples/circuits/teleport_qubit_circuit is self-explanatory.

// Qubit teleportation circuit simulator
// Source: ./examples/circuits/teleport_qubit_circuit.cpp

#include <iostream>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    std::cout << ">> Qubit teleportation quantum circuit simulation\n\n";

    // quantum circuit with 3 qubits and 2 classical bits
    QCircuit qc{3, 2};
    // set the qubit 0 to a random state
    cmat U = randU(2);
    // apply the gate U with name randU to qubit 0
    qc.gate(U, 0, "randU");

    // set the MES between qubits 1 and 2
    qc.gate(gt.H, 1);
    qc.CTRL(gt.X, 1, 2);

    // perform the Bell measurement between qubits 0 and 1
    qc.CTRL(gt.X, 0, 1);
    qc.gate(gt.H, 0);
    qc.measure({0, 1});

    // apply the classical controls
    qc.cCTRL(gt.X, 1, 2);
    qc.cCTRL(gt.Z, 0, 2);

    // initialize the quantum engine with a circuit
    QEngine engine{qc};

    // display the quantum circuit and its corresponding resources
    std::cout << qc << "\n\n" << qc.get_resources() << "\n\n";

    // execute the entire circuit
    engine.execute();

    // display the measurement statistics
    std::cout << engine << "\n\n";

    // verify that the teleportation was successful
    ket psi_in = U * 0_ket;
    ket psi_out = engine.get_state();
    std::cout << ">> Teleported state:\n";
    std::cout << disp(dirac(psi_out)) << '\n';
    std::cout << ">> Norm difference: " << norm(psi_out - psi_in) << '\n';
}

A possible output of this program is

>> Qubit teleportation quantum circuit simulation

[QCircuit nq: 3, nc: 2, d: 2]
0: SINGLE, target = [0], name = "randU"
1: SINGLE, target = [1], name = "H"
2: CTRL, ctrl = [1], target = [2], name = "CTRL-X"
3: CTRL, ctrl = [0], target = [1], name = "CTRL-X"
4: SINGLE, target = [0], name = "H"
5: |> MEASURE_MANY, target = [0, 1], c_reg = 0, name = "mZ"
6: cCTRL, c_ctrl = [1], target = [2], name = "cCTRL-X"
7: cCTRL, c_ctrl = [0], target = [2], name = "cCTRL-Z"

[Resources]
<QCircuit nq: 3, nc: 2, d: 2>
step count: 8
total gate count: 7
total gate depth: 4
total measurement count: 1
total measurement depth: 1
total depth: 5

[QEngine]
<QCircuit nq: 3, nc: 2, d: 2>
last probs: [0.5, 0.5]
last dits: [1, 1]
[Statistics]
	num_reps: 1
	num_outcomes: 1
	[1 1]: 1

>> Teleported state:
|0> * (0.760773 + 0.524461i) +
|1> * (-0.363785 + 0.117581i)
>> Norm difference: 2.77556e-16

Noisy quantum circuits

Quantum++ offers support for noisy quantum circuits (using noisy quantum engines). The following program examples/circuits/noisy_teleport_qubit_circuit is self-explanatory.

// Qubit noisy teleportation circuit simulator
// Source: ./examples/circuits/noisy_teleport_qubit_circuit.cpp

#include <iostream>

#include "qpp/qpp.h"

int main() {
    using namespace qpp;

    std::cout << ">> Qubit noisy teleportation quantum circuit simulation\n\n";

    // quantum circuit with 3 qubits and 2 classical bits
    QCircuit qc{3, 2};
    // set the qubit 0 to a random state
    cmat U = randU(2);
    // apply the gate U with name randU to qubit 0
    qc.gate(U, 0, "randU");

    // set the MES between qubits 1 and 2
    qc.gate(gt.H, 1);
    qc.CTRL(gt.X, 1, 2);

    // perform the Bell measurement between qubits 0 and 1
    qc.CTRL(gt.X, 0, 1);
    qc.gate(gt.H, 0);
    qc.measure({0, 1});

    // apply the classical controls
    qc.cCTRL(gt.X, 1, 2);
    qc.cCTRL(gt.Z, 0, 2);

    // display the quantum circuit and its corresponding resources
    std::cout << qc << "\n\n" << qc.get_resources() << "\n\n";

    // initialize the noisy quantum engine with an amplitude damping noise model
    // and a quantum circuit
    QNoisyEngine noisy_engine{qc, QubitAmplitudeDampingNoise{0.99}};

    // execute the entire circuit
    noisy_engine.execute();

    // display the measurement statistics
    std::cout << noisy_engine << "\n\n";

    // verify how successful the teleportation was
    ket psi_in = U * 0_ket;
    ket psi_out = noisy_engine.get_state();
    std::cout << ">> Initial state:\n";
    std::cout << disp(psi_in) << '\n';
    std::cout << ">> Teleported state:\n";
    std::cout << disp(psi_out) << '\n';
    std::cout << ">> Norm difference: " << norm(psi_out - psi_in) << '\n';
}

A possible output of this program is

>> Qubit noisy teleportation quantum circuit simulation

[QCircuit nq: 3, nc: 2, d: 2]
0: SINGLE, target = [0], name = "randU"
1: SINGLE, target = [1], name = "H"
2: CTRL, ctrl = [1], target = [2], name = "CTRL-X"
3: CTRL, ctrl = [0], target = [1], name = "CTRL-X"
4: SINGLE, target = [0], name = "H"
5: |> MEASURE_MANY, target = [0, 1], c_reg = 0, name = "mZ"
6: cCTRL, c_ctrl = [1], target = [2], name = "cCTRL-X"
7: cCTRL, c_ctrl = [0], target = [2], name = "cCTRL-Z"

[Resources]
<QCircuit nq: 3, nc: 2, d: 2>
step count: 8
total gate count: 7
total gate depth: 4
total measurement count: 1
total measurement depth: 1
total depth: 5

[QNoisyEngine]
<QCircuit nq: 3, nc: 2, d: 2>
last probs: [0.502513, 0.484558]
last dits: [0, 1]
[Statistics]
	num_reps: 1
	num_outcomes: 1
	[0 1]: 1

>> Initial state:
0.236511 - 0.925033i
0.0193838 + 0.29665i
>> Teleported state:
  0.236406 - 0.92462i
0.0194728 + 0.298012i
>> Norm difference: 0.00142931

OpenQASM

Quantum++ supports OpenQASM (i.e., can generate qpp::QCircuits from OpenQASM, but not the other way around, as qpp::QCircuit language description supports significantly more features than OpenQASM). The program examples/qasm/qpp_qasm is self-explanatory and executes (repeatedly, if the number of repetitions is passed as an argument) an OpenQASM file read from the standard input. You can use it to test your favourite OpenQASM programs.

// Executes an OpenQASM program read from the input stream, repeatedly if the
// number of repetitions is passed as the first argument. If there is a second
// argument (i.e., argc == 3), then the final quantum state is displayed. If
// there are three or more arguments (i.e., argc  > 3), the projector onto the
// final state is displayed.
// Source: ./examples/qasm/qpp_qasm.cpp

#include <iostream>

#include "qpp/qpp.h"

int main(int argc, char** argv) {
    using namespace qpp;

    // read the circuit from the input stream
    QCircuit qc = qasm::read(std::cin);

    // initialize the quantum engine with a circuit
    QEngine q_engine{qc};

    // display the quantum circuit and its corresponding resources
    std::cout << qc << "\n\n" << qc.get_resources() << "\n\n";

    // execute the quantum circuit
    idx reps = argc > 1 ? std::stoi(argv[1]) : 1; // repetitions
    q_engine.execute(reps);

    // display the measurement statistics
    std::cout << q_engine << '\n';

    // display the final state on demand
    if (argc == 3) {
        std::cout << ">> Final state (transpose):\n";
        std::cout << disp(transpose(q_engine.get_state())) << '\n';
    } else if (argc > 3) {
        std::cout << ">> Final density matrix:\n";
        std::cout << disp(prj(q_engine.get_state())) << '\n';
    }
}

If run with an input redirected from an OpenQASM file, such as this OpenQASM teleportation example

qpp_qasm <./teleport_minimal.qasm

the program will output (up to the intrinsic random measurement results)

[QCircuit nq: 3, nc: 2, d: 2]
0 : SINGLE, target = [0], name = "U"
1 : JOINT, target = [1], name = "h"
2 : JOINT, target = [1, 2], name = "cx"
3 : JOINT, target = [0, 1], name = "cx"
4 : JOINT, target = [0], name = "h"
5 : |] MEASURE_ND, target = [0], c_reg = 0, name = "mZ"
6 : |] MEASURE_ND, target = [1], c_reg = 1, name = "mZ"
7 : cCTRL, c_ctrl = [1], target = [2], shift = [0], name = "x"
8 : cCTRL, c_ctrl = [0], target = [2], shift = [0], name = "z"
9 : SINGLE, target = [2], name = "U"

[Resources]
<QCircuit nq: 3, nc: 2, d: 2>
step count: 10
total gate count: 8
total gate depth: 5
total measurement count: 2
total measurement depth: 1
total depth: 6

[QEngine]
<QCircuit nq: 3, nc: 2, d: 2>
last probs: [0.5, 0.5]
last dits: [0, 1]
[Statistics]
	num_reps: 1
	num_outcomes: 1
	[0 1]: 1