summaryrefslogtreecommitdiff
path: root/lib
diff options
context:
space:
mode:
Diffstat (limited to 'lib')
-rw-r--r--lib/convex.c8
-rw-r--r--lib/convex.h6
-rw-r--r--lib/graph.c67
-rw-r--r--lib/graph.h18
-rw-r--r--lib/rand.c20
-rw-r--r--lib/rand.h7
-rw-r--r--lib/stack.c (renamed from lib/queue.c)17
-rw-r--r--lib/stack.h (renamed from lib/queue.h)9
-rw-r--r--lib/tree.c81
-rw-r--r--lib/tree.h22
-rw-r--r--lib/types.h16
-rw-r--r--lib/wolff.h44
-rw-r--r--lib/wolff_tools.c300
13 files changed, 321 insertions, 294 deletions
diff --git a/lib/convex.c b/lib/convex.c
index 096504d..df506f9 100644
--- a/lib/convex.c
+++ b/lib/convex.c
@@ -5,7 +5,7 @@ double slope(point_t *P, point_t *Q) {
return (Q->y - P->y) / ((double)(Q->x) - (double)(P->x));
}
-double *get_convex_minorant(uint64_t n, double *Gammas) {
+double *get_convex_minorant(count_t n, double *Gammas) {
if (n < 2) {
return Gammas;
}
@@ -17,7 +17,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) {
list_t *pos = L;
- for (uint64_t i = 1; i < n; i++) {
+ for (count_t i = 1; i < n; i++) {
pos->next = (list_t *)calloc(1, sizeof(list_t));
pos->next->p = (point_t *)calloc(1, sizeof(point_t));
pos->next->p->x = i;
@@ -69,7 +69,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) {
double *g = (double *)calloc(n + 1, sizeof(double));
double rho = 0;
- for (uint64_t i = 0; i < n + 1; i++) {
+ for (count_t i = 0; i < n + 1; i++) {
if (i > pos->next->p->x) {
pos = pos->next;
}
@@ -87,7 +87,7 @@ double *get_convex_minorant(uint64_t n, double *Gammas) {
}
}
- for (uint64_t i = 0; i < n + 1; i++) {
+ for (count_t i = 0; i < n + 1; i++) {
g[i] += rho / 2;
}
diff --git a/lib/convex.h b/lib/convex.h
index 29f4ae0..5a405d4 100644
--- a/lib/convex.h
+++ b/lib/convex.h
@@ -6,8 +6,10 @@
#include <stdlib.h>
#include <string.h>
+#include "types.h"
+
typedef struct {
- uint64_t x;
+ count_t x;
double y;
} point_t;
@@ -17,5 +19,5 @@ typedef struct list_tag {
point_t *p;
} list_t;
-double *get_convex_minorant(uint64_t n, double *Gammas);
+double *get_convex_minorant(count_t n, double *Gammas);
diff --git a/lib/graph.c b/lib/graph.c
new file mode 100644
index 0000000..a1edf13
--- /dev/null
+++ b/lib/graph.c
@@ -0,0 +1,67 @@
+
+#include "graph.h"
+
+graph_t *graph_create_square(D_t D, L_t L) {
+ v_t nv = pow(L, D);
+ v_t ne = D * nv;
+
+ v_t *v_i = (v_t *)malloc((nv + 1) * sizeof(v_t));
+
+ for (v_t i = 0; i < nv + 1; i++) {
+ v_i[i] = 2 * D * i;
+ }
+
+ v_t *v_adj = (v_t *)malloc(2 * D * nv * sizeof(v_t));
+
+ for (v_t i = 0; i < nv; i++) {
+ for (D_t j = 0; j < D; j++) {
+ v_adj[v_i[i] + 2 * j] = pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(i + pow(L, j), pow(L, j + 1));
+ v_adj[v_i[i] + 2 * j + 1] = pow(L, j + 1) * (i / ((v_t)pow(L, j + 1))) + fmod(pow(L, j+1) + i - pow(L, j), pow(L, j + 1));
+ }
+ }
+
+ graph_t *g = (graph_t *)malloc(sizeof(graph_t));
+
+ g->ne = ne;
+ g->nv = nv;
+ g->v_i = v_i;
+ g->v_adj = v_adj;
+
+ return g;
+}
+
+graph_t *graph_add_ext(const graph_t *G) {
+ graph_t *tG = (graph_t *)calloc(1, sizeof(graph_t));
+
+ tG->nv = G->nv + 1;
+ tG->ne = G->ne + G->nv;
+
+ tG->v_i = (v_t *)malloc((tG->nv + 1) * sizeof(v_t));
+ tG->v_adj = (v_t *)malloc(2 * tG->ne * sizeof(v_t));
+
+ for (v_t i = 0; i < G->nv + 1; i++) {
+ tG->v_i[i] = G->v_i[i] + i;
+ }
+
+ tG->v_i[tG->nv] = 2 * tG->ne;
+
+ for (v_t i = 0; i < G->nv; i++) {
+ v_t nn = G->v_i[i + 1] - G->v_i[i];
+
+ for (v_t j = 0; j < nn; j++) {
+ tG->v_adj[tG->v_i[i] + j] = G->v_adj[G->v_i[i] + j];
+ }
+
+ tG->v_adj[tG->v_i[i] + nn] = G->nv;
+ tG->v_adj[tG->v_i[G->nv] + i] = i;
+ }
+
+ return tG;
+}
+
+void graph_free(graph_t *g) {
+ free(g->v_i);
+ free(g->v_adj);
+ free(g);
+}
+
diff --git a/lib/graph.h b/lib/graph.h
new file mode 100644
index 0000000..eda7630
--- /dev/null
+++ b/lib/graph.h
@@ -0,0 +1,18 @@
+
+#include <inttypes.h>
+#include <math.h>
+#include <stdlib.h>
+
+#include "types.h"
+
+typedef struct {
+ v_t ne;
+ v_t nv;
+ v_t *v_i;
+ v_t *v_adj;
+} graph_t;
+
+graph_t *graph_create_square(D_t D, L_t L);
+
+graph_t *graph_add_ext(const graph_t *G);
+
diff --git a/lib/rand.c b/lib/rand.c
new file mode 100644
index 0000000..76f537d
--- /dev/null
+++ b/lib/rand.c
@@ -0,0 +1,20 @@
+
+#include <assert.h>
+#include <stdio.h>
+
+unsigned long int rand_seed() {
+ /*
+ returns a random unsigned long integer read from the standard unix
+ pseudorandom device /dev/urandom
+ */
+
+ FILE *f = fopen("/dev/urandom", "r");
+ assert(f != NULL);
+
+ unsigned long int seed;
+ fread(&seed, sizeof(unsigned long int), 1, f);
+
+ fclose(f);
+
+ return seed;
+}
diff --git a/lib/rand.h b/lib/rand.h
new file mode 100644
index 0000000..2354f6a
--- /dev/null
+++ b/lib/rand.h
@@ -0,0 +1,7 @@
+
+#pragma once
+
+#include <assert.h>
+#include <stdio.h>
+
+unsigned long int rand_seed();
diff --git a/lib/queue.c b/lib/stack.c
index 0823518..9fde1e6 100644
--- a/lib/queue.c
+++ b/lib/stack.c
@@ -1,7 +1,7 @@
-#include "queue.h"
+#include "stack.h"
-void stack_push(ll_t **q, uint32_t x) {
+void stack_push(ll_t **q, v_t x) {
ll_t *nq = malloc(sizeof(ll_t));
nq->x = x;
nq->next = *q;
@@ -17,11 +17,11 @@ void stack_push_d(dll_t **q, double x) {
*q = nq;
}
-uint32_t stack_pop(ll_t **q) {
+v_t stack_pop(ll_t **q) {
ll_t *old_q = *q;
*q = old_q->next;
- uint32_t x = old_q->x;
+ v_t x = old_q->x;
free(old_q);
@@ -39,12 +39,3 @@ double stack_pop_d(dll_t **q) {
return x;
}
-bool stack_contains(const ll_t *q, uint32_t x) {
- if (q == NULL) {
- return false;
- } else if (q->x == x) {
- return true;
- } else {
- return stack_contains(q->next, x);
- }
-}
diff --git a/lib/queue.h b/lib/stack.h
index 23834ba..a354ab5 100644
--- a/lib/queue.h
+++ b/lib/stack.h
@@ -6,8 +6,10 @@
#include <stdlib.h>
#include <string.h>
+#include "types.h"
+
typedef struct ll_tag {
- uint32_t x;
+ v_t x;
struct ll_tag *next;
} ll_t;
@@ -16,10 +18,9 @@ typedef struct dll_tag {
struct dll_tag *next;
} dll_t;
-void stack_push(ll_t **q, uint32_t x);
+void stack_push(ll_t **q, v_t x);
void stack_push_d(dll_t **q, double x);
-uint32_t stack_pop(ll_t **q);
+v_t stack_pop(ll_t **q);
double stack_pop_d(dll_t **q);
-bool stack_contains(const ll_t *q, uint32_t x);
diff --git a/lib/tree.c b/lib/tree.c
new file mode 100644
index 0000000..081f1d1
--- /dev/null
+++ b/lib/tree.c
@@ -0,0 +1,81 @@
+
+#include "tree.h"
+
+node_t *tree_createNode(v_t value, v_t level, node_t *left, node_t *right) {
+ node_t *node = (node_t *)malloc(sizeof(node_t));
+
+ node->value = value;
+ node->level = level;
+ node->left = left;
+ node->right = right;
+
+ return node;
+}
+
+void tree_freeNode(node_t *node) {
+ if (node != NULL) {
+ tree_freeNode(node->left);
+ tree_freeNode(node->right);
+ free(node);
+ }
+}
+
+void tree_skew(node_t **T) {
+ if (*T != NULL) {
+ if ((*T)->left != NULL) {
+ if ((*T)->left->level == (*T)->level) {
+ node_t *L = (*T)->left;
+ (*T)->left = L->right;
+ L->right = (*T);
+
+ *T = L;
+ }
+ }
+ }
+}
+
+void tree_split(node_t **T) {
+ if (*T != NULL) {
+ if ((*T)->right != NULL) {
+ if ((*T)->right->right != NULL) {
+ if ((*T)->level == (*T)->right->right->level) {
+ node_t *R = (*T)->right;
+ (*T)->right = R->left;
+ R->left = *T;
+ R->level = R->level + 1;
+ *T = R;
+ }
+ }
+ }
+ }
+}
+
+void tree_insert(node_t **T, v_t x) {
+ if (*T == NULL) {
+ *T = tree_createNode(x, 1, NULL, NULL);
+ } else if (x < (*T)->value) {
+ node_t *L = (*T)->left;
+ tree_insert(&L, x);
+ (*T)->left = L;
+ } else if (x > (*T)->value) {
+ node_t *R = (*T)->right;
+ tree_insert(&R, x);
+ (*T)->right = R;
+ }
+
+ tree_skew(T);
+ tree_split(T);
+}
+
+bool tree_contains(node_t *T, v_t x) {
+ if (T == NULL) {
+ return false;
+ } else if (x < T->value) {
+ return tree_contains(T->left, x);
+ } else if (x > T->value) {
+ return tree_contains(T->right, x);
+ } else {
+ return true;
+ }
+}
+
diff --git a/lib/tree.h b/lib/tree.h
new file mode 100644
index 0000000..dc22c2d
--- /dev/null
+++ b/lib/tree.h
@@ -0,0 +1,22 @@
+
+#pragma once
+
+#include <inttypes.h>
+#include <stdlib.h>
+#include <stdbool.h>
+
+#include "types.h"
+
+typedef struct node_t {
+ v_t value;
+ v_t level;
+ struct node_t *left;
+ struct node_t *right;
+} node_t;
+
+void tree_insert(node_t **T, v_t x);
+
+bool tree_contains(node_t *T, v_t x);
+
+void tree_freeNode(node_t *T);
+
diff --git a/lib/types.h b/lib/types.h
new file mode 100644
index 0000000..982db47
--- /dev/null
+++ b/lib/types.h
@@ -0,0 +1,16 @@
+
+#include <inttypes.h>
+
+typedef uint_fast32_t v_t;
+typedef uint_fast8_t q_t;
+typedef uint_fast8_t D_t;
+typedef uint_fast16_t L_t;
+typedef uint_fast64_t count_t;
+
+#define MAX_Q 256
+
+#define PRID PRIu8
+#define PRIL PRIu16
+#define PRIcount PRIu64
+#define PRIq PRIu8
+
diff --git a/lib/wolff.h b/lib/wolff.h
index dfe7b1c..dae907a 100644
--- a/lib/wolff.h
+++ b/lib/wolff.h
@@ -13,34 +13,26 @@
#include <string.h>
#include <sys/types.h>
-#include <jst/graph.h>
-#include <jst/rand.h>
-
-#include "queue.h"
+#include "types.h"
+#include "rand.h"
+#include "stack.h"
#include "convex.h"
-
-typedef enum {
- WOLFF,
- WOLFF_GHOST,
- METROPOLIS
-} sim_t;
+#include "graph.h"
+#include "tree.h"
typedef struct {
graph_t *g;
- bool *spins;
- int32_t M;
- double H;
+ q_t *spins;
+ double T;
+ double *H;
+ double T_prob;
+ double *H_probs;
+ double E;
+ v_t *M;
+ q_t q;
} ising_state_t;
typedef struct {
- uint32_t nv;
- int32_t dJb;
- int32_t dHb;
- bool hit_ghost;
- ll_t *spins;
-} cluster_t;
-
-typedef struct {
uint64_t n;
double x;
double dx;
@@ -60,21 +52,13 @@ typedef struct {
double O2;
} autocorr_t;
-int8_t sign(double x);
-
-cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *x, bool stop_on_ghost,
- gsl_rng *r);
+v_t flip_cluster(ising_state_t *s, v_t v0, q_t s1, gsl_rng *r);
graph_t *graph_add_ext(const graph_t *g);
-uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
- double *ps);
-
void update_meas(meas_t *m, double x);
void update_autocorr(autocorr_t *OO, double O);
-double add_to_avg(double mx, double x, uint64_t n);
-
double rho(autocorr_t *o, uint64_t i);
diff --git a/lib/wolff_tools.c b/lib/wolff_tools.c
index a417d4e..416a498 100644
--- a/lib/wolff_tools.c
+++ b/lib/wolff_tools.c
@@ -1,257 +1,94 @@
-#include "queue.h"
#include "wolff.h"
-int8_t spin_to_sign(bool spin) {
- /* takes a spin (represented by a bool) and converts it to an integer (plus
- * or minus one). our convention takes true to be negative spin and false to
- * be positive spin
- */
- if (spin) {
- return -1;
- } else {
- return 1;
- }
-}
-
-int8_t sign(double x) {
- // the sign function. zero returns positive sign.
- return x >= 0 ? 1 : -1;
-}
-
-graph_t *graph_add_ext(const graph_t *G) {
- /* takes a graph object G and returns tG, the same graph with an extra vertex
- * that is connected to every other vertex.
- */
- graph_t *tG = (graph_t *)calloc(1, sizeof(graph_t));
-
- tG->nv = G->nv + 1;
- tG->ne = G->ne + G->nv;
-
- tG->ev = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t));
- tG->vei = (uint32_t *)malloc((tG->nv + 1) * sizeof(uint32_t));
- tG->ve = (uint32_t *)malloc(2 * tG->ne * sizeof(uint32_t));
- tG->vx = (double *)malloc(2 * tG->nv * sizeof(double));
- tG->bq = (bool *)malloc(tG->nv * sizeof(bool));
-
- memcpy(tG->ev, G->ev, 2 * G->ne * sizeof(uint32_t));
- memcpy(tG->vx, G->vx, 2 * G->nv * sizeof(double));
- memcpy(tG->bq, G->bq, G->nv * sizeof(bool));
-
- tG->vx[2 * G->nv] = -1;
- tG->vx[2 * G->nv + 1] = -0.5;
- tG->bq[G->nv] = false;
-
- for (uint32_t i = 0; i < G->nv; i++) {
- tG->ev[2 * G->ne + 2 * i] = i;
- tG->ev[2 * G->ne + 2 * i + 1] = G->nv;
- }
-
- for (uint32_t i = 0; i < G->nv; i++) {
- tG->vei[i] = G->vei[i] + i;
-
- for (uint32_t j = 0; j < G->vei[i + 1] - G->vei[i]; j++) {
- tG->ve[tG->vei[i] + j] = G->ve[G->vei[i] + j];
- }
-
- tG->ve[tG->vei[i] + G->vei[i + 1] - G->vei[i]] = G->ne + i;
- }
-
- tG->vei[G->nv] = G->vei[G->nv] + G->nv;
- tG->vei[G->nv + 1] = tG->vei[G->nv] + G->nv;
-
- for (uint32_t i = 0; i < G->nv; i++) {
- tG->ve[tG->vei[G->nv] + i] = G->ne + i;
- }
-
- return tG;
-}
-
-uint32_t get_neighbor(const graph_t *g, uint32_t v, uint32_t i) {
- // returns the index of the ith neighbor of vertex v on graph g.
- assert(i < g->vei[v + 1] - g->vei[v]); // don't request a neighbor over the total number of neighbors!
-
- uint32_t e, v1, v2;
-
- e = g->ve[g->vei[v] + i]; // select the ith bond connected to site
- v1 = g->ev[2 * e];
- v2 = g->ev[2 * e + 1];
-
- return v == v1 ? v2 : v1; // distinguish neighboring site from site itself
+q_t delta(q_t s0, q_t s1) {
+ return s0 == s1 ? 1 : 0;
}
-cluster_t *flip_cluster(const graph_t *g, const double *ps, bool *s, bool stop_on_ghost,
- gsl_rng *r) {
- /* flips a wolff cluster on the graph g with initial state s. s is always
- * changed by the routine. the probability of adding a normal bond to the
- * cluster is passed by ps[0], and the probability of adding a bond with the
- * external field spin to the cluster is passed by ps[1]. if stop_on_ghost is
- * true, adding the external field spin to the cluster stops execution of the
- * routine. flip_cluster returns an object of type cluster_t, which encodes
- * information about the flipped cluster.
- */
- uint32_t v0;
- int32_t n_h_bonds, n_bonds;
- bool s0;
- cluster_t *c;
-
- v0 = gsl_rng_uniform_int(r, g->nv); // pick a random vertex
- s0 = s[v0]; // record its orientation
+v_t flip_cluster(ising_state_t *s, v_t v0, q_t step, gsl_rng *r) {
+ q_t s0 = s->spins[v0];
+ v_t nv = 0;
ll_t *stack = NULL; // create a new stack
stack_push(&stack, v0); // push the initial vertex to the stack
- // initiate the data structure for returning flip information
- c = (cluster_t *)calloc(1, sizeof(cluster_t));
+ node_t *T = NULL;
while (stack != NULL) {
- uint32_t v;
- uint32_t nn;
-
- v = stack_pop(&stack);
- nn = g->vei[v + 1] - g->vei[v];
+ v_t v = stack_pop(&stack);
- if (s[v] == s0 && !(c->hit_ghost && stop_on_ghost)) { // if the vertex hasn't already been flipped
- s[v] = !s[v]; // flip the vertex
- if (stop_on_ghost) {
- stack_push(&(c->spins), v);
- }
-
- for (uint32_t i = 0; i < nn; i++) {
- bool is_ext;
- uint32_t vn;
- int32_t *bond_counter;
- double prob;
+ if (!tree_contains(T, v)) { // if the vertex hasn't already been flipped
+ q_t s_old = s->spins[v];
+ q_t s_new = (s->spins[v] + step) % s->q;
- vn = get_neighbor(g, v, i);
+ s->spins[v] = s_new; // flip the vertex
+ tree_insert(&T, v);
- is_ext = (v == g->nv - 1 || vn == g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph
+ v_t nn = s->g->v_i[v + 1] - s->g->v_i[v];
- bond_counter = is_ext ? &(c->dHb) : &(c->dJb);
- prob = is_ext ? ps[1] : ps[0];
+ for (v_t i = 0; i < nn; i++) {
+ v_t vn = s->g->v_adj[s->g->v_i[v] + i];
+ q_t sn = s->spins[vn];
+ double prob;
- if (s[vn] ==
- s0) { // if the neighboring site matches the flipping cluster...
- (*bond_counter)++;
+ bool is_ext = (v == s->g->nv - 1 || vn == s->g->nv - 1);
- if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
- if (is_ext && stop_on_ghost) {
- c->hit_ghost = true;
- }
- stack_push(&stack, vn); // push the neighboring vertex to the stack
+ if (is_ext) {
+ q_t M_ind_0;
+ q_t M_ind_1;
+ if (vn == s->g->nv - 1) {
+ M_ind_0 = (s_old + s->q - s->spins[s->g->nv - 1]) % s->q;
+ M_ind_1 = (s_new + s->q - s->spins[s->g->nv - 1]) % s->q;
+ } else {
+ M_ind_0 = (sn + s->q - s_old) % s->q;
+ M_ind_1 = (sn + s->q - s_new) % s->q;
}
+ prob = s->H_probs[M_ind_1 * s->q + M_ind_0];
+ s->M[M_ind_0]--;
+ s->M[M_ind_1]++;
+ s->E += - s->H[M_ind_1] + s->H[M_ind_0];
} else {
- (*bond_counter)--;
+ prob = s->T_prob * delta(s0, sn);
+ s->E += - ((double)delta(s->spins[v], sn)) + ((double)delta(s0, sn));
}
- }
-
- if (v != g->nv - 1) { // count the number of non-external sites that flip
- c->nv++;
- }
- }
- }
-
- return c;
-}
-
-uint32_t wolff_step(double T, double H, ising_state_t *s, sim_t sim, gsl_rng *r,
- double *ps) {
- uint32_t n_flips;
- bool no_r, no_ps;
- no_r = false;
- no_ps = false;
-
- if (r == NULL) {
- no_r = true;
- r = gsl_rng_alloc(gsl_rng_mt19937);
- gsl_rng_set(r, jst_rand_seed());
- }
-
- if (ps == NULL) { // computing exponentials is relatively expensive, so calling functions are given the option to supply values calculated once in the entire runtime
- no_ps = true;
- ps = (double *)malloc(2 * sizeof(double));
- ps[0] = 1 - exp(-2 / T);
- ps[1] = 1 - exp(-2 * fabs(H) / T);
- }
-
- switch (sim) {
- case METROPOLIS: {
- uint32_t v0 = gsl_rng_uniform_int(r, s->g->nv); // pick a random vertex
- uint32_t nn = s->g->vei[v0 + 1] - s->g->vei[v0];
-
- double dE = 0;
- for (uint32_t i = 0; i < nn; i++) {
- bool is_ext;
- uint32_t vn;
- int32_t *bond_counter;
- double prob;
- vn = get_neighbor(s->g, v0, i);
-
- is_ext = (v0 == s->g->nv - 1 || vn == s->g->nv - 1); // our edge contained the "ghost" spin if either of its vertices was the last in the graph
-
- if (is_ext) {
- dE += 2 * H * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]);
- } else {
- dE += 2 * spin_to_sign(s->spins[vn]) * spin_to_sign(s->spins[v0]);
+ if (gsl_rng_uniform(r) < prob) { // and with probability ps[e]...
+ stack_push(&stack, vn); // push the neighboring vertex to the stack
}
}
- double p = exp(-dE / T);
- if (gsl_rng_uniform(r) < p) {
- s->M += - sign(H) * 2 * spin_to_sign(s->spins[v0]) * spin_to_sign(s->spins[s->g->nv - 1]);
- s->H += dE;
- s->spins[v0] = !s->spins[v0];
+ if (v != s->g->nv - 1) { // count the number of non-external sites that flip
+ nv++;
}
+ }
+ }
- n_flips = 1;
- } break;
-
- case WOLFF: {
- cluster_t *c = flip_cluster(s->g, ps, s->spins, false, r);
- s->M += - sign(H) * 2 * c->dHb;
- s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
- n_flips = c->nv;
+ tree_freeNode(T);
- free(c);
- } break;
- case WOLFF_GHOST: {
- cluster_t *c = flip_cluster(s->g, ps, s->spins, true, r);
+ return nv;
+}
- if (c->hit_ghost) {
- while (c->spins != NULL) {
- uint32_t v = stack_pop(&(c->spins));
- s->spins[v] = !s->spins[v];
- // if we hit the external spin, undo the cluster flip
- }
- } else {
- while (c->spins != NULL) {
- stack_pop(&(c->spins));
- // we have to clear the memory on the stack anyway...
- }
- s->M += - sign(H) * 2 * c->dHb;
- s->H += 2 * (c->dJb + sign (H) * H * c->dHb);
- }
+double add_to_avg(double mx, double x, count_t n) {
+ return mx * (n / (n + 1.0)) + x * 1.0 / (n + 1.0);
+}
- n_flips = c->nv;
+void update_meas(meas_t *m, double x) {
+ count_t n = m->n;
- free(c);
- } break;
- }
+ m->x = add_to_avg(m->x, x, n);
+ m->x2 = add_to_avg(m->x2, pow(x, 2), n);
- if (no_ps) {
- free(ps);
- }
+ 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 (no_r) {
- gsl_rng_free(r);
+ 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);
}
- return n_flips;
-}
-
-double add_to_avg(double mx, double x, uint64_t n) {
- return mx * (n / (n + 1.)) + x * 1. / (n + 1.);
+ (m->n)++;
}
void update_autocorr(autocorr_t *OO, double O) {
@@ -260,7 +97,7 @@ void update_autocorr(autocorr_t *OO, double O) {
dll_t *Otmp = OO->Op;
dll_t *Osave;
- uint64_t t = 0;
+ count_t t = 0;
while (Otmp != NULL) {
OO->OO[t] = add_to_avg(OO->OO[t], O * (Otmp->x), OO->n - t - 1);
@@ -286,25 +123,6 @@ void update_autocorr(autocorr_t *OO, double O) {
OO->n++;
}
-double rho(autocorr_t *o, uint64_t i) {
+double rho(autocorr_t *o, count_t i) {
return (o->OO[i] - pow(o->O, 2)) / (o->O2 - pow(o->O, 2));
}
-
-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)++;
-}