all the R_t have been objectified

16 files changed, 247 insertions, 395 deletions

diff --git a/lib/cluster.h b/lib/cluster.h
index 427c3c8..c5f2be7 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -32,13 +32,13 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
       marks[v] = true;
 
       if (v == state->nv) {
-        R_new = act (r, R_old);
+        R_new = r.act(R_old);
       } else {
         si_old = state->spins[v];
-        si_new = act (r, si_old);
+        si_new = r.act(si_old);
       }
 
-      for (v_t vn : state->g.v_adj[v]) {
+      for (const v_t &vn : state->g.v_adj[v]) {
         X_t sj;
 
         if (vn != state->nv) {
@@ -53,17 +53,17 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
           X_t rs_old, rs_new;
           v_t non_ghost;
           if (vn == state->nv) {
-            rs_old = act_inverse (R_old, si_old);
-            rs_new = act_inverse (R_old, si_new);
+            rs_old = R_old.act_inverse(si_old);
+            rs_new = R_old.act_inverse(si_new);
             non_ghost = v;
           } else {
-            rs_old = act_inverse (R_old, sj);
-            rs_new = act_inverse (R_new, sj);
+            rs_old = R_old.act_inverse(sj);
+            rs_new = R_new.act_inverse(sj);
             non_ghost = vn;
           }
           double dE = state->H(rs_old) - state->H(rs_new);
 #ifdef FINITE_STATES
-          prob = H_probs[N_STATES * state_to_ind(rs_old) + state_to_ind(rs_new)];
+          prob = H_probs[state_to_ind(rs_old)][state_to_ind(rs_new)];
 #else
           prob = 1.0 - exp(-dE / state->T);
 #endif
@@ -85,7 +85,7 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
         } else {
           double dE = state->J(si_old, sj) - state->J(si_new, sj);
 #ifdef FINITE_STATES
-          prob = J_probs[N_STATES * N_STATES * state_to_ind(si_old) + N_STATES * state_to_ind(si_new) + state_to_ind(sj)];
+          prob = J_probs[state_to_ind(si_old)][state_to_ind(si_new)][state_to_ind(sj)];
 #else
           prob = 1.0 - exp(-dE / state->T);
 #endif
@@ -98,13 +98,9 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
       }
 
       if (v == state->nv) {
-        free_spin(state->R);
         state->R = R_new;
       } else {
         state->spins[v] = si_new;
-      }
-
-      if (v != state->nv) { // count the number of non-external sites that flip
         nv++;
       }
     }
diff --git a/lib/dihedral.h b/lib/dihedral.h
index 8238835..5588afa 100644
--- a/lib/dihedral.h
+++ b/lib/dihedral.h
@@ -4,76 +4,46 @@
 #include "types.h"
 #include "potts.h"
 
-template <class T, q_t q>
-struct dihedral_t { bool is_reflection; T x; };
-
-template <class T, q_t q>
-void init(dihedral_t<T, q> *ptr) {
-  ptr->is_reflection = false;
-  ptr->x = (T)0;
-}
-
-template <class T, q_t q>
-dihedral_t<T, q> copy(dihedral_t<T, q> r) {
-  return r;
-}
-
-template <class T, q_t q>
-void free_spin(dihedral_t<T, q> r) {
-  // do nothing!
-}
-
-template <q_t q>
-potts_t<q> act(dihedral_t<q_t, q> r, potts_t<q> s) {
-  potts_t<q> s2;
-  if (r.is_reflection) {
-    s2.x = ((q + r.x) - s.x) % q;
-  } else {
-    s2.x = (r.x + s.x) % q;
-  }
-
-  return s2;
-}
-
-template <q_t q>
-dihedral_t<q_t,q> act(dihedral_t<q_t,q> r1, dihedral_t<q_t,q> r2) {
-  dihedral_t<q_t,q> r3;
-
-  if (r1.is_reflection) {
-    r3.is_reflection = !(r2.is_reflection);
-    r3.x = ((q + r1.x) - r2.x) % q;
-  } else {
-    r3.is_reflection = r2.is_reflection;
-    r3.x = (r1.x + r2.x) % q;
-  }
-
-  return r3;
-}
-
-template <q_t q>
-potts_t<q> act_inverse(dihedral_t<q_t,q> r, potts_t<q> s) {
-  potts_t<q> s2;
-  if (r.is_reflection) {
-    s2.x = ((r.x + q) - s.x) % q;
-  } else {
-    s2.x = ((s.x + q) - r.x) % q;
-  }
-
-  return s2;
-}
-
 template <q_t q>
-dihedral_t<q_t, q> act_inverse(dihedral_t<q_t,q> r1, dihedral_t<q_t,q> r2) {
-  dihedral_t<q_t,q> r3;
-
-  if (r1.is_reflection) {
-    r3.is_reflection = !(r2.is_reflection);
-    r3.x = ((r1.x + q) - r2.x) % q;
-  } else {
-    r3.is_reflection = r2.is_reflection;
-    r3.x = ((r2.x + q) - r1.x) % q;
-  }
-
-  return r3;
-}
+class dihedral_t {
+  public:
+    bool is_reflection;
+    q_t x;
+
+    dihedral_t() : is_reflection(false), x(0) {}
+    dihedral_t(bool x, q_t y) : is_reflection(x), x(y) {}
+
+    potts_t<q> act(const potts_t<q>& s) const {
+      if (this->is_reflection) {
+        return potts_t<q>(((q + this->x) - s.x) % q);
+      } else {
+        return potts_t<q>((this->x + s.x) % q);
+      }
+    }
+
+
+    dihedral_t<q> act(dihedral_t<q> r) const {
+      if (this->is_reflection) {
+        return dihedral_t<q>(!(r.is_reflection), ((q + this->x) - r.x) % q); 
+      } else {
+        return dihedral_t<q>(r.is_reflection, (this->x + r.x) % q);
+      }
+    }
+
+    potts_t<q> act_inverse(potts_t<q> s) const {
+      if (this->is_reflection) {
+        return this->act(s);
+      } else {
+        return potts_t<q>(((s.x + q) - this->x) % q);
+      }
+    }
+
+    dihedral_t<q> act_inverse(dihedral_t<q> r) const {
+      if (this->is_reflection) {
+        return this->act(r);
+      } else {
+        return dihedral_t<q>(r.is_reflection, ((r.x + q) - this->x) % q);
+      }
+    }
+};
 
diff --git a/lib/dihedral_inf.h b/lib/dihedral_inf.h
index f7c4297..a064b48 100644
--- a/lib/dihedral_inf.h
+++ b/lib/dihedral_inf.h
@@ -4,75 +4,44 @@
 #include "height.h"
 
 template <class T>
-struct dihedral_inf_t { bool is_reflection; T x; };
-
-template <class T>
-void init(dihedral_inf_t<T> *ptr) {
-  ptr->is_reflection = false;
-  ptr->x = (T)0;
-}
-
-template <class T>
-dihedral_inf_t<T> copy(dihedral_inf_t<T> r) {
-  return r;
-}
-
-template <class T>
-void free_spin(dihedral_inf_t<T> r) {
-  // do nothing!
-}
-
-template <class T>
-height_t<T> act(dihedral_inf_t<T> r, height_t<T> h) {
-  height_t<T> h2;
-  if (r.is_reflection) {
-    h2.x = r.x - h.x;
-  } else {
-    h2.x = r.x + h.x;
+class dihedral_inf_t {
+  public:
+  bool is_reflection;
+  T x;
+
+  dihedral_inf_t() : is_reflection(false), x(0) {}
+  dihedral_inf_t(bool x, T y) : is_reflection(x), x(y) {}
+
+  height_t<T> act(const height_t<T>& h) const {
+    if (this->is_reflection) {
+      return height_t(this->x - h.x);
+    } else {
+      return height_t(this->x + h.x);
+    }
   }
 
-  return h2;
-}
-
-template <class T>
-dihedral_inf_t<T> act(dihedral_inf_t<T> r1, dihedral_inf_t<T> r2) {
-  dihedral_inf_t<T> r3;
-
-  if (r1.is_reflection) {
-    r3.is_reflection = !(r2.is_reflection);
-    r3.x = r1.x - r2.x;
-  } else {
-    r3.is_reflection = r2.is_reflection;
-    r3.x = r1.x + r2.x;
+  dihedral_inf_t<T> act(const dihedral_inf_t<T>& r) const {
+    if (this->is_reflection) {
+      return dihedral_inf_t<T>(!r.is_reflection, this->x - r.x);
+    } else {
+      return dihedral_inf_t<T>(r.is_reflection, this->x + r.x);
+    }
   }
 
-  return r3;
-}
-
-template <class T>
-height_t<T> act_inverse(dihedral_inf_t<T> r, height_t<T> h) {
-  height_t<T> h2;
-  if (r.is_reflection) {
-    h2.x = r.x - h.x;
-  } else {
-    h2.x = h.x - r.x;
+  height_t<T> act_inverse(const height_t<T>& h) const {
+    if (this->is_reflection) {
+      return this->act(h);
+    } else {
+      return height_t(h.x - this->x);
+    }
   }
 
-  return h2;
-}
-
-template <class T>
-dihedral_inf_t<T> act_inverse(dihedral_inf_t<T> r1, dihedral_inf_t<T> r2) {
-  dihedral_inf_t<T> r3;
-
-  if (r1.is_reflection) {
-    r3.is_reflection = !(r2.is_reflection);
-    r3.x = r1.x - r2.x;
-  } else {
-    r3.is_reflection = r2.is_reflection;
-    r3.x = r2.x - r1.x;
+  dihedral_inf_t<T> act_inverse(const dihedral_inf_t<T>& r) const {
+    if (this->is_reflection) {
+      return this->act(r);
+    } else {
+      return dihedral_inf_t<T>(r.is_reflection, r.x - this->x);
+    }
   }
-
-  return r3;
-}
+};
 
diff --git a/lib/finite_states.h b/lib/finite_states.h
index 1e8800c..08eff30 100644
--- a/lib/finite_states.h
+++ b/lib/finite_states.h
@@ -3,27 +3,28 @@
 
 #include <cmath>
 #include <functional>
+#include <array>
 
 #define FINITE_STATES
 
 // must have N_STATES, states[N_STATES], and state_to_ind defined before
 // invoking header
 
-double J_probs[N_STATES * N_STATES * N_STATES];
-double H_probs[N_STATES * N_STATES];
+std::array<std::array<std::array<double, N_STATES>, N_STATES>, N_STATES> J_probs;
+std::array<std::array<double, N_STATES>, N_STATES> H_probs;
 
 template <class X_t>
 void initialize_probs(std::function <double(X_t, X_t)> J, std::function <double(X_t)> H, double T) {
   for (q_t i = 0; i < N_STATES; i++) {
     for (q_t j = 0; j < N_STATES; j++) {
       for (q_t k = 0; k < N_STATES; k++) {
-        J_probs[i * N_STATES * N_STATES + j * N_STATES +k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T);
+        J_probs[i][j][k] = 1.0 - exp(-(J(states[i], states[k]) - J(states[j], states[k])) / T);
       }
     }
   }
   for (q_t i = 0; i < N_STATES; i++) {
     for (q_t j = 0; j < N_STATES; j++) {
-      H_probs[i * N_STATES + j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
+      H_probs[i][j] = 1.0 - exp(-(H(states[i]) - H(states[j])) / T);
     }
   }
 }
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index c0958f2..34cd44e 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -11,164 +11,134 @@
 #include "vector.h"
 
 template <q_t q, class T>
-struct orthogonal_t { bool is_reflection; T *x; };
+class orthogonal_t : public std::array<std::array<T, q>, q> {
+  public :
+  bool is_reflection;
 
-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];
+  orthogonal_t() : is_reflection(false) {
+    for (q_t i = 0; i < q; i++) {
+      (*this)[i].fill(0);
+      (*this)[i][i] = (T)1;
+    }
   }
 
-  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;
+  vector_t<q, T> act(const vector_t <q, T>& v) const {
+    vector_t <q, T> v_rot;
 
-  if (m.is_reflection) {
-    double prod = 0;
-    for (q_t i = 0; i < q; i++) {
-      prod += v[i] * m.x[i];
-    }
-    for (q_t i = 0; i < q; i++) {
-      v_rot[i] = v[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[i] += m.x[q * i + j] * v[j];
+    if (is_reflection) {
+      double prod = 0;
+      for (q_t i = 0; i < q; i++) {
+        prod += v[i] * (*this)[0][i];
+      }
+      for (q_t i = 0; i < q; i++) {
+        v_rot[i] = v[i] - 2 * prod * (*this)[0][i];
+      }
+    } else {
+      for (q_t i = 0; i < q; i++) {
+        for (q_t j = 0; j < q; j++) {
+          v_rot[i] += (*this)[i][j] * v[j];
+        }
       }
     }
-  }
 
-  return v_rot;
-}
+    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;
+  orthogonal_t<q, T> act(const orthogonal_t <q, T>& m) const {
+    orthogonal_t <q, T> m_rot;
 
-  m2_rot.is_reflection = false;
-  m2_rot.x = (T *)calloc(q * q, sizeof(T));
+    m_rot.is_reflection = false;
 
-  if (m1.is_reflection) {
-    for (q_t i = 0; i < q; i++) {
-      double akOki = 0;
+    if (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 k = 0; k < q; k++) {
+          akOki += (*this)[0][k] * m[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];
+        for (q_t j = 0; j < q; j++) {
+          m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][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];
+    } else {
+      for (q_t i = 0; i < q; i++) {
+        m_rot[i].fill(0);
+        for (q_t j = 0; j < q; j++) {
+          for (q_t k = 0; k < q; k++) {
+            m_rot[i][j] += (*this)[i][j] * m[j][k];
+          }
         }
       }
     }
-  }
 
-  return m2_rot;
-}
+    return m_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;
+  vector_t <q, T> act_inverse(const vector_t <q, T>& v) const {
+    if (is_reflection) {
+      return this->act(v); // reflections are their own inverse
+    } else {
+      vector_t <q, T> v_rot;
 
-    for (q_t i = 0; i < q; i++) {
-      for (q_t j = 0; j < q; j++) {
-        v_rot[i] += m.x[q * j + i] * v[j];
+      for (q_t i = 0; i < q; i++) {
+        for (q_t j = 0; j < q; j++) {
+          v_rot[i] += (*this)[j][i] * v[j];
+        }
       }
-    }
 
-    return v_rot;
+      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));
 
-    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];
+  vector_t <q, T> act_inverse(const orthogonal_t <q, T>& m) const {
+    if (is_reflection) {
+      return this->act(m); // reflections are their own inverse
+    } else {
+      orthogonal_t <q, T> m_rot;
+      m_rot.is_reflection = false;
+
+      for (q_t i = 0; i < q; i++) {
+        m_rot[i].fill(0);
+        for (q_t j = 0; j < q; j++) {
+          for (q_t k = 0; k < q; k++) {
+            m_rot[i][j] += (*this)[j][i] * m[j][k];
+          }
         }
       }
-    }
 
-    return m2_rot;
+      return m_rot;
+    }
   }
-}
+
+};
+
 
 template <q_t q>
-orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, vector_t <q, double> v) {
+orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, const vector_t <q, double>& v) {
   orthogonal_t <q, double> ptr;
   ptr.is_reflection = true;
-  ptr.x = (double *)calloc(q, sizeof(double));
 
   double v2 = 0;
 
   for (q_t i = 0; i < q; i++) {
-    ptr.x[i] = gsl_ran_ugaussian(r);
-    v2 += ptr.x[i] * ptr.x[i];
+    ptr[0][i] = gsl_ran_ugaussian(r);
+    v2 += ptr[0][i] * ptr[0][i];
   }
 
   double mag_v = sqrt(v2);
 
   for (q_t i = 0; i < q; i++) {
-    ptr.x[i] /= mag_v;
+    ptr[0][i] /= mag_v;
   }
 
   return ptr;
 }
 
 template <q_t q>
-orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q, double> v0, double epsilon, unsigned int n) {
+orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const vector_t <q, double>& v0, double epsilon, unsigned int n) {
   orthogonal_t <q, double> m;
   m.is_reflection = true;
-  m.x = (double *)malloc(q * sizeof(double));
 
   vector_t <q, double> v;
 
@@ -191,22 +161,22 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q
   double m_dot_v = 0;
 
   for (q_t i = 0; i < q; i++) {
-    m.x[i] = gsl_ran_ugaussian(r);
-    m_dot_v += m.x[i] * v[i];
+    m[0][i] = gsl_ran_ugaussian(r);
+    m_dot_v += m[0][i] * v[i];
   }
 
   double v2 = 0;
   double factor = epsilon * gsl_ran_ugaussian(r);
 
   for (q_t i = 0; i < q; i++) {
-    m.x[i] = m.x[i] - m_dot_v * v[i] + factor * v[i];
-    v2 += pow(m.x[i], 2);
+    m[0][i] = m[0][i] - m_dot_v * v[i] + factor * v[i];
+    v2 += pow(m[0][i], 2);
   }
 
   double mag_v = sqrt(v2);
 
   for (q_t i = 0; i < q; i++) {
-    m.x[i] /= mag_v;
+    m[0][i] /= mag_v;
   }
 
   return m;
diff --git a/lib/state.h b/lib/state.h
index b57a9e2..eefa0cb 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -30,12 +30,11 @@ class state_t {
     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), g(D, L), T(T), J(J), H(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), g(D, L), T(T), R(), J(J), H(H) {
       nv = g.nv;
       ne = g.ne;
       g.add_ext();
       spins.resize(nv);
-      init (&R);
       E = - (double)ne * J(spins[0], spins[0]) - (double)nv * H(spins[0]);
       M = spins[0] * nv;
       last_cluster_size = 0;
@@ -52,10 +51,6 @@ class state_t {
         precomputed_sin[i] = sin(2 * M_PI * (double)i / (double)L);
       }
     }
-
-    ~state_t() {
-      free_spin(R);
-    }
 };
 
 
diff --git a/lib/symmetric.h b/lib/symmetric.h
index 41c8fb6..9e9b6e4 100644
--- a/lib/symmetric.h
+++ b/lib/symmetric.h
@@ -2,84 +2,50 @@
 #pragma once
 
 #include <stdlib.h>
+#include <array>
 #include "types.h"
 #include "potts.h"
 
 template <q_t q>
-class symmetric_t {
+class symmetric_t : public std::array<q_t, q> {
   public:
-    q_t *perm;
-};
-
-template <q_t q>
-void init(symmetric_t<q> *p) {
-  p->perm = (q_t *)malloc(q * sizeof(q_t));
-
-  for (q_t i = 0; i < q; i++) {
-    p->perm[i] = i;
-  }
-}
-
-template <q_t q>
-void free_spin(symmetric_t<q> p) {
-  free(p.perm);
-}
 
-template <q_t q>
-symmetric_t<q> copy(symmetric_t<q> x) {
-  symmetric_t<q> x2;
-  x2.perm = (q_t *)malloc(q * sizeof(q_t));
-
-  for (q_t i = 0; i < q; i++) {
-    x2.perm[i] = x.perm[i];
-  }
-
-  return x2;
-}
-
-template <q_t q>
-potts_t<q> act(symmetric_t<q> r, potts_t<q> s) {
-  potts_t<q> s2;
-  s2.x = r.perm[s.x];
-  return s2;
-}
-
-template <q_t q>
-symmetric_t<q> act(symmetric_t<q> r1, symmetric_t<q> r2) {
-  symmetric_t<q> r3;
-  r3.perm = (q_t *)malloc(q * sizeof(q_t));
-  for (q_t i = 0; i < q; i++) {
-    r3.perm[i] = r1.perm[r2.perm[i]];
-  }
-
-  return r3;
-}
+    symmetric_t() {
+      for (q_t i = 0; i < q; i++) {
+        (*this)[i] = i;
+      }
+    }
 
-template <q_t q>
-potts_t<q> act_inverse(symmetric_t<q> r, potts_t<q> s) {
-  potts_t<q> s2;
+    potts_t<q> act(const potts_t<q> &s) const {
+      return potts_t<q>((*this)[s.x]);
+    }
 
-  q_t i;
+    symmetric_t<q> act(const symmetric_t<q>& r) const {
+      symmetric_t<q> r_rot;
+      for (q_t i = 0; i < q; i++) {
+        r_rot[i] = (*this)[r[i]];
+      }
 
-  for (i = 0; i < q; i++) {
-    if (r.perm[i] == s.x) {
-      break;
+      return r_rot;
     }
-  }
 
-  s2.x = i;
+    potts_t<q> act_inverse(const potts_t<q>& s) const {
+      for (q_t i = 0; i < q; i++) {
+        if ((*this)[i] == s.x) {
+          return potts_t<q>(i);
+        }
+      }
 
-  return s2;
-}
+      exit(EXIT_FAILURE);
+    }
 
-template <q_t q>
-symmetric_t<q> act_inverse(symmetric_t<q> r1, symmetric_t<q> r2) {
-  symmetric_t<q> r3;
-  r3.perm = (q_t *)malloc(q * sizeof(q_t));
-  for (q_t i = 0; i < q; i++) {
-    r3.perm[r1.perm[i]] = r2.perm[i];
-  }
+    symmetric_t<q> act_inverse(const symmetric_t<q>& r) const {
+      symmetric_t<q> r_rot;
+      for (q_t i = 0; i < q; i++) {
+        r_rot[(*this)[i]] = r[i];
+      }
 
-  return r3;
-}
+      return r_rot;
+    }
+};
 
diff --git a/lib/wolff.h b/lib/wolff.h
index da4b7b6..a4a663c 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -16,7 +16,6 @@ void wolff(count_t N, state_t <R_t, X_t> *s, std::function <R_t(gsl_rng *, X_t s
     v_t v0 = gsl_rng_uniform_int(r, s->nv);
     R_t step = gen_R(r, s->spins[v0]);
     flip_cluster <R_t, X_t> (s, v0, step, r);
-    free_spin(step);
 
     measurements(s);
   }
diff --git a/lib/z2.h b/lib/z2.h
index dae1eb7..fd4722e 100644
--- a/lib/z2.h
+++ b/lib/z2.h
@@ -17,49 +17,37 @@
  *
  */
 
-struct z2_t { bool x; };
+class z2_t {
+  public:
+  bool x;
 
-void init(z2_t *p) {
-  p->x = false;
-}
+  z2_t() : x(false) {}
 
-void free_spin(z2_t p) {
-  // do nothing!
-}
+  z2_t(bool x) : x(x) {}
 
-z2_t copy(z2_t x) {
-  return x;
-}
-
-ising_t act(z2_t r, ising_t s) {
-  ising_t rs;
-
-  if (r.x) {
-    rs.x = !s.x;
-    return rs;
-  } else {
-    rs.x = s.x;
-    return rs;
+  ising_t act(const ising_t& s) {
+    if (x) {
+      return ising_t(!s.x);
+    } else {
+      return ising_t(s.x);
+    }
   }
-}
 
-z2_t act(z2_t r1, z2_t r2) {
-  z2_t r3;
+  z2_t act(const z2_t& r) {
+    if (x) {
+      return z2_t(!r.x);
+    } else {
+      return z2_t(r.x);
+    }
+  }
 
-  if (r1.x) {
-    r3.x = !r2.x;
-    return r3;
-  } else {
-    r3.x = r2.x;
-    return r3;
+  ising_t act_inverse(const ising_t& s) {
+    return this->act(s);
   }
-}
 
-ising_t act_inverse(z2_t r, ising_t s) {
-  return act(r, s);
-}
+  z2_t act_inverse(const z2_t& r) {
+    return this->act(r);
+  }
+};
 
-z2_t act_inverse(z2_t r1, z2_t r2) {
-  return act(r1, r2);
-}
 
diff --git a/src/analyze_correlations.cpp b/src/analyze_correlations.cpp
index e108f98..5f3dda0 100644
--- a/src/analyze_correlations.cpp
+++ b/src/analyze_correlations.cpp
@@ -1,6 +1,7 @@
 
 #include <types.h>
-#include <measurement.h>
+#include <cmath>
+#include <cstring>
 #include <stdio.h>
 #include <stdlib.h>
 #include <getopt.h>
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index 15a9553..5458860 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -189,7 +189,7 @@ int main(int argc, char *argv[]) {
     other_f = [&] (const On_t *s) {
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(L, 2); i++) {
-        vector_R_t v_tmp = act_inverse(s->R, s->spins[i]);
+        vector_R_t v_tmp = s->R.act_inverse(s->spins[i]);
         double thetai = fmod(2 * M_PI + theta(v_tmp), 2 * M_PI);
         double saturation = 0.7;
         double value = 0.9;
diff --git a/src/wolff_cgm.cpp b/src/wolff_cgm.cpp
index 10477bd..7613691 100644
--- a/src/wolff_cgm.cpp
+++ b/src/wolff_cgm.cpp
@@ -128,7 +128,7 @@ int main(int argc, char *argv[]) {
       double max_h = INT64_MIN;
       double min_h = INT64_MAX;
       for (v_t i = 0; i < pow(L, 2); i++) {
-        double cur_h = act_inverse(s->R, s->spins[i]).x;
+        double cur_h = (s->R.act_inverse(s->spins[i])).x;
         if (cur_h < min_h) {
           min_h = cur_h;
         }
@@ -138,7 +138,7 @@ int main(int argc, char *argv[]) {
       }
 
       for (v_t i = 0; i < pow(L, 2); i++) {
-        double cur_h = act_inverse(s->R, s->spins[i]).x;
+        double cur_h = (s->R.act_inverse(s->spins[i])).x;
         double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h));
         glColor3f(mag, mag, mag);
         glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
diff --git a/src/wolff_clock.cpp b/src/wolff_clock.cpp
index e186c44..bc2c5d1 100644
--- a/src/wolff_clock.cpp
+++ b/src/wolff_clock.cpp
@@ -18,7 +18,7 @@
 #include <rand.h>
 #include <wolff.h>
 
-typedef state_t <dihedral_t<q_t,POTTSQ>, potts_t<POTTSQ>> sim_t;
+typedef state_t <dihedral_t<POTTSQ>, potts_t<POTTSQ>> sim_t;
 
 int main(int argc, char *argv[]) {
 
@@ -88,11 +88,11 @@ int main(int argc, char *argv[]) {
   };
 
   // initialize state object
-  state_t <dihedral_t<q_t,POTTSQ>, potts_t<POTTSQ>> s(D, L, T, Z, B);
+  state_t <dihedral_t<POTTSQ>, potts_t<POTTSQ>> s(D, L, T, Z, B);
 
   // define function that generates self-inverse rotations
-  std::function <dihedral_t<q_t,POTTSQ>(gsl_rng *, potts_t<POTTSQ>)> gen_R = [] (gsl_rng *r, potts_t<POTTSQ> v) -> dihedral_t<q_t,POTTSQ> {
-    dihedral_t<q_t,POTTSQ> rot;
+  std::function <dihedral_t<POTTSQ>(gsl_rng *, potts_t<POTTSQ>)> gen_R = [] (gsl_rng *r, potts_t<POTTSQ> v) -> dihedral_t<POTTSQ> {
+    dihedral_t<POTTSQ> rot;
     rot.is_reflection = true;
     q_t x = gsl_rng_uniform_int(r, POTTSQ - 1);
     rot.x = (2 * v.x + x + 1) % POTTSQ;
@@ -127,7 +127,7 @@ int main(int argc, char *argv[]) {
       average_M += (double)s->M[0] / (double)N / (double)s->nv;
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(L, 2); i++) {
-        potts_t<POTTSQ> tmp_s = act_inverse(s->R, s->spins[i]);
+        potts_t<POTTSQ> tmp_s = s->R.act_inverse(s->spins[i]);
         glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ));
         glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
       }
diff --git a/src/wolff_dgm.cpp b/src/wolff_dgm.cpp
index 6271b79..2583704 100644
--- a/src/wolff_dgm.cpp
+++ b/src/wolff_dgm.cpp
@@ -133,7 +133,7 @@ int main(int argc, char *argv[]) {
       int64_t max_h = INT64_MIN;
       int64_t min_h = INT64_MAX;
       for (v_t i = 0; i < pow(L, 2); i++) {
-        int64_t cur_h = act_inverse(s->R, s->spins[i]).x;
+        int64_t cur_h = (s->R.act_inverse(s->spins[i])).x;
         if (cur_h < min_h) {
           min_h = cur_h;
         }
@@ -143,7 +143,7 @@ int main(int argc, char *argv[]) {
       }
 
       for (v_t i = 0; i < pow(L, 2); i++) {
-        int64_t cur_h = act_inverse(s->R, s->spins[i]).x;
+        int64_t cur_h = (s->R.act_inverse(s->spins[i])).x;
         double mag = ((double)(cur_h - min_h)) / ((double)(max_h - min_h));
         glColor3f(mag, mag, mag);
         glRecti(i / L, i % L, (i / L) + 1, (i % L) + 1);
diff --git a/src/wolff_ising.cpp b/src/wolff_ising.cpp
index 93b7ea4..410e046 100644
--- a/src/wolff_ising.cpp
+++ b/src/wolff_ising.cpp
@@ -120,10 +120,8 @@ int main(int argc, char *argv[]) {
   state_t <z2_t, ising_t> s(D, L, T, Z, B);
 
   // define function that generates self-inverse rotations
-  std::function <z2_t(gsl_rng *, ising_t)> gen_R = [] (gsl_rng *, ising_t s) -> z2_t {
-    z2_t rot;
-    rot.x = true;
-    return rot;
+  std::function <z2_t(gsl_rng *, ising_t)> gen_R = [] (gsl_rng *, const ising_t& s) -> z2_t {
+    return z2_t(true);
   };
 
   FILE **outfiles = measure_setup_files(measurement_flags, timestamp);
diff --git a/src/wolff_potts.cpp b/src/wolff_potts.cpp
index cdc4c07..07663d5 100644
--- a/src/wolff_potts.cpp
+++ b/src/wolff_potts.cpp
@@ -118,7 +118,6 @@ int main(int argc, char *argv[]) {
   // define function that generates self-inverse rotations
   std::function <symmetric_t<POTTSQ>(gsl_rng *, potts_t<POTTSQ>)> gen_R = [] (gsl_rng *r, potts_t<POTTSQ> v) -> symmetric_t<POTTSQ> {
     symmetric_t<POTTSQ> rot;
-    init(&rot);
 
     q_t j = gsl_rng_uniform_int(r, POTTSQ - 1);
     q_t swap_v;
@@ -128,8 +127,8 @@ int main(int argc, char *argv[]) {
       swap_v = j + 1;
     }
 
-    rot.perm[v.x] = swap_v;
-    rot.perm[swap_v] = v.x;
+    rot[v.x] = swap_v;
+    rot[swap_v] = v.x;
 
     return rot;
   };
@@ -158,7 +157,7 @@ int main(int argc, char *argv[]) {
     other_f = [] (const sim_t *s) {
       glClear(GL_COLOR_BUFFER_BIT);
       for (v_t i = 0; i < pow(s->L, 2); i++) {
-        potts_t<POTTSQ> tmp_s = act_inverse(s->R, s->spins[i]);
+        potts_t<POTTSQ> tmp_s = s->R.act_inverse(s->spins[i]);
         glColor3f(hue_to_R(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_G(tmp_s.x * 2 * M_PI / POTTSQ), hue_to_B(tmp_s.x * 2 * M_PI / POTTSQ));
         glRecti(i / s->L, i % s->L, (i / s->L) + 1, (i % s->L) + 1);
       }