summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--lib/orthogonal.h51
-rw-r--r--src/wolff_On.cpp4
2 files changed, 30 insertions, 25 deletions
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index 523fd54..ce2d300 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -166,51 +166,56 @@ orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, vector_t <q, dou
}
template <q_t q>
-orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q, double> v, double epsilon) {
+orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q, double> v0, double epsilon, unsigned int n) {
orthogonal_t <q, double> m;
- vector_t <q, double> tmp_v;
m.is_reflection = true;
+ m.x = (double *)malloc(q * sizeof(double));
- tmp_v.x = (double *)malloc(q * sizeof(double));
+ vector_t <q, double> v;
- double tmp2 = 0;
- double M2 = 0;
- double tmpM = 0;
+ if (n > 1) {
+ unsigned int rotation = gsl_rng_uniform_int(r, n);
+ v.x = (double *)malloc(q * sizeof(double));
- for (q_t i = 0; i < q; i++) {
- tmp_v.x[i] = gsl_ran_ugaussian(r);
- M2 += pow(v.x[i], 2);
- tmpM += tmp_v.x[i] * v.x[i];
- }
+ double cosr = cos(2 * M_PI * rotation / (double)n / 2.0);
+ double sinr = sin(2 * M_PI * rotation / (double)n / 2.0);
- double v2 = 0;
- double factor = gsl_ran_ugaussian(r);
+ v.x[0] = v0.x[0] * cosr - v0.x[1] * sinr;
+ v.x[1] = v0.x[1] * cosr + v0.x[0] * sinr;
- for (q_t i = 0; i < q; i++) {
- tmp_v.x[i] = (tmp_v.x[i] - tmpM * v.x[i] / M2) + epsilon * factor * v.x[i] / sqrt(M2);
- v2 += tmp_v.x[i] * tmp_v.x[i];
+ for (q_t i = 2; i < q; i++) {
+ v.x[i] = v0.x[i];
+ }
+ } else {
+ v.x = v0.x;
}
- double mag_v = sqrt(v2);
+ double m2 = 0;
+ double m_dot_v = 0;
for (q_t i = 0; i < q; i++) {
- tmp_v.x[i] /= mag_v;
+ m.x[i] = gsl_ran_ugaussian(r);
+ m_dot_v += m.x[i] * v.x[i];
}
- m.x = tmp_v.x;
-
- v2 = 0;
+ double v2 = 0;
+ double factor = epsilon * gsl_ran_ugaussian(r);
for (q_t i = 0; i < q; i++) {
- v2 += m.x[i] * m.x[i];
+ m.x[i] = m.x[i] - m_dot_v * v.x[i] + factor * v.x[i];
+ v2 += pow(m.x[i], 2);
}
- mag_v = sqrt(v2);
+ double mag_v = sqrt(v2);
for (q_t i = 0; i < q; i++) {
m.x[i] /= mag_v;
}
+ if (n > 1) {
+ free(v.x);
+ }
+
return m;
}
diff --git a/src/wolff_On.cpp b/src/wolff_On.cpp
index ea80a28..997ec09 100644
--- a/src/wolff_On.cpp
+++ b/src/wolff_On.cpp
@@ -55,7 +55,7 @@ int main(int argc, char *argv[]) {
unsigned int window_size = 512;
bool modulated_field = false;
- int order = 2;
+ unsigned int order = 1;
int opt;
q_t H_ind = 0;
@@ -133,7 +133,7 @@ int main(int argc, char *argv[]) {
std::function <orthogonal_R_t(gsl_rng *, vector_R_t)> gen_R;
if (use_pert) {
- gen_R = std::bind(generate_rotation_perturbation <N_COMP>, std::placeholders::_1, std::placeholders::_2, epsilon);
+ gen_R = std::bind(generate_rotation_perturbation <N_COMP>, std::placeholders::_1, std::placeholders::_2, epsilon, order);
pert_type = "PERTURB";
} else {
gen_R = generate_rotation_uniform <N_COMP>;