summaryrefslogtreecommitdiff
path: root/lib/orthogonal.h
diff options
context:
space:
mode:
Diffstat (limited to 'lib/orthogonal.h')
-rw-r--r--lib/orthogonal.h200
1 files changed, 85 insertions, 115 deletions
diff --git a/lib/orthogonal.h b/lib/orthogonal.h
index c0958f2..34cd44e 100644
--- a/lib/orthogonal.h
+++ b/lib/orthogonal.h
@@ -11,164 +11,134 @@
#include "vector.h"
template <q_t q, class T>
-struct orthogonal_t { bool is_reflection; T *x; };
+class orthogonal_t : public std::array<std::array<T, q>, q> {
+ public :
+ bool is_reflection;
-template <q_t q, class T>
-void init(orthogonal_t <q, T> *ptr) {
- ptr->is_reflection = false;
- ptr->x = (T *)calloc(q * q, sizeof(T));
-
- for (q_t i = 0; i < q; i++) {
- ptr->x[q * i + i] = (T)1;
- }
-}
-
-template <q_t q, class T>
-orthogonal_t <q, T> copy (orthogonal_t <q, T> m) {
- orthogonal_t <q, T> m_copy;
- m_copy.is_reflection = m.is_reflection;
-
- q_t size;
-
- if (m.is_reflection) {
- size = q;
- } else {
- size = q * q;
- }
-
- m_copy.x = (T *)calloc(size, sizeof(T));
-
- for (q_t i = 0; i < size; i++) {
- m_copy.x[i] = m.x[i];
+ orthogonal_t() : is_reflection(false) {
+ for (q_t i = 0; i < q; i++) {
+ (*this)[i].fill(0);
+ (*this)[i][i] = (T)1;
+ }
}
- return m_copy;
-}
-
-template <q_t q, class T>
-void free_spin (orthogonal_t <q, T> m) {
- free(m.x);
-}
-
-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;
+ vector_t<q, T> act(const vector_t <q, T>& v) const {
+ vector_t <q, T> v_rot;
- if (m.is_reflection) {
- double prod = 0;
- for (q_t i = 0; i < q; i++) {
- prod += v[i] * m.x[i];
- }
- for (q_t i = 0; i < q; 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[i] += m.x[q * i + j] * v[j];
+ if (is_reflection) {
+ double prod = 0;
+ for (q_t i = 0; i < q; i++) {
+ prod += v[i] * (*this)[0][i];
+ }
+ for (q_t i = 0; i < q; i++) {
+ v_rot[i] = v[i] - 2 * prod * (*this)[0][i];
+ }
+ } else {
+ for (q_t i = 0; i < q; i++) {
+ for (q_t j = 0; j < q; j++) {
+ v_rot[i] += (*this)[i][j] * v[j];
+ }
}
}
- }
- return v_rot;
-}
+ return v_rot;
+ }
-template <q_t q, class T>
-orthogonal_t <q, T> act (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
- orthogonal_t <q, T> m2_rot;
+ orthogonal_t<q, T> act(const orthogonal_t <q, T>& m) const {
+ orthogonal_t <q, T> m_rot;
- m2_rot.is_reflection = false;
- m2_rot.x = (T *)calloc(q * q, sizeof(T));
+ m_rot.is_reflection = false;
- if (m1.is_reflection) {
- for (q_t i = 0; i < q; i++) {
- double akOki = 0;
+ if (is_reflection) {
+ for (q_t i = 0; i < q; i++) {
+ double akOki = 0;
- for (q_t k = 0; k < q; k++) {
- akOki += m1.x[k] * m2.x[q * k + i];
- }
+ for (q_t k = 0; k < q; k++) {
+ akOki += (*this)[0][k] * m[k][i];
+ }
- for (q_t j = 0; j < q; j++) {
- m2_rot.x[q * j + i] = m2.x[q * j + i] - 2 * akOki * m1.x[j];
+ for (q_t j = 0; j < q; j++) {
+ m_rot[j][i] = m[j][i] - 2 * akOki * (*this)[0][j];
+ }
}
- }
- } else {
- for (q_t i = 0; i < q; i++) {
- for (q_t j = 0; j < q; j++) {
- for (q_t k = 0; k < q; k++) {
- m2_rot.x[i * q + j] += m1.x[i * q + j] * m2.x[j * q + k];
+ } else {
+ for (q_t i = 0; i < q; i++) {
+ m_rot[i].fill(0);
+ for (q_t j = 0; j < q; j++) {
+ for (q_t k = 0; k < q; k++) {
+ m_rot[i][j] += (*this)[i][j] * m[j][k];
+ }
}
}
}
- }
- return m2_rot;
-}
+ return m_rot;
+ }
-template <q_t q, class T>
-vector_t <q, T> act_inverse (orthogonal_t <q, T> m, vector_t <q, T> v) {
- if (m.is_reflection) {
- return act(m, v); // reflections are their own inverse
- } else {
- vector_t <q, T> v_rot;
+ vector_t <q, T> act_inverse(const vector_t <q, T>& v) const {
+ if (is_reflection) {
+ return this->act(v); // reflections are their own inverse
+ } else {
+ vector_t <q, T> v_rot;
- for (q_t i = 0; i < q; i++) {
- for (q_t j = 0; j < q; j++) {
- v_rot[i] += m.x[q * j + i] * v[j];
+ for (q_t i = 0; i < q; i++) {
+ for (q_t j = 0; j < q; j++) {
+ v_rot[i] += (*this)[j][i] * v[j];
+ }
}
- }
- return v_rot;
+ return v_rot;
+ }
}
-}
-
-template <q_t q, class T>
-orthogonal_t <q, T> act_inverse (orthogonal_t <q, T> m1, orthogonal_t <q, T> m2) {
- if (m1.is_reflection) {
- return act(m1, m2); // reflections are their own inverse
- } else {
- orthogonal_t <q, T> m2_rot;
- m2_rot.x = (T *)calloc(q * q, sizeof(T));
- for (q_t i = 0; i < q; i++) {
- for (q_t j = 0; j < q; j++) {
- for (q_t k = 0; k < q; k++) {
- m2_rot.x[i * q + j] += m1.x[j * q + i] * m2.x[j * q + k];
+ vector_t <q, T> act_inverse(const orthogonal_t <q, T>& m) const {
+ if (is_reflection) {
+ return this->act(m); // reflections are their own inverse
+ } else {
+ orthogonal_t <q, T> m_rot;
+ m_rot.is_reflection = false;
+
+ for (q_t i = 0; i < q; i++) {
+ m_rot[i].fill(0);
+ for (q_t j = 0; j < q; j++) {
+ for (q_t k = 0; k < q; k++) {
+ m_rot[i][j] += (*this)[j][i] * m[j][k];
+ }
}
}
- }
- return m2_rot;
+ return m_rot;
+ }
}
-}
+
+};
+
template <q_t q>
-orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, vector_t <q, double> v) {
+orthogonal_t <q, double> generate_rotation_uniform (gsl_rng *r, const vector_t <q, double>& v) {
orthogonal_t <q, double> ptr;
ptr.is_reflection = true;
- ptr.x = (double *)calloc(q, sizeof(double));
double v2 = 0;
for (q_t i = 0; i < q; i++) {
- ptr.x[i] = gsl_ran_ugaussian(r);
- v2 += ptr.x[i] * ptr.x[i];
+ ptr[0][i] = gsl_ran_ugaussian(r);
+ v2 += ptr[0][i] * ptr[0][i];
}
double mag_v = sqrt(v2);
for (q_t i = 0; i < q; i++) {
- ptr.x[i] /= mag_v;
+ ptr[0][i] /= mag_v;
}
return ptr;
}
template <q_t q>
-orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q, double> v0, double epsilon, unsigned int n) {
+orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, const vector_t <q, double>& v0, double epsilon, unsigned int n) {
orthogonal_t <q, double> m;
m.is_reflection = true;
- m.x = (double *)malloc(q * sizeof(double));
vector_t <q, double> v;
@@ -191,22 +161,22 @@ orthogonal_t <q, double> generate_rotation_perturbation (gsl_rng *r, vector_t <q
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[i];
+ m[0][i] = gsl_ran_ugaussian(r);
+ m_dot_v += m[0][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[i] + factor * v[i];
- v2 += pow(m.x[i], 2);
+ m[0][i] = m[0][i] - m_dot_v * v[i] + factor * v[i];
+ v2 += pow(m[0][i], 2);
}
double mag_v = sqrt(v2);
for (q_t i = 0; i < q; i++) {
- m.x[i] /= mag_v;
+ m[0][i] /= mag_v;
}
return m;