1 // Shor's algorithm
2 // Source: ./examples/shor.cpp
3 #include <cmath>
4 #include <cstdlib>
5 #include <iostream>
6 #include <tuple>
7 #include <vector>
8 
9 #include "qpp.h"
10 
main()11 int main() {
12     using namespace qpp;
13 
14     bigint N = 15;                                  // number to factor
15     bigint a = rand(static_cast<bigint>(3), N - 1); // random co-prime with N
16     while (gcd(a, N) != 1) {
17         a = rand(static_cast<bigint>(3), N - 1);
18     }
19     // qubits required for half of the circuit, in total we need 2n qubits
20     // if you know the order 'r' of 'a', then you can take the smallest 'n' s.t.
21     // 2^n >= 2 * r^2, i.e., n = ceil(log2(2 * r^2))
22     auto n = static_cast<idx>(std::ceil(2 * std::log2(N)));
23     auto D = static_cast<idx>(std::llround(std::pow(2, n))); // dimension 2^n
24 
25     std::cout << ">> Factoring N = " << N << " with coprime a = " << a << '\n';
26     std::cout << ">> n = " << n << ", D = 2^n = " << D << " and 2^(2n) = ";
27     std::cout << D * D << '\n';
28 
29     // vector with labels of the first half of the qubits
30     std::vector<idx> first_subsys(n);
31     std::iota(std::begin(first_subsys), std::end(first_subsys), 0);
32 
33     // vector with labels of the second half of the qubits
34     std::vector<idx> second_subsys(n);
35     std::iota(std::begin(second_subsys), std::end(second_subsys), n);
36 
37     // QUANTUM STAGE
38     // prepare the initial state |0>^\otimes n \otimes |0...01>
39     ket psi = kron(st.zero(2 * n - 1), 1_ket);
40 
41     // apply Hadamards H^\otimes n on first half of the qubits
42     for (idx i = 0; i < n; ++i) {
43         psi = apply(psi, gt.H, {i});
44     }
45 
46     // perform the modular exponentiation as a sequence of
47     // modular multiplications
48     for (idx i = 0; i < n; ++i) {
49         // compute 2^(n-i-1) mod N
50         idx j = static_cast<idx>(std::llround(std::pow(2, n - i - 1)));
51         // compute the a^(2^(n-i-1)) mod N
52         idx aj = modpow(a, j, N);
53         // apply the controlled modular multiplication
54         psi = applyCTRL(psi, gt.MODMUL(aj, N, n), {i}, second_subsys);
55     }
56 
57     // apply inverse QFT on first half of the qubits
58     psi = applyTFQ(psi, first_subsys);
59     // END QUANTUM STAGE
60 
61     // FIRST MEASUREMENT STAGE
62     auto measured1 = measure_seq(psi, first_subsys); // measure first n qubits
63     std::vector<idx> vect_results1 = std::get<RES>(measured1); // results
64     double prob1 = std::get<PROB>(measured1); // probability of the result
65     idx n1 = multiidx2n(vect_results1, std::vector<idx>(n, 2)); // binary to int
66     auto x1 = static_cast<double>(n1) / D; // multiple of 1/r
67 
68     std::cout << ">> First measurement:  " << disp(vect_results1, " ") << " ";
69     std::cout << "i.e., j = " << n1 << " with probability " << prob1;
70     std::cout << '\n';
71 
72     bool failed = true;
73     idx r1 = 0, c1 = 0;
74     for (auto&& elem : convergents(x1, 10)) {
75         std::tie(c1, r1) = elem;
76         auto c1r1 = static_cast<double>(c1) / r1;
77         if (abs(x1 - c1r1) < 1. / std::pow(2, (n - 1.) / 2.)) {
78             failed = false;
79             break;
80         }
81     }
82     if (failed) {
83         std::cout << ">> Factoring failed at stage 1, please try again!\n";
84         std::exit(EXIT_FAILURE);
85     }
86     // END FIRST MEASUREMENT STAGE
87 
88     // SECOND MEASUREMENT STAGE
89     auto measured2 = measure_seq(psi, first_subsys); // measure first n qubits
90     std::vector<idx> vect_results2 = std::get<RES>(measured2); // results
91     double prob2 = std::get<PROB>(measured2); // probability of the result
92     idx n2 = multiidx2n(vect_results2, std::vector<idx>(n, 2)); // binary to int
93     auto x2 = static_cast<double>(n2) / D; // multiple of 1/r
94 
95     std::cout << ">> Second measurement: " << disp(vect_results2, " ") << " ";
96     std::cout << "i.e., j = " << n2 << " with probability " << prob2;
97     std::cout << '\n';
98 
99     failed = true;
100     idx r2 = 0, c2 = 0;
101     for (auto&& elem : convergents(x2, 10)) {
102         std::tie(c2, r2) = elem;
103         auto c2r2 = static_cast<double>(c2) / r2;
104         if (abs(x2 - c2r2) < 1. / std::pow(2, (n - 1.) / 2.)) {
105             failed = false;
106             break;
107         }
108     }
109     if (failed) {
110         std::cout << ">> Factoring failed at stage 2, please try again!\n";
111         std::exit(EXIT_FAILURE);
112     }
113     // END SECOND MEASUREMENT STAGE
114 
115     // THIRD POST-PROCESSING STAGE
116     idx r = lcm(r1, r2); // candidate order of a mod N
117     std::cout << ">> r = " << r << ", a^r mod N = " << modpow(a, r, N) << '\n';
118     if (r % 2 == 0 && modpow(a, r / 2, N) != static_cast<bigint>(N - 1)) {
119         std::cout << ">> Possible factors: ";
120         std::cout << gcd(modpow(a, r / 2, N) - 1, N) << " ";
121         std::cout << gcd(modpow(a, r / 2, N) + 1, N) << '\n';
122     } else {
123         std::cout << ">> Factoring failed at stage 3, please try again!\n";
124         std::exit(EXIT_FAILURE);
125     }
126     // END THIRD POST-PROCESSING STAGE
127 }
128