//  HOW TO COMPILE ON MAC
// export LDFLAGS="-L/opt/homebrew/opt/libomp/lib"
// export CPPFLAGS="-I/opt/homebrew/opt/libomp/include"
// clang++ -std=c++17 -O3 -Xpreprocessor -fopenmp $CPPFLAGS MainEx3.cpp -L/opt/homebrew/opt/libomp/lib -lomp -o Ex3
// ./Ex3

#include <iostream>
#include <vector>
#include <chrono>
#include <cmath>
#include <algorithm>
#include <utility>
#include <iomanip>
#include <omp.h>          
#include "mayer_primes.h" 

using namespace std;

int single_goldbach(int k, const vector<int>& primes, const vector<char>& is_prime)
{
    int count = 0;

    #pragma omp parallel for reduction(+:count)
    for (size_t i = 0; i < primes.size(); ++i)
    {
        int p = primes[i];
        if (p > k / 2) continue;
        int q = k - p;
        if (q >= 0 && q < (int)is_prime.size() && is_prime[q])
            count += 1;
    }
    return count;
}

vector<pair<int,int>> count_goldbach(int n, const vector<int>& primes, const vector<char>& is_prime)
{
    vector<pair<int,int>> results(n/2 - 1); // exact size

    #pragma omp parallel for
    for (int j = 0; j < (int)results.size(); ++j)
    {
        int k = 4 + 2*j;
        int c = single_goldbach(k, primes, is_prime);
        results[j] = {k, c};
    }
    return results;
}

vector<vector<pair<int,int>>> all_decompositions(int n, const vector<int>& primes, const vector<char>& is_prime)
{
    vector<vector<pair<int,int>>> out(n/2 - 1);

    #pragma omp parallel for
    for (int j = 0; j < (int)out.size(); ++j)
    {
        int k = 4 + 2*j;
        vector<pair<int,int>> local;

        for (size_t i = 0; i < primes.size(); ++i)
        {
            int p = primes[i];
            if (p > k / 2) break;
            int q = k - p;
            if (q >= 0 && q < (int)is_prime.size() && is_prime[q])
                local.push_back({p, q});
        }

        out[j] = std::move(local);
    }
    return out;
}

vector<char> make_is_prime(int n, const vector<int>& primes)
{
    vector<char> is_prime(n + 1, 0);
    for (int p : primes)
        if (p <= n) is_prime[p] = 1;
    return is_prime;
}

int main()
{
    vector<int> test_ns = {10000, 100000, 400000};

    for (int n : test_ns)
    {
        cout << "Computing Goldbach up to n = " << n << "\n";

        // CHECKS to see if I'm actually paralellizing 
        #pragma omp parallel
        {
            #pragma omp single
            cout << "Using " << omp_get_num_threads() 
                 << " OpenMP threads\n";
        }

        // Generate primes 
        auto t0 = chrono::system_clock::now();
        vector<int> primes = get_primes<int>(n);
        auto t1 = chrono::system_clock::now();
        chrono::duration<double> gen_time = t1 - t0;

        cout << "Generated " << primes.size()
             << " primes in " << gen_time.count() << " s\n";

        vector<char> is_prime = make_is_prime(n, primes);

        // Count Goldbach decompositions
        auto start = chrono::system_clock::now();
        vector<pair<int,int>> results = count_goldbach(n, primes, is_prime);
        auto end = chrono::system_clock::now();
        chrono::duration<double> elapsed = end - start;

        // Find k with maximum decompositions
        int max_k = 0;
        int max_count = 0;

        for (auto &pr : results)
        {
            if (pr.second > max_count)
            {
                max_k = pr.first;
                max_count = pr.second;
            }
        }

        cout << "Max decompositions: k = " << max_k
             << " with " << max_count << " decompositions\n";

        cout << "Time to count decompositions: "
             << elapsed.count() << " s\n";
    }

    return 0;
}