Written in QScript via Quantum Playground You can run the full simulation for this after compiling [here](https://www.quantumplayground.net/#/playground/5681717746597888). ## Shor's Algorithm ``` // QScript version of an example from libquantum library. num_regs = 4 proc swaptheleads width for i = 0; i < width; i++ CNot i, width + i CNot width + i, i CNot i, width + i endfor endproc proc swaptheleads_omuln_controlled control, width for i = 0; i < width; i++ Toffoli control, width + i, 2 * width + i + 2 Toffoli control, 2 * width + i + 2, width + i Toffoli control, width + i, 2 * width + i + 2 endfor endproc proc test_sum compare, width if compare & (1 << (width - 1)) CNot 2 * width - 1, width - 1 SigmaX 2 * width - 1 CNot 2 * width - 1, 0 else SigmaX 2 * width - 1 CNot 2 * width - 1, width - 1 endif for i = (width - 2) ; i > 0; i-- if compare & (1 << i) //is bit i set in compare? Toffoli i + 1, width + i, i SigmaX width + i Toffoli i + 1, width + i, 0 else SigmaX width + i Toffoli i + 1, width + i, i endif endfor if compare & 1 SigmaX width Toffoli width, 1, 0 endif Toffoli 2 * width + 1, 0, 2 * width //set output to 1 if enabled and b < compare if compare & 1 Toffoli width, 1, 0 SigmaX width endif for i = 1; i <= (width - 2) ; i++ if compare & (1 << i) //is bit i set in compare? Toffoli i + 1, width + i, 0 SigmaX width + i Toffoli i + 1, width + i, i else Toffoli i + 1, width + i, i SigmaX width + i endif endfor if compare & (1 << (width - 1)) CNot 2 * width - 1, 0 SigmaX 2 * width - 1 CNot 2 * width - 1, width - 1 else CNot 2 * width - 1, width - 1 SigmaX 2 * width - 1 endif endproc // This is a semi-quantum fulladder. It adds to b_in // a c-number. Carry-in bit is c_in and carry_out is // c_out. xlt-l and L are enablebits. See documentation // for further information proc muxfa a, b_in, c_in, c_out, xlt_l, L, total //a, if a == 0//00 Toffoli b_in, c_in, c_out CNot b_in, c_in endif if a == 3//11 Toffoli L, c_in, c_out CNot L, c_in Toffoli b_in, c_in, c_out CNot b_in, c_in endif if a == 1//01 Toffoli L, xlt_l, b_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, b_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, c_in Toffoli b_in, c_in, c_out CNot b_in, c_in endif if a == 2//10 SigmaX xlt_l Toffoli L, xlt_l, b_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, b_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, c_in Toffoli b_in, c_in, c_out CNot b_in, c_in SigmaX xlt_l endif endproc // This is just the inverse operation of the semi-quantum fulladder proc muxfa_inv a, b_in, c_in, c_out, xlt_l, L, total //a, if a == 0//00 CNot b_in, c_in Toffoli b_in, c_in, c_out endif if a == 3//11 CNot b_in, c_in Toffoli b_in, c_in, c_out CNot L, c_in Toffoli L, c_in, c_out endif if a == 1//01 CNot b_in, c_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, c_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, b_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, b_in endif if a == 2//10 SigmaX xlt_l CNot b_in, c_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, c_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, b_in Toffoli b_in, c_in, c_out Toffoli L, xlt_l, b_in SigmaX xlt_l endif endproc // This is a semi-quantum halfadder. It adds to b_in // a c-number. Carry-in bit is c_in and carry_out is // not necessary. xlt-l and L are enablebits. See // documentation for further information proc muxha a, b_in, c_in, xlt_l, L, total //a, if a == 0//00 CNot b_in, c_in endif if a == 3//11 CNot L, c_in CNot b_in, c_in endif if a == 1//01 Toffoli L, xlt_l, c_in CNot b_in, c_in endif if a == 2//10 SigmaX xlt_l Toffoli L, xlt_l, c_in CNot b_in, c_in SigmaX xlt_l endif endproc // just the inverse of the semi quantum-halfadder proc muxha_inv a, b_in, c_in, xlt_l, L, total //a, if a == 0//00 CNot b_in, c_in endif if a == 3//11 CNot b_in, c_in CNot L, c_in endif if a == 1//01 CNot b_in, c_in Toffoli L, xlt_l, c_in endif if a == 2//10 SigmaX xlt_l CNot b_in, c_in Toffoli L, xlt_l, c_in SigmaX xlt_l endif endproc // proc madd a, a_inv, width total = num_regs * width + 2 for i = 0; i < width - 1; i++ if (1 << i) & a j = 1 << 1 else j = 0 endif if (1 << i) & a_inv j += 1 endif muxfa j, width + i, i, i + 1, 2 * width, 2 * width + 1, total endfor j = 0 if (1 << (width - 1)) & a j = 2 endif if (1 << (width - 1)) & a_inv j += 1 endif muxha j, 2 * width - 1, width - 1, 2 * width, 2 * width + 1, total endproc proc madd_inv a, a_inv, width total = num_regs * width + 2 j = 0 if (1 << (width - 1)) & a j = 2; endif if (1 << (width - 1)) & a_inv j += 1; endif muxha_inv j, width - 1, 2 * width - 1, 2 * width, 2 * width + 1, total for i = width - 2; i >= 0; i-- if (1 << i) & a j = 1 << 1 else j = 0 endif if (1 << i) & a_inv j += 1 endif muxfa_inv j, i, width + i, width + 1 + i, 2 * width, 2 * width + 1, total endfor endproc proc addn N, a, width //add a to register reg (mod N) test_sum N - a, width //xlt N-a madd (1 << (width)) + a - N, a, width //madd 2^K+a-N endproc proc addn_inv N, a, width //inverse of add a to register reg (mod N) CNot 2 * width + 1, 2 * width //Attention! cnot gate instead of not, as in description madd_inv (1 << (width)) - a, N - a, width //madd 2^K+(N-a)-N = 2^K-a swaptheleads width test_sum a, width endproc proc add_mod_n N, a, width //add a to register reg (mod N) and clear the scratch bits addn N, a, width addn_inv N, a, width endproc proc emul a, L, width for i = width - 1; i >= 0; i-- if (a >> i) & 1 Toffoli 2 * width + 2, L, width + i endif endfor endproc proc muln N, a, ctl, width //ctl tells, which bit is the external enable bit L = 2 * width + 1 Toffoli ctl, 2 * width + 2, L emul a % N, L, width Toffoli ctl, 2 * width + 2, L for i = 1; i < width; i++ Toffoli ctl, 2 * width + 2 + i, L add_mod_n N, ((1 << i) * a) % N, width Toffoli ctl, 2 * width + 2 + i, L endfor endproc proc muln_inv N, a, ctl, width //ctl tells, which bit is the external enable bit L = 2 * width + 1 a = QMath.inverseMod(N, a) for i = width - 1; i > 0; i-- Toffoli ctl, 2 * width + 2 + i, L add_mod_n N, N - ((1 << i) * a) % N, width Toffoli ctl, 2 * width + 2 + i, L endfor Toffoli ctl, 2 * width + 2, L emul a % N, L, width Toffoli ctl, 2 * width + 2, L endproc proc mul_mod_n N, a, ctl, width muln N, a, ctl, width swaptheleads_omuln_controlled ctl, width muln_inv N, a, ctl, width endproc proc exp_mod_n N, x, width_input, width SigmaX 2 * width + 2 for i = 1; i <= width_input; i++ f = x % N //compute for j = 1; j < i; j++ f *= f //x^2^(i-1) f = f % N endfor mul_mod_n N, f, 3 * width + 1 + i, width endfor endproc proc FindFactors N x = 0 if N < 15 Print "Invalid number!" Breakpoint endif width = QMath.getWidth(N * N) swidth = QMath.getWidth(N) for x; (QMath.gcd(N, x) > 1) || (x < 2); x x = Math.floor(Math.random() * 10000) % N endfor Print "Random seed: " + x for i = 0; i < width; i++ Hadamard i endfor ShiftLeft 3 * swidth + 2 exp_mod_n N, x, width, swidth for i = 0; i < 3 * swidth + 2; i++ MeasureBit i endfor ShiftRight 3 * swidth + 2 InvQFT 0, width for i = 0; i < width / 2; i++ CNot i, width - i - 1 CNot width - i - 1, i CNot i, width - i - 1 endfor for trycnt = 100; trycnt >= 0; trycnt-- Measure c = measured_value if c == 0 Print "Measured zero, try again." continue endif q = 1 << width Print "Measured " + c + " (" + c / q + ")" tmp = QMath.fracApprox(c, q, width) c = tmp[0]; q = tmp[1]; Print "fractional approximation is " + c + "/" + q if (q % 2 == 1) && (2 * q < (1 << width)) Print "Odd denominator, trying to expand by 2." q *= 2 endif if q % 2 == 1 Print "Odd period, try again." continue endif Print "Possible period is " + q a = QMath.ipow(x, q / 2) + 1 % N b = QMath.ipow(x, q / 2) - 1 % N a = QMath.gcd(N, a) b = QMath.gcd(N, b) if a > b factor = a else factor = b endif if (factor < N) && (factor > 1) Display "