2 files changed, 39 insertions, 0 deletions

diff --git a/lib/wolff.h b/lib/wolff.h
index 428ce85..82ccb9e 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -31,6 +31,17 @@ typedef struct {
   int32_t dHb;
 } cluster_t;
 
+typedef struct {
+  uint64_t n;
+  double x;
+  double dx;
+  double x2;
+  double m2;
+  double m4;
+  double c;
+  double dc;
+} meas_t;
+
 int32_t sign(double x);
 
 cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x,
@@ -40,3 +51,8 @@ graph_t *graph_add_ext(const graph_t *g);
 
 uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
                     double *ps);
+
+void update_meas(meas_t *m, double x);
+
+double add_to_avg(double mx, double x, uint64_t n);
+
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index 7886bab..59d04fc 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -152,3 +152,26 @@ uint32_t wolff_step(double T, double H, ising_state_t *s, gsl_rng *r,
 
   return n_flips;
 }
+
+double add_to_avg(double mx, double x, uint64_t n) {
+  return mx * (n / (n + 1.)) + x * 1. / (n + 1.);
+}
+
+void update_meas(meas_t *m, double x) {
+  uint64_t n = m->n;
+
+  m->x = add_to_avg(m->x, x, n);
+  m->x2 = add_to_avg(m->x2, pow(x, 2), n);
+
+  m->m2 = add_to_avg(m->m2, pow(x - m->x, 2), n);
+  m->m4 = add_to_avg(m->m4, pow(x - m->x, 4), n);
+
+  if (n > 1) {
+    double s2 = n / (n - 1.) * (m->x2 - pow(m->x, 2));
+    m->dx = sqrt(s2 / n);
+    m->c = s2;
+    m->dc = sqrt((m->m4 - (n - 3.)/(n - 1.) * pow(m->m2, 2)) / n);
+  }
+
+  (m->n)++;
+}