Unverified Commit 753abc03 authored by Mccaskey, Alex's avatar Mccaskey, Alex Committed by GitHub
Browse files

Merge pull request #19 from tnguyen-ornl/tnguyen/devel

Added some arithmetic helper kernels
parents d6ba24a8 86372d2b
add_test(NAME qrt_bell_multi COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/bell_multi_qreg.cpp)
add_test(NAME qrt_bell_multi COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/bell/bell_multi_qreg.cpp)
find_library(STAQ_LIB xacc-staq-compiler HINTS ${CMAKE_INSTALL_PREFIX}/plugins)
if (STAQ_LIB)
add_test(NAME qrt_add_3_5 COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -v -c ${CMAKE_CURRENT_SOURCE_DIR}/add_3_5.cpp WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR})
add_test(NAME qrt_mixed_language COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/mixed_language.cpp)
add_test(NAME qrt_simple-demo COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/simple-demo.cpp)
add_test(NAME qrt_add_3_5 COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -v -c ${CMAKE_CURRENT_SOURCE_DIR}/adder/add_3_5.cpp WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/adder)
add_test(NAME qrt_mixed_language COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/simple/mixed_language.cpp)
add_test(NAME qrt_simple-demo COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/simple/simple-demo.cpp)
endif()
add_test(NAME qrt_deuteron_exp_inst COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/deuteron_exp_inst.cpp)
add_test(NAME qrt_deuteron_task_initiate COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/deuteron_task_initiate.cpp)
add_test(NAME qrt_qaoa_example COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/qaoa_example.cpp)
add_test(NAME qrt_simple-objective-function COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/simple-objective-function.cpp)
add_test(NAME qrt_qpe_example COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/qpe_example_qrt.cpp)
add_test(NAME qrt_kernel_include COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/multiple_kernels.cpp)
add_test(NAME qrt_deuteron_exp_inst COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/deuteron/deuteron_exp_inst.cpp)
add_test(NAME qrt_deuteron_task_initiate COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/deuteron/deuteron_task_initiate.cpp)
add_test(NAME qrt_qaoa_example COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/qaoa/qaoa_example.cpp)
add_test(NAME qrt_simple-objective-function COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c ${CMAKE_CURRENT_SOURCE_DIR}/simple/simple-objective-function.cpp)
add_test(NAME qrt_qpe_example COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c -I${CMAKE_CURRENT_SOURCE_DIR}/shared ${CMAKE_CURRENT_SOURCE_DIR}/qpe/qpe_example_qrt.cpp)
add_test(NAME qrt_kernel_include COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c -I${CMAKE_CURRENT_SOURCE_DIR}/shared ${CMAKE_CURRENT_SOURCE_DIR}/simple/multiple_kernels.cpp)
add_test(NAME qrt_period_finding COMMAND ${CMAKE_BINARY_DIR}/qcor -qrt -c -I${CMAKE_CURRENT_SOURCE_DIR}/shared ${CMAKE_CURRENT_SOURCE_DIR}/simple/period_finding.cpp)
// Example code structure whereby quantum kernels are defined
// in separate header files which can be included into a cpp file
// which is compiled by QCOR.
#pragma once
#include "qalloc.hpp"
// QFT kernel:
// Input: Qubit register and the max qubit index for the QFT,
// i.e. allow us to do QFT on a subset of the register [0, maxBitIdx)
__qpu__ void qft(qreg q, int maxBitIdx) {
// Local Declarations
const auto nQubits = maxBitIdx;
for (int qIdx = 0; qIdx < nQubits; ++qIdx) {
auto bitIdx = nQubits - qIdx - 1;
H(q[bitIdx]);
for (int j = 0; j < bitIdx; ++j) {
const double theta = M_PI/std::pow(2.0, bitIdx - j);
CPhase(q[j], q[bitIdx], theta);
}
}
// Swap qubits
for (int qIdx = 0; qIdx < nQubits / 2; ++qIdx) {
Swap(q[qIdx], q[nQubits - qIdx - 1]);
}
}
// Inverse QFT
__qpu__ void iqft(qreg q, int maxBitIdx) {
// Local Declarations
const auto nQubits = maxBitIdx;
// Swap qubits
for (int qIdx = 0; qIdx < nQubits / 2; ++qIdx) {
Swap(q[qIdx], q[nQubits - qIdx - 1]);
}
for (int qIdx = 0; qIdx < nQubits - 1; ++qIdx) {
H(q[qIdx]);
for (int j = 0; j < qIdx + 1; ++j) {
const double theta = -M_PI/std::pow(2.0, qIdx + 1 - j);
CPhase(q[j], q[qIdx + 1], theta);
}
}
H(q[nQubits - 1]);
}
\ No newline at end of file
......@@ -45,7 +45,9 @@ __qpu__ void QuantumPhaseEstimation(qreg q) {
}
// Inverse QFT on the counting qubits:
iqft(q, bitPrecision);
int startIdx = 0;
int shouldSwap = 1;
iqft(q, startIdx, bitPrecision, shouldSwap);
// Measure counting qubits
for (int qIdx = 0; qIdx < bitPrecision; ++qIdx) {
......
#pragma once
#include "qcor.hpp"
#include "qft.hpp"
// Note: We cannot put this in a namespace yet
// because we *do* want the XASM to recognize this as
// a quantum kernel hence invoking the token collector
// where we handle nested kernel calls (
// e.g. reset the __execute flag so that only the top-level
// kernel is submitted.
// Classical helper functions: wrapped it in a namespace to bypass XASM.
// These functions are used to construct circuit parameters.
namespace qcor { namespace util {
template <typename T>
T modpow(T base, T exp, T modulus) {
base %= modulus;
T result = 1;
while (exp > 0) {
if (exp & 1) result = (result * base) % modulus;
base = (base * base) % modulus;
exp >>= 1;
}
return result;
}
// Generates a list of angles to perform addition by a in the Fourier space.
void genAngles(std::vector<double>& io_angles, int a, int nbQubits) {
// Makes sure the vector appropriately sized
io_angles.resize(nbQubits);
std::fill(io_angles.begin(), io_angles.end(), 0.0);
for (int i = 0; i < nbQubits; ++i) {
int bitIdx = nbQubits - i - 1;
auto& angle = io_angles[bitIdx];
for (int j = i; j < nbQubits; ++j) {
// Check bit
int bitShift = nbQubits - j - 1;
if (((1 << bitShift) & a) != 0) {
angle += std::pow(2.0, -(j-i));
}
}
angle *= M_PI;
}
}
// Calculate 2^i * a % N;
inline void calcPowMultMod(int& result, int i, int a, int N) {
result = (1 << i) * a % N;
}
// Compute extended greatest common divisor of two integers
inline std::tuple<int, int, int> egcd(int a, int b) {
int m, n, c, q, r;
int m1 = 0, m2 = 1, n1 = 1, n2 = 0;
while (b) {
q = a / b, r = a - q * b;
m = m2 - q * m1, n = n2 - q * n1;
a = b, b = r;
m2 = m1, m1 = m, n2 = n1, n1 = n;
}
c = a, m = m2, n = n2;
// correct the signs
if (c < 0) {
m = -m;
n = -n;
c = -c;
}
return std::make_tuple(m, n, c);
}
// Modular inverse of a mod p
inline void modinv(int& result, int a, int p) {
int x, y;
int gcd_ap;
std::tie(x, y, gcd_ap) = egcd(p, a);
result = (y > 0) ? y : y + p;
}
// Compute the minimum number of bits required
// to represent an integer number N.
inline void calcNumBits(int& result, int N) {
int count = 0;
while (N) {
count++;
N = N >> 1;
}
result = count;
}
// Compute a^(2^i) mod N
inline void calcExpExp(int& result, int a, int i, int N) {
result = modpow(a, 1 << i, N);
}
}}
// Addition by a in Fourier Space
// If inverse != 0, run the reverse (minus)
// Note: the number to be added needs to be converted to an array of angles
// before calling this.
// Note: we supply the start Idx so that we can specify a subset of the register.
// Param: a -> the number to add
// startBitIdx and nbQubits specify a contiguous section of the qreg.
__qpu__ void phiAdd(qreg q, int a, int startBitIdx, int nbQubits, int inverse) {
std::vector<double> angles;
qcor::util::genAngles(angles, a, nbQubits);
for (int i = 0; i < nbQubits; ++i) {
double theta = (inverse == 0) ? angles[i] : -angles[i];
int idx = startBitIdx + i;
U1(q[idx], theta);
}
}
// Controlled add in Fourier Space
__qpu__ void cPhiAdd(qreg q, int ctrlBit, int a, int startBitIdx, int nbQubits, int inverse) {
std::vector<double> angles;
qcor::util::genAngles(angles, a, nbQubits);
for (int i = 0; i < nbQubits; ++i) {
double theta = (inverse == 0) ? angles[i] : -angles[i] ;
int idx = startBitIdx + i;
// Note: ctrlBit must be different from idx,
// i.e. ctrlBit must not be in the bitIdx vector
// We use CPhase which is equivalent with IBM cu1
CPhase(q[ctrlBit], q[idx], theta);
}
}
__qpu__ void ccPhase(qreg q, int ctl1, int ctl2, int tgt, double angle) {
CPhase(q[ctl1], q[tgt], angle/2);
CX(q[ctl2], q[ctl1]);
CPhase(q[ctl1], q[tgt], -angle/2);
CX(q[ctl2], q[ctl1]);
CPhase(q[ctl2], q[tgt], angle/2);
}
// Doubly-controlled Add operation in Fourier space
__qpu__ void ccPhiAdd(qreg q, int ctrlBit1, int ctrlBit2, int a, int startBitIdx, int nbQubits, int inverse) {
std::vector<double> angles;
qcor::util::genAngles(angles, a, nbQubits);
for (int i = 0; i < nbQubits; ++i) {
double theta = (inverse == 0) ? angles[i] : -angles[i] ;
int idx = startBitIdx + i;
ccPhase(q, ctrlBit1, ctrlBit2, idx, theta);
}
}
// Doubly-controlled *modular* addition by
// Inputs:
// - N the mod factor (represented as angles)
// - a the a addition operand (represented as angles)
__qpu__ void ccPhiAddModN(qreg q, int ctl1, int ctl2, int aux,
int a, int N,
int startBitIdx, int nbQubits) {
// FIXME: XASM should allow literal values in function calls
int inv0 = 0;
int inv1 = 1;
int withSwap = 0;
ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv0);
phiAdd(q, N, startBitIdx, nbQubits, inv1);
iqft(q, startBitIdx, nbQubits, withSwap);
int lastIdx = startBitIdx + nbQubits - 1;
CX(q[lastIdx], q[aux]);
qft(q, startBitIdx, nbQubits, withSwap);
cPhiAdd(q, aux, N, startBitIdx, nbQubits, inv0);
ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv1);
iqft(q, startBitIdx, nbQubits, withSwap);
X(q[lastIdx]);
CX(q[lastIdx], q[aux]);
X(q[lastIdx]);
qft(q, startBitIdx, nbQubits, withSwap);
ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv0);
}
// Inverse doubly-controlled Modular Addition (Fourier space)
__qpu__ void ccPhiAddModN_inv(qreg q, int ctl1, int ctl2, int aux,
int a, int N,
int startBitIdx, int nbQubits) {
// FIXME: XASM should allow literal values in function calls
int inv0 = 0;
int inv1 = 1;
int withSwap = 0;
ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv1);
iqft(q, startBitIdx, nbQubits, withSwap);
int lastIdx = startBitIdx + nbQubits - 1;
X(q[lastIdx]);
CX(q[lastIdx], q[aux]);
X(q[lastIdx]);
qft(q, startBitIdx, nbQubits, withSwap);
ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv0);
cPhiAdd(q, aux, N, startBitIdx, nbQubits, inv1);
iqft(q, startBitIdx, nbQubits, withSwap);
CX(q[lastIdx], q[aux]);
qft(q, startBitIdx, nbQubits, withSwap);
phiAdd(q, N, startBitIdx, nbQubits, inv0);
ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv1);
}
// Swap kernel: to be used to generate Controlled-Swap kernel
__qpu__ void SwapOp(qreg q, int idx1, int idx2) {
Swap(q[idx1], q[idx2]);
}
// Controlled-Swap
__qpu__ void cSwap(qreg q, int ctrlIdx, int idx1, int idx2) {
qcor::Controlled::Apply(ctrlIdx, SwapOp, q, idx1, idx2);
}
// Single controlled *modular* multiplication by a (mod N)
__qpu__ void cMultModN(qreg q, int ctrlIdx,
int a, int N,
int startBitIdx, int nbQubits,
// Auxilary qubit register
int auxStartBitIdx, int auxNbQubits) {
int inv0 = 0;
int inv1 = 1;
int withSwap = 0;
// QFT on the auxilary:
int qftSize = nbQubits + 1;
qft(q, auxStartBitIdx, qftSize, withSwap);
for (int i = 0; i < nbQubits; ++i) {
int tempA = 0;
qcor::util::calcPowMultMod(tempA, i, a, N);
int xIdx = startBitIdx + i;
int auxBit = auxStartBitIdx + auxNbQubits - 1;
// Doubly-controlled modular add on the Aux register.
ccPhiAddModN(q, xIdx, ctrlIdx, auxBit, tempA, N, auxStartBitIdx, qftSize);
}
iqft(q, auxStartBitIdx, qftSize, withSwap);
for (int i = 0; i < nbQubits; ++i) {
int idx1 = startBitIdx + i;
int idx2 = auxStartBitIdx + i;
cSwap(q, ctrlIdx, idx1, idx2);
}
qft(q, auxStartBitIdx, qftSize, withSwap);
int aInv = 0;
qcor::util::modinv(aInv, a, N);
for (int i = nbQubits - 1; i >= 0; --i) {
int tempA = 0;
qcor::util::calcPowMultMod(tempA, i, aInv, N);
int xIdx = startBitIdx + i;
int auxBit = auxStartBitIdx + auxNbQubits - 1;
ccPhiAddModN_inv(q, xIdx, ctrlIdx, auxBit, tempA, N, auxStartBitIdx, qftSize);
}
iqft(q, auxStartBitIdx, qftSize, withSwap);
}
// Period finding:
// Input a, N
// (1 < a < N)
// e.g. N = 15, a = 4
// Size of qreg must be at least (4n + 2)
// where n is the number of bits which can represent N,
// i.e. N = 15 -> n = 4 (4 bits) => qreg must have at lease 18 qubits.
__qpu__ void periodFinding(qreg q, int a, int N) {
// Down register: 0 -> (n-1): size = n
// Up register (QFT): n -> (3n - 1): size = 2n
// Aux register: 3n -> 4n + 1: size = n + 2
int n = 0;
qcor::util::calcNumBits(n, N);
int downStart = 0;
int downSize = n;
int upStart = n;
int upSize = 2*n;
int auxStart = 3*n;
int auxSize = n + 2;
// Hadamard up register:
for (int i = 0; i < upSize; ++i) {
int bitIdx = upStart + i;
H(q[bitIdx]);
}
// Init down register to 1
X(q[downStart]);
// Apply the multiplication gates
for (int i = 0; i < 2*n; ++i) {
int aTo2Toi = 0;
qcor::util::calcExpExp(aTo2Toi, a, i, N);
// Control bit is the upper register
int ctrlIdx = upStart + i;
cMultModN(q, ctrlIdx, aTo2Toi, N, downStart, downSize, auxStart, auxSize);
}
// Apply inverse QFT on the up register (with Swap)
int withSwap = 1;
iqft(q, upStart, upSize, withSwap);
// Measure the upper register
for (int i = 0; i < upSize; ++i) {
int bitIdx = upStart + i;
Measure(q[bitIdx]);
}
// Done: examine the shot count statistics
// to compute the period.
}
\ No newline at end of file
// Example code structure whereby quantum kernels are defined
// in separate header files which can be included into a cpp file
// which is compiled by QCOR.
#pragma once
#include "qalloc.hpp"
// Generic QFT and IQFT algorithms
// Input:
// (1) Qubit register (type qreg)
// (2) The start qubit index and the number of qubits (type integer)
// Note: these two parameters allow us to operate the QFT/IQFT on
// a contiguous subset of qubits of the register.
// If we want to use the entire register,
// just pass 0 and number of qubits in the register, respectively.
// (3) shouldSwap flag (as integer):
// If 0, no Swap gates will be added.
// Otherwise, Swap gates are added at the end (QFT) and beginning (IQFT).
__qpu__ void qft(qreg q, int startIdx, int nbQubits, int shouldSwap) {
for (int qIdx = nbQubits - 1; qIdx >= 0; --qIdx) {
auto shiftedBitIdx = qIdx + startIdx;
H(q[shiftedBitIdx]);
for (int j = qIdx - 1; j >= 0; --j) {
const double theta = M_PI/std::pow(2.0, qIdx - j);
auto targetIdx = j + startIdx;
CPhase(q[shiftedBitIdx], q[targetIdx], theta);
}
}
// A *hacky* way to do conditional (convert to a for loop)
int swapCount = (shouldSwap == 0) ? 0 : 1;
for (int count = 0; count < swapCount; ++count) {
for (int qIdx = 0; qIdx < nbQubits/2; ++qIdx) {
Swap(q[startIdx + qIdx], q[startIdx + nbQubits - qIdx - 1]);
}
}
}
__qpu__ void iqft(qreg q, int startIdx, int nbQubits, int shouldSwap) {
int swapCount = (shouldSwap == 0) ? 0 : 1;
for (int count = 0; count < swapCount; ++count) {
// Swap qubits
for (int qIdx = 0; qIdx < nbQubits/2; ++qIdx) {
Swap(q[startIdx + qIdx], q[startIdx + nbQubits - qIdx - 1]);
}
}
for (int qIdx = 0; qIdx < nbQubits - 1; ++qIdx) {
H(q[startIdx + qIdx]);
int j = qIdx + 1;
for (int y = qIdx; y >= 0; --y) {
const double theta = -M_PI/std::pow(2.0, j - y);
CPhase(q[startIdx + j], q[startIdx + y], theta);
}
}
H(q[startIdx + nbQubits - 1]);
}
\ No newline at end of file
......@@ -8,10 +8,12 @@ __qpu__ void f(qreg q) {
X(q[1]);
// Call qft kernel (defined in a separate header file)
qft(q, nQubits);
int startIdx = 0;
int shouldSwap = 1;
qft(q, startIdx, nQubits, shouldSwap);
// Inverse QFT:
iqft(q, nQubits);
iqft(q, startIdx, nQubits, shouldSwap);
// Measure all qubits
for (int qIdx = 0; qIdx < nQubits; ++qIdx) {
......
#include "arithmetic.hpp"
// Compile:
int main(int argc, char **argv) {
// Allocate 22 qubits required for a 5-bit number (4n*2)
auto q = qalloc(22);
int a = 11;
int N = 21;
// Call entry-point kernel
periodFinding(q, a, N);
q.print();
return 0;
}
......@@ -102,6 +102,10 @@ void rz(const qubit &qidx, const double theta) {
one_qubit_inst("Rz", qidx, {theta});
}
void u1(const qubit &qidx, const double theta) {
one_qubit_inst("U1", qidx, {theta});
}
void mz(const qubit &qidx) { one_qubit_inst("Measure", qidx); }
void cnot(const qubit &src_idx, const qubit &tgt_idx) {
......
......@@ -40,6 +40,8 @@ void tdg(const qubit& qidx);
void rx(const qubit &qidx, const double theta);
void ry(const qubit &qidx, const double theta);
void rz(const qubit &qidx, const double theta);
// U1(theta) gate
void u1(const qubit &qidx, const double theta);
void mz(const qubit &qidx);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment