summaryrefslogtreecommitdiff
path: root/lib/orthogonal.h
diff options
context:
space:
mode:
authorJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-26 00:32:38 -0400
committerJaron Kent-Dobias <jaron@kent-dobias.com>2018-07-26 00:32:38 -0400
commitd3b3e39a525d0c3aa9663f824ba96e0c429a8313 (patch)
tree942a902af332a532ab76fccbdb51c83e98338d5e /lib/orthogonal.h
parent225905d1ad11a0f10d63ab0c6c71e6a5265e0894 (diff)
downloadc++-d3b3e39a525d0c3aa9663f824ba96e0c429a8313.tar.gz
c++-d3b3e39a525d0c3aa9663f824ba96e0c429a8313.tar.bz2
c++-d3b3e39a525d0c3aa9663f824ba96e0c429a8313.zip
partially class-ified, ising and On work but potts and height do not
Diffstat (limited to 'lib/orthogonal.h')
-rw-r--r--lib/orthogonal.h27
1 files changed, 10 insertions, 17 deletions
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 2a391ae..c0958f2 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -53,20 +53,19 @@ void free_spin (orthogonal_t <q, T> m) {
template <q_t q, class T>
vector_t <q, T> act (orthogonal_t <q, T> m, vector_t <q, T> v) {
vector_t <q, T> v_rot;
- v_rot.x = (T *)calloc(q, sizeof(T));
if (m.is_reflection) {
double prod = 0;
for (q_t i = 0; i < q; i++) {
- prod += v.x[i] * m.x[i];
+ prod += v[i] * m.x[i];
}
for (q_t i = 0; i < q; i++) {
- v_rot.x[i] = v.x[i] - 2 * prod * m.x[i];
+ v_rot[i] = v[i] - 2 * prod * m.x[i];
}
} else {
for (q_t i = 0; i < q; i++) {
for (q_t j = 0; j < q; j++) {
- v_rot.x[i] += m.x[q * i + j] * v.x[j];
+ v_rot[i] += m.x[q * i + j] * v[j];
}
}
}
@@ -112,11 +111,10 @@ vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) {
return act(m, v); // reflections are their own inverse
} else {
vector_t <q, T> v_rot;
- v_rot.x = (T *)calloc(q, sizeof(T));
for (q_t i = 0; i < q; i++) {
for (q_t j = 0; j < q; j++) {
- v_rot.x[i] += m.x[q * j + i] * v.x[j];
+ v_rot[i] += m.x[q * j + i] * v[j];
}
}
@@ -176,33 +174,32 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q
if (n > 1) {
unsigned int rotation = gsl_rng_uniform_int(r, n);
- v.x = (double *)malloc(q * sizeof(double));
double cosr = cos(2 * M_PI * rotation / (double)n / 2.0);
double sinr = sin(2 * M_PI * rotation / (double)n / 2.0);
- v.x[0] = v0.x[0] * cosr - v0.x[1] * sinr;
- v.x[1] = v0.x[1] * cosr + v0.x[0] * sinr;
+ v[0] = v0[0] * cosr - v0[1] * sinr;
+ v[1] = v0[1] * cosr + v0[0] * sinr;
for (q_t i = 2; i < q; i++) {
- v.x[i] = v0.x[i];
+ v[i] = v0[i];
}
} else {
- v.x = v0.x;
+ v = v0;
}
double m_dot_v = 0;
for (q_t i = 0; i < q; i++) {
m.x[i] = gsl_ran_ugaussian(r);
- m_dot_v += m.x[i] * v.x[i];
+ m_dot_v += m.x[i] * v[i];
}
double v2 = 0;
double factor = epsilon * gsl_ran_ugaussian(r);
for (q_t i = 0; i < q; i++) {
- m.x[i] = m.x[i] - m_dot_v * v.x[i] + factor * v.x[i];
+ m.x[i] = m.x[i] - m_dot_v * v[i] + factor * v[i];
v2 += pow(m.x[i], 2);
}
@@ -212,10 +209,6 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q
m.x[i] /= mag_v;
}
- if (n > 1) {
- free(v.x);
- }
-
return m;
}