//--- Scaffold Code for Shor’s Algorithm ---//
//Input x -> H -> f(x) -> iQFT
//f(x) = a^x mod N
// Optimized implementation as described in
// http://arxiv.org/pdf/1207.0511v5.pdf
// Quantum Information and Computation, Vol. 14, No. 7&8 (2014) 0649–0682
#include <math.h>
#include <stdlib.h>
#define _n 512 //number of bits in N
#include "look_up_tables.h"
#define pi 3.141592
#define _N 10941738641570527421809707322040357612003732945449205990913842131476349984288934784717997257891267332497625752899781833797076537244027146743531593354333897 //number to be factorized
//char Nbits [] = "11010000111010100001101010111010100101111000110111110000000001100101101100100000000010011111011101011100100001000110111100101000101100000100111011010101000101000011101100100011011110110011111111000010010000100111001000100100010110101101111010000011011111101111111000000010011100011110000110100010100001010100111000001100100000011011101010011111011100001010100000111010110110000110110101000111101100001110101011001101000001100010101111000001010110111100011000011010100110011101110010000011000100100100111011001001"
//char Nbits [] = {1,1,0,1,0,0,0,0,1,1,1,0,1,0,1,0,0,0,0,1,1,0,1,0,1,0,1,1,1,0,1,0,1,0,0,1,0,1,1,1,1,0,0,0,1,1,0,1,1,1,1,1,0,0,0,0,0,0,0,0,0,1,1,0,0,1,0,1,1,0,1,1,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,1,1,1,1,0,1,1,1,0,1,0,1,1,1,0,0,1,0,0,0,0,1,0,0,0,1,1,0,1,1,1,1,0,0,1,0,1,0,0,0,1,0,1,1,0,0,0,0,0,1,0,0,1,1,1,0,1,1,0,1,0,1,0,1,0,0,0,1,0,1,0,0,0,0,1,1,1,0,1,1,0,0,1,0,0,0,1,1,0,1,1,1,1,0,1,1,0,0,1,1,1,1,1,1,1,1,0,0,0,0,1,0,0,1,0,0,0,0,1,0,0,1,1,1,0,0,1,0,0,0,1,0,0,1,0,0,0,1,0,1,1,0,1,0,1,1,0,1,1,1,1,0,1,0,0,0,0,0,1,1,0,1,1,1,1,1,1,0,1,1,1,1,1,1,1,0,0,0,0,0,0,0,1,0,0,1,1,1,0,0,0,1,1,1,1,0,0,0,0,1,1,0,1,0,0,0,1,0,1,0,0,0,0,1,0,1,0,1,0,0,1,1,1,0,0,0,0,0,1,1,0,0,1,0,0,0,0,0,0,1,1,0,1,1,1,0,1,0,1,0,0,1,1,1,1,1,0,1,1,1,0,0,0,0,1,0,1,0,1,0,0,0,0,0,1,1,1,0,1,0,1,1,0,1,1,0,0,0,0,1,1,0,1,1,0,1,0,1,0,0,0,1,1,1,1,0,1,1,0,0,0,0,1,1,1,0,1,0,1,0,1,1,0,0,1,1,0,1,0,0,0,0,0,1,1,0,0,0,1,0,1,0,1,1,1,1,0,0,0,0,0,1,0,1,0,1,1,0,1,1,1,1,0,0,0,1,1,0,0,0,0,1,1,0,1,0,1,0,0,1,1,0,0,1,1,1,0,1,1,1,0,0,1,0,0,0,0,0,1,1,0,0,0,1,0,0,1,0,0,1,0,0,1,1,1,0,1,1,0,0,1,0,0,1};
//_dc = _N - (1<<_n)
//char dcbits[] = {101111000101011110010101000101011010000111001000001111111110011010010011011111111101100000100010100011011110111001000011010111010011111011000100101010111010111100010011011100100001001100000000111101101111011000110111011011101001010010000101111100100000010000000111111101100011100001111001011101011110101011000111110011011111100100010101100000100011110101011111000101001001111001001010111000010011110001010100110010111110011101010000111110101001000011100111100101011001100010001101111100111011011011000100110111};
//FIXME: negative this
#define _a 2 //randomly chosen s.t. gcd(a,N) = 1
// pre-computed constants for the GMQDIV2 circuit
// we use the fact that this divider is only ever used
// to divide by _N
#define _l _n //1+floor(log2(_N))
#define _mdash floor(((1<<(_n+1))-1)/_N) - (1<<_N) //FIXME: definition not clear from paper
//AQFT parameters
#define _n_upper 2*_n
#define log_2(x) log(x)/log(2)
#define AQFTk (int)(3+log_2(_n_upper)) //value of k for approximate QFT // k>=2+log2(n)
module swap (qbit a, qbit b) {
CNOT(a, b);
CNOT(b, a);
CNOT(a, b);
}
module cRz (qbit ctrl, qbit target, double angle) {
Rz (target, angle/2);
CNOT (target, ctrl);
Rz (target, -1*angle/2);
CNOT (target, ctrl);
}
module ccRz (qbit ctrl1, qbit ctrl2, qbit target, double angle)
{
// doubly - controlled Rz circuit with angle
// decompose in terms of V = sqrt(U)
cRz (ctrl2, target, angle/2);
CNOT (ctrl2, ctrl1);
cRz (ctrl1, target, -angle/2);
CNOT (ctrl2, ctrl1);
cRz (ctrl1, target, angle/2);
}
module QFT(qbit x[]) {
int i, j;
double angle = 0.0;
for (i=0; i<_n; i++){
H(x[i]);
for (j=i+1; j<_n; j++) {
angle = angle_R[j];
cRz(x[j], x[i], angle);
}
}
// TODO: swaps for reordering
}
// inverse QFT
module iQFT (qbit x[]){
int i, j;
double angle = 0.0;
for (i=_n-1; i>=0; i--){
for (j=_n-1; j>i; j--) {
angle = -angle_R[j];
cRz(x[j], x[i], angle);
}
H(x[i]);
}
// TODO: swaps for reordering
}
// QFT with input broken into two equal parts, each size n
module QFT_2 (qbit x0[], qbit x1[]) {
int i, j;
double angle = 0.0;
// the first register
for (i=0; i<_n; i++){
H(x0[i]);
for (j=i+1; j<_n; j++) {
angle = angle_R[j];
cRz(x0[j], x0[i], angle);
}
for (j=0; j<_n; j++) {
angle = angle_R[j+_n];
cRz(x1[j], x0[i], angle);
}
}
// the second register
for (i=0; i<_n; i++){
H(x1[i]);
for (j=i+1; j<_n; j++) {
angle = angle_R[j+_n];
cRz(x1[j], x1[i], angle);
}
}
// TODO: swaps for reordering
}
// inverse QFT with input broken into two equal parts, each size n
module iQFT_2 (qbit x0[], qbit x1[]){
int i, j;
double angle = 0.0;
// the second register
for (i=_n-1; i>=0; i--){
for (j=_n-1; j>i; j--) {
angle = -angle_R[j+_n];
cRz(x1[j], x1[i], angle);
}
H(x1[i]);
}
// the first register
for (i=_n-1; i>=0; i--){
for (j=_n-1; j>=0; j--) {
angle = -angle_R[j+_n];
cRz(x1[j], x0[i], angle);
}
for (j=_n-1; j>i; j--) {
angle = -angle_R[j];
cRz(x0[j], x0[i], angle);
}
H(x0[i]);
}
// TODO: swaps for reordering
}
// add constant _N to qubit register which is broken down into 2 registers of length _n
module QADD_2 (qbit R0[], qbit R2[])
{
int j = 0;
double angle = 0.0;
for(j=0; j<_n; j++){
angle = angle_A_N[j];
Rz(R0[j],angle);
}
for(j=0; j<_n; j++){
angle = angle_A_N[j+_n];
Rz(R2[j],angle);
}
}
// add constant _N to qubit register which is broken down into 2 registers of length _n
module iQADD_2 (qbit R0[], qbit R2[])
{
int j = 0;
double angle = 0.0;
for(j=_n-1; j>=0; j--){
angle = -angle_A_N[j+_n];
Rz(R2[j],angle);
}
for(j=_n-1; j>=0; j--){
angle = -angle_A_N[j];
Rz(R0[j],angle);
}
}
module CQADD_1 (qbit ctrl, qbit target[])
{
int j = 0;
double angle = 0.0;
for(j=0; j<_n; j++){
angle = angle_A_1[j]; //compute angle of rotation
cRz(ctrl,target[j],angle);
}
}
module iCQADD_1 (qbit ctrl, qbit target[])
{
int j = 0;
double angle = 0.0;
for(j=_n-1; j>=0; j--){
angle = -angle_A_1[j]; //compute angle of rotation
cRz(ctrl,target[j],angle);
}
}
module CQADD_dc (qbit ctrl, qbit target[])
{
int j = 0;
double angle = 0.0;
for(j=0; j<_n; j++){
angle = angle_A_dc[j]; //compute angle of rotation
cRz(ctrl,target[j],angle);
}
}
module iCQADD_dc (qbit ctrl, qbit target[])
{
int j = 0;
double angle = 0.0;
for(j=_n-1; j>=0; j--){
angle = -angle_A_dc[j]; //compute angle of rotation
cRz(ctrl,target[j],angle);
}
}
// add two n-qubit registers (Fig.7 in paper)
// qbit a[n]
// qbit b[n]
// (left is lsb)
module QADDn( qbit a[], // n bits
qbit b[] // n bits
)
{
int i = 0;
int j = 0;
double angle = 0.0;
// is this going to be parallelized?
// lower bits
for (i=_n-1; i>=0; i--) {
for (j=i; j<_n; j++) {
angle = angle_R [j-i+1];
cRz (a[i], b[j], angle);
}
}
}
module iQADDn( qbit a[], // n bits
qbit b[] // n bits
)
{
int i = 0;
int j = 0;
double angle = 0.0;
// is this going to be parallelized?
// lower bits
for (i=0; i<_n; i++) {
for (j=_n-1; j>=i; j--) {
angle = -angle_R [j-i+1];
cRz (a[i], b[j], angle);
}
}
}
// multiplier-accumulator (refer to Fig.11b and Fig.12 in paper)
// parallel version with 8n depth
// Is this parallelized?
module QMAC_a(qbit ctrl,
qbit x[], // size n
// b is divided into two registers of size n
qbit R10[], // size n
qbit R20[] // size n
)
{
double angle = 0.0;
int i=0;
int j=0;
int l=0;
// V section
for (l=0; l<_n; l++){
for (i=0; i<_n; i++){
j = (l-i)%_n;
angle = angle_V_a[l][j];
cRz(x[j], R10[l], angle);
}
}
for (l=0; l<_n; l++){
for (i=0; i<_n; i++){
j = (l+_n-i)%_n;
angle = angle_V_a[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
//CNOTs
for(i=0; i<_n; i++)
CNOT(x[i], ctrl);
// V_dagger section
for(l=0; l<_n; l++){
for(i=0; i<_n; i++){
j = (l-i)%_n;
angle = -angle_V_a[l][j];
cRz(x[j], R10[l], angle);
}
}
for(l=0; l<_n; l++){
for(i=0; i<_n; i++){
j = (l+_n-i)%_n;
angle = -angle_V_a[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
//CNOTs
for(i=0; i<_n; i++)
CNOT(x[i], ctrl);
// W section
for(l=0; l<_n; l++) {
angle = angle_W_a[l];
cRz(ctrl, R10[l], angle);
}
for(l=0; l<_n; l++) {
angle = angle_W_a[l+_n];
cRz(ctrl, R20[l], angle);
}
}
module iQMAC_a(qbit ctrl,
qbit x[], // size n
// b is divided into two registers of size n
// which are again divided into two of size t and n-t
qbit R10[], // size n
qbit R20[] // size n
)
{
double angle = 0.0;
int i=0;
int j=0;
int l=0;
// More parallel version of above with O(n) depth
// Is this parallelized?
// W section
for(l=_n-1; l>=0; l--) {
angle = -angle_W_a[l+_n];
cRz(ctrl, R20[l], angle);
}
for(l=_n-1; l>=0; l--) {
angle = -angle_W_a[l];
cRz(ctrl, R10[l], angle);
}
//CNOTs
for(i=_n-1; i>=0; i--)
CNOT(x[i], ctrl);
// V_dagger section
for(l=_n-1; l>=0; l--){
for(i=_n-1; i>=0; i--){
j = (l+_n-i)%_n;
angle = angle_V_a[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
for(l=_n-1; l>=0; l--){
for(i=_n-1; i>=0; i--){
j = (l-i)%_n;
angle = angle_V_a[l][j];
cRz(x[j], R10[l], angle);
}
}
//CNOTs
for(i=_n-1; i>=0; i--)
CNOT(x[i], ctrl);
// V section
for (l=_n-1; l>=0; l--){
for (i=_n-1; i>=0; i--){
j = (l+_n-i)%_n;
angle = -angle_V_a[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
for (l=_n-1; l>=0; l--){
for (i=_n-1; i>=0; i--){
j = (l-i)%_n;
angle = -angle_V_a[l][j];
cRz(x[j], R10[l], angle);
}
}
}
// multiplier-accumulator (refer to Fig.11b and Fig.12 in paper)
// parallel version with 8n depth
// Is this parallelized?
module QMAC_N(qbit ctrl,
qbit x[], // size n
// b is divided into two registers of size n
qbit R10[], // size n
qbit R20[] // size n
)
{
double angle = 0.0;
int i=0;
int j=0;
int l=0;
// V section
for (l=0; l<_n; l++){
for (i=0; i<_n; i++){
j = (l-i)%_n;
angle = angle_V_N[l][j];
cRz(x[j], R10[l], angle);
}
}
for (l=0; l<_n; l++){
for (i=0; i<_n; i++){
j = (l+_n-i)%_n;
angle = angle_V_N[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
//CNOTs
for(i=0; i<_n; i++)
CNOT(x[i], ctrl);
// V_dagger section
for(l=0; l<_n; l++){
for(i=0; i<_n; i++){
j = (l-i)%_n;
angle = -angle_V_N[l][j];
cRz(x[j], R10[l], angle);
}
}
for(l=0; l<_n; l++){
for(i=0; i<_n; i++){
j = (l+_n-i)%_n;
angle = -angle_V_N[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
//CNOTs
for(i=0; i<_n; i++)
CNOT(x[i], ctrl);
// W section
for(l=0; l<_n; l++) {
angle = angle_W_N[l];
cRz(ctrl, R10[l], angle);
}
for(l=0; l<_n; l++) {
angle = angle_W_N[l+_n];
cRz(ctrl, R20[l], angle);
}
}
module iQMAC_N(qbit ctrl,
qbit x[], // size n
// b is divided into two registers of size n
// which are again divided into two of size t and n-t
qbit R10[], // size n
qbit R20[] // size n
)
{
double angle = 0.0;
int i=0;
int j=0;
int l=0;
// More parallel version of above with O(n) depth
// Is this parallelized?
// W section
for(l=_n-1; l>=0; l--) {
angle = -angle_W_N[l+_n];
cRz(ctrl, R20[l], angle);
}
for(l=_n-1; l>=0; l--) {
angle = -angle_W_N[l];
cRz(ctrl, R10[l], angle);
}
//CNOTs
for(i=_n-1; i>=0; i--)
CNOT(x[i], ctrl);
// V_dagger section
for(l=_n-1; l>=0; l--){
for(i=_n-1; i>=0; i--){
j = (l+_n-i)%_n;
angle = angle_V_N[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
for(l=_n-1; l>=0; l--){
for(i=_n-1; i>=0; i--){
j = (l-i)%_n;
angle = angle_V_N[l][j];
cRz(x[j], R10[l], angle);
}
}
//CNOTs
for(i=_n-1; i>=0; i--)
CNOT(x[i], ctrl);
// V section
for (l=_n-1; l>=0; l--){
for (i=_n-1; i>=0; i--){
j = (l+_n-i)%_n;
angle = -angle_V_N[l+_n][j];
cRz(x[j], R20[l], angle);
}
}
for (l=_n-1; l>=0; l--){
for (i=_n-1; i>=0; i--){
j = (l-i)%_n;
angle = -angle_V_N[l][j];
cRz(x[j], R10[l], angle);
}
}
}
// Granlund-Montgomery division by constant circuit, giving n-bit quotient and n-bit remainder
// Circuit Depth = 244n timesteps
// Circuit Width = 7n+1 qubits
// Figures 13 & 14 & 18
// In this case, divisor d = _N always.
module GMQDIV2(
// dividend z at the start of computation (2*_n qubits)
qbit Reg00[_n], // lower bits: l bits
qbit Reg10[_n] // lower bits: l bits
)
{
// quotient q at the end of computation (n qubits)
qbit Reg2[_n]; // n bits
// ancilla qubits
qbit Reg3[_n]; // ancilla -- q1
qbit Reg4[_n]; // ancilla -- n2, n2+n1, n2, 0, n1, 0
qbit Reg51[_n];
qbit Reg6[_n]; // ancilla -- LOW(m'(n2+n1)+nadj)
qbit Aqbit[1];// ancilla - sign(n1)
qbit anc[1];
// loop variable
int i;
// -------------------- Compute (fig. 13) ------------------------/
// SLL(HIGH(z),n-l) --> Reg0[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg4[0 .. n-1]
QADDn(Reg00, Reg4);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
QADDn(Reg10, Reg6);
// SRL(LOW(z), l) --> Reg1[l .. n-1] | Reg5[n-l .. n-1]
// Reg4[0 .. n-1]
QADDn(Reg51, Reg4);
iQFT(Reg6);
// CNOT(target, control)
CNOT(Aqbit[0], Reg6[_n-1]);
CQADD_1(Aqbit[0], Reg4);
QFT(Reg6);
CQADD_dc(Aqbit[0], Reg6);
iQFT(Reg6);
// Reg6[n]|Reg51[l]
QFT_2(Reg6, Reg51);
iQFT(Reg4);
//FIXME: code me (algorithm does not show control signal)
//QMAC_mdash();
iQFT_2(Reg6, Reg51);
for (i=0; i<_n ; i++) {
CNOT(Reg2[_n-_n+i], Reg51[i]);
}
QFT(Reg2);
QADDn(Reg2, Reg4);
iCQADD_1(Aqbit[0], Reg2);
iQFT(Reg2);
QFT_2(Reg10, Reg00);
for (i=0; i<_n; i++)
CNOT(Reg3[i], Reg2[i]);
// FIXME: code me
//QMAC(, -_N);
iQADD_2(Reg10, Reg00);
QFT(Reg2);
iQFT_2(Reg10, Reg00);
iQFT (Reg2);
QFT_2(Reg10, Reg00);
QADD_2(Reg10, Reg00);
// FIXME: incomplete document specification -- create |1> ctrl signal
X(anc[0]);
//QMAC_N(anc[0], Reg2, Reg10, Reg00);
iQFT_2(Reg10, Reg00);
//------------------------- UnCompute (fig. 14) ------------------------/
QFT(Reg3);
iQADDn(Reg3, Reg4);
CQADD_1(Aqbit[0], Reg3);
iQFT(Reg3);
for (i=0; i<_n; i++) {
CNOT(Reg3[i], Reg51[i]);
}
QFT_2(Reg6, Reg51);
//FIXME: code me
//QMAC_mdash();
iQFT_2(Reg6, Reg51);
QFT(Reg4);
iCQADD_1(Aqbit[0], Reg4);
QFT(Reg6);
// SLL(HIGH(z),n-l) --> Reg0[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg4[0 .. n-1]
iQADDn(Reg00, Reg4);
iCQADD_1(Aqbit[0], Reg6);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
QADDn(Reg10, Reg6);
iQFT(Reg6);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
iQADDn(Reg10, Reg4);
iQFT(Reg4);
QFT(Reg4);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
iQADDn(Reg10, Reg4);
iQFT(Reg4);
QFT_2(Reg10, Reg00);
// FIXME: code me
//QMAC(, -_N);
iQFT_2(Reg10, Reg00);
}
module iGMQDIV2(
// dividend z at the start of computation (2*_n qubits)
qbit Reg00[_n], // lower bits: l bits
qbit Reg10[_n] // lower bits: l bits
)
{
// quotient q at the end of computation (n qubits)
qbit Reg2[_n]; // n bits
// ancilla qubits
qbit Reg3[_n]; // ancilla -- q1
qbit Reg4[_n]; // ancilla -- n2, n2+n1, n2, 0, n1, 0
qbit Reg51[_n];
qbit Reg6[_n]; // ancilla -- LOW(m'(n2+n1)+nadj)
qbit Aqbit[1]; // ancilla - sign(n1)
qbit anc[1];
// loop variable
int i;
//------------------------- UnCompute (fig. 14) ------------------------/
QFT_2(Reg10, Reg00);
// FIXME: code me
//iQMAC(, -_N);
iQFT_2(Reg10, Reg00);
QFT(Reg4);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
QADDn(Reg10, Reg4);
iQFT(Reg4);
QFT(Reg4);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
QADDn(Reg10, Reg4);
QFT(Reg6);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
iQADDn(Reg10, Reg6);
CQADD_dc(Aqbit[0], Reg6);
// SLL(HIGH(z),n-l) --> Reg0[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg4[0 .. n-1]
QADDn(Reg00, Reg4);
iQFT(Reg6);
CQADD_1(Aqbit[0], Reg4);
iQFT(Reg4);
QFT_2(Reg6, Reg51);
//FIXME: code me
//iQMAC_mdash();
iQFT_2(Reg6, Reg51);
for (i=_n-1; i>=0; i--) {
CNOT(Reg3[i], Reg51[i]);
}
QFT(Reg3);
iCQADD_1(Aqbit[0], Reg3);
iQADDn(Reg3, Reg4);
iQFT(Reg3);
// -------------------- Compute (fig. 13) ------------------------/
QFT_2(Reg10, Reg00);
// FIXME: incomplete document specification
X(anc[0]);
//iQMAC_N(anc[0], Reg2, Reg10, Reg00);
iQADD_2(Reg10, Reg00);
iQFT_2(Reg10, Reg00);
QFT (Reg2);
QFT_2(Reg10, Reg00);
iQFT(Reg2);
QADD_2(Reg10, Reg00);
// FIXME: code me
//iQMAC(, -_N);
for (i=_n-1; i>=0; i--)
CNOT(Reg3[i], Reg2[i]);
iQFT_2(Reg10, Reg00);
QFT(Reg2);
CQADD_1(Aqbit[0], Reg2);
iQADDn(Reg2, Reg4);
iQFT(Reg2);
for (i=_n-1; i>=0 ; i--) {
CNOT(Reg2[_n-_n+i], Reg51[i]);
}
QFT_2(Reg6, Reg51);
//FIXME: code me (algorithm does not show control signal)
//iQMAC_mdash();
QFT(Reg4);
// Reg6[n]|Reg51[l]
iQFT_2(Reg6, Reg51);
QFT(Reg6);
iCQADD_dc(Aqbit[0], Reg6);
iQFT(Reg6);
iCQADD_1(Aqbit[0], Reg4);
// CNOT(target, control)
CNOT(Aqbit[0], Reg6[_n-1]);
QFT(Reg6);
// SRL(LOW(z), l) --> Reg1[l .. n-1] | Reg5[n-l .. n-1]
// Reg4[0 .. n-1]
iQADDn(Reg51, Reg4);
// SLL(LOW(z),n-l) --> Reg1[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg6[0 .. n-1]
iQADDn(Reg10, Reg6);
// SLL(HIGH(z),n-l) --> Reg0[0 .. l-1] | Reg5[0 .. n-l-1]
// Reg4[0 .. n-1]
iQADDn(Reg00, Reg4);
}
// Optimized controlled modular multiplier/accumulator (Fig. 19)
// 4n+1 qubits + 5n+1 ancillas
module QMAC_MOD2 (qbit ctrl,
qbit y[_n],
qbit R10[_n],
qbit R20[_n],
qbit R30[_n]
)
{
int i = 0;
QFT_2(R20, R10);
QMAC_a(ctrl, y, R10, R20);
iQFT_2(R20, R10);
GMQDIV2(R20, R10); // only GMQDIV2 with parameter _N is used
for (i=0; i<_n; i++)
CNOT(R30[i], R10[i]);
iGMQDIV2(R20, R10); // only iGMQDIV2 with parameter _N is used
QFT_2(R20, R30);
iQMAC_a(ctrl, y, R20, R30);
iQFT_2(R20, R30);
}
module iQMAC_MOD2 (qbit ctrl,
qbit y[_n],
qbit R10[_n],
qbit R20[_n],
qbit R30[_n])
{
int i = 0;
QFT_2(R20, R30);
QMAC_a(ctrl, y, R20, R30);
iQFT_2(R20, R30);
GMQDIV2(R20, R10);
for (i=_n-1; i>=0; i--)
CNOT(R30[i], R10[i]);
iGMQDIV2(R20, R10);
QFT_2(R20, R10);
iQMAC_a(ctrl, y, R10, R20);
iQFT_2(R20, R10);
}
// Optimized controlled modular multiplier (Fig. 20) -- specialized a=_a and N=_N
// 4n+1 qubits + 5n+1 ancillas
// In total: depth = 1045n-38
module CU(qbit ctrl, qbit x[_n], int k)
{
//Compute U^(2^k) = a^(2^k) mod N
//Uf = a^x mod N.
//When x is expressed in binary as x[n-1]:x[n-2]:...:x[1]:x[0],
//a^x mod N = PRODUCT{a^(2^(k)*x[k]}, where k goes from 0 to (n-1)
// iteration indices
int i, j;
// ancilla qubits
qbit R10[_n];
qbit R20[_n];
qbit R30[_n];
for (i=0; i<_n; i++){
PrepZ(R10[i], 0);
PrepZ(R20[i], 0);
PrepZ(R30[i], 0);
}
// QMAC_MOD2 is only used in the overall design with (_a, _N) as parameters
QMAC_MOD2(ctrl, x, R10, R20, R30);
for (i=0; i<_n; i++) {
//Fredkin(x[i], R10[i], ctrl);
// Fredkin decomposition (J. A. Smolin and D. P. DiVincenzo, Phys. Rev. A, 53, 2855 (1996).)
CNOT(R10[i], ctrl);
Toffoli(x[i], R10[i], ctrl);
CNOT(R10[i], ctrl);
}
iQMAC_MOD2(ctrl, x, R10, R20, R30);
}
module RotBy(qbit q, int k, int neg)
{
//neg = +1 for positive angle
//neg = -1 for negative angle
//Rk = [1 0; 0 e^i(2*pi/2^k) ]
//Rk(theta) performed using Rz(theta/2), and ignoring global phase shift
//angle is negative for inverse QFT
double angle = neg*pi/pow(2,k); //theta/2
Rz(q, angle);
}
module CRinv(qbit ctrl, qbit q, int k) //Controlled RotBy
{
RotBy(q,k-1, -1); //SP: negative angle
CNOT(q, ctrl); // ajavadia: changed order of ctrl and target
RotBy(q,k-1, 1); //positive angle
CNOT(q, ctrl);
}
module AQFT(qbit x[_n_upper], qbit y[_n])
{
cbit measBits[_n_upper];
int i = 0, j = 0;
for(i=0; i<_n_upper; i++)
H(x[i]);
for(i=0; i<_n_upper; i++)
CU(x[i], y, i);
for(i=_n_upper-1; i>=0; i--)
{
for(j=i+AQFTk; j>i; j--)
if(j<_n_upper)
CRinv(x[j], x[i], j-i+1);
H(x[i]);
}
for(i=0; i<_n_upper; i++)
measBits[i] = MeasZ(x[i]);
}
int main()
{
qbit x[_n_upper];
qbit y[_n];
AQFT(x,y); //Implement AQFT
return 0;
}
