Added ex11 solutions

exercises/ex11/solutions/factorial_binomial/compiletime_version/main_clever.cpp

+/*
+ * Programming Techniques for Scientific Simulations I
+ * HS 2020
+ * Week 11
+ */
+
+#include <iostream>
+
+// we can't change the factorial, since loops are done with recursions in
+// the functional paradigm
+template<size_t N>
+struct factorial {
+    enum{ value =  N * factorial<N-1>::value };
+};
+
+template<>
+struct factorial<0> {
+    enum { value = 1 };
+};
+
+
+// we are using the recursive formula (N, k) = (N-1, k-1) + (N-1, k)
+// with (N, 0) = (N, N) = 1
+template<size_t N, size_t k>
+struct binomial {
+    enum{ value = binomial<N-1, k-1>::value + binomial<N-1, k>::value };
+};
+
+template<size_t N>
+struct binomial<N, 0> {
+    enum{ value =  1 };
+};
+
+template<size_t N>
+struct binomial<N, N> {
+    enum{ value =  1 };
+};
+
+
+int main()
+{
+    std::cout << "Factorial:" << std::endl;
+    std::cout << "    5!=" << factorial<5>::value << std::endl;
+    std::cout << "    20!=" << factorial<20>::value << std::endl;
+    std::cout << "    21!=" << factorial<21>::value << std::endl;
+    std::cout << "    70!=" << factorial<70>::value << std::endl;
+
+    std::cout << "Binomial:" << std::endl;
+    std::cout << "    N=20, k=1 ==> " << binomial<20, 1>::value << std::endl;
+    std::cout << "    N=21, k=1 ==> " << binomial<21, 1>::value << std::endl;
+    std::cout << "    N=70, k=1 ==> " << binomial<70, 1>::value << std::endl;
+    return 0;
+}

exercises/ex11/solutions/factorial_binomial/compiletime_version/main_simple.cpp

+/*
+ * Programming Techniques for Scientific Simulations I
+ * HS 2020
+ * Week 11
+ */
+
+#include <iostream>
+
+
+template<size_t N>
+struct factorial {
+    enum{ value =  N * factorial<N-1>::value };
+};
+
+template<>
+struct factorial<0> {
+    enum { value = 1 };
+};
+
+template<size_t N, size_t k>
+struct binomial {
+    enum{ value =  factorial<N>::value / (factorial<N-k>::value * factorial<k>::value) };
+};
+
+
+int main()
+{
+    std::cout << "Factorial:" << std::endl;
+    std::cout << "    5!=" << factorial<5>::value << std::endl;
+    std::cout << "    20!=" << factorial<20>::value << std::endl;
+    std::cout << "    21!=" << factorial<21>::value << std::endl;
+    std::cout << "    70!=" << factorial<70>::value << std::endl;
+
+    std::cout << "Binomial:" << std::endl;
+    std::cout << "    N=20, k=1 ==> " << binomial<20, 1>::value << std::endl;
+    std::cout << "    N=21, k=1 ==> " << binomial<21, 1>::value << std::endl;
+    // will trigger a compiletime error
+    // std::cout << "    N=70, k=1 ==> " << binomial<70, 1>::value << std::endl;
+    return 0;
+}

exercises/ex11/solutions/factorial_binomial/constexpr_version_bonus/main_clever.cpp

+/*
+ * Programming Techniques for Scientific Simulations I
+ * HS 2020
+ * Week 11
+ */
+
+#include <iostream>
+#include <iomanip>
+#include <assert.h>
+
+// note the constexpr
+constexpr size_t factorial(const size_t & value)
+{
+    // in a constexpr function, only a very limited set of c++ works
+    // forget arrays, containers (in c++14), pointers, modifying globals, ...
+    // If you can't do the same with templates only, it's not allowed.
+    size_t res = 1;
+    for(size_t i = 1; i <= value; ++i) {
+        res *= i;
+    }
+    return res;
+}
+
+// note the constexpr
+constexpr size_t binomial(const size_t & N, const size_t & k)
+{
+    size_t res = 1;
+    for(size_t i = 1; i <= k; ++i) {
+        res *= (N+1-i);
+        assert(res % i == 0); // make sure there is no remainder (int-arithmetics)
+        res /= i; // I divide after the mult, so we dont need double-arithmetics
+    }
+    return res;
+}
+
+// we only need this echo struct to demonstrate that the constexpr
+// funtions can and will be executed during compiletime
+// (alternative: look at assembler code)
+template<size_t N>
+struct echo {
+    // this replaces the enum { value = value } trick
+    static constexpr size_t value = N;
+};
+
+
+int main()
+{
+    // our normal runtime code
+    // (no change needed, a constexpr fct is also a normal function)
+    std::cout << "Factorial:" << std::endl;
+    for(size_t N = 0; N < 10; ++N) {
+        std::cout << "    " << N << "! = " << factorial(N) << std::endl;
+    }
+
+    std::cout << "Binomial:" << std::endl;
+    for(size_t N = 1; N < 10; ++N) {
+        std::cout << "    " << std::setw(5) << "N=" << N ;
+        for(size_t k = 0; k <= N; ++k) {
+            std::cout << std::setw(5) << binomial(N, k);
+        }
+        std::cout << std::endl;
+    }
+
+    std::cout << "Limits:" << std::endl;
+    std::cout << "    20! = " << factorial(20) << std::endl;
+    std::cout << "    21! = " << factorial(21) << std::endl;
+    std::cout << "    70! = " << factorial(70) << std::endl;
+
+    std::cout << "Limits:" << std::endl;
+    std::cout << "    N=20, k=1 ==> " << binomial(20,1) << std::endl;
+    std::cout << "    N=21, k=1 ==> " << binomial(21,1) << " <- right " << std::endl;
+    std::cout << "    N=70, k=1 ==> " << binomial(70,1) << std::endl;
+
+    // using our functions in a context where the result needs to be knows
+    // at compiletime
+
+    std::cout << "Compiletime Factorial:" << std::endl;
+    std::cout << "    20! = " << echo<factorial(20)>::value << std::endl;
+    std::cout << "    21! = " << echo<factorial(21)>::value << std::endl;
+    std::cout << "    70! = " << echo<factorial(70)>::value << std::endl;
+
+    std::cout << "Compiletime Binomial:" << std::endl;
+    std::cout << "    N=20, k=1 ==> " << echo<binomial(20,1)>::value << std::endl;
+    std::cout << "    N=21, k=1 ==> " << echo<binomial(21,1)>::value << std::endl;
+    std::cout << "    N=70, k=1 ==> " << echo<binomial(70,1)>::value << std::endl;
+
+    return 0;
+}

exercises/ex11/solutions/factorial_binomial/runtime_version/main_clever.cpp

+/*
+ * Programming Techniques for Scientific Simulations I
+ * HS 2020
+ * Week 11
+ */
+
+#include <iostream>
+#include <iomanip>
+#include <assert.h>
+
+size_t factorial(const size_t & value)
+{
+    size_t res = 1;
+    for(size_t i = 1; i <= value; ++i) {
+        res *= i;
+    }
+
+    return res;
+}
+
+size_t binomial(const size_t & N, const size_t & k)
+{
+    size_t res = 1;
+    for(size_t i = 1; i <= k; ++i) {
+        res *= (N+1-i);
+        assert(res % i == 0); // make sure there is no remainder (int-arithmetics)
+        res /= i; // I divide after the mult, so we dont need double-arithmetics
+    }
+    return res;
+}
+
+
+int main()
+{
+    std::cout << "Factorial:" << std::endl;
+    for(size_t N = 0; N < 10; ++N) {
+        std::cout << "    " << N << "! = " << factorial(N) << std::endl;
+    }
+
+    std::cout << "Binomial:" << std::endl;
+    for(size_t N = 1; N < 10; ++N) {
+        std::cout << "    " << std::setw(5) << "N=" << N ;
+        for(size_t k = 0; k <= N; ++k) {
+            std::cout << std::setw(5) << binomial(N, k);
+        }
+        std::cout << std::endl;
+    }
+
+    std::cout << "Limits:" << std::endl;
+    std::cout << "    20! = " << factorial(20) << std::endl;
+    std::cout << "    21! = " << factorial(21) << std::endl;
+    std::cout << "    70! = " << factorial(70) << std::endl;
+
+    std::cout << "Limits:" << std::endl;
+    std::cout << "    N=20, k=1 ==> " << binomial(20,1) << std::endl;
+    std::cout << "    N=21, k=1 ==> " << binomial(21,1) << " <- right " << std::endl;
+    std::cout << "    N=70, k=1 ==> " << binomial(70,1) << std::endl;
+
+    return 0;
+}

exercises/ex11/solutions/factorial_binomial/runtime_version/main_simple.cpp

+/*
+ * Programming Techniques for Scientific Simulations I
+ * HS 2020
+ * Week 11
+ */
+
+#include <iostream>
+#include <iomanip>
+
+size_t factorial(const size_t & value)
+{
+    if(value == 0) {
+        return 1;
+    }
+    return value * factorial(value - 1);
+}
+
+size_t binomial(const size_t & N, const size_t & k)
+{
+    return factorial(N) / (factorial(N-k) * factorial(k));
+}
+
+
+int main()
+{
+    std::cout << "Factorial:" << std::endl;
+    for(size_t N = 0; N < 10; ++N) {
+        std::cout << "    " << N << "! = " << factorial(N) << std::endl;
+    }
+
+    std::cout << "Binomial:" << std::endl;
+    for(size_t N = 1; N < 10; ++N) {
+        std::cout << std::setw(5) << "N=" << N ;
+        for(size_t k = 0; k <= N; ++k) {
+            std::cout << "    " << std::setw(5) << binomial(N, k);
+        }
+        std::cout << std::endl;
+    }
+
+    std::cout << "Limits:" << std::endl;
+    std::cout << "    20! = " << factorial(20) << std::endl;
+    std::cout << "    21! = " << factorial(21) << std::endl;
+    std::cout << "    70! = " << factorial(70) << std::endl;
+
+    std::cout << "Limits:" << std::endl;
+    std::cout << "    N=20, k=1 ==> " << binomial(20,1) << std::endl;
+    std::cout << "    N=21, k=1 ==> " << binomial(21,1) << " <- wrong " << std::endl;
+    // will trigger floating point exception
+    // std::cout << "    N=70, k=1 ==> " << binomial(70,1) << std::endl;
+
+    return 0;
+}

exercises/ex11/solutions/harmonic_chain/CMakeLists.txt

+cmake_minimum_required (VERSION 3.15)
+project (ex11-harmonic-oscillator)
+
+set (CMAKE_CXX_STANDARD 17)
+set (CMAKE_CXX_STANDARD_REQUIRED TRUE)
+set (CMAKE_CXX_EXTENSIONS FALSE)
+
+add_compile_options (-Wall -Wextra -Wpedantic -march=native)
+
+find_package (LAPACK REQUIRED)
+message(INFO "${LAPACK_LINKER_FLAGS} ${LAPACK_LIBRARIES}")
+
+add_executable (harmonic_chain harmonic_chain.cpp)
+target_link_options (harmonic_chain PRIVATE ${LAPACK_LINKER_FLAGS})
+target_link_libraries (harmonic_chain ${LAPACK_LIBRARIES})

exercises/ex11/solutions/harmonic_chain/harmonic_chain.cpp

+/*
+ * Programming Techniques for Scientific Simulations I
+ * HS 2020
+ * Week 11
+ */
+
+#include <cmath>
+#include <stdexcept>
+#include <vector>
+#include <iostream>
+
+extern "C" void dsyev_(
+	char   const & JOBZ,
+	char   const & UPLO,
+	int    const & N,
+	double       * A,
+	int    const & LDA,
+	double       * W,
+	double       * WORK,
+	int    const & LWORK,
+	int          & INFO
+);
+
+std::vector<double> hamiltonian(double const K, double const m, std::size_t const N) {
+	std::vector<double> M(N * N);
+        //Note that the matrix is symmetric, we don't need to worry about majors
+        //But usually, Fortran expects column major ordering
+	for (std::size_t i = 0; i < N; ++i) {
+		M[i + i * N] =  2 * K / m; //diagonal
+		if (i != 0) {
+			M[(i - 1) + i * N] = -1 * K / m; //to the left
+		}
+		if (i != N - 1) {
+			M[(i + 1) + i * N] = -1 * K / m; //to the right
+		}
+	}
+
+	return M;
+}
+
+std::vector<double> solve(std::vector<double> & M, std::size_t const N) {
+	std::vector<double> omega(N);
+	int info;
+	double dwork;
+
+	dsyev_('V', 'L', N, M.data(), N, omega.data(), &dwork, -1, info);
+	if (info != 0) {
+		throw std::runtime_error("Something went wrong!");
+	}
+	int lwork = static_cast<int>(dwork);
+
+	double * const work = new double[lwork];
+	dsyev_('V', 'L', N, M.data(), N, omega.data(), work, lwork, info);
+	delete[] work;
+	if (info != 0) {
+		throw std::runtime_error("Something went wrong!");
+	}
+
+	for (auto & w : omega) {
+		w = std::sqrt(w);
+	}
+	
+	return omega;
+}
+
+int main() {
+	double const K = 1;
+	double const m = 1;
+	std::size_t const N = 16;
+
+	auto M = hamiltonian(K, m, N);
+	auto omega = solve(M, N);
+
+        std::cout << "Computed eigenvalues\n";
+	for (auto && w : omega) {
+		std::cout << w << "\t ";
+	}
+	std::cout << '\n';
+
+        std::cout << "Exact eigenvalues\n";
+	for (std::size_t i = 1; i <= N; ++i) {
+		std::cout <<  std::sqrt(K / m * (2. - 2. * std::cos(M_PI * static_cast<double>(i) / (N + 1))))<< "\t ";
+	}
+	std::cout << '\n';
+
+        std::cout << "Eigenvectors (stored column-wise)\n";
+        for(std::size_t i = 0; i < N; ++i){
+                for(std::size_t j = 0; j < N; ++j){
+                        std::cout << M[j + N*i] << "\t ";
+                }
+                std::cout << std::endl;
+        }
+        std::cout << '\n';
+}