1 files changed, 31 insertions, 82 deletions

diff --git a/euclidean.hpp b/euclidean.hpp
index fb2dbcd..f0e54d4 100644
--- a/euclidean.hpp
+++ b/euclidean.hpp
@@ -1,6 +1,7 @@
 
 #pragma once
 
+#include "mod.hpp"
 #include "matrix.hpp"
 #include "spin.hpp"
 #include "vector.hpp"
@@ -14,12 +15,16 @@ public:
   Vector<T, D> relativePosition;
 };
 
-template <class T, int D> class Euclidean {
+template <class T, int D> class Affine {
+protected:
+  virtual Vector<T, D> actV(const Vector<T, D>& x) const { return t + r * x; }
+  Affine actA(const Affine& x) const { return Affine(act(x.t), r * x.r); }
+
 public:
   Vector<T, D> t;
   Matrix<T, D> r;
 
-  Euclidean(T L) {
+  Affine() {
     for (unsigned i = 0; i < D; i++) {
       t(i) = 0;
       r(i, i) = 1;
@@ -29,21 +34,13 @@ public:
     }
   }
 
-  Euclidean(Vector<T, D> t0, Matrix<T, D> r0) {
-    t = t0;
-    r = r0;
-  }
-
-  Vector<T, D> act(const Vector<T, D>& x) const { return t + r * x; }
+  Affine(const Vector<T, D>& t0, const Matrix<T, D>& r0) : t(t0), r(r0) {}
 
-  template <typename U> Vector<U, D> act(const Vector<U, D>& x) const {
-    return t.template cast<U>() + r.template cast<U>() * x;
-  }
+  Affine inverse() const { return Affine(-r.transpose() * t, r.transpose()); }
 
+  Vector<T, D> act(const Vector<T, D>& x) const { return actV(x); }
   Radius act(Radius r) const { return r; }
-
   IsingSpin act(IsingSpin s) const { return s; }
-
   Dimer<T, D> act(const Dimer<T, D>& d) const {
     return {.radius = d.radius, .relativePosition = r * d.relativePosition};
   }
@@ -51,90 +48,42 @@ public:
   template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
     return {.x = act(s.x), .s = act(s.s)};
   }
+};
 
-  Euclidean act(const Euclidean& x) const { return Euclidean(r * x.t + t, r * x.r); }
+template <class T, int D> class Euclidean : public Affine<T, D> {
+public:
+  Euclidean(T L = 0) : Affine<T, D>() {}
+  Euclidean(const Affine<T, D>& a) : Affine<T, D>(a) {}
+  Euclidean(Vector<T, D> t0, Matrix<T, D> r0) : Affine<T, D>(t0, r0) {}
+  
+  Euclidean act(const Euclidean& t) const { return Affine<T, D>::actA(t); };
 
-  Euclidean inverse() const { return Euclidean(-r.transpose() * t, r.transpose()); }
+  using Affine<T, D>::act;
 };
 
-template <class T, int D> class TorusGroup {
+template <class T, int D> class TorusGroup : public Affine<T, D>{
 private:
   T L;
 
-public:
-  Vector<T, D> t;
-  Matrix<T, D> r;
-
-  /** brief TorusGroup - default constructor, constructs the identity
-   */
-  TorusGroup(T L) : L(L) {
-    for (unsigned i = 0; i < D; i++) {
-      t(i) = 0;
-      r(i, i) = 1;
-      for (unsigned j = 1; j < D; j++) {
-        r(i, (i + j) % D) = 0;
-      }
-    }
-  }
-
-  TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : L(L) {
-    t = t0;
-    r = r0;
-  }
-
-  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
-    Spin<T, D, S> s_new;
-
-    s_new.x = t + r * s.x;
-    s_new.s = s.s;
-
-    for (unsigned i = 0; i < D; i++) {
-      s_new.x(i) = fmod(L + s_new.x(i), L);
-    }
-
-    return s_new;
-  }
-
-  Vector<T, D> act(const Vector<T, D>& s) const {
-    Vector<T, D> s_new = t + r * s;
-
-    for (unsigned i = 0; i < D; i++) {
-      s_new(i) = fmod(L + s_new(i), L);
-    }
-
-    return s_new;
-  }
-
-  template <typename U> Vector<U, D> act(const Vector<U, D>& s) const {
-    Vector<U, D> s_new = t.template cast<U>() + r.template cast<U>() * s;
+protected:
+  Vector<T, D> actV(const Vector<T, D>& s) const override {
+    Vector<T, D> s_new = Affine<T, D>::actV(s);
 
     for (unsigned i = 0; i < D; i++) {
-      s_new(i) = fmod(L + s_new(i), L);
+      s_new(i) = mod(s_new(i), L);
     }
 
     return s_new;
   }
 
-  TorusGroup act(const TorusGroup& x) const {
-    Vector<T, D> tnew = r * x.t + t;
-    Matrix<T, D> rnew = r * x.r;
-
-    for (unsigned i = 0; i < D; i++) {
-      tnew(i) = fmod(L + tnew(i), L);
-    }
-
-    TorusGroup pnew(this->L, tnew, rnew);
-
-    return pnew;
-  }
-
-  TorusGroup inverse() const {
-    Vector<T, D> tnew = -r.transpose() * t;
-    Matrix<T, D> rnew = r.transpose();
+public:
+  TorusGroup(T L) : Affine<T, D>(), L(L) {}
+  TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : Affine<T, D>(), L(L) {}
+  TorusGroup(T L, const Affine<T, D>& t) : Affine<T, D>(t), L(L) {}
 
-    TorusGroup pnew(this->L, tnew, rnew);
+  TorusGroup act(const TorusGroup& t) const { return TorusGroup(L, Affine<T, D>::actA(t)); }
+  TorusGroup inverse() const { return TorusGroup(L, Affine<T, D>::inverse()); }
 
-    return pnew;
-  }
+  using Affine<T, D>::act;
 };