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

Merge pull request #139 from tnguyen-ornl/tnguyen/arithmethic-updates 

Update arithmethic standard library and add unit testing
parents dbc5c59f c49b8e7e

examples/CMakeLists.txt

+6 −1
Original line number Diff line number Diff line
@@ -20,7 +20,6 @@ add_test(NAME qrt_qaoa_example COMMAND ${CMAKE_BINARY_DIR}/qcor -c ${CMAKE_CURRE 
add_test(NAME qrt_simple-objective-function COMMAND ${CMAKE_BINARY_DIR}/qcor -c ${CMAKE_CURRENT_SOURCE_DIR}/simple/simple-objective-function.cpp)

add_test(NAME qrt_qpe_example COMMAND ${CMAKE_BINARY_DIR}/qcor -c ${CMAKE_CURRENT_SOURCE_DIR}/qpe/qpe_example_qrt.cpp)

add_test(NAME qrt_kernel_include COMMAND ${CMAKE_BINARY_DIR}/qcor -c ${CMAKE_CURRENT_SOURCE_DIR}/simple/multiple_kernels.cpp)

add_test(NAME qrt_period_finding COMMAND ${CMAKE_BINARY_DIR}/qcor -c ${CMAKE_CURRENT_SOURCE_DIR}/simple/period_finding.cpp)

add_test(NAME qrt_grover COMMAND ${CMAKE_BINARY_DIR}/qcor -c ${CMAKE_CURRENT_SOURCE_DIR}/grover/grover.cpp)

add_test(NAME qrt_adapt COMMAND ${CMAKE_BINARY_DIR}/qcor -c ${CMAKE_CURRENT_SOURCE_DIR}/hybrid/adapt_h2.cpp)

add_test(NAME qrt_ftqc_simple COMMAND ${CMAKE_BINARY_DIR}/qcor -c -qrt ftqc ${CMAKE_CURRENT_SOURCE_DIR}/ftqc_qrt/simple-demo.cpp)
@@ -42,6 +41,12 @@ add_test(NAME hadamard_ctrl_test COMMAND ${CMAKE_BINARY_DIR}/qcor ${CMAKE_CURREN 
add_test(NAME multi_ctrl_test COMMAND ${CMAKE_BINARY_DIR}/qcor ${CMAKE_CURRENT_SOURCE_DIR}/ctrl-gates/multiple_controls.cpp)



add_qcor_compile_and_exe_test(qrt_bell_ctrl bell/bell_control.cpp)



# Lambda tests

add_qcor_compile_and_exe_test(qrt_qpu_lambda_simple qpu_lambda/lambda_test.cpp)

add_qcor_compile_and_exe_test(qrt_qpu_lambda_bell qpu_lambda/lambda_test_bell.cpp)

add_qcor_compile_and_exe_test(qrt_qpu_lambda_grover qpu_lambda/grover_lambda_oracle.cpp)



# Arithmetic tests

add_qcor_compile_and_exe_test(qrt_qpu_arith_adder arithmetic/simple.cpp)

add_qcor_compile_and_exe_test(qrt_qpu_arith_integer_add arithmetic/integer_add.cpp)
 
 No newline at end of file

examples/arithmetic/integer_add.cpp

0 → 100644
+103 −0
Original line number Diff line number Diff line 
#include <qcor_arithmetic>



__qpu__ void test_add_integer(qreg q) {

  H(q[2]); // |000> + |001>

  add_integer(q, 3); // ==> |111> + |110>

  Measure(q);

}



__qpu__ void test_add_integer_mod(qreg q, qubit anc) {

  H(q[2]); // |000> + |001>

  add_integer_mod(q, anc, 3, 5); // |3 mod 5> and |7 mod 5>

  Measure(q);

  Measure(anc);

}



__qpu__ void test_mul_integer(qreg x, qreg b, qubit anc, int a, int N) {

  // b = 1;

  X(b[0]); 

  // x = |1> + |3>

  X(x[0]);

  H(x[1]);



  // ==> |a*x + b> 

  

  mul_integer_mod(x, b, anc, a, N);

  Measure(b);

  Measure(x);

}



__qpu__ void test_mul_integer_inline(qreg x, qreg anc, int a, int N) {

  // x = |1> + |3>

  X(x[0]);

  H(x[1]);

  // ==> |a*x> inline (save in x)

  // anc register is just for scratch pad.

  mul_integer_mod_in_place(x, anc.head(x.size() + 1), anc[x.size() + 1], a, N);

  Measure(x);

}



int main(int argc, char **argv) {

  set_shots(1024);

  auto a = qalloc(3);

  test_add_integer::print_kernel(a);

  test_add_integer(a);

  a.print();

  // Add 3 to a superposition of 0 and 4

  // => superposition of 3 and 7

  qcor_expect(a.counts().size() == 2);

  qcor_expect(a.counts()["111"] > 400);

  qcor_expect(a.counts()["110"] > 400);



  // Test modular add

  auto b = qalloc(3);

  auto anc = qalloc(1);

  test_add_integer_mod::print_kernel(b, anc[0]);

  test_add_integer_mod(b, anc[0]);

  // |3 mod 5> and |7 mod 5> == |2> + |3>

  b.print();

  qcor_expect(b.counts().size() == 2);

  qcor_expect(b.counts()["110"] > 400);

  qcor_expect(b.counts()["010"] > 400);

  anc.print();

  // anc returns to 0

  qcor_expect(anc.counts()["0"] == 1024);



  // Test modular multiply 

  int a_val = 3;

  // Large modulo => exact answer

  int N_val = 16;

  // More qubits to save the result

  auto b_reg = qalloc(5);

  auto x_reg = qalloc(2);

  auto anc_reg = qalloc(1);

  test_mul_integer::print_kernel(x_reg, b_reg, anc_reg[0], a_val, N_val);

  test_mul_integer(x_reg, b_reg, anc_reg[0], a_val, N_val);

  // x = |1> + |3>; |b> = 1

  // |a*x + b> = |4> + |10>

  b_reg.print();

  qcor_expect(b_reg.counts().size() == 2);

  // 4

  qcor_expect(b_reg.counts()["00100"] > 400);

  // 10 = 8 + 2

  qcor_expect(b_reg.counts()["01010"] > 400);

  x_reg.print();



  // Test in-place modular multiplication:

  // |x> ==> |ax mod N> on the same register.

  // x = |1> + |3>; a = 2

  // --> |2> + |6>

  // Simple test:

  a_val = 2;

  N_val = 8;

  auto x_reg2 = qalloc(3);

  auto anc_reg2 = qalloc(5);

  test_mul_integer_inline(x_reg2, anc_reg2, a_val, N_val);

  x_reg2.print();

  qcor_expect(x_reg2.counts().size() == 2);

  // 2

  qcor_expect(x_reg2.counts()["010"] > 400);

  // 6

  qcor_expect(x_reg2.counts()["011"] > 400);

  return 0;

}
 
 No newline at end of file

examples/arithmetic/simple.cpp

0 → 100644
+30 −0
Original line number Diff line number Diff line 
#include <qcor_arithmetic>



__qpu__ void test_adder(qreg a, qreg b, qubit cin, qubit cout) {

  integer_init(a, 1); // Set input a = 01

  integer_init(b, 3); // Set input b = 11

  // Apply the adder

  ripple_add(a, b, cin, cout);

  Measure(b);

  Measure(cout);

}



int main(int argc, char **argv) {

  set_shots(1024);

  // Set the inputs to the adder

  auto a = qalloc(2);

  auto b = qalloc(2);

  auto carry_in = qalloc(1);

  auto carry_out = qalloc(1);

  // Execute:

  test_adder::print_kernel(a, b, carry_in[0], carry_out[0]);

  test_adder(b, a, carry_in[0], carry_out[0]);

  b.print(); // 00

  qcor_expect(b.counts().size() == 1);

  qcor_expect(b.counts()["00"] == 1024);



  carry_out.print(); // 1

  qcor_expect(carry_out.counts().size() == 1);

  qcor_expect(carry_out.counts()["1"] == 1024);

  return 0;

}
 
 No newline at end of file

examples/simple/period_finding.cpp

deleted100644 → 0
+0 −14
Original line number Diff line number Diff line 
#include <qcor_arithmetic>



// 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;

}

lib/impl/arithmetic.hpp

deleted100644 → 0
+0 −306
Original line number Diff line number Diff line 
#pragma once



#include "qcor.hpp"

#include <qcor_qft>

// 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_range_opt_swap(q, startBitIdx, nbQubits, withSwap);

  int lastIdx = startBitIdx + nbQubits - 1;

  CX(q[lastIdx], q[aux]);

  qft_range_opt_swap(q, startBitIdx, nbQubits, withSwap);

  cPhiAdd(q, aux, N, startBitIdx, nbQubits, inv0);

  ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv1);

  iqft_range_opt_swap(q, startBitIdx, nbQubits, withSwap);

  X(q[lastIdx]);

  CX(q[lastIdx], q[aux]);

  X(q[lastIdx]);

  qft_range_opt_swap(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_range_opt_swap(q, startBitIdx, nbQubits, withSwap);

  int lastIdx = startBitIdx + nbQubits - 1;

  X(q[lastIdx]);

  CX(q[lastIdx], q[aux]);

  X(q[lastIdx]);

  qft_range_opt_swap(q, startBitIdx, nbQubits, withSwap);

  ccPhiAdd(q, ctl1, ctl2, a, startBitIdx, nbQubits, inv0);

  cPhiAdd(q, aux, N, startBitIdx, nbQubits, inv1);

  iqft_range_opt_swap(q, startBitIdx, nbQubits, withSwap);

  CX(q[lastIdx], q[aux]);

  qft_range_opt_swap(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) {

    SwapOp::ctrl(ctrlIdx, 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_range_opt_swap(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_range_opt_swap(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_range_opt_swap(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_range_opt_swap(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_range_opt_swap(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