Initial work on Period Finding

#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 function: wrapped it in a namespace to bypass XASM.

// These functions are used to construct circuit parameters.

namespace qcor { namespace util {

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

      if ((1 << j & 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;

}

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

}

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;

}

inline void calcNumBits(int& result, int N) {

  int count = 0;

  while (N) {

    count++;

    N = N >> 1;

  }

  result = count;

}

// a^(2^i)

inline void calcExpExp(int& result, int a, int i) {

  result = std::pow(result, std::pow(2, i));

}

}}

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

    // If theta is zero, just ignored.

    double eps = 1e-12;

    // FIXME: XASM handles if conditional

    int opCount = (abs(theta) < eps) ? 0 : 1;

    for (int ii = 0; ii < opCount; ++ii) {

      int idx = startBitIdx + i;

      // We don't have U1 op, hence use Rz instead

      Rz(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] ;

    // If theta is zero, just ignored.

    double eps = 1e-12;

    // FIXME: XASM handles if conditional

    int opCount = (abs(theta) < eps) ? 0 : 1;

    for (int ii = 0; ii < opCount; ++ii) {

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

    }

  }

}

// Doubly-controlled Add operation in Fourier space

__qpu__ void ccPhiAdd(qreg q, int ctrlBit1, int ctrlBit2, int a, int startBitIdx, int nbQubits, int inverse) {

  qcor::Controlled::Apply(ctrlBit2, cPhiAdd, q, ctrlBit1, a, startBitIdx, nbQubits, inverse);

}

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

}

__qpu__ void SwapOp(qreg q, int idx1, int idx2) {

  Swap(q[idx1], q[idx2]);

}

__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, and must be coprime with 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);

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

}

 No newline at end of file

#include "arithmetic.hpp"

// Compile:

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

  // Allocate 18 qubits required for a 4-bit number (4n*2)

  auto q = qalloc(18);

  int a = 4;

  int N = 15;

  // Call entry-point kernel

  periodFinding(q, a, N);

  q.print();

  return 0;

}

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

__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);

    }

  }

  // Swap qubits

  for (int qIdx = 0; qIdx < nQubits / 2; ++qIdx) {

    Swap(q[qIdx], q[nQubits - qIdx - 1]);

  // 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 - 1)/2; ++qIdx) {

      Swap(q[startIdx + qIdx], q[startIdx + nbQubits - qIdx - 1]);

    }

  }

}

// Inverse QFT

__qpu__ void iqft(qreg q, int maxBitIdx) {

  // Local Declarations

  const auto nQubits = maxBitIdx;

__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 < nQubits / 2; ++qIdx) {

    Swap(q[qIdx], q[nQubits - qIdx - 1]);

    for (int qIdx = 0; qIdx < (nbQubits - 1)/2; ++qIdx) {

      Swap(q[startIdx + qIdx], q[startIdx + nbQubits - 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);

  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[nQubits - 1]);

  H(q[startIdx + nbQubits - 1]);

}

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