summaryrefslogtreecommitdiff
path: root/euclidean.hpp
blob: acf65acbad58434ddc6f7f3255fba39149d1ff9a (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90

#pragma once

#include "mod.hpp"
#include "matrix.hpp"
#include "spin.hpp"
#include "vector.hpp"

typedef double Radius;
typedef signed IsingSpin;

template <class T, int D> class Dimer {
public:
  double radius;
  Vector<T, D> relativePosition;
};

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;

  Affine() {
    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;
      }
    }
  }

  Affine(const Vector<T, D>& t0, const Matrix<T, D>& r0) : t(t0), r(r0) {}

  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};
  }

  template <class S> Spin<T, D, S> act(const Spin<T, D, S>& s) const {
    return {.x = act(s.x), .s = act(s.s)};
  }
};

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 Euclidean(Affine<T, D>::actA(t)); };
  Euclidean inverse() const { return Euclidean(Affine<T, D>::inverse()); }

  using Affine<T, D>::act;
};

template <class T, int D> class TorusGroup : public Affine<T, D>{
private:
  T L;

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) = mod(s_new(i), L);
    }

    return s_new;
  }

public:
  TorusGroup(T L) : Affine<T, D>(), L(L) {}
  TorusGroup(T L, const Affine<T, D>& t) : Affine<T, D>(t), L(L) {}
  TorusGroup(T L, Vector<T, D> t0, Matrix<T, D> r0) : Affine<T, D>(t0, r0), L(L) {}

  TorusGroup act(const TorusGroup& t) const { return TorusGroup(L, Affine<T, D>::actA(t)); }
  TorusGroup inverse() const { return TorusGroup(L, Affine<T, D>::inverse()); }

  using Affine<T, D>::act;
};