fieldpattern:: an option to use 'quotient group' fieldpatterns in 3D for more efficient representation

fieldpattern.cpp

@@ -11,6 +11,8 @@ namespace hr {
 
 EX namespace fieldpattern {
 
+EX bool use_quotient_fp = false;
+
 int limitsq = 10;
 int limitp = 10000;
 int limitv = 100000;
@@ -311,6 +313,8 @@ struct fpattern {
   unsigned compute_hash();
 
   void set_field(int p, int sq);
+  
+  unsigned hashv;
 
   #if MAXMDIM >= 4
   // general 4D
@@ -328,6 +332,7 @@ struct fpattern {
   #endif
 
   vector<triplet_info> find_triplets();
+  void generate_quotientgroup();
   };
 
 #if CAP_THREAD && MAXMDIM >= 4
@@ -505,11 +510,78 @@ bool fpattern::generate_all3() {
       for(int j=0; j<local_group; j++) add1(mmul(E, matrices[j]));
     if(isize(matrices) >= limitv) { println(hlog, "limitv exceeded"); return false; }
     }
-  unsigned hashv = compute_hash();
+  hashv = compute_hash();
   DEBB(DF_FIELD, ("all = ", isize(matrices), "/", local_group, " = ", isize(matrices) / local_group, " hash = ", hashv, " count = ", ++hash_found[hashv]));
+  
+  if(use_quotient_fp) 
+    generate_quotientgroup();
   return true;
   }
 
+void fpattern::generate_quotientgroup() {
+  int MS = isize(matrices);
+  int best_p = 0, best_i = 0;
+  for(int i=0; i<MS; i++) {
+    int j = i, p = 1;
+    while(j >= local_group) 
+      j = gmul(j, i), p++;
+    if(j == 0 && p > best_p) {
+      bool okay = true;
+
+      vector<bool> visited(MS, false);
+      for(int ii=0; ii<MS; ii++) if(!visited[ii]) {
+        int jj = ii;
+        for(int k=0; k<p; k++) {
+          if(k && jj/local_group == ii/local_group) okay = false;
+          visited[jj] = true;
+          jj = gmul(i, jj);
+          }
+        }
+      
+      if(okay) {
+        bool chk = (MS/p) % local_group;
+        println(hlog, "quotient by ", i, " : ", p, " times less, ", (MS/p/local_group), " tiles, check ", chk);
+        best_p = p; best_i = i;
+        if(chk) {
+          exit(1);
+          }
+        }
+      }
+    }
+  
+  if(best_p > 1) {
+    vector<int> new_id(MS, -1);
+    vector<int> orig_id(MS, -1);
+    vector<matrix> new_matrices;
+    int nv = 0;
+    for(int i=0; i<MS; i++) if(new_id[i] == -1) {
+      int prode = i;
+      for(int l=0; l<local_group; l++) {
+        new_matrices.push_back(matrices[i+l]);
+        }
+      for(int k=0; k<best_p; k++) {
+        for(int l=0; l<local_group; l++) {
+          new_id[gmul(prode, l)] = nv + l;
+          }
+        prode = gmul(best_i, prode);
+        }
+      nv += local_group;
+      }
+    println(hlog, "got nv = ", nv, " / ", local_group);
+    
+    for(int i=0; i<MS; i++)
+      matcode[matrices[i]] = new_id[i];
+    matrices = std::move(new_matrices);
+    println(hlog, "size matrices = ", isize(matrices), " size matcode = ", isize(matcode));
+    println(hlog, tie(P, R, X));
+    
+    /*println(hlog, "TRY AGAIN");
+    generate_quotientgroup();
+    exit(1);*/
+    }  
+  
+  }
+
 int fpattern::solve3() {
   reg3::construct_relations();
   
@@ -558,7 +630,7 @@ int fpattern::solve3() {
     #if CAP_THREAD && MAXMDIM >= 4
     if(dis) { dis->discovered(); continue; }
     #endif
-    if(force_hash && compute_hash() != force_hash) continue;
+    if(force_hash && hashv != force_hash) continue;
     cmb++;
     goto ok;
     bad: ;
@@ -1415,7 +1487,7 @@ void discovery::activate() {
 void discovery::discovered() {
   std::unique_lock<std::mutex> lk(lock);
   auto& e = experiment;
-  hashes_found[e.compute_hash()] = make_tuple(e.Prime, e.wsquare, e.R, e.P, e.X, isize(e.matrices) / e.local_group);
+  hashes_found[e.hashv] = make_tuple(e.Prime, e.wsquare, e.R, e.P, e.X, isize(e.matrices) / e.local_group);
   }
 
 void discovery::suspend() { is_suspended = true; }

reg3.cpp

@@ -386,7 +386,7 @@ EX namespace reg3 {
         tmatrices[a].resize(S7);
         for(int b=0; b<S7; b++) {
           int k = lgr*a;
-          k = f->gmul(f->gmul(k, moveid[b]), lgr);
+          k = f->matcode[ f->mmul(f->mmul(f->matrices[k], f->matrices[moveid[b]]), f->P) ];
           for(int l=0; l<lgr; l++) if(f->gmul(k, l) % lgr == 0) {
             tmatrices[a][b] = cgi.adjmoves[b] * f->fullv[l];
             allh[a]->c.connect(b, allh[k/lgr], movedir[l], false);