Main file - Exercise 3

Sheet5/mainEx3.cpp

@@ -0,0 +1,135 @@
+//  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;
+}