quantum_simulator.cpp

/**
 * High-Performance Quantum Circuit Simulator in C++
 *
 * This implements a basic quantum circuit simulator using state vector simulation.
 * Designed for performance comparison with Python and Rust implementations.
 *
 * Features:
 * - State vector representation
 * - Basic quantum gates (RY, Rot, CNOT)
 * - Circuit execution benchmarking
 * - Optimized with Eigen library
 */

#include <iostream>
#include <vector>
#include <complex>
#include <cmath>
#include <chrono>
#include <random>
#include <iomanip>
#include <fstream>
#include <unistd.h>
#include <sys/resource.h>

using namespace std;
using namespace std::chrono;

// Type aliases
using Complex = complex<double>;
using StateVector = vector<Complex>;

const double PI = 3.14159265358979323846;

/**
 * Get current RSS (Resident Set Size) memory usage in MB
 */
double getMemoryUsageMB() {
#ifdef __APPLE__
    // macOS implementation
    struct rusage usage;
    getrusage(RUSAGE_SELF, &usage);
    return usage.ru_maxrss / 1024.0 / 1024.0;  // Convert bytes to MB on macOS
#elif __linux__
    // Linux implementation
    ifstream stat_stream("/proc/self/stat", ios_base::in);
    string pid, comm, state, ppid, pgrp, session, tty_nr;
    string tpgid, flags, minflt, cminflt, majflt, cmajflt;
    string utime, stime, cutime, cstime, priority, nice;
    string O, itrealvalue, starttime;
    unsigned long vsize;
    long rss;

    stat_stream >> pid >> comm >> state >> ppid >> pgrp >> session >> tty_nr
                >> tpgid >> flags >> minflt >> cminflt >> majflt >> cmajflt
                >> utime >> stime >> cutime >> cstime >> priority >> nice
                >> O >> itrealvalue >> starttime >> vsize >> rss;
    stat_stream.close();

    long page_size_kb = sysconf(_SC_PAGE_SIZE) / 1024;
    return rss * page_size_kb / 1024.0;  // Convert to MB
#else
    return 0.0;  // Unknown platform
#endif
}

/**
 * Quantum Simulator Class
 */
class QuantumSimulator {
private:
    int n_qubits;
    int state_size;
    StateVector state;

    // Get the index after applying a gate to specific qubits
    inline int getIndex(int idx, int qubit, int bit) const {
        return (idx & ~(1 << qubit)) | (bit << qubit);
    }

public:
    QuantumSimulator(int qubits) : n_qubits(qubits), state_size(1 << qubits) {
        state.resize(state_size);
        reset();
    }

    // Reset to |0...0> state
    void reset() {
        fill(state.begin(), state.end(), Complex(0.0, 0.0));
        state[0] = Complex(1.0, 0.0);
    }

    // Get current state
    const StateVector& getState() const {
        return state;
    }

    // Apply RY gate (rotation around Y-axis)
    void applyRY(int qubit, double theta) {
        double cos_half = cos(theta / 2.0);
        double sin_half = sin(theta / 2.0);

        StateVector new_state(state_size);

        for (int i = 0; i < state_size; ++i) {
            int idx0 = getIndex(i, qubit, 0);
            int idx1 = getIndex(i, qubit, 1);

            if (i & (1 << qubit)) {
                // Qubit is in |1> state
                new_state[i] = -sin_half * state[idx0] + cos_half * state[idx1];
            } else {
                // Qubit is in |0> state
                new_state[i] = cos_half * state[idx0] + sin_half * state[idx1];
            }
        }

        state = new_state;
    }

    // Apply RZ gate (rotation around Z-axis)
    void applyRZ(int qubit, double phi) {
        Complex phase0 = exp(Complex(0.0, -phi / 2.0));
        Complex phase1 = exp(Complex(0.0, phi / 2.0));

        for (int i = 0; i < state_size; ++i) {
            if (i & (1 << qubit)) {
                state[i] *= phase1;
            } else {
                state[i] *= phase0;
            }
        }
    }

    // Apply RX gate (rotation around X-axis)
    void applyRX(int qubit, double theta) {
        double cos_half = cos(theta / 2.0);
        double sin_half = sin(theta / 2.0);
        Complex i_unit(0.0, 1.0);

        StateVector new_state(state_size);

        for (int i = 0; i < state_size; ++i) {
            int idx0 = getIndex(i, qubit, 0);
            int idx1 = getIndex(i, qubit, 1);

            if (i & (1 << qubit)) {
                new_state[i] = cos_half * state[idx1] - i_unit * sin_half * state[idx0];
            } else {
                new_state[i] = cos_half * state[idx0] - i_unit * sin_half * state[idx1];
            }
        }

        state = new_state;
    }

    // Apply general rotation gate (Rot = RZ(phi) RY(theta) RZ(omega))
    void applyRot(int qubit, double phi, double theta, double omega) {
        applyRZ(qubit, omega);
        applyRY(qubit, theta);
        applyRZ(qubit, phi);
    }

    // Apply CNOT gate
    void applyCNOT(int control, int target) {
        for (int i = 0; i < state_size; ++i) {
            if (i & (1 << control)) {
                // Control qubit is |1>, flip target
                int flipped = i ^ (1 << target);
                if (i < flipped) {
                    swap(state[i], state[flipped]);
                }
            }
        }
    }

    // Measure expectation value of Pauli-Z on a qubit
    double measurePauliZ(int qubit) const {
        double expectation = 0.0;

        for (int i = 0; i < state_size; ++i) {
            double prob = norm(state[i]);
            if (i & (1 << qubit)) {
                expectation -= prob;  // |1> state contributes -1
            } else {
                expectation += prob;  // |0> state contributes +1
            }
        }

        return expectation;
    }

    // Print state vector (for debugging)
    void printState() const {
        cout << "State vector:" << endl;
        for (int i = 0; i < min(8, state_size); ++i) {
            cout << "|" << i << ">: " << state[i] << endl;
        }
        if (state_size > 8) {
            cout << "..." << endl;
        }
    }
};

/**
 * Quantum Circuit for benchmarking
 */
class QuantumCircuit {
private:
    int n_qubits;
    int n_layers;
    vector<vector<vector<double>>> params;
    vector<double> input_data;

public:
    QuantumCircuit(int qubits, int layers) : n_qubits(qubits), n_layers(layers) {
        // Initialize random parameters
        random_device rd;
        mt19937 gen(rd());
        uniform_real_distribution<> dis(0.0, PI);

        params.resize(n_layers);
        for (int l = 0; l < n_layers; ++l) {
            params[l].resize(n_qubits);
            for (int q = 0; q < n_qubits; ++q) {
                params[l][q] = {dis(gen), dis(gen), dis(gen)};
            }
        }

        // Initialize random input
        uniform_real_distribution<> input_dis(-PI, PI);
        input_data.resize(n_qubits);
        for (int i = 0; i < n_qubits; ++i) {
            input_data[i] = input_dis(gen);
        }
    }

    // Execute the circuit
    double execute(QuantumSimulator& sim) {
        sim.reset();

        // Encoding layer
        for (int i = 0; i < n_qubits; ++i) {
            sim.applyRY(i, input_data[i]);
        }

        // Variational layers
        for (int layer = 0; layer < n_layers; ++layer) {
            // Parameterized rotations
            for (int i = 0; i < n_qubits; ++i) {
                sim.applyRot(i, params[layer][i][0], params[layer][i][1], params[layer][i][2]);
            }

            // Entangling layer
            for (int i = 0; i < n_qubits; ++i) {
                sim.applyCNOT(i, (i + 1) % n_qubits);
            }
        }

        // Measurement
        return sim.measurePauliZ(0);
    }
};

/**
 * Benchmark function
 */
void runBenchmark(int n_qubits, int n_layers, int n_executions) {
    cout << "\n" << string(60, '=') << endl;
    cout << "C++ Quantum Circuit Benchmark" << endl;
    cout << string(60, '=') << endl;
    cout << "Qubits: " << n_qubits << endl;
    cout << "Layers: " << n_layers << endl;
    cout << "Executions: " << n_executions << endl;
    cout << string(60, '-') << endl;

    // Get initial memory
    double mem_before = getMemoryUsageMB();

    // Create simulator and circuit
    QuantumSimulator sim(n_qubits);
    QuantumCircuit circuit(n_qubits, n_layers);

    // Warm-up run
    circuit.execute(sim);

    // Benchmark execution
    auto start = high_resolution_clock::now();

    for (int i = 0; i < n_executions; ++i) {
        double result = circuit.execute(sim);
        (void)result;  // Suppress unused variable warning
    }

    auto end = high_resolution_clock::now();
    auto duration = duration_cast<microseconds>(end - start);

    // Get final memory
    double mem_after = getMemoryUsageMB();
    double mem_used = mem_after - mem_before;

    // Calculate metrics
    double total_time = duration.count() / 1e6;  // Convert to seconds
    double avg_time = total_time / n_executions;
    double exec_per_sec = n_executions / total_time;

    cout << fixed << setprecision(6);
    cout << "Total time: " << total_time << "s" << endl;
    cout << "Average time per execution: " << avg_time << "s" << endl;
    cout << setprecision(2);
    cout << "Executions per second: " << exec_per_sec << endl;
    cout << "Memory used: " << max(mem_used, 0.0) << " MB" << endl;
    cout << string(60, '=') << endl;
}

/**
 * Main function
 */
int main() {
    cout << "\n" << string(60, '=') << endl;
    cout << "HIGH-PERFORMANCE QUANTUM SIMULATOR (C++)" << endl;
    cout << string(60, '=') << endl;

    // Run benchmarks with different configurations
    runBenchmark(4, 3, 100);
    runBenchmark(6, 3, 50);
    runBenchmark(8, 3, 20);

    // Test simulator functionality
    cout << "\nFunctionality Test:" << endl;
    cout << string(60, '-') << endl;
    QuantumSimulator test_sim(2);
    test_sim.applyRY(0, PI / 4);
    test_sim.applyCNOT(0, 1);

    cout << "2-qubit circuit test:" << endl;
    cout << "Pauli-Z expectation on qubit 0: " << test_sim.measurePauliZ(0) << endl;
    cout << "Pauli-Z expectation on qubit 1: " << test_sim.measurePauliZ(1) << endl;

    cout << "\nBenchmark completed!" << endl;

    return 0;
}