added visualization, and started potts

5 files changed, 209 insertions, 2 deletions

diff --git a/lib/ising.h b/lib/ising.h
index 4ad88f4..b4856c3 100644
--- a/lib/ising.h
+++ b/lib/ising.h
@@ -76,6 +76,14 @@ int scalar_multiple(int factor, ising_t s) {
   }
 }
 
+double scalar_multiple(double factor, ising_t s) {
+  if (s.x) {
+    return -factor;
+  } else {
+    return factor;
+  }
+}
+
 double norm_squared(double s) {
   return pow(s, 2);
 }
diff --git a/lib/potts.h b/lib/potts.h
new file mode 100644
index 0000000..e7f0899
--- /dev/null
+++ b/lib/potts.h
@@ -0,0 +1,100 @@
+#pragma once
+
+#include <cmath>
+#include <stdio.h>
+
+#include "types.h"
+
+/* The following is the minimum definition of a spin class.
+ *
+ * The class must contain an M_t and an F_t for holding the sum of an
+ * integer number of spins and a double-weighted number of spins,
+ * respectively.
+ *
+ * void init(X_t *p);
+ * void free_spin(X_t p);
+ * void free_spin(M_t p);
+ * void free_spin(F_t p);
+ * X_t copy(X_t x);
+ * void add(M_t *x1, int factor, X_t x2);
+ * void add(F_t *x1, double factor, X_t x2);
+ * M_t scalar_multiple(int factor, X_t x);
+ * F_t scalar_multiple(double factor, X_t x);
+ * double norm_squared(F_t x);
+ * void write_magnetization(M_t M, FILE *outfile);
+ *
+ */
+
+template <q_t q>
+class potts_t {
+  public:
+    q_t x;
+
+    typedef int *M_t;
+    typedef double *F_t;
+};
+
+template <q_t q>
+void init(potts_t <q> *p) {
+  p->x = 0;
+}
+
+template <q_t q>
+void free_spin(potts_t <q> s) {
+  // do nothing!
+}
+
+template <q_t q>
+void free_spin(typename potts_t<q>::M_t s) {
+  free(s);
+}
+
+template <q_t q>
+void free_spin(typename potts_t<q>::F_t 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;
+}
+
+template <q_t q>
+void add(typename potts_t<q>::F_t s1, double a, potts_t <q> s2) {
+  s1[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;
+  return M;
+}
+
+template <q_t q>
+typename potts_t<q>::F_t scalar_multiple(double factor, potts_t <q> s) {
+  int *F = (double *)calloc(q, sizeof(double));
+  M[s.x] += factor;
+  return M;
+}
+
+template <q_t q>
+double norm_squared(typename potts<q>::F_t s) {
+  double total = 0;
+  for (q_t i = 0; i < q; i++) {
+    total += pow(s[i], 2);
+  }
+
+  return total * (double)q / ((double)q - 1.0);
+}
+
+template <q_t q>
+void write_magnetization(typename potts_t<q>::M_t M, FILE *outfile) {
+  fwrite(&M, sizeof(int), q, outfile);
+}
+
diff --git a/lib/state.h b/lib/state.h
index 76d3e5a..8630810 100644
--- a/lib/state.h
+++ b/lib/state.h
@@ -47,8 +47,8 @@ class state_t {
       ReF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t));
       ImF = (typename X_t::F_t *)malloc(D * sizeof(typename X_t::F_t));
       for (D_t i = 0; i < D; i++) {
-        ReF[i] = scalar_multiple(0, spins[0]);
-        ImF[i] = scalar_multiple(0, spins[0]);
+        ReF[i] = scalar_multiple(0.0, spins[0]);
+        ImF[i] = scalar_multiple(0.0, spins[0]);
       }
       precomputed_cos = (double *)malloc(L * sizeof(double));
       precomputed_sin = (double *)malloc(L * sizeof(double));
diff --git a/lib/symmetric.h b/lib/symmetric.h
index c71521d..0b292a6 100644
--- a/lib/symmetric.h
+++ b/lib/symmetric.h
@@ -5,6 +5,10 @@
 
 #include "types.h"
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 q_t *symmetric_compose(q_t q, const q_t *g1, const q_t *g2);
 
 q_t symmetric_act(const q_t *g, q_t s);
@@ -13,3 +17,87 @@ q_t *symmetric_invert(q_t q, const q_t *g);
 
 q_t *symmetric_gen_transformations(q_t q);
 
+#ifdef __cplusplus
+}
+#endif
+
+#ifdef __cplusplus
+template <q_t q>
+class symmetric_t {
+  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_t> copy(symmetric_t<q_t> 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;
+}
+
+template <q_t q>
+potts_t<q> act_inverse(symmetric_t<q> r, potts_t<q> s) {
+  potts_t<q> s2;
+
+  q_t i;
+
+  for (i = 0; i < q; i++) {
+    if (r.perm[i] == s.x) {
+      break;
+    }
+  }
+
+  s2.x = i;
+
+  return s2;
+}
+
+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];
+  }
+
+  return r3;
+}
+#endif
+
diff --git a/lib/vector.h b/lib/vector.h
index c478618..2fe6ab8 100644
--- a/lib/vector.h
+++ b/lib/vector.h
@@ -80,6 +80,17 @@ vector_t <q, T> scalar_multiple(int a, vector_t <q, T> v) {
 }
 
 template <q_t q, class T>
+vector_t <q, T> scalar_multiple(double 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>
 double norm_squared (vector_t <q, T> v) {
   double tmp = 0;