Quantum Computing Theory

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct QuantumState {
5
    int n; // number of qubits
6
    vector<complex<double>> amp; // state vector of size 2^n
7
8
    QuantumState(int n_) : n(n_) {
9
        size_t N = 1ull << n;
10
        amp.assign(N, {0.0, 0.0});
11
        amp[0] = {1.0, 0.0}; // start in |0...0>
12
    }
13
14
    // Normalize the state (useful after numerical noise)
15
    void normalize() {
16
        double s = 0.0;
17
        for (auto &a : amp) s += norm(a);
18
        double inv = 1.0 / sqrt(s);
19
        for (auto &a : amp) a *= inv;
20
    }
21
22
    // Apply a single-qubit gate (2x2 complex matrix) to qubit q (0 = least significant bit)
23
    void applySingleQubitGate(int q, array<complex<double>,4> U) {
24
        size_t N = amp.size();
25
        size_t mask = 1ull << q;
26
        for (size_t i = 0; i < N; ++i) {
27
            if ((i & mask) == 0) {
28
                size_t j = i | mask; // partner index with bit q flipped
29
                complex<double> a = amp[i];
30
                complex<double> b = amp[j];
31
                amp[i] = U[0]*a + U[1]*b;  // m00*a + m01*b
32
                amp[j] = U[2]*a + U[3]*b;  // m10*a + m11*b
33
            }
34
        }
35
    }
36
37
    // Apply a CNOT with control c and target t
38
    void applyCNOT(int c, int t) {
39
        if (c == t) return;
40
        size_t N = amp.size();
41
        size_t cm = 1ull << c;
42
        size_t tm = 1ull << t;
43
        for (size_t i = 0; i < N; ++i) {
44
            if ((i & cm) && ((i & tm) == 0)) {
45
                size_t j = i | tm; // flip target bit
46
                swap(amp[i], amp[j]);
47
            }
48
        }
49
    }
50
51
    // Compute probabilities over all computational basis states
52
    vector<double> probabilities() const {
53
        vector<double> p(amp.size());
54
        for (size_t i = 0; i < amp.size(); ++i) p[i] = norm(amp[i]);
55
        return p;
56
    }
57
58
    // Measure all qubits in the computational basis (collapses state)
59
    uint64_t measureAll(mt19937_64 &rng) {
60
        vector<double> p = probabilities();
61
        // build CDF
62
        vector<double> cdf(p.size());
63
        partial_sum(p.begin(), p.end(), cdf.begin());
64
        uniform_real_distribution<double> dist(0.0, 1.0);
65
        double r = dist(rng);
66
        size_t idx = lower_bound(cdf.begin(), cdf.end(), r) - cdf.begin();
67
        // collapse
68
        for (size_t i = 0; i < amp.size(); ++i) amp[i] = {0.0, 0.0};
69
        amp[idx] = {1.0, 0.0};
70
        return idx;
71
    }
72
};
73
74
int main() {
75
    ios::sync_with_stdio(false);
76
    cin.tie(nullptr);
77
78
    // Build a 2-qubit Bell state: (|00> + |11>)/sqrt(2)
79
    QuantumState qs(2);
80
81
    // Hadamard gate matrix
82
    const double is2 = 1.0 / sqrt(2.0);
83
    array<complex<double>,4> H = {complex<double>(is2,0), complex<double>(is2,0),
84
                                  complex<double>(is2,0), complex<double>(-is2,0)};
85
86
    // Apply H on qubit 0 (LSB), then CNOT(0->1)
87
    qs.applySingleQubitGate(0, H);
88
    qs.applyCNOT(0, 1);
89
90
    // Show probabilities (should be ~0.5 on |00> and |11>)
91
    auto p = qs.probabilities();
92
    cout << fixed << setprecision(6);
93
    for (size_t i = 0; i < p.size(); ++i) {
94
        cout << "P(|" << bitset<2>(i) << ">) = " << p[i] << "\n";
95
    }
96
97
    // Sample measurements to highlight perfect correlation
98
    mt19937_64 rng(42);
99
    int same = 0, trials = 20;
100
    for (int t = 0; t < trials; ++t) {
101
        QuantumState trial(2);
102
        trial.applySingleQubitGate(0, H);
103
        trial.applyCNOT(0, 1);
104
        uint64_t outcome = trial.measureAll(rng);
105
        int b0 = outcome & 1;
106
        int b1 = (outcome >> 1) & 1;
107
        if (b0 == b1) same++;
108
    }
109
    cout << "Same bits count in " << trials << " trials: " << same << "\n";
110
    return 0;
111
}
112

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
struct State {
5
    int n; vector<complex<double>> a; // amplitudes
6
    State(int n_) : n(n_) { a.assign(1ull<<n, {0,0}); a[0] = {1,0}; }
7
    void normalize(){ double s=0; for(auto &x:a) s+=norm(x); double inv=1.0/sqrt(s); for(auto &x:a) x*=inv; }
8
};
9
10
// Initialize uniform superposition |s> = (1/sqrt(N)) sum_x |x>
11
void make_uniform(State &st){
12
    size_t N = st.a.size();
13
    complex<double> c = {1.0/sqrt((double)N), 0.0};
14
    for(size_t i=0;i<N;++i) st.a[i]=c;
15
}
16
17
// Oracle for single marked index m: phase flip |m> -> -|m>
18
void oracle_phase_flip(State &st, size_t m){ st.a[m] *= complex<double>{-1.0, 0.0}; }
19
20
// Diffusion operator: reflect about the mean: a_i -> 2*avg - a_i
21
void diffusion(State &st){
22
    size_t N = st.a.size();
23
    complex<double> sum = {0.0, 0.0};
24
    for(auto &x: st.a) sum += x;
25
    complex<double> avg = sum / (double)N;
26
    for(auto &x: st.a) x = 2.0*avg - x;
27
}
28
29
// Measure most probable outcome (argmax probability) without collapsing randomness
30
size_t argmax_prob(const State &st){
31
    size_t idx=0; double best=-1.0;
32
    for(size_t i=0;i<st.a.size();++i){ double p=norm(st.a[i]); if(p>best){best=p; idx=i;} }
33
    return idx;
34
}
35
36
int main(){
37
    ios::sync_with_stdio(false);
38
    cin.tie(nullptr);
39
40
    int n = 3; // try 3 or 4
41
    size_t N = 1ull<<n;
42
    size_t marked = 5; // the hidden solution index in [0, N)
43
44
    State st(n);
45
    make_uniform(st); // prepare |s>
46
47
    // Number of Grover iterations r ≈ floor( (pi/4) * sqrt(N) ) for a single marked item
48
    int r = (int)floor((M_PI/4.0)*sqrt((double)N));
49
50
    for(int i=0;i<r;++i){
51
        oracle_phase_flip(st, marked);
52
        diffusion(st);
53
    }
54
55
    // Report the most likely state and its probability
56
    size_t guess = argmax_prob(st);
57
    double p = norm(st.a[guess]);
58
59
    cout << "N=" << N << ", marked=" << marked << ", iterations=" << r << "\n";
60
    cout << "Most likely outcome: " << guess << ", probability=" << fixed << setprecision(6) << p << "\n";
61
62
    // Optional: print top probabilities
63
    vector<pair<double,size_t>> v; v.reserve(N);
64
    for(size_t i=0;i<N;++i) v.push_back({norm(st.a[i]), i});
65
    sort(v.rbegin(), v.rend());
66
    cout << "Top 4 states by probability:" << "\n";
67
    for(int i=0;i<min((size_t)4,N);++i){
68
        cout << v[i].second << ": " << v[i].first << "\n";
69
    }
70
    return 0;
71
}
72

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
using cd = complex<double>;
5
6
struct State { int n; vector<cd> a; State(int n_):n(n_){ a.assign(1ull<<n, {0,0}); a[0]={1,0}; } };
7
8
// Build a periodic state with period r: |psi> ∝ sum_{k} |k*r mod N>
9
void make_periodic(State &st, int r){
10
    size_t N = st.a.size();
11
    for(auto &x: st.a) x = {0,0};
12
    for(size_t x=0; x<N; x+=r) st.a[x] = {1.0, 0.0};
13
    // normalize
14
    double cnt = (double)( (N + r - 1) / r );
15
    cd c = {1.0/sqrt(cnt), 0.0};
16
    for(size_t x=0; x<N; x+=r) st.a[x] = c;
17
}
18
19
// Naive dense QFT: b[y] = (1/sqrt(N)) * sum_x a[x] * exp(2pi i x y / N)
20
void qft_dense(State &st){
21
    size_t N = st.a.size();
22
    vector<cd> b(N, {0,0});
23
    const double tau = 2.0 * M_PI;
24
    double invsqrtN = 1.0 / sqrt((double)N);
25
    for(size_t y=0; y<N; ++y){
26
        cd sum = {0.0, 0.0};
27
        for(size_t x=0; x<N; ++x){
28
            double ang = tau * ((double)(x*y) / (double)N);
29
            cd w = { cos(ang), sin(ang) }; // e^{i ang}
30
            sum += st.a[x] * w;
31
        }
32
        b[y] = sum * invsqrtN;
33
    }
34
    st.a.swap(b);
35
}
36
37
int main(){
38
    ios::sync_with_stdio(false);
39
    cin.tie(nullptr);
40
41
    int n = 3; // N=8
42
    State st(n);
43
44
    int r = 2; // period
45
    make_periodic(st, r); // state on |0>,|2>,|4>,|6>
46
47
    qft_dense(st);
48
49
    // Print probabilities; peaks should occur at multiples of N/r = 4
50
    cout << fixed << setprecision(6);
51
    for(size_t y=0; y<st.a.size(); ++y){
52
        cout << y << ": " << norm(st.a[y]) << "\n";
53
    }
54
    return 0;
55
}
56