Merge branch 'master' of m5:/srv/git/wolff

20 files changed, 1088 insertions, 630 deletions

diff --git a/lib/cluster.c b/lib/cluster.c
deleted file mode 100644
index 7274eb9..0000000
--- a/lib/cluster.c
+++ /dev/null
@@ -1,368 +0,0 @@
-
-#include "cluster.h"
-
-v_t flip_cluster(ising_state_t *s, v_t v0, q_t rot, gsl_rng *r) {
-  v_t nv = 0;
-
-  ll_t *stack = NULL;     // create a new stack
-  stack_push(&stack, v0); // push the initial vertex to the stack
-
-  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
-
-  while (stack != NULL) {
-    v_t v = stack_pop(&stack);
-
-    if (!marks[v]) {
-      q_t s_old, s_new;
-      dihedral_t *R_new; 
-      bool external_flipped;
-
-      marks[v] = true;
-
-      if (v == s->g->nv - 1) {
-        R_new = dihedral_compose(s->q, rot, s->R);
-        external_flipped = true;
-      } else {
-        s_old = s->spins[v];
-        s_new = dihedral_act(s->q, rot, s_old);
-        external_flipped = false;
-      }
-
-      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
-
-      for (v_t i = 0; i < nn; i++) {
-        q_t sn;
-        double prob;
-        bool external_neighbor = false;
-
-        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
-
-        if (vn == s->g->nv - 1) {
-          external_neighbor = true;
-        } else {
-          sn = s->spins[vn];
-        }
-
-        if (external_flipped || external_neighbor) {
-          q_t rot_s_old, rot_s_new;
-
-          if (external_neighbor) {
-            rot_s_old = dihedral_inverse_act(s->q, s->R, s_old);
-            rot_s_new = dihedral_inverse_act(s->q, s->R, s_new);
-          } else {
-            rot_s_old = dihedral_inverse_act(s->q, s->R, sn);
-            rot_s_new = dihedral_inverse_act(s->q, R_new, sn);
-          }
-
-          prob = s->H_probs[rot_s_new * s->q + rot_s_old];
-
-          s->M[rot_s_old]--;
-          s->M[rot_s_new]++;
-
-          s->E += - s->H[rot_s_new] + s->H[rot_s_old];
-        } else {
-          q_t diff_old = (s_old + s->q - sn) % s->q;
-          q_t diff_new = (s_new + s->q - sn) % s->q;
-
-          prob = s->J_probs[diff_new * s->q + diff_old];
-
-          s->E += - s->J[diff_new] + s->J[diff_old];
-        }
-
-        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
-          stack_push(&stack, vn); // push the neighboring vertex to the stack
-        }
-      }
-
-      if (external_flipped) {
-        free(s->R);
-        s->R = R_new;
-      } else {
-        s->spins[v] = s_new;
-      }
-
-      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
-        nv++;
-      }
-    }
-  }
-
-  free(marks);
-
-  return nv;
-}
-
-v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t rot, gsl_rng *r) {
-  v_t nv = 0;
-
-  ll_t *stack = NULL;     // create a new stack
-  stack_push(&stack, v0); // push the initial vertex to the stack
-
-  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
-
-  while (stack != NULL) {
-    v_t v = stack_pop(&stack);
-
-    if (!marks[v]) {
-      h_t s_old, s_new;
-      dihinf_t *R_new; 
-      bool external_flipped;
-
-      marks[v] = true;
-
-      if (v == s->g->nv - 1) {
-        R_new = dihinf_compose(rot, s->R);
-        external_flipped = true;
-      } else {
-        s_old = s->spins[v];
-        s_new = dihinf_act(rot, s_old);
-        external_flipped = false;
-      }
-
-      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
-
-      for (v_t i = 0; i < nn; i++) {
-        h_t sn;
-        double prob;
-        bool external_neighbor = false;
-
-        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
-
-        if (vn == s->g->nv - 1) {
-          external_neighbor = true;
-        } else {
-          sn = s->spins[vn];
-        }
-
-        if (external_flipped || external_neighbor) {
-          h_t rot_s_old, rot_s_new;
-
-          if (external_neighbor) {
-            rot_s_old = dihinf_inverse_act(s->R, s_old);
-            rot_s_new = dihinf_inverse_act(s->R, s_new);
-          } else {
-            rot_s_old = dihinf_inverse_act(s->R, sn);
-            rot_s_new = dihinf_inverse_act(R_new, sn);
-          }
-
-          double dE = s->H(s->H_info, rot_s_old) - s->H(s->H_info, rot_s_new);
-          prob = 1.0 - exp(-dE / s->T);
-
-          s->M += rot_s_new - rot_s_old;
-          s->E += dE;
-        } else {
-          double dE = (s->J)(s_old - sn) - (s->J)(s_new - sn);
-          prob = 1.0 - exp(-dE / s->T);
-          s->E += dE;
-        }
-
-        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
-          stack_push(&stack, vn); // push the neighboring vertex to the stack
-        }
-      }
-
-      if (external_flipped) {
-        free(s->R);
-        s->R = R_new;
-      } else {
-        s->spins[v] = s_new;
-      }
-
-      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
-        nv++;
-      }
-    }
-  }
-
-  free(marks);
-
-  return nv;
-}
-
-v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r) {
-  v_t nv = 0;
-
-  ll_t *stack = NULL;     // create a new stack
-  stack_push(&stack, v0); // push the initial vertex to the stack
-
-  //node_t *T = NULL;
-  bool *marks = (bool *)calloc(s->g->nv, sizeof(bool));
-
-  while (stack != NULL) {
-    v_t v = stack_pop(&stack);
-
-//    if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
-    if (!marks[v]) {
-      bool v_is_external = false;
-      double *s_old, *s_new, *R_tmp; 
-
-      if (v == s->g->nv - 1) {
-        v_is_external = true;
-      }
-
-      //tree_insert(&T, v);
-      marks[v] = true;
-
-      if (v == s->g->nv - 1) {
-        R_tmp = orthogonal_rotate(s->n, rot, s->R);
-      } else {
-        s_old = &(s->spins[s->n * v]); // don't free me! I'm a pointer within array s->spins
-        s_new = vector_rotate(s->n, rot, s_old); // free me! I'm a new vector
-      }
-
-      v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
-
-      for (v_t i = 0; i < nn; i++) {
-        v_t vn = s->g->v_adj[s->g->v_i[v] + i];
-
-        bool vn_is_external = false;
-
-        if (vn == s->g->nv - 1) {
-          vn_is_external = true;
-        }
-
-        double *sn;
-
-        if (!vn_is_external) {
-          sn = &(s->spins[s->n * vn]);
-        }
-
-        double prob;
-
-        if (v_is_external || vn_is_external) {
-          double *rs_old, *rs_new;
-          if (vn_is_external) {
-            rs_old = vector_rotate_inverse(s->n, s->R, s_old);
-            rs_new = vector_rotate_inverse(s->n, s->R, s_new);
-          } else {
-            rs_old = vector_rotate_inverse(s->n, s->R, sn);
-            rs_new = vector_rotate_inverse(s->n, R_tmp, sn);
-          }
-          double dE = s->H(s->n, s->H_info, rs_old) - s->H(s->n, s->H_info, rs_new);
-          prob = 1.0 - exp(-dE / s->T);
-          vector_subtract(s->n, s->M, rs_old);
-          vector_add(s->n, s->M, rs_new);
-          s->E += dE;
-
-          free(rs_old);
-          free(rs_new);
-        } else {
-          double dE = (s->J)(vector_dot(s->n, sn, s_old)) - (s->J)(vector_dot(s->n, sn, s_new));
-          prob = 1.0 - exp(-dE / s->T);
-          s->E += dE;
-        }
-
-        if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
-          stack_push(&stack, vn); // push the neighboring vertex to the stack
-        }
-      }
-
-      if (v == s->g->nv - 1) {
-        free(s->R);
-        s->R = R_tmp;
-      } else {
-        vector_replace(s->n, s_old, s_new);
-        free(s_new);
-      }
-
-      if (v != s->g->nv - 1) { // count the number of non-external sites that flip
-        nv++;
-      }
-    }
-  }
-
-  //tree_freeNode(T);
-  free(marks);
-
-  return nv;
-}
-
-/*G
-template <class R_t, class X_t>
-v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t *r, gsl_rng *rand) {
-  v_t nv = 0;
-
-  ll_t *stack = NULL;     // create a new stack
-  stack_push(&stack, v0); // push the initial vertex to the stack
-
-  bool *marks = (bool *)calloc(state->g->nv, sizeof(bool));
-
-  while (stack != NULL) {
-    v_t v = stack_pop(&stack);
-
-    if (!marks[v]) {
-      X_t *si_old, *si_new;
-      R_t *s0_old, *s0_new;
-
-      si_old = state->s[v];
-      s0_old = state->s0;
-
-      marks[v] = true;
-
-      if (v == state->g->nv - 1) {
-        s0_new = act <R_t, R_t> (r, s0_old);
-      } else {
-        si_new = act <R_t, X_t> (r, si_old);
-      }
-
-      v_t nn = state->g->v_i[v + 1] - state->g->v_i[v];
-
-      for (v_t i = 0; i < nn; i++) {
-        v_t vn = state->g->v_adj[state->g->v_i[v] + i];
-
-        X_t *sj;
-
-        if (vn != state->g->nv - 1) {
-          sj = state->s[vn];
-        }
-
-        double prob;
-
-        bool is_ext = (v == state->g->nv - 1 || vn == state->g->nv - 1);
-
-        if (is_ext) {
-          X_t *rs_old, *rs_new;
-          if (vn == state->g->nv - 1) {
-            rs_old = inverse_act <class R_t, class X_t> (s0_old, si_old);
-            rs_new = inverse_act <class R_t, class X_t> (s0_old, si_new);
-          } else {
-            rs_old = inverse_act <class R_t, class X_t> (s0_old, sj);
-            rs_new = inverse_act <class R_t, class X_t> (s0_new, sj);
-          }
-          double dE = state->B(rs_old) - state->B(rs_new);
-          prob = 1.0 - exp(-dE / state->T);
-          update_magnetization <X_t> (state->M, rs_old, rs_new);
-          state->E += dE;
-
-          free_X <X_t> (rs_old);
-          free_X <X_t> (rs_new);
-        } else {
-          double dE = state->Z(si_old, sj) - state->Z(si_new, sj);
-          prob = 1.0 - exp(-dE / state->T);
-          state->E += dE;
-        }
-
-        if (gsl_rng_uniform(rand) < prob) { // and with probability...
-          stack_push(&stack, vn); // push the neighboring vertex to the stack
-        }
-      }
-
-      if (v == state->g->nv - 1) {
-        free_R <R_t> (state->s0);
-        state->s0 = s0_new;
-      } else {
-        free_X <X_t> (state->s[v]);
-        state->s[v] = si_new;
-      }
-
-      if (v != state->g->nv - 1) { // count the number of non-external sites that flip
-        nv++;
-      }
-    }
-  }
-
-  free(marks);
-
-  return nv;
-}
-*/
diff --git a/lib/cluster.h b/lib/cluster.h
index d118735..a20912e 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -1,13 +1,14 @@
 
 #pragma once
 
+#include <functional>
 #include <assert.h>
 #include <fftw3.h>
 #include <float.h>
 #include <gsl/gsl_randist.h>
 #include <gsl/gsl_rng.h>
 #include <inttypes.h>
-#include <math.h>
+#include <cmath>
 #include <stdbool.h>
 #include <string.h>
 #include <sys/types.h>
@@ -19,62 +20,169 @@
 #include "graph.h"
 #include "tree.h"
 #include "measurement.h"
+#include "vector.h"
 #include "orthogonal.h"
 #include "dihedral.h"
 #include "dihinf.h"
 #include "yule_walker.h"
 
-typedef struct {
-  graph_t *g;
-  q_t *spins;
-  double T;
-  double *J;
-  double *H;
-  double *J_probs;
-  double *H_probs;
-  dihedral_t *R;
-  double E;
-  v_t *M;
-  q_t q;
-} ising_state_t;
-
-typedef struct {
-  graph_t *g;
-  h_t *spins;
-  double T;
-  double (*J)(h_t);
-  double (*H)(double *, h_t);
-  double *H_info;
-  dihinf_t *R;
-  double E;
-  h_t M;
-} dgm_state_t;
-
-typedef struct {
-  graph_t *g;
-  double *spins;
-  double T;
-  double (*J)(double);
-  double (*H)(q_t, double *, double *);
-  double *H_info;
-  double *R;
-  double E;
-  double *M;
-  q_t n;
-} vector_state_t;
-
-typedef enum {
-  VECTOR,
-  MODULATED,
-  CUBIC,
-  QUADRATIC
-} vector_field_t;
-
-v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r);
-
-v_t flip_cluster_vector(vector_state_t *s, v_t v0, double *rot, gsl_rng *r);
-
-v_t flip_cluster_dgm(dgm_state_t *s, v_t v0, h_t rot, gsl_rng *r);
-
-graph_t *graph_add_ext(const graph_t *g);
+template <class T>
+void init(T*);
+
+template <class T>
+T scalar_multiple(v_t a, T b);
+
+template <class R_t, class X_t>
+X_t act(R_t a, X_t b);
+
+template <class R_t, class X_t>
+X_t act_inverse(R_t a, X_t b);
+
+template <class T>
+T copy(T a);
+
+template <class T>
+void free_spin(T a);
+
+template <class T>
+T add(T, T);
+
+template <class T>
+T subtract(T, T);
+
+template <class T>
+T gen_rot(gsl_rng *r);
+
+template <class R_t, class X_t>
+class state_t {
+  public:
+    D_t D;
+    L_t L;
+    v_t nv;
+    v_t ne;
+    graph_t *g;
+    double T;
+    X_t *spins;
+    R_t R;
+    double E;
+    X_t M; // the "sum" of the spins, like the total magnetization
+
+    std::function <double(X_t, X_t)> J;
+    std::function <double(X_t)> H;
+
+    state_t(D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H) : D(D), L(L), T(T), J(J), H(H) {
+      graph_t *h = graph_create_square(D, L);
+      nv = h->nv;
+      ne = h->ne;
+      g = graph_add_ext(h);
+      graph_free(h);
+      spins = (X_t *)malloc(nv * sizeof(X_t));
+      for (v_t i = 0; i < nv; i++) {
+        init (&(spins[i]));
+      }
+      init (&R);
+      E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
+      M = scalar_multiple (nv, spins[0]);
+    }
+
+    ~state_t() {
+      graph_free(g);
+      for (v_t i = 0; i < nv; i++) {
+        free_spin(spins[i]);
+      }
+      free(spins);
+      free_spin(R);
+      free_spin(M);
+    }
+};
+
+template <class R_t, class X_t>
+v_t flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
+  v_t nv = 0;
+
+  ll_t *stack = NULL;     // create a new stack
+  stack_push(&stack, v0); // push the initial vertex to the stack
+
+  bool *marks = (bool *)calloc(state->g->nv, sizeof(bool));
+
+  while (stack != NULL) {
+    v_t v = stack_pop(&stack);
+
+    if (!marks[v]) {
+      X_t si_old, si_new;
+      R_t R_old, R_new;
+
+      si_old = state->spins[v];
+      R_old = state->R;
+
+      marks[v] = true;
+
+      if (v == state->g->nv - 1) {
+        R_new = act (r, R_old);
+      } else {
+        si_new = act (r, si_old);
+      }
+
+      v_t nn = state->g->v_i[v + 1] - state->g->v_i[v];
+
+      for (v_t i = 0; i < nn; i++) {
+        v_t vn = state->g->v_adj[state->g->v_i[v] + i];
+
+        X_t sj;
+
+        if (vn != state->g->nv - 1) {
+          sj = state->spins[vn];
+        }
+
+        double prob;
+
+        bool is_ext = (v == state->g->nv - 1 || vn == state->g->nv - 1);
+
+        if (is_ext) {
+          X_t rs_old, rs_new;
+          if (vn == state->g->nv - 1) {
+            rs_old = act_inverse (R_old, si_old);
+            rs_new = act_inverse (R_old, si_new);
+          } else {
+            rs_old = act_inverse (R_old, sj);
+            rs_new = act_inverse (R_new, sj);
+          }
+          double dE = state->H(rs_old) - state->H(rs_new);
+          prob = 1.0 - exp(-dE / state->T);
+
+          subtract (state->M, rs_old);
+          add (state->M, rs_new);
+          state->E += dE;
+
+          free_spin (rs_old);
+          free_spin (rs_new);
+        } else {
+          double dE = state->J(si_old, sj) - state->J(si_new, sj);
+          prob = 1.0 - exp(-dE / state->T);
+          state->E += dE;
+        }
+
+        if (gsl_rng_uniform(rand) < prob) { // and with probability...
+          stack_push(&stack, vn); // push the neighboring vertex to the stack
+        }
+      }
+
+      if (v == state->g->nv - 1) {
+        free_spin(state->R);
+        state->R = R_new;
+      } else {
+        free_spin(state->spins[v]);
+        state->spins[v] = si_new;
+      }
+
+      if (v != state->g->nv - 1) { // count the number of non-external sites that flip
+        nv++;
+      }
+    }
+  }
+
+  free(marks);
+
+  return nv;
+}
 
diff --git a/lib/cluster_finite.c b/lib/cluster_finite.c
index f11a3ea..9392cf8 100644
--- a/lib/cluster_finite.c
+++ b/lib/cluster_finite.c
@@ -1,8 +1,8 @@
 
 #include "cluster_finite.h"
 
-v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) {
-  q_t *rot = s->transformations + s->q * rot_ind;
+v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot_ind, gsl_rng *r) {
+  q_t *rot = s->transformations + s->q * s->involutions[rot_ind];
   q_t *R_inv = symmetric_invert(s->q, s->R);
   v_t nv = 0;
 
@@ -62,14 +62,14 @@ v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot_ind, gsl_rng *r) {
           s->M[rot_s_old]--;
           s->M[rot_s_new]++;
 
-          s->E += - s->H[rot_s_new] + s->H[rot_s_old];
         } else {
-          q_t diff_old = (s_old + s->q - sn) % s->q;
-          q_t diff_new = (s_new + s->q - sn) % s->q;
+          q_t diff_old = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_old]];
+          q_t diff_new = s->bond_with_zero_type[s->transformations[s->q * s->transform_site_to_zero[sn] + s_new]];
 
-          prob = s->J_probs[diff_new * s->q + diff_old];
+          prob = s->J_probs[diff_new * s->n_bond_types + diff_old];
 
-          s->E += - s->J[diff_new] + s->J[diff_old];
+          s->B[diff_old]--;
+          s->B[diff_new]++;
         }
 
         if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
diff --git a/lib/cluster_finite.h b/lib/cluster_finite.h
index abdc8fc..b2d764e 100644
--- a/lib/cluster_finite.h
+++ b/lib/cluster_finite.h
@@ -23,10 +23,19 @@
 #include "yule_walker.h"
 
 typedef struct {
+  D_t D;
+  L_t L;
+  v_t nv;
+  v_t ne;
   graph_t *g;
   q_t q;
   R_t n_transformations;
   q_t *transformations;
+  R_t n_involutions;
+  R_t *involutions;
+  R_t *transform_site_to_zero;
+  q_t n_bond_types;
+  q_t *bond_with_zero_type;
   double T;
   double *J;
   double *H;
@@ -34,9 +43,9 @@ typedef struct {
   double *H_probs;
   q_t *spins;
   q_t *R;
-  double E;
+  v_t *B;
   v_t *M;
 } state_finite_t;
 
-v_t flip_cluster_finite(state_finite_t *s, v_t v0, q_t rot, gsl_rng *r);
+v_t flip_cluster_finite(state_finite_t *s, v_t v0, R_t rot, gsl_rng *r);
 
diff --git a/lib/dihedral.c b/lib/dihedral.c
index ac74a23..8158b43 100644
--- a/lib/dihedral.c
+++ b/lib/dihedral.c
@@ -11,10 +11,14 @@ dihedral_t *dihedral_compose(q_t q, q_t g1i, const dihedral_t *g2) {
   return g3;
 }
 
-q_t dihedral_act(q_t q, q_t gi, q_t s) {
+q_t dihedral_act(q_t q, q_t gi, bool r, q_t s) {
   // we only need to consider the action of reflections
 
-  return (gi + q - s) % q;
+  if (r) {
+    return (gi + q - s) % q;
+  } else {
+    return (gi + s) % q;
+  }
 }
 
 q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) {
@@ -26,15 +30,26 @@ q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s) {
 }
 
 q_t *dihedral_gen_transformations(q_t q) {
-  q_t *transformations = (q_t *)malloc(q * q * sizeof(q_t));
+  q_t *transformations = (q_t *)malloc(2 * q * q * sizeof(q_t));
 
   for (q_t i = 0; i < q; i++) {
     for (q_t j = 0; j < q; j++) {
-      transformations[q * i + j] = dihedral_act(q, i, j);
+      transformations[q * i + j] = dihedral_act(q, i, false, j);
+      transformations[q * q + q * i + j] = dihedral_act(q, i, true, j);
     }
   }
 
   return transformations;
 }
 
+R_t *dihedral_gen_involutions(q_t q) {
+  R_t *transformations = (R_t *)malloc(q * sizeof(R_t));
+
+  for (q_t i = 0; i < q; i++) {
+    transformations[i] = q + i;
+  }
+
+  return transformations;
+}
+
 
diff --git a/lib/dihedral.h b/lib/dihedral.h
index e5e4cbd..c95b23a 100644
--- a/lib/dihedral.h
+++ b/lib/dihedral.h
@@ -11,9 +11,11 @@ typedef struct {
 
 dihedral_t *dihedral_compose(q_t q, q_t gti, const dihedral_t *g2);
 
-q_t dihedral_act(q_t q, q_t gi, q_t s);
+q_t dihedral_act(q_t q, q_t gi, bool r, q_t s);
 
 q_t dihedral_inverse_act(q_t q, const dihedral_t *g, q_t s);
 
 q_t *dihedral_gen_transformations(q_t q);
+R_t *dihedral_gen_involutions(q_t q);
 
+R_t factorial(q_t);
diff --git a/lib/graph.h b/lib/graph.h
index 9c80dd6..beb7f4c 100644
--- a/lib/graph.h
+++ b/lib/graph.h
@@ -1,10 +1,16 @@
 
+#pragma once
+
 #include <inttypes.h>
 #include <math.h>
 #include <stdlib.h>
 
 #include "types.h"
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 typedef struct {
   v_t ne;
   v_t nv;
@@ -13,8 +19,10 @@ typedef struct {
 } graph_t;
 
 graph_t *graph_create_square(D_t D, L_t L);
-
 graph_t *graph_add_ext(const graph_t *G);
-
 void graph_free(graph_t *h);
 
+#ifdef __cplusplus
+}
+#endif
+
diff --git a/lib/initial_finite.c b/lib/initial_finite.c
new file mode 100644
index 0000000..6ea76ef
--- /dev/null
+++ b/lib/initial_finite.c
@@ -0,0 +1,326 @@
+
+#include "initial_finite.h"
+
+double *Jprobs_from_J(q_t q, double T, double *J) {
+  double *J_probs = (double *)calloc(pow(q, 2), sizeof(double));
+
+  for (q_t i = 0; i < q; i++) {
+    for (q_t j = 0; j < q; j++) {
+      J_probs[q * i + j] = 1.0 - exp((J[i] - J[j]) / T);
+    }
+  }
+
+  return J_probs;
+}
+
+q_t *initialize_R(q_t q) {
+  q_t *R = (q_t *)malloc(q * sizeof(q_t));
+
+  for (q_t i = 0; i < q; i++) {
+    R[i] = i;
+  }
+
+  return R;
+}
+
+R_t *transformation_bringing_to_zero(q_t q, R_t n_transformations, q_t *transformations) {
+  R_t *destination = (R_t *)malloc(q * sizeof(R_t));
+
+  for (q_t i = 0; i < q; i++) {
+    for (R_t j = 0; j < n_transformations; j++) {
+      if (transformations[q * j + i] == 0) {
+        destination[i] = j;
+      }
+    }
+  }
+
+  return destination;
+}
+
+R_t find_involutions(R_t *destination, q_t q, R_t n_transformations, q_t *transformations) {
+  R_t n_involutions = 0;
+
+  for (R_t i = 1; i < n_transformations; i++) {
+    bool is_involution = true;
+    for (q_t j = 0; j < q; j++) {
+      if (j != transformations[q * i + transformations[q * i + j]]) {
+        is_involution = false;
+        break;
+      }
+    }
+    if (is_involution) {
+      destination[n_involutions] = i;
+      n_involutions++;
+    }
+  }
+
+  return n_involutions;
+}
+
+state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H) {
+  state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t));
+
+  s->D = D;
+  s->L = L;
+
+  {
+    graph_t *g = graph_create_square(D, L);
+    s->nv = g->nv;
+    s->ne = g->ne;
+    s->g = graph_add_ext(g);
+    graph_free(g);
+  }
+
+  s->q = 2;
+
+  s->n_transformations = 2;
+  s->transformations = (q_t *)malloc(2 * 2 * sizeof(q_t));
+  s->transformations[0] = 0;
+  s->transformations[1] = 1;
+  s->transformations[2] = 1;
+  s->transformations[3] = 0;
+
+  s->n_involutions = 1;
+  s->involutions = (R_t *)malloc(1 * sizeof(R_t));
+  s->involutions[0] = 1;
+
+  s->transform_site_to_zero = (R_t *)malloc(2 * sizeof(R_t));
+  s->transform_site_to_zero[0] = 0;
+  s->transform_site_to_zero[1] = 1;
+
+  s->n_bond_types = 2;
+  s->bond_with_zero_type = (q_t *)malloc(2 * sizeof(q_t));
+  s->bond_with_zero_type[0] = 0;
+  s->bond_with_zero_type[1] = 1;
+
+  s->T = T;
+  s->J = (double *)malloc(2 * sizeof(double)); 
+  s->J[0] = 1.0;
+  s->J[1] = -1.0;
+  s->H = (double *)malloc(2 * sizeof(double)); 
+  s->H[0] = H[0];
+  s->H[1] = -H[0];
+
+  s->J_probs = Jprobs_from_J(2, T, s->J);
+  s->H_probs = Jprobs_from_J(2, T, s->H);
+
+  s->spins = (q_t *)calloc(s->nv, sizeof(q_t));
+  s->R = initialize_R(2);
+
+  s->M = (v_t *)calloc(2, sizeof(v_t));
+  s->M[0] = s->nv; // everyone starts in state 0, remember?
+  s->B = (v_t *)calloc(2, sizeof(v_t));
+  s->B[0] = s->ne;
+
+  return s;
+}
+
+state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, double *H) {
+  state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t));
+
+  s->D = D;
+  s->L = L;
+
+  {
+    graph_t *g = graph_create_square(D, L);
+    s->nv = g->nv;
+    s->ne = g->ne;
+    s->g = graph_add_ext(g);
+    graph_free(g);
+  }
+
+  s->q = q;
+  s->n_transformations = factorial(q);
+  s->transformations = symmetric_gen_transformations(q);
+  s->involutions = (R_t *)malloc(s->n_transformations * sizeof(R_t));
+  s->n_involutions = find_involutions(s->involutions, q, s->n_transformations, s->transformations);
+
+  s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations);
+
+  s->n_bond_types = 2;
+
+  s->bond_with_zero_type = (q_t *)malloc(q * sizeof(q_t));
+
+  s->bond_with_zero_type[0] = 0;
+
+  for (q_t i = 1; i < q; i++) {
+    s->bond_with_zero_type[i] = 1;
+  }
+
+  s->T = T;
+  s->J = (double *)calloc(2, sizeof(double)); 
+  s->J[0] = 1.0;
+  s->J[1] = 0.0;
+
+  s->H = (double *)malloc(q * sizeof(double)); 
+  for (q_t i = 0; i < q; i++) {
+    s->H[i] = H[i];
+  }
+
+  s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J);
+  s->H_probs = Jprobs_from_J(q, T, s->H);
+
+  s->spins = (q_t *)calloc(s->nv, sizeof(q_t));
+  s->R = initialize_R(q);
+
+  s->M = (v_t *)calloc(q, sizeof(v_t));
+  s->M[0] = s->nv; // everyone starts in state 0, remember?
+  s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t));
+  s->B[0] = s->ne; // everyone starts in state 0, remember?
+
+  return s;
+}
+
+state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, double *H) {
+  state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t));
+
+  s->D = D;
+  s->L = L;
+
+  {
+    graph_t *g = graph_create_square(D, L);
+    s->nv = g->nv;
+    s->ne = g->ne;
+    s->g = graph_add_ext(g);
+    graph_free(g);
+  }
+
+  s->q = q;
+
+  s->n_transformations = 2 * q;
+  s->transformations = dihedral_gen_transformations(q);
+  s->n_involutions = q;
+  s->involutions = dihedral_gen_involutions(q);
+
+  s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations);
+  s->bond_with_zero_type = malloc(q * sizeof(q_t));
+
+  s->n_bond_types = q / 2 + 1;
+
+  for (q_t i = 0; i < q / 2 + 1; i++) {
+    s->bond_with_zero_type[i] = i;
+  }
+
+  for (q_t i = 1; i < (q + 1) / 2; i++) { 
+    s->bond_with_zero_type[q - i] = i;
+  }
+
+
+  s->T = T;
+  s->J = (double *)malloc(s->n_bond_types * sizeof(double)); 
+
+  for (q_t i = 0; i < s->n_bond_types; i++) {
+    s->J[i] = cos(2 * M_PI * i / ((double)q));
+  }
+
+
+  s->H = (double *)malloc(q * sizeof(double)); 
+  for (q_t i = 0; i < q; i++) {
+    s->H[i] = H[i];
+  }
+
+  s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J);
+  s->H_probs = Jprobs_from_J(q, T, s->H);
+
+  s->spins = (q_t *)calloc(s->nv, sizeof(q_t));
+  s->R = initialize_R(q);
+
+  s->M = (v_t *)calloc(q, sizeof(v_t));
+  s->M[0] = s->nv; // everyone starts in state 0, remember?
+  s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t));
+  s->B[0] = s->ne; // everyone starts in state 0, remember?
+
+  return s;
+}
+
+
+state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double *H) {
+  state_finite_t *s = (state_finite_t *)calloc(1, sizeof(state_finite_t));
+
+  s->D = D;
+  s->L = L;
+
+  {
+    graph_t *g = graph_create_square(D, L);
+    s->nv = g->nv;
+    s->ne = g->ne;
+    s->g = graph_add_ext(g);
+    graph_free(g);
+  }
+
+  s->q = q;
+
+  s->n_transformations = 2 * q;
+  s->transformations = dihedral_gen_transformations(q);
+  s->n_involutions = q;
+  s->involutions = dihedral_gen_involutions(q);
+
+  s->transform_site_to_zero = transformation_bringing_to_zero(q, s->n_transformations, s->transformations);
+  s->bond_with_zero_type = malloc(q * sizeof(q_t));
+
+  s->n_bond_types = q / 2 + 1;
+
+  for (q_t i = 0; i < q / 2 + 1; i++) {
+    s->bond_with_zero_type[i] = i;
+  }
+
+  for (q_t i = 1; i < (q + 1) / 2; i++) { 
+    s->bond_with_zero_type[(int)q - (int)i] = i;
+  }
+
+  s->T = T;
+  s->J = (double *)malloc(s->n_bond_types * sizeof(double)); 
+
+  for (q_t i = 0; i < s->n_bond_types; i++) {
+    s->J[i] = -pow(i, 2);
+  }
+
+  s->H = (double *)malloc(q * sizeof(double)); 
+  for (q_t i = 0; i < q; i++) {
+    s->H[i] = H[i];
+  }
+
+  s->J_probs = Jprobs_from_J(s->n_bond_types, T, s->J);
+  s->H_probs = Jprobs_from_J(q, T, s->H);
+
+  s->spins = (q_t *)calloc(s->nv, sizeof(q_t));
+  s->R = initialize_R(q);
+
+  s->M = (v_t *)calloc(q, sizeof(v_t));
+  s->M[0] = s->nv; // everyone starts in state 0, remember?
+  s->B = (v_t *)calloc(s->n_bond_types, sizeof(v_t));
+  s->B[0] = s->nv; // everyone starts in state 0, remember?
+
+  return s;
+}
+
+double state_finite_energy(state_finite_t *s) {
+  double E = 0;
+
+  for (q_t i = 0; i < s->n_bond_types; i++) {
+    E += s->J[i] * s->B[i];
+  }
+  for (q_t i = 0; i < s->q; i++) {
+    E += s->H[i] * s->M[i];
+  }
+
+  return -E;
+}
+
+void state_finite_free(state_finite_t *s) {
+  graph_free(s->g);
+  free(s->J);
+  free(s->H);
+  free(s->J_probs);
+  free(s->H_probs);
+  free(s->spins);
+  free(s->R);
+  free(s->M);
+  free(s->B);
+  free(s->transformations);
+  free(s->involutions);
+  free(s->transform_site_to_zero);
+  free(s->bond_with_zero_type);
+  free(s);
+}
+
diff --git a/lib/initial_finite.h b/lib/initial_finite.h
new file mode 100644
index 0000000..542f923
--- /dev/null
+++ b/lib/initial_finite.h
@@ -0,0 +1,27 @@
+
+#pragma once
+
+#include <stdbool.h>
+
+#include "types.h"
+#include "dihedral.h"
+#include "cluster_finite.h"
+
+static char *finite_model_t_strings[] = {"ISING", "POTTS", "CLOCK", "DGM"};
+
+typedef enum {
+  ISING,
+  POTTS,
+  CLOCK,
+  DGM
+} finite_model_t;
+
+state_finite_t *initial_finite_prepare_ising(D_t D, L_t L, double T, double *H);
+state_finite_t *initial_finite_prepare_potts(D_t D, L_t L, q_t q, double T, double *H);
+state_finite_t *initial_finite_prepare_clock(D_t D, L_t L, q_t q, double T, double *H);
+state_finite_t *initial_finite_prepare_dgm(D_t D, L_t L, q_t q, double T, double *H);
+
+void state_finite_free(state_finite_t *s);
+
+double state_finite_energy(state_finite_t *s);
+
diff --git a/lib/measurement.c b/lib/measurement.c
index ad824f6..b30cf6b 100644
--- a/lib/measurement.c
+++ b/lib/measurement.c
@@ -1,6 +1,15 @@
 
+#include "convex.h"
 #include "measurement.h"
 
+meas_t *meas_initialize(count_t W) {
+  meas_t *m = (meas_t *)calloc(1, sizeof(meas_t));
+  m->W = W;
+  m->xx = (double *)calloc(2 * W + 1, sizeof(double));
+
+  return m;
+}
+
 double add_to_avg(double mx, double x, count_t n) {
   return mx * (n / (n + 1.0)) + x / (n + 1.0);
 }
@@ -10,24 +19,42 @@ void meas_update(meas_t *m, double x) {
 
   m->x = add_to_avg(m->x, x, n);
   m->x2 = add_to_avg(m->x2, pow(x, 2), n);
+  m->x4 = add_to_avg(m->x4, pow(x, 4), n);
 
   m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n);
   m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n);
 
-  /*
-  if (n > 1) {
-    double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2));
-    m->dx = sqrt(s2 / n);
-    m->c = s2;
-    m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n);
+  dll_t *tmp_window = m->x_window;
+  dll_t *pos_save;
+  count_t t = 0;
+
+  while (tmp_window != NULL) {
+    m->xx[t] = add_to_avg(m->xx[t], x * (tmp_window->x), m->n - t - 1);
+    t++;
+    if (t == 2 * (m->W)) {
+      pos_save = tmp_window;
+    }
+    tmp_window = tmp_window->next;
   }
-  */
+
+  if (t == 2 * (m->W) + 1) {
+    if (2 * (m->W) + 1 == 1) {
+      free(m->x_window);
+      m->x_window = NULL;
+    } else {
+      free(pos_save->next);
+      pos_save->next = NULL;
+    }
+  }
+
+  stack_push_d(&(m->x_window), x);
 
   (m->n)++;
 }
 
 double meas_dx(const meas_t *m) {
-  return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2)));
+  return 2 * get_tau(m) * Cxx(m, 0) / m->n;
+//  return sqrt(1. / (m->n - 1.) * (m->x2 - pow(m->x, 2)));
 }
 
 double meas_c(const meas_t *m) {
@@ -74,3 +101,105 @@ double rho(const autocorr_t *o, count_t i) {
   return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2));
 }
 
+double Cxx(const meas_t *m, count_t t) {
+  return m->xx[t] - pow(m->x, 2);
+}
+
+double rho_m(const meas_t *m, count_t t) {
+  return Cxx(m, t) / Cxx(m, 0);
+}
+
+double get_tau(const meas_t *m) {
+  double *Gammas = (double *)malloc((m->W + 1) * sizeof(double));
+
+  Gammas[0] = 1 + rho_m(m, 0);
+  for (uint64_t i = 0; i < m->W; i++) {
+    Gammas[1 + i] = rho_m(m, 2 * i + 1) + rho_m(m, 2 * i + 2);
+  }
+
+  uint64_t n;
+  for (n = 0; n < m->W + 1; n++) {
+    if (Gammas[n] <= 0) {
+      break;
+    }
+  }
+
+  double *conv_Gamma = get_convex_minorant(n, Gammas);
+
+  double tau = - 0.5;
+
+  for (uint64_t i = 0; i < n + 1; i++) {
+    tau += conv_Gamma[i];
+  }
+
+  free(Gammas);
+
+  return tau;
+}
+
+void print_meas(const meas_t *m, const char *sym, FILE *outfile) {
+  fprintf(outfile, "%s-><|n->%" PRIcount ",x->%.15f,x^2->%.15f,x^4->%.15f,xx->{", sym, m->n, m->x, m->x2, m->x4);
+  for (count_t i = 0; i < 2 * (m->W) + 1; i++) {
+    fprintf(outfile, "%.15f", m->xx[i]);
+    if (i < 2 * (m->W)) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "}|>");
+}
+
+void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile) {
+  fprintf(outfile, "%s-><|n->{", sym);
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%" PRIcount, m[i]->n);
+    if (i < q - 1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},x->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", m[i]->x);
+    if (i < q - 1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},x^2->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", m[i]->x2);
+    if (i < q - 1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},x^4->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "%.15f", m[i]->x4);
+    if (i < q - 1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "},xx->{");
+  for (q_t i = 0; i < q; i++) {
+    fprintf(outfile, "{");
+    for (count_t j = 0; j < 2 * (m[i]->W) + 1; j++) {
+      fprintf(outfile, "%.15f", m[i]->xx[j]);
+      if (j < 2 * (m[i]->W)) {
+        fprintf(outfile, ",");
+      }
+    }
+    fprintf(outfile, "}");
+    if (i < q - 1) {
+      fprintf(outfile, ",");
+    }
+  }
+  fprintf(outfile, "}|>");
+}
+
+void free_meas(meas_t *m) {
+  free(m->xx);
+  while (m->x_window != NULL) {
+    stack_pop_d(&(m->x_window));
+  }
+  free(m);
+}
+
+
diff --git a/lib/measurement.h b/lib/measurement.h
index 46c034f..d9bd52e 100644
--- a/lib/measurement.h
+++ b/lib/measurement.h
@@ -3,16 +3,21 @@
 
 #include <math.h>
 #include <stdlib.h>
+#include <stdio.h>
 
 #include "types.h"
 #include "stack.h"
 
 typedef struct {
-  uint64_t n;
+  count_t n;
   double x;
   double x2;
+  double x4;
   double m2;
   double m4;
+  count_t W;
+  double *xx;
+  dll_t *x_window;
 } meas_t;
 
 typedef struct {
@@ -24,11 +29,6 @@ typedef struct {
   double O2;
 } autocorr_t;
 
-typedef struct {
-  void (*f)(state_finite_t *, void *);
-  void *data;
-} measurement_t;
-
 void meas_update(meas_t *m, double x);
 
 double meas_dx(const meas_t *m);
@@ -41,3 +41,14 @@ void update_autocorr(autocorr_t *OO, double O);
 
 double rho(const autocorr_t *o, uint64_t i);
 
+void print_meas(const meas_t *m, const char *sym, FILE *outfile);
+void print_vec_meas(q_t q, const meas_t **m, const char *sym, FILE *outfile);
+
+void free_meas(meas_t *m);
+
+meas_t *meas_initialize(count_t W);
+
+double get_tau(const meas_t *m);
+
+double Cxx(const meas_t *m, count_t t);
+
diff --git a/lib/orthogonal.c b/lib/orthogonal.c
deleted file mode 100644
index 87569ae..0000000
--- a/lib/orthogonal.c
+++ /dev/null
@@ -1,99 +0,0 @@
-
-#include "orthogonal.h"
-
-void vector_replace(q_t n, double *v1, const double *v2) {
-  // writes vector v2 of length n to memory located at v1
-  for (q_t i = 0; i < n; i++) {
-    v1[i] = v2[i];
-  }
-}
-
-void vector_add(q_t n, double *v1, const double *v2) {
-  // adds vector v2 of length n to vector v1
-  for (q_t i = 0; i < n; i++) {
-    v1[i] += v2[i];
-  }
-}
-
-void vector_subtract(q_t n, double *v1, const double *v2) {
-  // subtracts vector v2 of length n from vector v1
-  for (q_t i = 0; i < n; i++) {
-    v1[i] -= v2[i];
-  }
-}
-
-double *vector_rotate(q_t n, const double *rot, const double *vec) {
-  // multiplies n by n rotation matrix rot to vector vec
-  double *rot_vec = (double *)malloc(n * sizeof(double));
-
-  double prod = 0.0;
-  for (q_t i = 0; i < n; i++) {
-    prod += rot[i] * vec[i];
-  }
-
-  for (q_t i = 0; i < n; i++) {
-    rot_vec[i] = vec[i] - 2 * prod * rot[i];
-  }
-
-  return rot_vec;
-}
-
-double *vector_rotate_inverse(q_t n, const double *rot, const double *vec) {
-  double *rot_vec = (double *)calloc(n, sizeof(double));
-
-  for (q_t i = 0; i < n; i++) {
-    for (q_t j = 0; j < n; j++) {
-      rot_vec[i] += rot[n * j + i] * vec[j];
-    }
-  }
-
-  return rot_vec;
-}
-
-double vector_dot(q_t n, const double *v1, const double *v2) {
-  double dot = 0;
-
-  for (q_t i = 0; i < n; i++) {
-    dot += v1[i] * v2[i];
-  }
-
-  return dot;
-}
-
-double *orthogonal_rotate(q_t n, const double *r, const double *m) {
-  double *mul = (double *)calloc(n * n, sizeof(double));
-
-  for (q_t i = 0; i < n; i++) {
-    double akOki = 0;
-
-    for (q_t k = 0; k < n; k++) {
-      akOki += r[k] * m[n * k + i];
-    }
-
-    for (q_t j = 0; j < n; j++) {
-      mul[n * j + i] = m[n * j + i] - 2 * akOki * r[j];
-    }
-  }
-
-  return mul;
-}
-
-double *gen_rot(gsl_rng *r, q_t n) {
-  double *v = (double *)malloc(n * sizeof(double));
-
-  double v2 = 0;
-
-  for (q_t i = 0; i < n; i++) {
-    v[i] = gsl_ran_ugaussian(r);
-    v2 += v[i] * v[i];
-  }
-
-  double magv = sqrt(v2);
-
-  for (q_t i = 0; i < n; i++) {
-    v[i] /= magv;
-  }
-
-  return v;
-}
-
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 60d5f49..85f7a20 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -1,24 +1,163 @@
 
+#pragma once
+
 #include <stdlib.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
-#include <math.h>
+#include <cmath>
 
 #include "types.h"
 
-void vector_replace(q_t n, double *v1, const double *v2);
+template <q_t q, class T>
+struct orthogonal_t { bool is_reflection; T *x; };
+
+template <q_t q, class T>
+void init(orthogonal_t <q, T> *ptr) {
+  ptr->is_reflection = false;
+  ptr->x = (T *)calloc(q * q, sizeof(T));
+
+  for (q_t i = 0; i < q; i++) {
+    ptr->x[q * i + i] = (T)1;
+  }
+}
+
+template <q_t q, class T>
+orthogonal_t <q, T> copy (orthogonal_t <q, T> m) {
+  orthogonal_t <q, T> m_copy;
+  m_copy.is_reflection = m.is_reflection;
+
+  q_t size;
+
+  if (m.is_reflection) {
+    size = q;
+  } else {
+    size = q * q;
+  }
+
+  m_copy.x = (T *)calloc(size, sizeof(T));
+
+  for (q_t i = 0; i < size; i++) {
+    m_copy.x[i] = m.x[i];
+  }
+
+  return m_copy;
+}
+
+template <q_t q, class T>
+void free_spin (orthogonal_t <q, T> m) {
+  free(m.x);
+}
+
+template <q_t q, class T>
+vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) {
+  vector_t <q, T> v_rot;
+  v_rot.x = (T *)calloc(q, sizeof(T));
+
+  if (m.is_reflection) {
+    double prod = 0;
+    for (q_t i = 0; i < q; i++) {
+      prod += v.x[i] * m.x[i];
+    }
+    for (q_t i = 0; i < q; i++) {
+      v_rot.x[i] = v.x[i] - 2 * prod * m.x[i];
+    }
+  } else {
+    for (q_t i = 0; i < q; i++) {
+      for (q_t j = 0; j < q; j++) {
+        v_rot.x[i] += m.x[q * i + j] * v.x[j];
+      }
+    }
+  }
+
+  return v_rot;
+}
+
+template <q_t q, class T>
+orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
+  orthogonal_t <q, T> m2_rot;
+
+  m2_rot.is_reflection = false;
+  m2_rot.x = (T *)calloc(q * q, sizeof(T));
+
+  if (m1.is_reflection) {
+    for (q_t i = 0; i < q; i++) {
+      double akOki = 0;
+
+      for (q_t k = 0; k < q; k++) {
+        akOki += m1.x[k] * m2.x[q * k + i];
+      }
+
+      for (q_t j = 0; j < q; j++) {
+        m2_rot.x[q * j + i] = m2.x[q * j + i] - 2 * akOki * m1.x[j];
+      }
+    }
+  } else {
+    for (q_t i = 0; i < q; i++) {
+      for (q_t j = 0; j < q; j++) {
+        for (q_t k = 0; k < q; k++) {
+          m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k];
+        }
+      }
+    }
+  }
+
+  return m2_rot;
+}
+
+template <q_t q, class T>
+vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) {
+  if (m.is_reflection) {
+    return act(m, v); // reflections are their own inverse
+  } else {
+    vector_t <q, T> v_rot;
+    v_rot.x = (T *)calloc(q, sizeof(T));
+
+    for (q_t i = 0; i < q; i++) {
+      for (q_t j = 0; j < q; j++) {
+        v_rot.x[i] += m.x[q * j + i] * v.x[j];
+      }
+    }
+
+    return v_rot;
+  }
+}
+
+template <q_t q, class T>
+orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
+  if (m1.is_reflection) {
+    return act(m1, m2); // reflections are their own inverse
+  } else {
+    orthogonal_t <q, T> m2_rot;
+    m2_rot.x = (T *)calloc(q * q, sizeof(T));
 
-void vector_add(q_t n, double *v1, const double *v2);
+    for (q_t i = 0; i < q; i++) {
+      for (q_t j = 0; j < q; j++) {
+        for (q_t k = 0; k < q; k++) {
+          m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k];
+        }
+      }
+    }
 
-void vector_subtract(q_t n, double *v1, const double *v2);
+    return m2_rot;
+  }
+}
 
-double *vector_rotate(q_t n, const double *rot, const double *vec);
+template <q_t q>
+void generate_rotation (gsl_rng *r, orthogonal_t <q, double> *ptr) {
+  ptr->is_reflection = true;
+  ptr->x = (double *)calloc(q, sizeof(double));
 
-double *vector_rotate_inverse(q_t n, const double *rot, const double *vec);
+  double v2 = 0;
 
-double vector_dot(q_t n, const double *v1, const double *v2);
+  for (q_t i = 0; i < q; i++) {
+    ptr->x[i] = gsl_ran_ugaussian(r);
+    v2 += ptr->x[i] * ptr->x[i];
+  }
 
-double *orthogonal_rotate(q_t n, const double *m1, const double *m2);
+  double mag_v = sqrt(v2);
 
-double *gen_rot(gsl_rng *r, q_t n);
+  for (q_t i = 0; i < q; i++) {
+    ptr->x[i] /= mag_v;
+  }
+}
 
diff --git a/lib/rand.h b/lib/rand.h
index 2354f6a..7bb5354 100644
--- a/lib/rand.h
+++ b/lib/rand.h
@@ -4,4 +4,13 @@
 #include <assert.h>
 #include <stdio.h>
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 unsigned long int rand_seed();
+
+#ifdef __cplusplus
+}
+#endif
+
diff --git a/lib/stack.h b/lib/stack.h
index a354ab5..8d25aff 100644
--- a/lib/stack.h
+++ b/lib/stack.h
@@ -8,6 +8,11 @@
 
 #include "types.h"
 
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 typedef struct ll_tag {
   v_t x;
   struct ll_tag *next;
@@ -24,3 +29,7 @@ void stack_push_d(dll_t **q, double x);
 v_t stack_pop(ll_t **q);
 double stack_pop_d(dll_t **q);
 
+#ifdef __cplusplus
+}
+#endif
+
diff --git a/lib/symmetric.c b/lib/symmetric.c
index 729b38c..4487538 100644
--- a/lib/symmetric.c
+++ b/lib/symmetric.c
@@ -25,3 +25,46 @@ q_t *symmetric_invert(q_t q, const q_t *g) {
   return g_inv;
 }
 
+void swap(q_t *q1, q_t *q2) {
+  q_t temp = *q1;
+  *q1 = *q2;
+  *q2 = temp;
+}
+
+R_t factorial(q_t q) {
+  if (q == 0) {
+    return 1;
+  } else {
+    return q * factorial(q - 1);
+  }
+}
+
+void permute(q_t *a, q_t l, q_t r, R_t pos, q_t *transformations) {
+  if (l == r - 1) {
+    for (q_t i = 0; i < r; i++) {
+      transformations[r * pos + i] = a[i];
+    }
+  } else {
+    for (q_t i = l; i < r; i++) {
+      swap((a+l), (a+i));
+      permute(a, l+1, r, pos + (i - l) * factorial(r - l - 1), transformations);
+      swap((a+l), (a+i));
+    }
+  }
+}
+
+q_t *symmetric_gen_transformations(q_t q) {
+  q_t *transformations = (q_t *)malloc(q * factorial(q) * sizeof(q_t));
+  q_t *tmp = (q_t *)malloc(q * sizeof(q_t));
+
+  for (q_t i = 0; i < q; i++) {
+    tmp[i] = i;
+  }
+
+  permute(tmp, 0, q, 0, transformations);
+
+  free(tmp);
+
+  return transformations;
+}
+
diff --git a/lib/symmetric.h b/lib/symmetric.h
index 6e00f52..c71521d 100644
--- a/lib/symmetric.h
+++ b/lib/symmetric.h
@@ -11,3 +11,5 @@ q_t symmetric_act(const q_t *g, q_t s);
 
 q_t *symmetric_invert(q_t q, const q_t *g);
 
+q_t *symmetric_gen_transformations(q_t q);
+
diff --git a/lib/vector.h b/lib/vector.h
new file mode 100644
index 0000000..62ce59e
--- /dev/null
+++ b/lib/vector.h
@@ -0,0 +1,85 @@
+
+#pragma once
+
+#include <stdlib.h>
+#include <cmath>
+
+#include "types.h"
+
+template <q_t q, class T>
+struct vector_t { T *x; };
+
+template <q_t q, class T>
+void init(vector_t <q, T> *ptr) {
+  ptr->x = (T *)calloc(q, sizeof(T));
+
+  ptr->x[0] = (T)1;
+}
+
+template <q_t q, class T>
+vector_t <q, T> copy (vector_t <q, T> v) {
+  vector_t <q, T> v_copy;
+ 
+ v_copy.x = (T *)calloc(q, sizeof(T));
+
+  for (q_t i = 0; i < q; i++) {
+    v_copy.x[i] = v.x[i];
+  }
+
+  return v_copy;
+}
+
+template <q_t q, class T>
+void add (vector_t <q, T> v1, vector_t <q, T> v2) {
+  for (q_t i = 0; i < q; i++) {
+    v1.x[i] += v2.x[i];
+  }
+}
+
+template <q_t q, class T>
+void subtract (vector_t <q, T> v1, vector_t <q, T> v2) {
+  for (q_t i = 0; i < q; i++) {
+    v1.x[i] -= v2.x[i];
+  }
+}
+
+template <q_t q, class T>
+vector_t <q, T> scalar_multiple(v_t a, vector_t <q, T> v) {
+  vector_t <q, T> multiple;
+  multiple.x = (T *)malloc(q * sizeof(T));
+  for (q_t i = 0; i < q; i++) {
+    multiple.x[i] = a * v.x[i];
+  }
+
+  return multiple;
+}
+
+template <q_t q, class T>
+T dot(vector_t <q, T> v1, vector_t <q, T> v2) {
+  T prod = 0;
+
+  for (q_t i = 0; i < q; i++) {
+    prod += v1.x[i] * v2.x[i];
+  }
+
+  return prod;
+}
+
+template <q_t q, class T>
+void free_spin (vector_t <q, T> v) {
+  free(v.x);
+}
+
+template <q_t q, class T>
+void write_magnetization(vector_t <q, T> M, FILE *outfile) {
+  fwrite(M.x, sizeof(double), q, outfile);
+}
+
+template <q_t q> // save some space and don't write whole doubles
+void write_magnetization(vector_t <q, double> M, FILE *outfile) {
+  for (q_t i = 0; i < q; i++) {
+    float M_tmp = (float)M.x[i];
+    fwrite(&M_tmp, sizeof(float), 1, outfile);
+  }
+}
+
diff --git a/lib/wolff.h b/lib/wolff.h
new file mode 100644
index 0000000..caf413b
--- /dev/null
+++ b/lib/wolff.h
@@ -0,0 +1,73 @@
+
+#include <time.h>
+#include <getopt.h>
+
+#include <cluster.h>
+
+template <q_t q, class T>
+double H_vector(vector_t <q, T> v1, T *H) {
+  vector_t <q, T> H_vec;
+  H_vec.x = H;
+  return (double)(dot <q, T> (v1, H_vec));
+}
+
+template <class R_t, class X_t>
+void wolff(count_t N, D_t D, L_t L, double T, std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, unsigned long timestamp, bool silent) {
+
+  state_t <R_t, X_t> s(D, L, T, J, H);
+
+  // initialize random number generator
+  gsl_rng *r = gsl_rng_alloc(gsl_rng_mt19937);
+  gsl_rng_set(r, rand_seed());
+
+  char *filename_M = (char *)malloc(255 * sizeof(char));
+  char *filename_E = (char *)malloc(255 * sizeof(char));
+  char *filename_S = (char *)malloc(255 * sizeof(char));
+
+  sprintf(filename_M, "wolff_%lu_M.dat", timestamp);
+  sprintf(filename_E, "wolff_%lu_E.dat", timestamp);
+  sprintf(filename_S, "wolff_%lu_S.dat", timestamp);
+
+  FILE *outfile_M = fopen(filename_M, "wb");
+  FILE *outfile_E = fopen(filename_E, "wb");
+  FILE *outfile_S = fopen(filename_S, "wb");
+
+  free(filename_M);
+  free(filename_E);
+  free(filename_S);
+
+  v_t cluster_size = 0;
+
+  if (!silent) printf("\n");
+  for (count_t steps = 0; steps < N; steps++) {
+    if (!silent) printf("\033[F\033[JWOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", steps, N, s.E, cluster_size);
+
+    v_t v0 = gsl_rng_uniform_int(r, s.nv);
+     
+    R_t step;
+    generate_rotation(r, &step);
+
+    cluster_size = flip_cluster <R_t, X_t> (&s, v0, step, r);
+
+    free_spin(step);
+
+    {
+      float smaller_E = (float)s.E;
+      fwrite(&smaller_E, sizeof(float), 1, outfile_E);
+    }
+    write_magnetization(s.M, outfile_M);
+    fwrite(&cluster_size, sizeof(uint32_t), 1, outfile_S);
+
+  }
+  if (!silent) {
+    printf("\033[F\033[J");
+  }
+  printf("WOLFF: sweep %" PRIu64 " / %" PRIu64 ": E = %.2f, S = %" PRIv "\n", N, N, s.E, cluster_size);
+
+  fclose(outfile_M);
+  fclose(outfile_E);
+  fclose(outfile_S);
+
+  gsl_rng_free(r);
+}
+
diff --git a/lib/wolff_finite.c b/lib/wolff_finite.c
deleted file mode 100644
index 64de9ba..0000000
--- a/lib/wolff_finite.c
+++ /dev/null
@@ -1,70 +0,0 @@
-
-#include "cluster_finite.h"
-
-void wolff_finite(state_finite_t *s, count_t sweeps, count_t sweeps_per_measurement, count_t n_measurements, measurement_t *measurements) {
-  for (count_t i = 0; i < sweeps; i++) {
-
-    count_t n_flips = 0;
-
-    while (n_flips / h->nv < sweeps_per_measurement) {
-      v_t v0 = gsl_rng_uniform_int(r, h->nv);
-      R_t step;
-     
-      bool changed = false;
-      while (!changed) {
-        step = gsl_rng_uniform_int(r, s->n_transformations);
-        if v(symmetric_act(s->transformations + q * step, s->spins[v0]) != s->spins[v0]) {
-          changed = true;
-        }
-      }
-
-      v_t tmp_flips = flip_cluster_finite(s, v0, step, r);
-      n_flips += tmp_flips;
-
-      if (n_runs > 0) {
-        n_steps++;
-        meas_update(clust, tmp_flips);
-
-        if (record_autocorrelation && n_steps % ac_skip == 0) {
-          update_autocorr(autocorr, s->E);
-        }
-
-      }
-
-    }
-
-    for (q_t i = 0; i < q; i++) {
-      meas_update(M[i], s->M[i]);
-    }
-    meas_update(E, s->E);
-
-    q_t n_at_max = 0;
-    q_t max_M_i = 0;
-    v_t max_M = 0;
-
-    for (q_t i = 0; i < q; i++) {
-      if (s->M[i] > max_M) {
-        n_at_max = 1;
-        max_M_i = i;
-        max_M = s->M[i];
-      } else if (s->M[i] == max_M) {
-        n_at_max++;
-      }
-    }
-
-    if (record_distribution) {
-      mag_dist[s->M[0]]++;
-    }
-
-    if (n_at_max == 1) {
-      for (q_t i = 0; i < q; i++) {
-        meas_update(sM[max_M_i][i], s->M[i]);
-      }
-      meas_update(sE[max_M_i], s->E);
-      freqs[max_M_i]++;
-    }
-
-    diff = fabs(meas_dx(clust) / clust->x);
-  }
-}
-