Hi,

I want to report a bug in qpp::applyCTRL function, that there is a discrepancy in the behavior of the function regarding vectors and density matrices. Following is the code, that reproduces the bug.

qpp::ket psi (8) ; //some random 3 qubit pure state
// initialize psi
qpp::cmat rho = psi*qpp::adjoint(psi); //density matrix

qpp::cmat U (2,2) ; // some random unitary
// initialize U

auto A = qpp::applyCTRL(psi,U,{0,1},{2},{2,2,2}) ; // apply channel on pure state
auto B = qpp::applyCTRL(rho,U,{0,1},{2},{2,2,2}) ; // apply channel on density matrix

// Now A*adjoint(A) and B have to be equal. But from the code, they are not. Or may be, I am missing something much trivial.

I suspect the bug is in the density matrix part, and the pure state or vector part is absolutely fine.

When I look at this file, I notice you have got the extra r and cswap gates. That is cool. So, I use that as my template to make my own qelib1.inc, which includes more of my own custom gates.

That does not work if I do it verbatim. The trouble is on these lines:

    "gate cswap a,b,c { cx c,b; ccx a,b,c; cx c,b; }\n"
    "gate ch a,b { s b; h b; t b; cx a,b; tdg b; h b; sdg b; }\n"
    "gate ccx a,b,c { h c; cx b,c; tdg c; cx a,c; t c; cx b,c; tdg c; cx a,c; "

The ccx is referenced in the first line before it is defined later. Well, I can fix my qelib1.inc easily.

I gather you may be interested in assessing whether those lines affect other parts of the system.

This question is motivated by using staq to map a small circuit to Tokyo. That results in this in the output QASM:

qreg q[20];

The small circuit (for educational purposes) only needs first 4 qubits, so I discard() the rest of the qubits (q[4]...q[19]). That drastically reduces the size of the ket returned by the QEngine.get_psi(). That is very good.

The trouble is that the time to run QEngine.execute() is almost same with or without the discard().

The 16 qubits do not involve any quantum operations. I should expect a huge performance boost after discarding them, shouldn't I?

Say I have two qubits psi_1 and psi_2. I entangle them and wish to measure them separately. Please consider the code below taken from example, slightly edited to illustrate my question:

// Measurements
// Source: ./examples/measurements.cpp
#include <iostream>
#include <tuple>

#include "qpp.h"

int main() {
    using namespace qpp;
    ket psi_1 = 0_ket;
    ket psi_2 = 0_ket;
    ket psi = kron(psi_1, psi_2);
    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';

    // measure the first qubit in the Z basis
    std::cout << "Measuring the first ket:\n";
    auto measured_one = measure(result, gt.Z, {0});
    std::cout << ">> Measurement result: " << std::get<0>(measured_one) << '\n';

    std::cout << "\n\nMeasuring the second ket:\n";
    auto measured_two = measure(result, gt.Z, {1});
    std::cout << ">> Measurement result: " << std::get<0>(measured_two) << '\n';
}

If the psi_1 measurement result is <0>, psi_2 measurement result must also be zero for a Bell-00 state.
But if one runs the code above a few times, it will produce results were psi_1 results in <0> while psi_2 results in <1>.
I am aware that if I use measure_seq I can measure both at the same time and get consistent results but things get muddled quickly.

Is there a way to maintain consistency of results without having to measure all the subsystems at once using measure_seq?

I tried with cmake -DBUILD_TESTING:BOOL=ON && gmake test but got:

-- Build files have been written to: /usr/ports/misc/quantum++/work/.build
ninja: no work to do.
[0/1] cd /usr/ports/misc/quantum++/work/.build && /usr/local/bin/ctest --force-new-ctest-process
Test project /usr/ports/misc/quantum++/work/.build
    Start 1: unit_tests_NOT_BUILT
Could not find executable unit_tests_NOT_BUILT
Looked in the following places:
unit_tests_NOT_BUILT
unit_tests_NOT_BUILT
Release/unit_tests_NOT_BUILT
Release/unit_tests_NOT_BUILT
Debug/unit_tests_NOT_BUILT
Debug/unit_tests_NOT_BUILT
MinSizeRel/unit_tests_NOT_BUILT
MinSizeRel/unit_tests_NOT_BUILT
RelWithDebInfo/unit_tests_NOT_BUILT
RelWithDebInfo/unit_tests_NOT_BUILT
Deployment/unit_tests_NOT_BUILT
Deployment/unit_tests_NOT_BUILT
Development/unit_tests_NOT_BUILT
Development/unit_tests_NOT_BUILT
Unable to find executable: unit_tests_NOT_BUILT
1/1 Test #1: unit_tests_NOT_BUILT .............***Not Run   0.00 sec

Here is a snippet. The bell_state is (|00> + |11>) / sqrt(2). We measure them like this:

std::tie(a, b, c) = measure(bell_state, kron(gt.H, gt.H), {0, 1});

If we print out (disp) the bell_state, a, b, and c , they look like:

0.707107 
       0 
       0 
0.707107 

Result: 3

Prob: [0.5, 0, 0, 0.5]

Post-measurement:
1 

0 

0 

1 

It is easy to understand. The post-measurement state vector has 4 elements. Now, if we measure only one bit:

std::tie(a, b, c) = measure(bell_state, gt.H, {0});

We print them out, and they look like:

0.707107 
       0 
       0 
0.707107 

Result: 1

Prob: [0.5, 0.5]

0.707107 
0.707107 

 0.707107 
-0.707107 

It is difficult for me to figure out the post-measurement vector, which has two elements. Each element has a pair of values. Are they real and imaginary parts of two complex numbers like this?

1/sqrt(2) + i * 1/sqrt(2)

1/sqrt(2) - i * 1/sqrt(2)

Could you please explain how the post-measurement vector is produced in the latter (1 bit measurement) case?

Thanks.

Why bother to multiply the randH(D) by 10? Did I miss something? Since all the elements in the matrix are of double data type, I seem unable to see the merits of the adjustment of the factor of 10.

I did a test. It produces the same final results with or without 10.

inline cmat randrho(idx D = 2) {
    //Exception handling...

    cmat result = 10 * randH(D);
    result = result * result.adjoint();

    return result / result.trace();
}

I run the Quantum++ example qasm2 on this QASM file:

OPENQASM 2.0;
include "qelib1.inc";

qreg q[1];
rz(0) q[0];

The output is like this:

>> BEGIN CIRCUIT
nq = 1, nc = 0, d = 2, name = ""
0 CUSTOM, target = [0], name = "rz"
step count: 1
total gate count: 1
total gate depth: 1
total measurement count: 0
total measurement depth: 0
total depth: 1
>> END CIRCUIT

>> BEGIN ENGINE STATISTICS
last probs: []
last dits: []
>> END ENGINE STATISTICS

>> Final state:
0   0 

The final state is not valid, is it?

The example qpp_qasm.cpp takes this snippet as an input:

OPENQASM 2.0;
include "qelib1.inc";

qreg q[2];
qreg qqq[1];
x q[0];

creg c[2];
measure q -> c;

// if (c==1) x qqq[0]; // This is to prove the c is equal to one here.

creg ccc[1];
measure qqq[0] -> ccc[0];

We will get this:

last probs: [1, 1, 1]
last dits: [1, 0, 0]
stats:
	reps: 100
	outcomes: 1
	[1 0 0]: 100
>> Final state (transpose):
0  0  0  0  1  0  0  0

That looks intuitive because the "last dits" can be evaluated to 4 if the 1 is the most significant bit, and the final state matches that.

But, if we look at theif (c==1)in the QASM code above, the classical register c is evaluated to 1 at that point. That seems to suggest the 1 in the "last dits" is the least significant bit.

I cannot find anything specific to bit numbering of the classical registers in the QASM spec. If the particular evaluation of classical registers is a Quantum++ implementation decision, would it be nice to have an engine function to print out the values of the classical registers in addition to the existing statistics? That would clearly indicate the presentation of the "last dits" and the classical register values are different things.

More important, I hope the bit numbering of other quantum platforms is compatible with that of the Quantum++.

Thanks.

In singleton.h file, get_thread_local_instance() contains thread_local keyword in return type, like...

 thread_local static T& get_thread_local_instance()
    noexcept(std::is_nothrow_constructible<T>::value)
    {
        // Guaranteed to be destroyed.
        // Instantiated on first use.
        // Thread safe in C++11
        thread_local static T instance;

        return instance;
    }

I think thread_local keyword should be removed from the return type declaration, as the instance is already constructed as thread_local and thus it has static duration, and the function returns by its static ref. So the returned object is automatically becomes thread_local. Also clang complains about this keyword in the return type. I have removed the keyword, and it is working as it should be, like the following,

 static T& get_thread_local_instance() // no thread_local needed here
    noexcept(std::is_nothrow_constructible<T>::value)
    {
        // Guaranteed to be destroyed.
        // Instantiated on first use.
        // Thread safe in C++11
        thread_local static T instance;

        return instance;
    }

Kindly consider the following example:

qpp::ket ket_1 = 1_ket;
qpp::ket ket_2 = 0_ket;
qpp::ket ket_3 = 0_ket;

// take the tensor product of ket_1 and ket_2
qpp::ket result_1 = qpp::kron(ket_1, ket_2);

// take the tensor product again of result_1 and ket_3
qpp::ket result_2 = qpp::kron(result_1, ket_3);

// take the tensor product of ket_3 and result_1
qpp::ket result_3 = qpp::kron(ket_3, result_1);

What I wish to know are the indices of the different subsystems after the application of each tensor product:

• In result_1, what is the index of ket_1 and of ket_2:
• In result_2, what is the index of ket_1, ket_2 and ket_3?
• In result_3, what is the index of ket_1, ket_2 and ket_3?

What is the general rule to follow to find the indices?
Once those indices are found, what is the right way to use them? For instance, consider the following:

qpp::ket ket_1 = 1_ket;
qpp::ket ket_2 = 0_ket;
qpp::ket ket_3 = 0_ket;

// take the tensor product of ket_1 and ket_2
qpp::ket result_1 = qpp::kron(ket_1, ket_2);
// is this the correct way to apply the hadamard gate to ket_2 in result_1? 
result_1 = qpp::apply(result_1, gt.H, {ket_2_index});

// take the tensor product again of result_1 and ket_3
qpp::ket result_2 = qpp::kron(result_1, ket_3);
// Is the line below the correct way to apply controlled NOT with ket_1 as control and ket_3 as target?
qpp::ket result_2 = qpp::applyCTRL(result_2, gt.X, {ket_1_index}, {ket_3_index});
// Or this the correct way?
qpp::ket result_2 = qpp::applyCTRL(result_2, gt.X, {0}, {1});

This doesn't seem to be documented in the reference manual or the wiki.

Help again greatly appreciated.

I am looking for the doc folder with the reference manual, and I cannot find it

I think the example has a minor issue. The qpp convention is that the q0 is the most significant bit, but the example uses it as the least significant bit.

The developer may be given the freedom to interpret the endianness in his own way as long as he arranges the code subsequently to treat q0 as the least significant bit. But, that may confuse some users down the road, since the examples serve the purpose of documentation.

I'd like to report some documentation errors. Well, I don't know if that is a doxygen problem or your problem.

A lot of images are missing. For example, I access the applyCTRL() API html doc and see the missing image pointing to qpp/doc/html/form_0.png, which does not exist.

Please see the attached screenshot.

Here is a test case.

    auto a = gt.T;
    std::cout<<disp(a);
    std::cout<<"\n--------------------------------\n\n";

    std::cout<<disp(hevals(a));
    std::cout<<"\n--------------------------------\n\n";

    std::cout<<disp(evals(a));
    std::cout<<"\n--------------------------------\n\n";

The same matrix produces two different sets of eigenvalues:

1                     0 
0   0.707107 + 0.707107i
--------------------------------

0.707107 
       1 
--------------------------------

                  1 
0.707107 + 0.707107i
--------------------------------

The evals() is correct and hevals() is not. The T gate is not even Hermitian, is it? Should hevals() throw an exception or something when its input is not Hermitian?

The problem seems to boil down to the the Eigen::SelfAdjointEigenSolver<cmat>, which does not check whether the input is Hermitian or not.

The regetti / grove API has this very convenient and useful function:

create_arbitrary_state()

I used to think it would be difficult to achieve that, but I am amazed a function of about a hundred lines of Python code based on a 5 page paper could transform 0 to any states with only pow(2, n+2) - 4n - 4 CNOT gates and pow(2, n+2) - 5 one-qubit elementary rotation gates.

I cannot find a similar function in the Quantum++. If it doesn't exist, I think it would be nice to implement it.

This issue may be related to the staq Issue 34, but I am not sure. So, I open this one to track the problem.

On qpp:

OPENQASM 2.0;
include "qelib1.inc";
qreg q[2];
rz (1.1544377492147073) q[0];
h q[1];
u1 (-1.1544377492147073) q[0];
h q[0];
ry (0.501899064626972) q[0];
rz (-2.664259129025202) q[0];
swap q[0], q[1];
cx q[1], q[0];
ry (-0.354312360364173) q[0];
cx q[1], q[0];
ry (-0.1912917581366851) q[0];
cx q[1], q[0];
rz (0.9715277168265499) q[0];
cx q[1], q[0];
rz (0.23743868408440671) q[0];

The final state:

0.4439      
-0.0677419-0.552801i     
0.0884309+0.233639i     
-0.495458-0.429988i

That can be compiled to run on the Rigetti simulator:

RZ(1.1544377492147073) 0
PHASE(-1.1544377492147073) 0
H 0
H 1
RY(0.501899064626972) 0
RZ(-2.664259129025202) 0
SWAP 0 1
CNOT 1 0
RY(-0.354312360364173) 0
CNOT 1 0
RY(-0.1912917581366851) 0
CNOT 1 0
RZ(0.9715277168265499) 0
CNOT 1 0
RZ(0.23743868408440671) 0

The final state:

[ 
0.43888697+0.0665232j      
0.0524189+0.24425298j      
0.0158663-0.55670996j      
-0.42542397-0.49938197j
]

Would you please consider setting that to YES in the Doxyfile?

I think the search engine costs very little, but offers a lot of convenience. After all, the staq project sets that to YES by default.

Otherwise, the git would have to deal with that whenever I pull the master, after I make a local change to YES.

Thanks.

You published this nearly a year ago. Any updates on benchmarks?

You pointed out:

The most likely explanation is that the former [QuTiP] uses sparse matrices during computation, whereas the latter [Quantum++] does not.

The fundamental data types in Quantum++ are non-sparse vectors and matrices.

My questions:

Do you know by now the QuTiP outperforming Quantum++ in that test is indeed attributed to the sparse matrices?

Could you please elaborate a little more on why your design decision favors non-sparse vectors and matrices considering the apparent performance benefits of using their sparse counterparts?

Thanks.

Have you read this paper?

I understand Quantum++ has the MODMUL to facilitate Shor's algorithm, but an oracle built from fundamental quantum gates would be much more intuitive and compelling. I have googled and found that a lot of papers are trying to achieve that with little success, let alone implementing that in a C++ library.

If Quantum++ implements a function that could return a string that represents an oracle in QASM format for , would that be a first in the industry?

Right now, I have to make that oracle piece by piece following the recipe of that paper, but have not succeeded. It is quite a hassle! Waiting for the beacon of hope from Quantum++! LOL!

OPENQASM 2.0;
include "qelib1.inc";

qreg q[1];
rx(3e-1) q[0];

That will crash because of the parsing error of the 3e-1.

That was fixed at staq, but not here.

Here are some simple test cases:

OPENQASM 2.0;
include "qelib1.inc";

qreg q[2];

h q;

// Case 1: Crashes
// U(+ 0.1 - 0.4  , 0, 0) q[1];
// U(+0.1 - 0.4 , 0, 0) q[1];

// Case 2: Wrong final states 
// U(- 0.4 + 0.1 , 0, 0) q[1];
// U(-0.4 + 0.1 , 0, 0) q[1];

// Case 3: Correct
// U(0.1 - 0.4 , 0, 0) q[1];

// Case 4: Correct
U(-0.3 , 0, 0) q[1];

Essentially, all the cases should produce the same final state as the Case 4 does.

Just a friendly reminder: a similar issue was fixed yesterday at staq.

Suppose we do this:

    cmat rho = randrho(4);
    cmat U = randU(4); 
    std::tie(result, probs, states) = measure(rho, U, {0,1});

It will produce the information similar to what we do it by hand like this:

adjoint(U) * rho * U

This lecture and this lecture (slide 8), however, seem to like the following convention. By the way, the latter is from your very own Michele Mosca.

U * rho * adjoint(U)

The Qiskit API has the DensityMatrix.evolve(U) that also follows

U * rho * adjoint(U)

When I exchange data between Quantum++ and Qiskit, I have to flip the U and adjoint(U) to suit both APIs.

Well, I understand it is technically correct for Quantum++ to place the adjoint(U) first because it is still unitary by definition.

Could you please comment on why Quantum++ API chooses to do business that way?

Here is a test case:

OPENQASM 2.0;
include "qelib1.inc";

gate Aa a
{
    x a;
}

qreg q[1];

That will fail because of the name Aa.

According to the Open QASM spec, the regular expression of the <id> is [a-z][A-Za-z0-9_]*

The Aa is a match, isn't it?

Here are two QASM files. I expect them to produce the same final states, but they don't. Could you please comment on them?

This file produces the correct state:

OPENQASM 2.0;
include "qelib1.inc";

qreg ctrl[1];
qreg q[1];

x ctrl[0];

h q[0];
rz((pi/4)+(pi/2)) q[0];

Its result is:

>> Final state:
0   0   0.707107   -0.5 + 0.5i

But, this file doesn't.

OPENQASM 2.0;
include "qelib1.inc";

qreg ctrl[1];
qreg q[1];

x ctrl[0];

ch ctrl[0], q[0];
crz((pi/4)+(pi/2)) ctrl[0], q[0];

Its result is:

>> Final state:
0   0   0.270598 - 0.653281i  0.270598 + 0.653281i

Well, it is almost correct -- the probabilities are correct. The phases of the state seem to be wrong.

The second file basically uses the ctrl[0] to control the other qubit. Since the ctrl[0] is 1, the two circuits should be deemed identical.

This cannot be compiled:

#include "qpp.h"
int main()
{
    return 0;
}

This is fine:

#include <iostream>
#include "qpp.h"
int main()
{
    return 0;
}

So, the qpp.h is not self-contained because of the external dependency on iostream.

Not sure if that was recently introduced or not.

This function is very convenient.

QCircuit qpp::qasm::read_from_file(const std::string& fname)

Intuitively, however, one would expect its counterpart function -- something like this:

void qpp::qasm::write_to_file(const std::string& fname, QCircuit qc)

or

std::string QCircuit::to_QASM()

After all, it has already got a to_JSON(). Translation from JSON to QASM looks easy, but it is not trivial to get it done thoroughly. Would you please take the issue into consideration in the future releases?

One of my use cases is to have a distributed graphical interface to edit and test many QCircuit objects and save them to a QASM file when the testing is done. Then, I could run the QASM file in a physical quantum computer such as the IBM Q Experience.

As long as I can get the QASM, I can use this tool to translate it to pyQuil and run it on Rigetti Quantum Cloud too.

I notice the base of the logarithm in the entropy() is 2. It is also documented in the API as such.

result -= ev(i) * std::log2(ev(i));

Many papers, however, tend to use the natural logarithm. So does the Qiskit.

I understand base 2, e, or 10 can be used to calculate entropy. Is your choice of the base 2 just a personal preference, or is it conceptually significant in the context of Quantum++?

Here is a QASM file, representing the value 1.

OPENQASM 2.0;
include "qelib1.inc";

qreg a[3];

x a[0]; // a = 001

If we run it, we will get the final state:

0 0 0 0 1 0 0 0

The final state indicates the value 4. That suggests each follows a different endianness. I guess the latter is from a tensor product of |a[0]> |a[1]> |a[2]>. Well, is there anything wrong to produce the final state from a tensor product of |a[2]> |a[1]> |a[0]>? One would think the final state should be:

0 1 0 0 0 0 0 0

You explained the endianness of the classical registers clearly in the Issue 75. Would you please comment on this issue?

The staq issue has been fixed.

It turns out that the qpp has a similar issue. If you use the QEngine to execute this simple file, i.e., a file with a comment line without a newline/linefeed at the end, the engine will hang.

OPENQASM 2.0;
include "qelib1.inc";

qreg a[1];
//

Suppose I have this:

auto m_16x16 = gt.MODMUL(a, N, 4);

I hope there is another function to achieve this:

std::vector<cmat> d= gt.MODMUL(a, N, 4);
auto a_2x2 = d[0];
auto b_2x2 = d[1];
auto c_2x2 = d[2];
auto d_2x2 = d[3];

auto diff = norm(m_16x16  - kron(a_2x2, b_2x2, c_2x2, d_2x2));

I want to get diff=0. That would be so nice to work on those 2x2 matrices on each of the 4 qubits.

I understand it would be very difficult to process an arbitrary unitary. But, for this particular gt.MODMUL(a, N, n), is it easy for you to implement it?

when i execute
g++ -I ../include/ minimal.cpp -o minimal In file included from minimal.cpp:5: ../include/qpp.h:81:10: fatal error: Eigen/Dense: No such file or directory 81 | #include <Eigen/Dense> | ^~~~~~~~~~~~~

I constructed the following circuit in qpp:

 >> probability q0: [0.5, 0.5]
 >> probability q1: [0.75, 0.25]

However, the actual measurement probabilities were

 >> probability q0: [0.5, 0.5]
 >> probability q1: [1, 0]

The qpp code I used is here:

using namespace qpp;

std::vector<idx> qubits{0, 0};

ket psi = mket(qubits);
ket result = psi;

cmat U(2,2);
U << 1, 0, 0, omega(-4); // e^-i*pi/2

result = apply(result, gt.H, {0});
result = apply(result, gt.H, {1});
result = applyCTRL(result, U, {0}, {1});
result = apply(result, gt.H, {1});

auto measured0 = measure(result, gt.Z, {0});
auto measured1 = measure(result, gt.Z, {1});

std::cout << ">> probability q0: ";
std::cout << disp(std::get<PROB>(measured0), ", ") << '\n';

std::cout << ">> probability q1: ";
std::cout << disp(std::get<PROB>(measured1), ", ") << '\n';

/*
 
 EXPECTED (CORRECT) OUTPUT
 >> probability q0: [0.5, 0.5]
 >> probability q1: [0.75, 0.25]

 ACTUAL (INCORRECT) OUTPUT
 >> probability q0: [0.5, 0.5]
 >> probability q1: [1, 0]
 
 */

In both the Quirk simulation and the qpp program I use

for λ = -π/2.

The wiki article here states:

Repeating the same measurement without any evolution of the quantum state will lead to the same result. If the preparation is repeated, subsequent measurements will likely lead to different results.

This simple test shows the 10 measurements do often produce different results (a).

    for(int i=0; i<10; i++)
    {
        std::tie(a, b, c) = measure(bell_state, gt.H, {0});
    }

Does that imply the experiments simulated in your Quantum++ assume the preparation is indeed repeated?

In a physical quantum computer, that would represent a scenario where the input device keeps on producing the bell state. Is that correct?

Here is a code snippet from the shor example.

    for (idx i = 0; i < n; ++i) {
        // compute 2^(n-i-1) mod N
        idx j = static_cast<idx>(std::llround(std::pow(2, n - i - 1)));
        // compute the a^(2^(n-i-1)) mod N
        idx aj = modpow(a, j, N);
        // apply the controlled modular multiplication
        psi = applyCTRL(psi, gt.MODMUL(aj, N, n), {i}, second_subsys);
    }

That shows how to use the gt.MODMUL() properly.

The doc says something like this:

Note
... a and N should be co-prime.

Parameters
a Positive integer less than N

The a in the documentation basically means the aj in the code, doesn't it? Since the gate does not check co-primality, it is the user's responsibility to make sure the argument such as modpow(a, j, N) is a co-prime of N. Is that correct?

I am not sure if this should be filed against qpp or staq. Anyway, here is a simple test case.

OPENQASM 2.0;
include "qelib1.inc";

qreg h[1];
qreg q[1];
h q[0];

That cannot be run by qpp because the h is declared as a gate and a register. The trouble is that the script can be parsed by the staq perfectly.

I think they should succeed or fail consistently.

Two of the signatures of the function have this:

const std::vector< idx > & | target,

But the textual description of the function says:

Applies the single qudit controlled gate U...

Parameters

U Single qudit quantum gate

When the target has more than one qudit in it, the U won't be a single qudit gate, will it?

Reallocations of the state matrix happen at each

const dyn_mat<typename Derived1::Scalar>& rstate = state;

line in operations.h and each default copy thereof. I suspect this is because the reference is upcasted when passed as const Eigen::MatrixBase<Derived1>& state in the function argument and then needs to be converted back to a dyn_mat. For example, calling

psi = qpp::apply(psi, qpp::gt.X, {1});

from my code, a total of nine copies of the state are created within operations.h:

1. at line 578 (apply overload with idx d)
2. at line 492 (apply overload with dims)
3. at line 61 (applyCTRL)
4. for a copy reference in coeff_idx_ket (the copy constructor of Matrix is demonstrably called at the moment the coeff_idx_ket is defined) – is there a reason for it to be [=] anyway?
5. at line 358 (result = rstate),

along with temporary copies for the purposes of internal::_check_cvector and internal::_check_dims_match_cvect (at lines 532, 535, 349, 352). (Line numbers referring to f0765d7.) Obviously only number 5 is necessary.

This is no problem for small state vectors or density matrices but when the data reaches 1 GB it becomes a real issue.

It appears I am unable to either understand the manual or be imaginative enough.
I have two separate kets and I wish to create an entanglement between them.
What I seek to accomplish is something akin to the snippet below:

qpp::ket ket1 = 0_ket;
qpp::ket ket2 = 0_ket;

// we apply the hadamard gate to the first ket
ket1 = qpp::gt.H * ket1;

// Now here, I want to apply CNOT between ket1 and ket2.
// In fact, in my use case, it can be any controlled gate between two kets.
qpp::ket result = applyCTRL(gt.CNOT, ket1, ket2); // I seek something similar to this
// Later, if I measure ket1, the state of ket2 should be automatically determined as expected of a Bell pair.

I'm using QPP in a project of mine at Avalon PL where users declare qubits that are represented internally as kets so I'm trying to figure out how QPP is to handle things like entanglement between two separate variables.
Any and all suggestions will be extremely useful.

A simple Google search shows many papers about the universal quantum gates.

Can gt.MODMUL be constructed with some universal quantum gates? If impossible, constructing it with even simpler gates of the existing API would be better than the current implementation of the gt.MODMUL, which starts from scratch.

More general question is what criteria you use to include certain gates in the API. They don't seem to be built up from any universal quantum gates, do they?

Hi!
The following code:

typedef std::complex<double> scalar;
void pretty_print(const qpp::ket& psi)
{	
    for (auto i = 0; i < psi.rows(); i++) {
        const scalar & v = psi(i);
        std::cout << std::bitset<8>(i++).to_string() << ": " <<  v << std::endl;
    }
}

gives a differently-ordered psi's as compared to qpp::disp, and it's not a simple reversed order, but some reshuffling I cannot decipher.
Why is that? What's the internal ordering qpp::disp prefers? I couldn't figure it out looking at your code.

Thanks!

After my evaluation of the design of your library and its API, I highly appreciate it, especially these design decisions, as you pointed out in your publication.

...library written in C++11 and composed solely of header files.

The design is inspired from functional programming...

You may extend it to incorporate graphical input, ... or build a more sophisticated library on top of it.

When I think about the header-only approach, I believe all the open source C/C++ library projects ought to be header-only. I am afraid those library developers who reject the header-only doctrine have forgotten why the header-source separation concept was introduced in the first place. When C was invented in 1972, there was no open source concept!

A lot of libraries have been contaminated heavily with GUI (Qt, etc., even worse with the Windows SDK) inadvertently or intentionally! Either way, that is a sin.

I am glad you offer a clean and sharp library. Equipped with the generous MIT license, that would enable the community to do a lot of good things in the years to come.

Thank you.

Although the official QASM spec does not say much about the ancilla, the staq does produce the code that requires ancilla. Consequently, as you know, the code cannot be run by the qpp QEngine.

Would it be nice for the qpp to support the ancilla?

How do you think of the subject?

Thanks.

Header-only libraries should install headers into include/.

Also it's best to install many headers into a dedicated directory under include/, for example include/qpp.

The current implementation of the gate RZ is simply to call Rn(theta, {0, 0, 1}). That is good. The current definition of the rz in the QASM is gate rz(phi) a { u1(phi) a; }. They will end up with two different matrices. Both are correct, but they just look different. Let me illustrate them in this image:

The Rn(theta, {0, 0, 1}) will be in the form of "A"; gate rz(phi) a { u1(phi) a; } "B".

Since the latter looks prettier, would it be nice to shift the phase by θ/2, as suggested in the image, before the cmat RZ(double theta) const returns?

How do you think?

I'd like to suggest you put more references into the API doc or source code comments, preferably as what SciPy people do in their "References" section. I don't mean that you spend a lot of time to put a "References" section into every piece of code. But, some references would be very helpful for understanding your particular implementation of those components of quintessential nature of quantum computing.

Thanks.

I have two simple QASM files:

//test1.qasm
OPENQASM 2.0;
include "qelib1.inc";

qreg q[1];
h q;

//test2.qasm
OPENQASM 2.0;
include "qelib1.inc";

qreg q[1];
u2(0,pi) q;

Let's pass test1.qasm to this program to run first, then test2.qasm.

#include <iostream>
#include <cassert>
#include "qpp.h"

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

    QCircuit qc = qasm::read_from_file(argv[1]);
    QEngine q_engine{qc};
    q_engine.execute();

    std::cout << disp(q_engine.get_psi()) << '\n';
}

Running test1.qasm produces this result:

0.707107
0.707107

and test2.qasm this result:

-0.707107i
-0.707107i

Logically, however, we should expect the same results because both the QASM specification and the Quantum++ implementation, which has faithfully followed the spec, define this:

// Clifford gate: Hadamard
gate h a { u2(0,pi) a; }

I did some analysis on the issue. Using Hadamard gate directly like the one in the test1.qasm is always correct because the H gate is hard-coded correctly in the Quantum++, bypassing the definition in the qelib1.inc.

So, the problem seems to be twofold. The QASM spec is wrong. Then, the implementation leads the API users to an illusion that the H gate is equivalent to gate h a { u2(0,pi) a; }. But, in fact, it does not use that definition at all.

By the way, I use this QASM paper as my reference.

consider

03:42 PM noon@ext>git clone [email protected]:vsoftco/qpp.git                                         △▲▲ 
Cloning into 'qpp'...
remote: Counting objects: 17774, done.
remote: Compressing objects: 100% (348/348), done.
remote: Total 17774 (delta 303), reused 0 (delta 0)
Receiving objects: 100% (17774/17774), 89.52 MiB | 88.00 KiB/s, done.
Resolving deltas: 100% (14175/14175), done.
Checking connectivity... done.
03:47 PM noon@ext>cd qpp                                                                           △△▲ 
03:47 PM noon@master qpp>ls                                                                        △△▲ 
CMakeLists.txt  COPYING  doc  examples  include  README.md  RELEASE.md  run_OSX_MATLAB  VERSION
03:47 PM noon@master qpp>mkdir build                                                               △△▲ 
03:48 PM noon@master qpp>cd build                                                                  △△▲ 
03:48 PM noon@master build>cmake ..                                                                △△▲ 
-- The C compiler identification is GNU 4.8.2
-- The CXX compiler identification is GNU 4.8.2
-- Check for working C compiler: /usr/bin/cc
-- Check for working C compiler: /usr/bin/cc -- works
-- Detecting C compiler ABI info
-- Detecting C compiler ABI info - done
-- Check for working CXX compiler: /usr/bin/c++
-- Check for working CXX compiler: /usr/bin/c++ -- works
-- Detecting CXX compiler ABI info
-- Detecting CXX compiler ABI info - done
-- Configuring done
-- Generating done
-- Build files have been written to: /home/noon/dev/ext/qpp/build
03:48 PM noon@master build>make                                                                    △△▲ 
Scanning dependencies of target qpp
[100%] Building CXX object CMakeFiles/qpp.dir/examples/example.cpp.o
In file included from /home/noon/dev/ext/qpp/examples/example.cpp:22:0:
/home/noon/dev/ext/qpp/include/qpp.h:63:23: fatal error: Eigen/Dense: No such file or directory
 #include <Eigen/Dense>
                       ^
compilation terminated.
make[2]: *** [CMakeFiles/qpp.dir/examples/example.cpp.o] Error 1
make[1]: *** [CMakeFiles/qpp.dir/all] Error 2
make: *** [all] Error 2

system information

03:51 PM noon@master qpp>uname -a                                                                             Tue Dec 16 △△▲ 
Linux qma 3.11.0-26-generic #45-Ubuntu SMP Tue Jul 15 04:02:06 UTC 2014 x86_64 x86_64 x86_64 GNU/Linux

Here is a test case:

OPENQASM 2.0;
include "qelib1.inc";

qreg q[2];
creg c[2];

h q;
measure q[1] -> c[0];
if(c==1) measure q[0] -> c[1];

The if statement does not take effect no matter what the c value is. In other words, the measure q[0] -> c[1]; is executed always.

The Open QASM spec seems vague on whether conditional measurements should be allowed. I can't see why not. The qpp does compile them, but just doesn't execute them correctly.