From 802b63ddf121b520db7942fe330cce6004fbeb6d Mon Sep 17 00:00:00 2001
From: Jaron Kent-Dobias <jaron@kent-dobias.com>
Date: Wed, 25 Jul 2018 16:22:50 -0400
Subject: got everyone recording data, and fixed huge bug in the updating of
 ReF and ImF

---
 lib/cluster.h    |  8 ++++----
 lib/orthogonal.h |  1 +
 lib/potts.h      | 47 +++++++++++++++++++++++------------------------
 3 files changed, 28 insertions(+), 28 deletions(-)

(limited to 'lib')

diff --git a/lib/cluster.h b/lib/cluster.h
index b8c98e5..3261969 100644
--- a/lib/cluster.h
+++ b/lib/cluster.h
@@ -77,11 +77,11 @@ void flip_cluster(state_t <R_t, X_t> *state, v_t v0, R_t r, gsl_rng *rand) {
           for (D_t i = 0; i < state->D; i++) {
             L_t x = (non_ghost / (v_t)pow(state->L, state->D - i - 1)) % state->L;
 
-            add(&(state->ReF[i]), -state->precomputed_cos[i], rs_old);
-            add(&(state->ReF[i]),  state->precomputed_cos[i], rs_new);
+            add(&(state->ReF[i]), -state->precomputed_cos[x], rs_old);
+            add(&(state->ReF[i]),  state->precomputed_cos[x], rs_new);
 
-            add(&(state->ImF[i]), -state->precomputed_sin[i], rs_old);
-            add(&(state->ImF[i]),  state->precomputed_sin[i], rs_new);
+            add(&(state->ImF[i]), -state->precomputed_sin[x], rs_old);
+            add(&(state->ImF[i]),  state->precomputed_sin[x], rs_new);
           }
 
           free_spin (rs_old);
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index ce2d300..0ba5eee 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -8,6 +8,7 @@
 
 #include "state.h"
 #include "types.h"
+#include "vector.h"
 
 template <q_t q, class T>
 struct orthogonal_t { bool is_reflection; T *x; };
diff --git a/lib/potts.h b/lib/potts.h
index 732a76f..e61e4e1 100644
--- a/lib/potts.h
+++ b/lib/potts.h
@@ -4,6 +4,7 @@
 #include <stdio.h>
 
 #include "types.h"
+#include "vector.h"
 
 /* The following is the minimum definition of a spin class.
  *
@@ -30,8 +31,8 @@ class potts_t {
   public:
     q_t x;
 
-    typedef int *M_t;
-    typedef double *F_t;
+    typedef vector_t<q, int> M_t;
+    typedef vector_t<q, double> F_t;
 };
 
 template <q_t q>
@@ -44,56 +45,54 @@ void free_spin(potts_t <q> s) {
   // do nothing!
 }
 
-void free_spin(int *s) {
-  free(s);
-}
-
-void free_spin(double *s) {
-  free(s);
-}
-
 template <q_t q>
 potts_t <q> copy(potts_t <q> s) {
   return s;
 }
 
 template <q_t q>
-void add(typename potts_t<q>::M_t *s1, int a, potts_t <q> s2) {
-  (*s1)[s2.x] += a;
+void add(vector_t<q, int> *s1, int a, potts_t <q> s2) {
+  (s1->x)[s2.x] += a;
 }
 
 template <q_t q>
-void add(typename potts_t<q>::F_t *s1, double a, potts_t <q> s2) {
-  (*s1)[s2.x] += a;
+void add(vector_t<q, double> *s1, double a, potts_t <q> s2) {
+  (s1->x)[s2.x] += a;
 }
 
 template <q_t q>
-typename potts_t<q>::M_t scalar_multiple(int factor, potts_t <q> s) {
-  int *M = (int *)calloc(q, sizeof(int));
-  M[s.x] += factor;
+vector_t<q, int> scalar_multiple(int factor, potts_t <q> s) {
+  vector_t<q, int> M;
+  M.x = (int *)calloc(q, sizeof(int));
+  M.x[s.x] += factor;
   return M;
 }
 
 template <q_t q>
-typename potts_t<q>::F_t scalar_multiple(double factor, potts_t <q> s) {
-  double *F = (double *)calloc(q, sizeof(double));
-  F[s.x] += factor;
+vector_t<q, double> scalar_multiple(double factor, potts_t <q> s) {
+  vector_t<q, double> F;
+  F.x = (double *)calloc(q, sizeof(double));
+  F.x[s.x] += factor;
   return F;
 }
 
+// we could inherit norm_squared from vector.h, but convention dictates that
+// potts norms be changed by a constant factor
 template <q_t q>
-double norm_squared(typename potts_t<q>::F_t s) {
+double norm_squared(vector_t<q, double> s) {
   double total = 0;
   for (q_t i = 0; i < q; i++) {
-    total += pow(s[i], 2);
+    total += pow(s.x[i], 2);
   }
 
   return total * (double)q / ((double)q - 1.0);
 }
 
+// we could inherit write_magnetization from vector.h, but since M.x must sum
+// to nv we don't need to write the last element
 template <q_t q>
-void write_magnetization(typename potts_t<q>::M_t M, FILE *outfile) {
-  fwrite(&M, sizeof(int), q, outfile);
+void write_magnetization(vector_t<q, int> M, FILE *outfile) {
+  fwrite(M.x, sizeof(int), q - 1, outfile);
 }
 
 // knock yourself out
-- 
cgit v1.2.3-70-g09d2