Solutions

sheet1/F/.vscode/settings.json

@@ -0,0 +1,5 @@
+{
+    "files.associations": {
+        "ostream": "cpp"
+    }
+}

sheet1/F/Makefile

@@ -0,0 +1,30 @@
+#
+# use GNU-Compiler tools
+COMPILER=GCC_
+# alternatively from the shell
+# export COMPILER=GCC_
+# or, alternatively from the shell
+# make COMPILER=GCC_
+
+# use Intel compilers
+#COMPILER=ICC_
+
+# use PGI compilers
+# COMPILER=PGI_
+
+
+SOURCES = main.cpp 
+OBJECTS = $(SOURCES:.cpp=.o)
+
+PROGRAM	= main.${COMPILER}
+
+# uncomment the next to lines for debugging and detailed performance analysis
+CXXFLAGS += -g
+LINKFLAGS += -g
+# do not use -pg with PGI compilers
+
+ifndef COMPILER
+  COMPILER=GCC_
+endif
+
+include ../${COMPILER}default.mk

sheet1/F/main.cpp

@@ -0,0 +1,130 @@
+#include <iostream>
+#include "mayer_primes.h"
+#include <vector>
+#include <algorithm>
+#include "timing.h"
+using namespace std;
+
+unsigned int single_goldbach(unsigned int k)
+{
+    unsigned int decomp = 0;
+    unsigned int p,q;
+    vector<unsigned int> primes = get_primes(k);
+    for(unsigned i=0; i< primes.size(); i++)
+    {
+        p = primes.at(i);
+        if(p > k/2)
+        {
+            break;
+        }
+        
+        q = k-p;
+        vector<unsigned int>::iterator it = lower_bound(primes.begin(), primes.end(), q);
+        if (it != primes.end() && *it == q)
+        {
+            decomp++;
+        }
+
+    }
+
+    return decomp;
+}
+
+//one can call single_goldbach but this way it is faster cause you dont have to regenerate primes and search in list
+vector<unsigned int> count_goldbach(unsigned int n)
+{
+     vector<unsigned int> counts(n+1,0); 
+     unsigned int p,q;
+     vector<unsigned int> primes = get_primes(n);
+     for(unsigned int j=0; j<primes.size(); j++)
+     {
+         for(unsigned i=j; i< primes.size(); i++) //start a j elsewise double counting of pairs
+        {
+           p = primes.at(j);
+           q = primes.at(i);
+
+           if(p+q > n)
+           {
+            continue;
+           }
+           counts[p+q] += 1;
+
+        }
+     }
+     return counts;
+}
+
+vector<vector<vector<unsigned int>>> count_goldbach_all(unsigned int n)
+{
+     vector<vector<vector<unsigned int>>> counts(n+1,{{0,0}}); 
+     unsigned int p,q;
+     vector<unsigned int> primes = get_primes(n);
+     for(unsigned int j=0; j<primes.size(); j++)
+     {
+         for(unsigned i=j; i< primes.size(); i++)
+        {
+           p = primes.at(j);
+           q = primes.at(i);
+
+           if(p+q > n)
+           {
+            continue;
+           }
+           if (counts[p+q].at(0).at(0) == 0) //first pair found
+           {
+            counts[p+q]={{p,q}};
+           }
+           else{
+            counts[p+q].push_back({p,q});
+           }
+
+        }
+     }
+     return counts;
+}
+
+int main()
+{   
+
+    //2
+    cout << single_goldbach(694) << endl;
+
+    //3
+    vector<unsigned int> counts = count_goldbach(100000);
+    cout << (max_element(counts.begin(), counts.end()) - counts.begin()) << endl; 
+
+    //4
+    vector<unsigned int> nvalues = {10000, 100000, 400000, 1000000, 2000000, 10000000};
+    double time;
+    unsigned int n;
+    for(unsigned int i = 0; i< nvalues.size(); i++)
+    {
+        n = nvalues.at(i);
+        tic();
+        count_goldbach(n);
+        time = toc();
+        cout << "Time for n=" << n << ": " << time << endl;
+    }
+    /*Time for n=10000: 0.0006853
+    Time for n=100000: 0.0371858
+    Time for n=400000: 0.505129
+    Time for n=1000000: 2.85873
+    Time for n=2000000: 15.0026
+    Time for n=10000000: 549.658*/
+
+
+    //*)
+    unsigned int n2 = 694;
+    vector<vector<vector<unsigned int>>> counts2 = count_goldbach_all(n2);
+    for(unsigned int i=4; i<n2; i+=2)
+    {
+        cout << "Decomp for " << i << ":"<<endl;
+        for(unsigned int j=0; j<counts2.at(i).size(); j++)
+        {
+            //cout << counts2.at(i).at(j).at(0) << "+" << counts2.at(i).at(j).at(1) << "=" << i << endl;
+        }
+    }
+
+
+    return 0;
+}

sheet1/F/mayer_primes.h

@@ -0,0 +1,73 @@
+#pragma once
+
+#include <cstring> //memset
+#include <vector>
+//using namespace std;
+
+/** \brief Determines all prime numbers in interval [2, @p max].
+ *
+ *  The sieve of Eratosthenes is used.
+ *
+ *  The implementation originates from <a href="http://code.activestate.com/recipes/576559-fast-prime-generator/">Florian Mayer</a>.
+ *
+ * \param[in]   max end of interval for the prime number search.
+ * \return  vector of prime numbers @f$2,3,5, ..., p<=max @f$.
+ *
+ * \copyright
+ * Copyright (c) 2008 Florian Mayer          (adapted by Gundolf Haase 2018)
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ *  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED,
+ * INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ *
+ */
+template <class T>
+std::vector<T> get_primes(T max)
+{
+    std::vector<T> primes;
+    char *sieve;
+    sieve = new char[max / 8 + 1];
+    // Fill sieve with 1
+    memset(sieve, 0xFF, (max / 8 + 1) * sizeof(char));
+    for (T x = 2; x <= max; x++)
+    {
+        if (sieve[x / 8] & (0x01 << (x % 8))) {
+            primes.push_back(x);
+            // Is prime. Mark multiplicates.
+            for (T j = 2 * x; j <= max; j += x)
+            {
+                sieve[j / 8] &= ~(0x01 << (j % 8));
+			}
+        }
+	}
+    delete[] sieve;
+    return primes;
+}
+
+//---------------------------------------------------------------
+//int main()        // by  Florian Mayer
+//{g++ -O3  -std=c++14 -fopenmp main.cpp && ./a.out
+//    vector<unsigned long> primes;
+//    primes = get_primes(10000000);
+//    // return 0;
+//    // Print out result.
+//    vector<unsigned long>::iterator it;
+//    for(it=primes.begin(); it < primes.end(); it++)
+//        cout << *it << " ";
+//
+//    cout << endl;
+//    return 0;
+//}

sheet1/F/mayer_primes.h:Zone.Identifier

@@ -0,0 +1,3 @@
+[ZoneTransfer]
+ZoneId=3
+HostUrl=https://imsc.uni-graz.at/haasegu/Lectures/Math2CPP/Examples/goldbach/mayer_primes.h

sheet1/F/timing.h

@@ -0,0 +1,51 @@
+//
+//    Gundolf Haase, Oct 18 2024
+//
+#pragma once
+#include <chrono>                  // timing
+#include <stack>
+
+//using Clock = std::chrono::system_clock;   //!< The wall clock timer chosen
+using Clock = std::chrono::high_resolution_clock;
+using TPoint= std::chrono::time_point<Clock>;
+
+// [Galowicz, C++17 STL Cookbook, p. 29]
+
+std::stack<TPoint> MyStopWatch; //!< starting time of stopwatch
+
+/** Starts stopwatch timer.
+ *  Use as @code tic(); myfunction(...) ; double tsec = toc();  @endcode
+ * 
+ *  The timining can be nested and the recent time point is stored on top of the stack.
+ * 
+ *  @return recent time point
+ *  @see toc
+ */
+ auto tic()
+{
+    MyStopWatch.push(Clock::now());
+    return MyStopWatch.top();
+}
+ 
+/** Returns the elapsed time from stopwatch. 
+ * 
+ * The time point from top of the stack is used
+ * if time point @p t_b is not passed as input parameter.
+ * Use as @code tic(); myfunction(...) ; double tsec = toc();  @endcode
+ * or as @code auto t_b = tic(); myfunction(...) ; double tsec = toc(t_b);  @endcode
+ * The last option is to be used in the case of 
+ * non-nested but overlapping time measurements.
+ * 
+ * @param[in]  t_b start time of some stop watch
+ * @return elapsed time in seconds.
+ *
+*/ 
+ double toc(TPoint const &t_b = MyStopWatch.top())
+{
+    // https://en.cppreference.com/w/cpp/chrono/treat_as_floating_point
+    using Unit      = std::chrono::seconds;
+    using FpSeconds = std::chrono::duration<double, Unit::period>;        
+    auto t_e = Clock::now();
+    MyStopWatch.pop();
+    return FpSeconds(t_e-t_b).count();
+}