//--- 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 #include #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(b, a); CNOT(a, b); CNOT(b, a); } module cRz (qbit ctrl, qbit target, double angle) { Rz (target, angle/2); CNOT (ctrl, target); Rz (target, -1*angle/2); CNOT (ctrl, target); } 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 (ctrl1, ctrl2); cRz (ctrl1, target, -angle/2); CNOT (ctrl1, ctrl2); 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(ctrl, x[i]); // 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(ctrl, x[i]); // 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(ctrl, x[i]); // 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(ctrl, x[i]); // 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(ctrl, x[i]); // 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(ctrl, x[i]); // 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(ctrl, x[i]); // 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(ctrl, x[i]); // 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(control, target) CNOT(Reg6[_n-1], Aqbit[0]); 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(Reg51[i], Reg2[_n-_n+i]); } QFT(Reg2); QADDn(Reg2, Reg4); iCQADD_1(Aqbit[0], Reg2); iQFT(Reg2); QFT_2(Reg10, Reg00); for (i=0; i<_n; i++) CNOT(Reg2[i], Reg3[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(Reg51[i], Reg3[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(Reg51[i], Reg3[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(Reg2[i], Reg3[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(Reg51[i], Reg2[_n-_n+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(control, target) CNOT(Reg6[_n-1], Aqbit[0]); 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(R10[i], R30[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(R10[i], R30[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(ctrl, R10[i]); Toffoli(ctrl, x[i], R10[i]); CNOT(ctrl, R10[i]); } 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(ctrl, q); RotBy(q,k-1, 1); //positive angle CNOT(ctrl, q); } 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; }