diff --git a/exercises/ex12_skeleton/Makefile b/exercises/ex12_skeleton/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..e2200a6231aaa278b0cf17c254d784414983d561
--- /dev/null
+++ b/exercises/ex12_skeleton/Makefile
@@ -0,0 +1,20 @@
+CXX_FLAGS ?= -DNDEBUG -O2
+
+all: sim_thermal_annealing sim_quantum_annealing
+
+run_quantum:
+	./sim_quantum_annealing 2> succ_rates-quantum.txt
+	./plot_succ_rates.py
+
+run_thermal:
+	./sim_thermal_annealing > output-thermal.txt
+
+sim_thermal_annealing: sim_thermal_annealer.hpp sim_thermal.cpp
+	$(CXX) $(CXX_FLAGS) sim_thermal.cpp -o sim_thermal_annealing
+
+sim_quantum_annealing: sim_quantum_annealer.hpp sim_quantum.cpp
+	$(CXX) $(CXX_FLAGS) sim_quantum.cpp -o sim_quantum_annealing
+
+clean:
+	-rm sim_thermal_annealing
+	-rm sim_quantum_annealing
diff --git a/exercises/ex12_skeleton/instance.txt b/exercises/ex12_skeleton/instance.txt
new file mode 100644
index 0000000000000000000000000000000000000000..db9957158fce43b815483718eb532fbf1d40d301
--- /dev/null
+++ b/exercises/ex12_skeleton/instance.txt
@@ -0,0 +1 @@
+-136 10 1 1 -1 1 1 -1 -1 1 1 1 -1 1 1 -1 1 -1 -1 1 1 1 1 -1 1 1 1 1 1 -1 1 -1 1 -1 -1 -1 -1 1 1 -1 1 -1 -1 -1 1 1 -1 -1 -1 1 1 1 -1 1 1 -1 -1 -1 1 -1 1 -1 1 -1 1 1 1 1 1 -1 -1 -1 1 -1 1 -1 1 -1 -1 -1 1 -1 -1 1 1 1 1 -1 1 1 -1 1 -1 -1 1 1 1 -1 1 -1 -1 -1 -1 1 1 -1 -1 1 -1 1 1 1 1 1 -1 1 -1 -1 -1 1 1 1 1 -1 -1 1 -1 1 -1 -1 -1 -1 1 -1 1 -1 1 1 -1 1 -1 1 -1 -1 -1 -1 -1 -1 -1 1 1 1 1 -1 1 1 -1 1 1 1 1 -1 1 1 1 -1 -1 1 -1 1 -1 1 1 -1 -1 -1 -1 1 1 1 1 1 -1 1 1 1 1 1 -1 1 -1 1 -1 -1 1 -1 1 -1 -1 -1 1 -1 
diff --git a/exercises/ex12_skeleton/sim_quantum.cpp b/exercises/ex12_skeleton/sim_quantum.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c0c81e92a7064292e5ac440734561fc1e846787a
--- /dev/null
+++ b/exercises/ex12_skeleton/sim_quantum.cpp
@@ -0,0 +1,50 @@
+#include <iostream>
+#include <fstream>
+#include <stdexcept>
+#include "sim_quantum_annealer.hpp"
+
+
+void load(std::string const& filename, std::size_t& N, int& exact_energy, std::vector<short>& row_couplings, std::vector<short>& col_couplings) {
+    std::ifstream in(filename.c_str());
+    in >> exact_energy;
+    in >> N;
+    row_couplings.resize(N*N);
+    col_couplings.resize(N*N);
+    for(std::size_t i=0; i < N*N; ++i)
+        in >> row_couplings[i];
+    for(std::size_t i=0; i < N*N; ++i)
+        in >> col_couplings[i];
+}
+
+int main() {
+    int const copies = 1000;
+    int exact_energy = 0.0;
+    std::size_t N = 10;
+    std::vector<short> row_couplings;
+    std::vector<short> col_couplings;
+    load("instance.txt", N, exact_energy, row_couplings, col_couplings);
+
+
+    for(double beta=4e-2; beta <= 11; beta *= 2)
+    for(double dt=1e-5; dt <= 1e-1 ; dt *= 10) {
+        int minimum = 0;
+        int gs_hits = 0;
+        for(int i = 0 ; i < copies; ++i) {
+            cqp::quantum_annealer annealer(N,row_couplings, col_couplings, 10, 23*i);
+            annealer.run(beta, dt); {
+                std::cout <<"Run " << i <<": E = " <<annealer.energy() << std::endl;
+                minimum = (std::min)(minimum,annealer.energy());
+                if(annealer.energy() == exact_energy)
+                    ++gs_hits;
+            }
+        }
+
+        double rate = static_cast<double>(gs_hits) / copies;
+        double err  = std::sqrt((static_cast<double>(gs_hits) / copies - rate*rate)/(copies-1));
+        std::cout<<"exact         E = "<< exact_energy <<std::endl;
+        std::cout<<"annealing min E = "<< minimum << std::endl;
+        std::cout<<"%:                "<< rate << "+/-" << err << std::endl;
+        std::cerr<<"beta = " << beta << " dt = " << dt << " rate = "<< rate << " +/- " << err<< std::endl;
+    }
+    return 0;
+}
diff --git a/exercises/ex12_skeleton/sim_quantum_annealer.hpp b/exercises/ex12_skeleton/sim_quantum_annealer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..48e0b29e199194d7bf3ddb8535a8d956a5356e4b
--- /dev/null
+++ b/exercises/ex12_skeleton/sim_quantum_annealer.hpp
@@ -0,0 +1,85 @@
+#ifndef SIM_QUANTUM_ANNEALER_HPP
+#define SIM_QUANTUM_ANNEALER_HPP
+#include <vector>
+#include <limits>
+#include <random>
+#include <assert.h>
+
+
+namespace cqp {
+
+class quantum_annealer {
+  public:
+    typedef int energy_t;
+
+    quantum_annealer(std::size_t N, std::vector<short> const& row_c, std::vector<short> const& col_c, std::size_t M, unsigned int seed = 3000)
+    : rng(seed), random_site(0,N-1), random_time_slice(0,M-1), row_c(row_c), col_c(col_c), spin_rows(N*M), N(N), M(M) {
+        // We also assume row_c and col_c have only +1,-1 elements
+        assert(N < 64);
+        init();
+    }
+
+    void init() {
+        // Generate random spin config
+        for(std::size_t i = 0; i < spin_rows.size(); ++i) {
+            spin_rows[i] = rng();
+            spin_rows[i] <<= 32;
+            spin_rows[i] |= rng();
+        }
+    }
+
+    // +1 if spins are antialigned, -1 if spins are aligned
+    inline int64_t get_sign( /* TODO */ ) const {
+
+
+        /* TODO */
+        return 0;
+
+
+    }
+
+    void sweep(double beta, double alpha, double field) {
+
+        double delta = beta / M;
+        double j_t = std::numeric_limits<double>::max(); // limit of the coupling strength as gamma ->0 (all the replicas are perfectly coupled)
+
+
+        /* TODO */
+
+
+    }
+
+    energy_t energy() const {
+        energy_t e = 0;
+
+
+        /* TODO */
+
+
+        return e;
+    }
+
+    void run(double beta, double delta_t) {
+        init();
+        for(double t=0; t <= 1.0; t += delta_t) {
+
+
+            /* TODO */
+
+
+        }
+    }
+  private:
+    std::mt19937 rng;
+    std::uniform_real_distribution<> rnd;
+    std::uniform_int_distribution<> random_site;
+    std::uniform_int_distribution<> random_time_slice;
+    std::vector<short> const row_c;
+    std::vector<short> const col_c;
+    std::vector<uint64_t> spin_rows;
+    std::size_t const N;
+    std::size_t const M;
+};
+
+}
+#endif //SIM_QUANTUM_ANNEALER_HPP
diff --git a/exercises/ex12_skeleton/sim_thermal.cpp b/exercises/ex12_skeleton/sim_thermal.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..3ee9785011cba14ab4665e1792bc8ce8fec361db
--- /dev/null
+++ b/exercises/ex12_skeleton/sim_thermal.cpp
@@ -0,0 +1,45 @@
+#include <iostream>
+#include <fstream>
+#include <stdexcept>
+#include "sim_thermal_annealer.hpp"
+
+
+void load(std::string const& filename, std::size_t& N, int& exact_energy, std::vector<short>& row_couplings, std::vector<short>& col_couplings) {
+    std::ifstream in(filename.c_str());
+    in >> exact_energy;
+    in >> N;
+    row_couplings.resize(N*N);
+    col_couplings.resize(N*N);
+    for(std::size_t i=0; i < N*N; ++i)
+        in >> row_couplings[i];
+    for(std::size_t i=0; i < N*N; ++i)
+        in >> col_couplings[i];
+}
+
+int main() {
+    int minimum = 0;
+    int const copies = 100;
+    int gs_hits = 0;
+
+
+    int exact_energy = 0.0;
+    std::size_t N = 10;
+    std::vector<short> row_couplings;
+    std::vector<short> col_couplings;
+    load("instance.txt", N, exact_energy, row_couplings, col_couplings);
+
+    cqp::thermal_annealer annealer(N,row_couplings, col_couplings);
+
+    for(int i = 0 ; i < copies; ++i) {
+        annealer.run(1e-5, 0.1, 5.0);
+        std::cout <<"Run " << i <<": E = " <<annealer.energy() << std::endl;
+        minimum = (std::min)(minimum,annealer.energy());
+        if(annealer.energy() == exact_energy)
+            ++gs_hits;
+    }
+
+    std::cout<<"exact         E = "<< exact_energy <<std::endl;
+    std::cout<<"annealing min E = "<< minimum << std::endl;
+    std::cout<<"%:                "<< static_cast<double>(gs_hits)/copies << std::endl;
+    return 0;
+}
diff --git a/exercises/ex12_skeleton/sim_thermal_annealer b/exercises/ex12_skeleton/sim_thermal_annealer
new file mode 100644
index 0000000000000000000000000000000000000000..7a6a60d2b1ba6a1bf28803e23327713345af5720
Binary files /dev/null and b/exercises/ex12_skeleton/sim_thermal_annealer differ
diff --git a/exercises/ex12_skeleton/sim_thermal_annealer.hpp b/exercises/ex12_skeleton/sim_thermal_annealer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..22dd3447a9d14983181f70305aadae0ec9cd7b1c
--- /dev/null
+++ b/exercises/ex12_skeleton/sim_thermal_annealer.hpp
@@ -0,0 +1,72 @@
+#ifndef SIM_THERMAL_ANNEALER_HPP
+#define SIM_THERMAL_ANNEALER_HPP
+#include <vector>
+#include <random>
+#include <assert.h>
+
+namespace cqp {
+
+class thermal_annealer {
+  public:
+    typedef int energy_t;
+
+    thermal_annealer(std::size_t N, std::vector<short> const& row_c, std::vector<short> const& col_c, unsigned int seed = 5000)
+    : rng(seed), random_site(0,N-1), row_c(row_c), col_c(col_c), spin_rows(N), N(N) {
+        // We also assume row_c and col_c have only +1,-1 elements
+        assert(N < 64);
+        init();
+    }
+
+    void init() {
+        // Generate random spin config
+        for(std::size_t i = 0; i < N; ++i) {
+            spin_rows[i] = rng();
+            spin_rows[i] <<= 32;
+            spin_rows[i] |= rng();
+        }
+    }
+
+    inline int64_t get_sign( /* TODO */ ) const {
+
+
+        /* TODO */
+        return 0;
+
+
+    }
+
+    void sweep(double beta) {
+
+
+        /* TODO */
+
+
+    }
+
+    energy_t energy() const {
+        energy_t e = 0;
+
+
+        /* TODO */
+
+
+        return e;
+    }
+
+    void run(double delta_beta, double beta = 1e-5, double beta_max = 5.0) {
+        init();
+        for(; beta < beta_max; beta += delta_beta)
+            sweep(beta);
+    }
+  private:
+    std::mt19937 rng;
+    std::uniform_real_distribution<> rnd;
+    std::uniform_int_distribution<> random_site;
+    std::vector<short> const row_c;
+    std::vector<short> const col_c;
+    std::vector<uint64_t> spin_rows;
+    std::size_t const N;
+};
+
+}
+#endif //SIM_THERMAL_ANNEALER_HPP
diff --git a/exercises/exercise12.pdf b/exercises/exercise12.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..b5e2d75fcce902c42cd2ac5a9bc716f878da0344
Binary files /dev/null and b/exercises/exercise12.pdf differ