1 files changed, 90 insertions, 17 deletions

diff --git a/space_wolff.hpp b/space_wolff.hpp
index b60e274..9a4b4b4 100644
--- a/space_wolff.hpp
+++ b/space_wolff.hpp
@@ -237,6 +237,10 @@ class quantity {
     double serr() const {
       return sqrt(this->σ());
     }
+
+    unsigned num_added() const {
+      return n - wait;
+    }
 };
 
 template <class U, int D, class state>
@@ -249,16 +253,87 @@ class model {
     std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> neighbors;
     std::function<double(spin<U, D, state>, spin<U, D, state>)> Z;
     std::function<double(spin<U, D, state>)> B;
+    std::vector<matrix<U, D>> mats;
+    std::vector<vector<U, D>> steps;
     long double E;
     quantity Eq;
     quantity Cq;
 
+    void one_sequences(std::list<std::array<double, D>>& sequences, unsigned level) {
+      if (level > 0) {
+        unsigned new_level = level - 1;
+        unsigned old_length = sequences.size();
+        for (std::array<double, D>& sequence : sequences) {
+          std::array<double, D> new_sequence = sequence;
+          new_sequence[new_level] = -1;
+          sequences.push_front(new_sequence);
+        }
+        one_sequences(sequences, new_level);
+      }
+    }
+
     model(unsigned L, std::function<double(spin<U, D, state>, spin<U, D, state>)> Z,
         std::function<double(spin<U, D, state>)> B,
         std::function<std::set<unsigned>(model<U, D, state>&, unsigned, spin<U, D, state>)> ns) : 
       L(L), s0(L), dict(L), neighbors(ns), Z(Z), B(B), Eq(1000, 1000), Cq(1000, 1000) {
+        std::array<double, D> ini_sequence;
+        ini_sequence.fill(1);
+        std::list<std::array<double, D>> sequences;
+        sequences.push_back(ini_sequence);
+
+        one_sequences(sequences, D);
+
+        sequences.pop_front(); // don't want the identity matrix!
+
+        for (std::array<double, D> sequence : sequences) {
+          matrix<U, D> m;
+          for (unsigned i = 0; i < D; i++) {
+            for (unsigned j = 0; j < D; j++) {
+              if (i == j) {
+                m(i, j) = sequence[i];
+              } else {
+                m(i, j) = 0;
+              }
+            }
+          }
+
+          mats.push_back(m);
+
+          vector<U, D> v;
+          for (unsigned i = 0; i < D; i++) {
+            if (sequence[i] == 1) {
+              v(i) = 0;
+            } else {
+              v(i) = L / 2;
+            }
+          }
+
+          steps.push_back(v);
         }
 
+        for (unsigned i = 0; i < D; i++) {
+          for (unsigned j = 0; j < D; j++) {
+            if (i != j) {
+              matrix<U, D> m;
+              for (unsigned k = 0; k < D; k++) {
+                for (unsigned l = 0; l < D; l++) {
+                  if ((k == i && l == j) || (k == j && l == i)) {
+                    if (i < j) {
+                      m(k, l) = 1;
+                    } else {
+                      m(k, l) = -1;
+                    }
+                  } else {
+                    m(k, l) = 0;
+                  }
+                }
+              }
+              mats.push_back(m);
+            }
+          }
+        }
+      }
+
     void update_energy() {
       E = 0;
       for (unsigned i = 0; i < s.size(); i++) {
@@ -347,31 +422,29 @@ class model {
       std::uniform_real_distribution<double> t_dist(0, L);
       std::uniform_int_distribution<unsigned> r_dist(0, D - 1);
       std::uniform_int_distribution<unsigned> ind_dist(0, s.size() - 1);
+      std::uniform_int_distribution<unsigned> coin(0, mats.size() + steps.size() - 1);
 
       for (unsigned i = 0; i < N; i++) {
         vector<U, D> t;
         matrix<U, D> m;
-        for (unsigned j = 0; j < D; j++) {
-          t(j) = (U)t_dist(rng);
-        }
-
-        unsigned flip_D1 = r_dist(rng);
-        unsigned flip_D2 = r_dist(rng);
-
-        for (unsigned j = 0; j < D; j++) {
-          for (unsigned k = 0; k < D; k++) {
-            if ((j == flip_D1 && k == flip_D2) || (j == flip_D2 && k == flip_D1)) {
-              if (flip_D1 <= flip_D2) {
-                m(j, k) = -1;
-              } else {
+        unsigned flip = coin(rng);
+        if (flip < mats.size()) {
+          for (unsigned j = 0; j < D; j++) {
+            t(j) = (U)t_dist(rng);
+          }
+          m = mats[flip];
+        } else {
+          for (unsigned j = 0; j < D; j++) {
+            for (unsigned k = 0; k < D; k++) {
+              if (j == k) {
                 m(j, k) = 1;
+              } else {
+                m(j, k) = 0;
               }
-            } else if ((j == k && j != flip_D1) && j != flip_D2) {
-              m(j, k) = 1;
-            } else {
-              m(j, k) = 0;
             }
           }
+
+          t = steps[flip - mats.size()];
         }
 
         euclidean<U, D> g(L, t, m);