From 7038df605cc036b36c85c027baeac726e17392de Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Wed, 15 Jan 2020 19:18:17 +0100
Subject: [PATCH] field:: quick distance computations in 3D

---
 cell.cpp         |  2 +-
 fieldpattern.cpp | 33 ++++++++++++++++++++++++++-------
 pattern2.cpp     |  3 ++-
 3 files changed, 29 insertions(+), 9 deletions(-)

diff --git a/cell.cpp b/cell.cpp
index 5ee179a4..a76ff566 100644
--- a/cell.cpp
+++ b/cell.cpp
@@ -983,7 +983,7 @@ EX int celldistance(cell *c1, cell *c2) {
   if(hybri) return hybrid::celldistance(c1, c2);
   
   #if CAP_FIELD
-  if(geometry == gFieldQuotient && !GOLDBERG && WDIM == 2)
+  if(geometry == gFieldQuotient && !GOLDBERG)
     return currfp.getdist(fieldpattern::fieldval(c1), fieldpattern::fieldval(c2));
   #endif
   
diff --git a/fieldpattern.cpp b/fieldpattern.cpp
index 783335cb..80be1d1e 100644
--- a/fieldpattern.cpp
+++ b/fieldpattern.cpp
@@ -59,8 +59,12 @@ struct matrix : array<array<int, MAXMDIM>, MAXMDIM> {
   };
 #endif
 
+EX int groupspin(int id, int d, int group) {
+  return group*(id/group) + (id + d) % group;
+  }
+
 EX int btspin(int id, int d) {
-  return S7*(id/S7) + (id + d) % S7;
+  return groupspin(id, d, S7);
   }
 
 #if HDR
@@ -688,7 +692,7 @@ int fpattern::getdist(pair<int,bool> a, pair<int,bool> b) {
   }
  
 int fpattern::dijkstra(vector<char>& dists, vector<int> indist[MAXDIST]) {
-  int N = connections.size();
+  int N = isize(matrices);
   dists.resize(N);
   for(int i=0; i<N; i++) dists[i] = MAXDIST-1;
   int maxd = 0;
@@ -698,15 +702,18 @@ int fpattern::dijkstra(vector<char>& dists, vector<int> indist[MAXDIST]) {
     if(dists[at] <= i) continue;
     maxd = i;
     dists[at] = i;
-    for(int q=0; q<S7; q++) {
+    int lg = MWDIM == 4 ? local_group : S7;
+    for(int q=0; q<lg; q++) {
       dists[at] = i;
-      if(PURE) // todo-variation: PURE here?
+      if(WDIM == 3)
+        indist[i+1].push_back(gmul(at, local_group));
+      else if(PURE) // todo-variation: PURE here?
         indist[i+1].push_back(connections[at]);
       else {
         indist[i+2].push_back(connections[at]);
         indist[i+3].push_back(connections[btspin(connections[at], 2)]);
         }
-      at = btspin(at, 1);
+      at = groupspin(at, 1, lg);
       }
     }
   return maxd;
@@ -714,8 +721,18 @@ int fpattern::dijkstra(vector<char>& dists, vector<int> indist[MAXDIST]) {
 
 void fpattern::analyze() {
 
-  if(WDIM == 3) return;
-
+  if(MWDIM == 4) {
+    /* we need to compute inverses */
+    int N = isize(matrices);
+    inverses.resize(N);
+    for(int i=0; i<N; i++) {
+      matrix M = matrices[i];
+      matrix M2 = mpow(M, N-1);
+      inverses[i] = matcode[M2];
+      println(hlog, mmul(M, M2));
+      }
+    }
+    
   DEBB(DF_FIELD, ("variation = %d\n", int(variation)));
   int N = connections.size();
   
@@ -725,6 +742,8 @@ void fpattern::analyze() {
 
   indist[0].push_back(0);
   int md0 = dijkstra(disthep, indist);
+  
+  if(MWDIM == 4) return;
 
   indist[1].push_back(0);
   indist[1].push_back(connections[3]);
diff --git a/pattern2.cpp b/pattern2.cpp
index fd308532..0f1c2635 100644
--- a/pattern2.cpp
+++ b/pattern2.cpp
@@ -366,7 +366,8 @@ EX int zebra3(cell *c) {
 EX namespace fieldpattern {
 
 EX pair<int, bool> fieldval(cell *c) {
-  if(ctof(c)) return make_pair(int(c->master->fieldval), false);
+  if(WDIM == 3) return make_pair(int(c->master->fieldval) * currfp.local_group, false);
+  else if(ctof(c)) return make_pair(int(c->master->fieldval), false);
   else return make_pair(btspin(c->master->fieldval, c->c.spin(0)), true);
   }