5 files changed, 45 insertions, 59 deletions
diff --git a/lib/include/wolff/cluster.hpp b/lib/include/wolff/cluster.hpp
index b66f367..3912446 100644
--- a/lib/include/wolff/cluster.hpp
+++ b/lib/include/wolff/cluster.hpp
@@ -68,9 +68,8 @@ void system<R_t, X_t>::flip_cluster(v_t i0, const R_t& r,
           dE = B(s0s_old) - B(s0s_new);
 #endif
 
-#ifdef WOLFF_FINITE_STATES
-          p = finite_states_Bp[finite_states_enum(s0s_old)]
-                              [finite_states_enum(s0s_new)];
+#ifdef WOLFF_USE_FINITE_STATES
+          p = Bp[s0s_old.enumerate()][s0s_new.enumerate()];
 #endif
 
           // run measurement hooks for encountering a ghost bond
@@ -84,17 +83,15 @@ void system<R_t, X_t>::flip_cluster(v_t i0, const R_t& r,
           dE = Z(s[i], s[j]) - Z(si_new, s[j]);
 #endif
 
-#ifdef WOLFF_FINITE_STATES
-          p = finite_states_Zp[finite_states_enum(s[i])]
-                              [finite_states_enum(si_new)]
-                              [finite_states_enum(s[j])];
+#ifdef WOLFF_USE_FINITE_STATES
+          p = Zp[s[i].enumerate()][si_new.enumerate()][s[j].enumerate()];
 #endif
 
           // run measurement hooks for encountering a plain bond
           A.plain_bond_visited(*this, i, si_new, j, dE);
         }
 
-#ifndef FINITE_STATES
+#ifndef WOLFF_USE_FINITE_STATES
         p = 1.0 - exp(-dE / T);
 #endif
 
diff --git a/lib/include/wolff/finite_states.hpp b/lib/include/wolff/finite_states.hpp
deleted file mode 100644
index 7e54eff..0000000
--- a/lib/include/wolff/finite_states.hpp
+++ /dev/null
@@ -1,39 +0,0 @@
-
-
-#ifndef WOLFF_FINITE_STATES
-#define WOLFF_FINITE_STATES
-
-#include <cmath>
-#include <array>
-
-#include "system.hpp"
-
-namespace wolff {
-
-std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Zp;
-#ifndef WOLFF_NO_FIELD
-std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> finite_states_Bp;
-#endif
-
-template <class R_t, class X_t>
-void finite_states_init(const system<R_t, X_t>& S) {
-#ifndef WOLFF_NO_FIELD
-  for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
-    for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
-      finite_states_Bp[i][j] = 1.0 - exp(-(S.B(finite_states_possible[i]) - S.B(finite_states_possible[j])) / S.T);
-    }
-  }
-#endif
-  for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
-    for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
-      for (q_t k = 0; k < WOLFF_FINITE_STATES_N; k++) {
-        finite_states_Zp[i][j][k] = 1.0 - exp(-(S.Z(finite_states_possible[i], finite_states_possible[k]) - S.Z(finite_states_possible[j], finite_states_possible[k])) / S.T);
-      }
-    }
-  }
-}
-
-}
-
-#endif
-
diff --git a/lib/include/wolff/models/ising.hpp b/lib/include/wolff/models/ising.hpp
index 824c063..39d0ad6 100644
--- a/lib/include/wolff/models/ising.hpp
+++ b/lib/include/wolff/models/ising.hpp
@@ -2,6 +2,8 @@
 #ifndef WOLFF_MODELS_ISING
 #define WOLFF_MODELS_ISING
 
+#define WOLFF_FINITE_STATES_N 2
+
 #include "../types.h"
 #include "../system.hpp"
 
@@ -14,7 +16,7 @@ class ising_t {
     ising_t() : x(false) {}
 
     ising_t(bool x) : x(x) {} 
-    ising_t(int x) : x((bool)x) {}
+    ising_t(q_t x) : x((bool)x) {}
 
     ising_t act(const ising_t& s) const {
       if (x) {
@@ -66,6 +68,10 @@ class ising_t {
         return -1;
       }
     }
+
+    q_t enumerate() const {
+      return (q_t)x;
+    }
 };
 
 inline ising_t::M_t operator*(v_t a, const ising_t& s) {
@@ -80,10 +86,6 @@ ising_t gen_ising(std::mt19937&, const system<ising_t, ising_t>&, v_t) {
   return ising_t(true);
 };
 
-#define WOLFF_FINITE_STATES_N 2
-const ising_t finite_states_possible[2] = {ising_t(0), ising_t(1)};
-q_t finite_states_enum(const ising_t& state) { return (q_t)state.x; }
-
 }
 
 #endif
diff --git a/lib/include/wolff/models/potts.hpp b/lib/include/wolff/models/potts.hpp
index ce419ce..3d0a0e9 100644
--- a/lib/include/wolff/models/potts.hpp
+++ b/lib/include/wolff/models/potts.hpp
@@ -44,6 +44,10 @@ class potts_t {
 
       return result;
     }
+
+    q_t enumerate() const {
+      return x;
+    }
 };
 
 template<q_t q>
@@ -56,10 +60,5 @@ inline typename potts_t<q>::F_t operator*(double a, const potts_t<q>& s) {
   return s * a;
 }
 
-#define WOLFF_FINITE_STATES_N WOLFF_POTTSQ
-const potts_t<WOLFF_POTTSQ> finite_states_possible[256] = {{0}, {1}, {2}, {3}, {4}, {5}, {6}, {7}, {8}, {9}, {11}, {12}, {13}, {14}, {15}, {16}, {17}, {18}, {19}, {20}, {21}, {22}, {23}, {24}, {25}, {26}, {27}, {28}, {29}, {30}, {31}, {32}, {33}, {34}, {35}, {36}, {37}, {38}, {39}, {40}, {41}, {42}, {43}, {44}, {45}, {46}, {47}, {48}, {49}, {50}, {51}, {52}, {53}, {54}, {55}, {56}, {57}, {58}, {59}, {60}, {61}, {62}, {63}, {64}, {65}, {66}, {67}, {68}, {69}, {70}, {71}, {72}, {73}, {74}, {75}, {76}, {77}, {78}, {79}, {80}, {81}, {82}, {83}, {84}, {85}, {86}, {87}, {88}, {89}, {90}, {91}, {92}, {93}, {94}, {95}, {96}, {97}, {98}, {99}, {100}, {101}, {102}, {103}, {104}, {105}, {106}, {107}, {108}, {109}, {110}, {111}, {112}, {113}, {114}, {115}, {116}, {117}, {118}, {119}, {120}, {121}, {122}, {123}, {124}, {125}, {126}, {127}, {128}, {129}, {130}, {131}, {132}, {133}, {134}, {135}, {136}, {137}, {138}, {139}, {140}, {141}, {142}, {143}, {144}, {145}, {146}, {147}, {148}, {149}, {150}, {151}, {152}, {153}, {154}, {155}, {156}, {157}, {158}, {159}, {160}, {161}, {162}, {163}, {164}, {165}, {166}, {167}, {168}, {169}, {170}, {171}, {172}, {173}, {174}, {175}, {176}, {177}, {178}, {179}, {180}, {181}, {182}, {183}, {184}, {185}, {186}, {187}, {188}, {189}, {190}, {191}, {192}, {193}, {194}, {195}, {196}, {197}, {198}, {199}, {200}, {201}, {202}, {203}, {204}, {205}, {206}, {207}, {208}, {209}, {210}, {211}, {212}, {213}, {214}, {215}, {216}, {217}, {218}, {219}, {220}, {221}, {222}, {223}, {224}, {225}, {226}, {227}, {228}, {229}, {230}, {231}, {232}, {233}, {234}, {235}, {236}, {237}, {238}, {239}, {240}, {241}, {242}, {243}, {244}, {245}, {246}, {247}, {248}, {249}, {250}, {251}, {252}, {253}, {254}, {255}};
-template <q_t q>
-q_t finite_states_enum(const potts_t<q>& state) { return (q_t)state.x; }
-
 }
 
diff --git a/lib/include/wolff/system.hpp b/lib/include/wolff/system.hpp
index 32ad38e..4aa8ee4 100644
--- a/lib/include/wolff/system.hpp
+++ b/lib/include/wolff/system.hpp
@@ -39,6 +39,13 @@ class system {
 #endif
 #endif
 
+#ifdef WOLFF_USE_FINITE_STATES
+    std::array<std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Zp;
+#ifndef WOLFF_NO_FIELD
+    std::array<std::array<double, WOLFF_FINITE_STATES_N>, WOLFF_FINITE_STATES_N> Bp;
+#endif
+#endif
+
     system(graph g, double T,
 #ifdef WOLFF_BOND_DEPENDENCE
         std::function <double(v_t, const X_t&, v_t, const X_t&)> Z
@@ -63,13 +70,33 @@ class system {
 #ifndef WOLFF_NO_FIELD
       G.add_ghost();
 #endif
-#ifdef WOLFF_FINITE_STATES
-      finite_states_init(*this);
+#ifdef WOLFF_USE_FINITE_STATES
+      this->finite_states_init();
 #endif
     }
 
     void flip_cluster(v_t, const R_t&, std::mt19937&, measurement<R_t, X_t>&);
     void run_wolff(N_t, std::function <R_t(std::mt19937&, const system<R_t, X_t>&, v_t)> r_gen, measurement<R_t, X_t>& A, std::mt19937& rng);
+
+#ifdef WOLFF_USE_FINITE_STATES
+    void finite_states_init() {
+#ifndef WOLFF_NO_FIELD
+      for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+        for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+          Bp[i][j] = 1.0 - exp(-(B(X_t(i)) - B(X_t(j))) / T);
+        }
+      }
+#endif
+      for (q_t i = 0; i < WOLFF_FINITE_STATES_N; i++) {
+        for (q_t j = 0; j < WOLFF_FINITE_STATES_N; j++) {
+          for (q_t k = 0; k < WOLFF_FINITE_STATES_N; k++) {
+            Zp[i][j][k] = 1.0 - exp(-(Z(X_t(i), X_t(k)) - Z(X_t(j), X_t(k))) / T);
+          }
+        }
+      }
+    }
+#endif
+
 };
 
 }