1 // Hyperbolic Rogue -- Field Quotient geometry
2 // Copyright (C) 2011-2018 Zeno Rogue, see 'hyper.cpp' for details
3 
4 /** \file fieldpattern.cpp
5  *  \brief Field Quotient geometry
6  */
7 
8 #include "hyper.h"
9 #if CAP_FIELD
10 namespace hr {
11 
12 EX namespace fieldpattern {
13 
14 EX bool use_rule_fp = false;
15 
16 EX bool use_quotient_fp = false;
17 
18 int limitsq = 10;
19 int limitp = 10000;
20 int limitv = 100000;
21 
22 #if HDR
23 #define currfp fieldpattern::getcurrfp()
24 
25 struct primeinfo {
26   int p;
27   int cells;
28   bool squared;
29   };
30 
31 struct fgeomextra {
32   eGeometry base;
33   vector<primeinfo> primes;
34   vector<int> dualval;
35   int current_prime_id;
fgeomextrahr::fieldpattern::fgeomextra36   fgeomextra(eGeometry b, int i) : base(b), current_prime_id(i) {}
37   };
38 #endif
39 
isprime(int n)40 EX bool isprime(int n) {
41   for(int k=2; k<n; k++) if(n%k == 0) return false;
42   return true;
43   }
44 
45 #if HDR
46 #define MWDIM (prod ? 3 : WDIM+1)
47 
48 struct matrix : array<array<int, MAXMDIM>, MAXMDIM> {
operator ==hr::fieldpattern::matrix49   bool operator == (const matrix& B) const {
50     for(int i=0; i<MWDIM; i++) for(int j=0; j<MWDIM; j++)
51       if(self[i][j] != B[i][j]) return false;
52     return true;
53     }
54 
operator !=hr::fieldpattern::matrix55   bool operator != (const matrix& B) const {
56     for(int i=0; i<MWDIM; i++) for(int j=0; j<MWDIM; j++)
57       if(self[i][j] != B[i][j]) return true;
58     return false;
59     }
60 
operator <hr::fieldpattern::matrix61   bool operator < (const matrix& B) const {
62     for(int i=0; i<MWDIM; i++) for(int j=0; j<MWDIM; j++)
63       if(self[i][j] != B[i][j]) return self[i][j] < B[i][j];
64     return false;
65     }
66 
67   };
68 #endif
69 
groupspin(int id,int d,int group)70 EX int groupspin(int id, int d, int group) {
71   return group*(id/group) + (id + d) % group;
72   }
73 
btspin(int id,int d)74 EX int btspin(int id, int d) {
75   return groupspin(id, d, S7);
76   }
77 
78 #if HDR
79 
80 static const int ERR = -99;
81 
82 struct triplet_info {
83   int i, j, size;
84   };
85 
86 struct fpattern {
87 
88   unsigned force_hash;
89 
90   int Prime, wsquare, Field, dual;
91   // we perform our computations in the field Z_Prime[w] where w^2 equals wsquare
92   // (or simply Z_Prime for wsquare == 0)
93 
94   #define EASY
95   // 'easy' assumes that all elements of the field actually used
96   // are of form n or mw (not n+mw), and cs and ch are both of form n
97   // by experimentation, such cs and ch always exist
98   // many computations are much simpler under that assumption
99 
100   #ifndef EASY
101   static int neasy;
102 
mhr::fieldpattern::fpattern103   int m(int x) { x %= Prime; if(x<0) x+= Prime; return x; }
104   #endif
105 
subhr::fieldpattern::fpattern106   int sub(int a, int b) {
107     #ifdef EASY
108     return (a + b * (Prime-1)) % Prime;
109     #else
110     return m(a%Prime-b%Prime) + Prime * m(a/Prime-b/Prime);
111     #endif
112     }
113 
addhr::fieldpattern::fpattern114   int add(int a, int b) {
115     #ifdef EASY
116     if(a == ERR || b == ERR || a*b<0) return ERR;
117     return (a+b)%Prime;
118     #else
119     return m(a%Prime+b%Prime) + Prime * m(a/Prime+b/Prime);
120     #endif
121     }
122 
mulhr::fieldpattern::fpattern123   int mul(int tx, int ty) {
124     #ifdef EASY
125     return (tx*ty*((tx<0&&ty<0)?wsquare:1)) % Prime;
126     #else
127     if(tx >= Prime && tx % Prime) neasy++;
128     if(ty >= Prime && ty % Prime) neasy++;
129     int x[2], y[2], z[3];
130     for(int i=0; i<3; i++) z[i] = 0;
131     for(int i=0; i<2; i++)
132       x[i] = tx%Prime, tx /= Prime;
133     for(int i=0; i<2; i++)
134       y[i] = ty%Prime, ty /= Prime;
135     for(int i=0; i<2; i++)
136     for(int j=0; j<2; j++)
137       z[i+j] = (z[i+j] + x[i] * y[j]) % Prime;
138     z[0] += z[2] * wsquare;
139 
140     return m(z[0]) + Prime * m(z[1]);
141     #endif
142     }
143 
sqrhr::fieldpattern::fpattern144   int sqr(int x) { return mul(x,x); }
145 
mmulhr::fieldpattern::fpattern146   matrix mmul(const matrix& A, const matrix& B) {
147     matrix res;
148     for(int i=0; i<MWDIM; i++) for(int k=0; k<MWDIM; k++) {
149       int t = 0;
150   #ifdef EASY
151       for(int j=0; j<MWDIM; j++) t += mul(A[i][j], B[j][k]);
152       t %= Prime;
153   #else
154       for(int j=0; j<MWDIM; j++) t = add(t, mul(A[i][j], B[j][k]));
155   #endif
156       res[i][k] = t;
157       }
158     return res;
159     }
160 
161   map<matrix, int> matcode;
162   vector<matrix> matrices;
163 
164   vector<string> qpaths;
165 
166   vector<matrix> qcoords;
167 
168   // S7 in 2D, but e.g. 4 for a 3D cube
169   int rotations;
170 
171   // S7 in 2D, but e.g. 24 for a 3D cube
172   int local_group;
173 
174   // Id: Identity
175   // R : rotate by 1/rotations of the full circle
176   // P : make a step and turn backwards
177   // X : in 3-dim, turn by 90 degrees
178 
179   matrix Id, R, P, X;
180 
strtomatrixhr::fieldpattern::fpattern181   matrix strtomatrix(string s) {
182     matrix res = Id;
183     matrix m = Id;
184     for(int i=isize(s)-1; i>=0; i--)
185       if(s[i] == 'R') res = mmul(R, res);
186       else if (s[i] == 'P') res = mmul(P, res);
187       else if (s[i] == 'x') { m[0][0] = -1; res = mmul(m, res); m[0][0] = +1; }
188       else if (s[i] == 'y') { m[1][1] = -1; res = mmul(m, res); m[1][1] = +1; }
189       else if (s[i] == 'z') { m[2][2] = -1; res = mmul(m, res); m[2][2] = +1; }
190     return res;
191     }
192 
addashr::fieldpattern::fpattern193   void addas(const matrix& M, int i) {
194     if(!matcode.count(M)) {
195       matcode[M] = i;
196       for(int j=0; j<isize(qcoords); j++)
197         addas(mmul(M, qcoords[j]), i);
198       }
199     }
200 
addhr::fieldpattern::fpattern201   void add(const matrix& M) {
202     if(!matcode.count(M)) {
203       int i = isize(matrices);
204       matcode[M] = i, matrices.push_back(M);
205       for(int j=0; j<isize(qcoords); j++)
206         addas(mmul(M, qcoords[j]), i);
207       if(WDIM == 3) add(mmul(X, M));
208       add(mmul(R, M));
209       }
210     }
211 
212   #define MXF 1000000
213 
214   vector<int> connections;
215 
216   vector<int> inverses; // NYI in 3D
217 
218   // 2D only
219   vector<int> rrf; // rrf[i] equals gmul(i, rotations-1)
220   vector<int> rpf; // rpf[i] equals gmul(i, rotations)
221 
mpowhr::fieldpattern::fpattern222   matrix mpow(matrix M, int N) {
223     while((N&1) == 0) N >>= 1, M = mmul(M, M);
224     matrix res = M;
225     N >>= 1;
226     while(N) {
227       M = mmul(M,M); if(N&1) res = mmul(res, M);
228       N >>= 1;
229       }
230     return res;
231     }
232 
gmulhr::fieldpattern::fpattern233   int gmul(int a, int b) { return matcode[mmul(matrices[a], matrices[b])]; }
gpowhr::fieldpattern::fpattern234   int gpow(int a, int N) { return matcode[mpow(matrices[a], N)]; }
235 
gorderhr::fieldpattern::fpattern236   int gorder(int a) {
237     int b = a;
238     int q = 1;
239     while(b) b = gmul(b, a), q++;
240     return q;
241     }
242 
gmulhr::fieldpattern::fpattern243   pair<int,bool> gmul(pair<int, bool> a, int b) {
244     return make_pair(gmul(a.first,b), a.second);
245     }
246 
247   int order(const matrix& M);
248 
decodepathhr::fieldpattern::fpattern249   string decodepath(int i) {
250     string s;
251     while(i) {
252       if(i % S7) i--, s += 'R';
253       else i = connections[i], s += 'P';
254       }
255     return s;
256     }
257 
258   int orderstats();
259 
260   int cs, sn, ch, sh;
261 
262   int solve();
263 
264   void build();
265 
266   static const int MAXDIST = 120;
267 
268   vector<char> disthep;
269   vector<char> disthex;
270 
271   vector<char> distwall, distriver, distwall2, distriverleft, distriverright, distflower;
272   int distflower0;
273 
274   vector<eItem> markers;
275 
276   int getdist(pair<int,bool> a, vector<char>& dists);
277   int getdist(pair<int,bool> a, pair<int,bool> b);
278   int dijkstra(vector<char>& dists, vector<int> indist[MAXDIST]);
279   void analyze();
280 
281   int maxdist, otherpole, circrad, wallid, wallorder, riverid;
282 
easyhr::fieldpattern::fpattern283   bool easy(int i) {
284     return i < Prime || !(i % Prime);
285     }
286 
287   // 11 * 25
288   // (1+z+z^3) * (1+z^3+z^4) ==
289   // 1+z+z^7 == 1+z+z^2(z^5) == 1+z+z^2(1+z^2) = 1+z+z^2+z^4
290 
inithr::fieldpattern::fpattern291   void init(int p) {
292     Prime = p;
293     if(solve()) {
294       printf("error: could not solve the fieldpattern\n");
295       exit(1);
296       }
297     build();
298     }
299 
fpatternhr::fieldpattern::fpattern300   fpattern(int p) {
301     force_hash = 0;
302     #if CAP_THREAD && MAXMDIM >= 4
303     dis = nullptr;
304     #endif
305     if(!p) return;
306     init(p);
307     }
308 
309   void findsubpath();
310 
311   vector<matrix> generate_isometries();
312 
313   bool check_order(matrix M, int req);
314 
315   unsigned compute_hash();
316 
317   void set_field(int p, int sq);
318 
319   unsigned hashv;
320 
321   #if MAXMDIM >= 4
322   // general 4D
323   vector<transmatrix> fullv;
324 
325   void add1(const matrix& M);
326   void add1(const matrix& M, const transmatrix& Full);
327   vector<matrix> generate_isometries3();
328   int solve3();
329   bool generate_all3();
330 
331   #if CAP_THREAD
332   struct discovery *dis;
333   #endif
334   #endif
335 
336   vector<triplet_info> find_triplets();
337   void generate_quotientgroup();
338   };
339 
340 #if CAP_THREAD && MAXMDIM >= 4
341 struct discovery {
342   fpattern experiment;
343   std::shared_ptr<std::thread> discoverer;
344   std::mutex lock;
345   std::condition_variable cv;
346   bool is_suspended;
347   bool stop_it;
348 
349   map<unsigned, tuple<int, int, matrix, matrix, matrix, int> > hashes_found;
discoveryhr::fieldpattern::discovery350   discovery() : experiment(0) { is_suspended = false; stop_it = false; experiment.dis = this; experiment.Prime = experiment.Field = experiment.wsquare = 0; }
351 
352   void activate();
353   void suspend();
354   void check_suspend();
355   void schedule_destruction();
356   void discovered();
357   ~discovery();
358   };
359 #endif
360 
361 #endif
362 
check_order(matrix M,int req)363 bool fpattern::check_order(matrix M, int req) {
364   matrix P = M;
365   for(int i=1; i<req; i++) {
366     if(P == Id) return false;
367     P = mmul(P, M);
368     }
369   return P == Id;
370   }
371 
generate_isometries()372 vector<matrix> fpattern::generate_isometries() {
373   matrix T = Id;
374   int low = wsquare ? 1-Prime : 0;
375   vector<matrix> res;
376 
377   auto colprod = [&] (int a, int b) {
378     return add(add(mul(T[0][a], T[0][b]), mul(T[1][a], T[1][b])), mul(T[2][a], T[2][b]));
379     };
380 
381   for(T[0][0]=low; T[0][0]<Prime; T[0][0]++)
382   for(T[1][0]=low; T[1][0]<Prime; T[1][0]++)
383   for(T[2][0]=low; T[2][0]<Prime; T[2][0]++)
384   if(colprod(0, 0) == 1)
385   for(T[0][1]=low; T[0][1]<Prime; T[0][1]++)
386   for(T[1][1]=low; T[1][1]<Prime; T[1][1]++)
387   for(T[2][1]=low; T[2][1]<Prime; T[2][1]++)
388   if(colprod(1, 1) == 1)
389   if(colprod(1, 0) == 0)
390   for(T[0][2]=low; T[0][2]<Prime; T[0][2]++)
391   for(T[1][2]=low; T[1][2]<Prime; T[1][2]++)
392   for(T[2][2]=low; T[2][2]<Prime; T[2][2]++)
393   if(colprod(2, 2) == 1)
394   if(colprod(2, 0) == 0)
395   if(colprod(2, 1) == 0)
396     res.push_back(T);
397 
398   return res;
399   }
400 
401 #if MAXMDIM >= 4
generate_isometries3()402 vector<matrix> fpattern::generate_isometries3() {
403 
404   matrix T = Id;
405   int low = wsquare ? 1-Prime : 0;
406   vector<matrix> res;
407 
408   auto colprod = [&] (int a, int b) {
409     return add(add(mul(T[0][a], T[0][b]), mul(T[1][a], T[1][b])), sub(mul(T[2][a], T[2][b]), mul(T[3][a], T[3][b])));
410     };
411 
412   auto rowcol = [&] (int a, int b) {
413     return add(add(mul(T[a][0], T[0][b]), mul(T[a][1], T[1][b])), add(mul(T[a][2], T[2][b]), mul(T[a][3], T[3][b])));
414     };
415 
416   for(T[0][0]=low; T[0][0]<Prime; T[0][0]++)
417   for(T[1][0]=low; T[1][0]<Prime; T[1][0]++)
418   for(T[2][0]=low; T[2][0]<Prime; T[2][0]++)
419   for(T[3][0]=low; T[3][0]<Prime; T[3][0]++)
420   if(colprod(0, 0) == 1)
421   for(T[0][1]=low; T[0][1]<Prime; T[0][1]++)
422   for(T[1][1]=low; T[1][1]<Prime; T[1][1]++)
423   for(T[2][1]=low; T[2][1]<Prime; T[2][1]++)
424   for(T[3][1]=low; T[3][1]<Prime; T[3][1]++)
425   if(colprod(1, 1) == 1)
426   if(colprod(1, 0) == 0) {
427     #if CAP_THREAD && MAXMDIM >= 4
428     if(dis) dis->check_suspend();
429     if(dis && dis->stop_it) return res;
430     #endif
431 
432     for(T[0][2]=low; T[0][2]<Prime; T[0][2]++)
433     for(T[0][3]=low; T[0][3]<Prime; T[0][3]++)
434     if(rowcol(0, 0) == 1)
435     if(rowcol(0, 1) == 0)
436     for(T[1][2]=low; T[1][2]<Prime; T[1][2]++)
437     for(T[1][3]=low; T[1][3]<Prime; T[1][3]++)
438     if(rowcol(1, 0) == 0)
439     if(rowcol(1, 1) == 1)
440     for(T[2][2]=low; T[2][2]<Prime; T[2][2]++)
441     for(T[3][2]=low; T[3][2]<Prime; T[3][2]++)
442     if(colprod(2, 2) == 1)
443     if(colprod(2, 0) == 0)
444     if(colprod(2, 1) == 0)
445     for(T[2][3]=low; T[2][3]<Prime; T[2][3]++)
446     for(T[3][3]=low; T[3][3]<Prime; T[3][3]++)
447     if(rowcol(2, 0) == 0)
448     if(rowcol(2, 1) == 0)
449     if(rowcol(2, 2) == 1)
450     // if(colprod(3, 3) == 1)
451     if(add(colprod(3, 3), 1) == 0)
452     if(colprod(3, 0) == 0)
453     if(colprod(3, 1) == 0)
454     if(colprod(3, 2) == 0)
455     if(rowcol(3, 3) == 1)
456     if(rowcol(3, 0) == 0)
457     if(rowcol(3, 1) == 0)
458     if(rowcol(3, 2) == 0)
459       res.push_back(T);
460     if(isize(res) > limitp) return res;
461     }
462 
463   return res;
464   }
465 
add1(const matrix & M)466 void fpattern::add1(const matrix& M) {
467   if(!matcode.count(M)) {
468     int i = isize(matrices);
469     matcode[M] = i, matrices.push_back(M);
470     }
471   }
472 
add1(const matrix & M,const transmatrix & Full)473 void fpattern::add1(const matrix& M, const transmatrix& Full) {
474   if(!matcode.count(M)) {
475     int i = isize(matrices);
476     matcode[M] = i, matrices.push_back(M), fullv.push_back(Full);
477     }
478   }
479 #endif
480 
481 map<unsigned,int> hash_found;
482 
compute_hash()483 unsigned fpattern::compute_hash() {
484   unsigned hashv = 0;
485   int iR = matcode[R];
486   int iP = matcode[P];
487   int iX = matcode[X];
488   for(int i=0; i<isize(matrices); i++) {
489     hashv = 3 * hashv + gmul(i, iP) + 7 * gmul(i, iR);
490     if(MWDIM == 4) hashv += 11 * gmul(i, iX);
491     }
492   return hashv;
493   }
494 
495 #if MAXMDIM >= 4
generate_all3()496 bool fpattern::generate_all3() {
497 
498   reg3::generate_fulls();
499 
500   matrices.clear();
501   matcode.clear();
502   add1(Id);
503   fullv = {hr::Id};
504   for(int i=0; i<isize(matrices); i++) {
505     add1(mmul(matrices[i], R), fullv[i] * cgi.full_R);
506     add1(mmul(matrices[i], X), fullv[i] * cgi.full_X);
507     }
508   local_group = isize(matrices);
509   for(int i=0; i<(int)matrices.size(); i++) {
510     matrix E = mmul(matrices[i], P);
511     if(!matcode.count(E))
512       for(int j=0; j<local_group; j++) add1(mmul(E, matrices[j]));
513     if(isize(matrices) >= limitv) { println(hlog, "limitv exceeded"); return false; }
514     }
515   hashv = compute_hash();
516   DEBB(DF_FIELD, ("all = ", isize(matrices), "/", local_group, " = ", isize(matrices) / local_group, " hash = ", hashv, " count = ", ++hash_found[hashv]));
517 
518   if(use_quotient_fp)
519     generate_quotientgroup();
520   return true;
521   }
522 
generate_quotientgroup()523 void fpattern::generate_quotientgroup() {
524   int MS = isize(matrices);
525   int best_p = 0, best_i = 0;
526   for(int i=0; i<MS; i++) {
527     int j = i, p = 1;
528     while(j >= local_group)
529       j = gmul(j, i), p++;
530     if(j == 0 && p > best_p) {
531       bool okay = true;
532 
533       vector<bool> visited(MS, false);
534       for(int ii=0; ii<MS; ii++) if(!visited[ii]) {
535         int jj = ii;
536         for(int k=0; k<p; k++) {
537           if(k && jj/local_group == ii/local_group) okay = false;
538           visited[jj] = true;
539           jj = gmul(i, jj);
540           }
541         }
542 
543       if(okay) {
544         bool chk = (MS/p) % local_group;
545         println(hlog, "quotient by ", i, " : ", p, " times less, ", (MS/p/local_group), " tiles, check ", chk);
546         best_p = p; best_i = i;
547         if(chk) {
548           exit(1);
549           }
550         }
551       }
552     }
553 
554   if(best_p > 1) {
555     vector<int> new_id(MS, -1);
556     vector<int> orig_id(MS, -1);
557     vector<matrix> new_matrices;
558     int nv = 0;
559     for(int i=0; i<MS; i++) if(new_id[i] == -1) {
560       int prode = i;
561       for(int l=0; l<local_group; l++) {
562         new_matrices.push_back(matrices[i+l]);
563         }
564       for(int k=0; k<best_p; k++) {
565         for(int l=0; l<local_group; l++) {
566           new_id[gmul(prode, l)] = nv + l;
567           }
568         prode = gmul(best_i, prode);
569         }
570       nv += local_group;
571       }
572     println(hlog, "got nv = ", nv, " / ", local_group);
573 
574     for(int i=0; i<MS; i++)
575       matcode[matrices[i]] = new_id[i];
576     matrices = std::move(new_matrices);
577     println(hlog, "size matrices = ", isize(matrices), " size matcode = ", isize(matcode));
578     println(hlog, tie(P, R, X));
579 
580     /*println(hlog, "TRY AGAIN");
581     generate_quotientgroup();
582     exit(1);*/
583     }
584 
585   }
586 
solve3()587 int fpattern::solve3() {
588   reg3::construct_relations();
589 
590   DEBB(DF_FIELD, ("generating isometries for ", Field));
591 
592   auto iso3 = generate_isometries();
593   auto iso4 = generate_isometries3();
594 
595   int cmb = 0;
596 
597   int N = isize(cgi.rels);
598 
599   vector<int> fails(N);
600 
601   vector<matrix> possible_P, possible_X, possible_R;
602 
603   for(auto& M: iso3) {
604     if(check_order(M, 2))
605       possible_X.push_back(M);
606     if(check_order(M, cgi.r_order))
607       possible_R.push_back(M);
608     }
609   for(auto& M: iso4)
610     if(check_order(M, 2))
611       possible_P.push_back(M);
612 
613   DEBB(DF_FIELD, ("field = ", Field, " #P = ", isize(possible_P), " #X = ", isize(possible_X), " #R = ", isize(possible_R), " r_order = ", cgi.r_order, " xp_order = ", cgi.xp_order));
614 
615   for(auto& xX: possible_X)
616   for(auto& xP: possible_P) if(check_order(mmul(xP, xX), cgi.xp_order))
617   for(auto& xR: possible_R) if(check_order(mmul(xR, xX), cgi.rx_order)) { // if(xR[0][0] == 1 && xR[0][1] == 0)
618     #if CAP_THREAD && MAXMDIM >+ 4
619     if(dis) dis->check_suspend();
620     if(dis && dis->stop_it) return 0;
621     #endif
622     auto by = [&] (char ch) -> matrix& { return ch == 'X' ? xX : ch == 'R' ? xR : xP; };
623     for(int i=0; i<N; i++) {
624       matrix ml = Id;
625       for(char c: cgi.rels[i].first) { ml = mmul(ml, by(c)); if(ml == Id) { fails[i]++; goto bad; }}
626       matrix mr = Id;
627       for(char c: cgi.rels[i].second) { mr = mmul(mr, by(c)); if(mr == Id) { fails[i]++; goto bad; }}
628       if(ml != mr) { fails[i]++; goto bad;}
629       }
630     P = xP; R = xR; X = xX;
631     if(!generate_all3()) continue;
632     #if CAP_THREAD && MAXMDIM >= 4
633     if(dis) { dis->discovered(); continue; }
634     #endif
635     if(force_hash && hashv != force_hash) continue;
636     cmb++;
637     goto ok;
638     bad: ;
639     }
640 
641   ok:
642 
643   DEBB(DF_FIELD, ("cmb = ", cmb, " for field = ", Field));
644   for(int i=0; i<N; i++) if(fails[i]) DEBB(DF_FIELD, (cgi.rels[i], " fails = ", fails[i]));
645 
646   return cmb;
647   }
648 #endif
649 
set_field(int p,int sq)650 void fpattern::set_field(int p, int sq) {
651   Prime = p;
652   Field = sq ? Prime*Prime : Prime;
653   wsquare = sq;
654   for(int a=0; a<MWDIM; a++) for(int b=0; b<MWDIM; b++) Id[a][b] = a==b?1:0;
655   }
656 
solve()657 int fpattern::solve() {
658 
659   for(int a=0; a<MWDIM; a++) for(int b=0; b<MWDIM; b++) Id[a][b] = a==b?1:0;
660 
661   if(!isprime(Prime)) {
662     return 1;
663     }
664 
665   rotations = WDIM == 2 ? S7 : 4;
666   local_group = WDIM == 2 ? S7 : 24;
667 
668   for(dual=0; dual<3; dual++) {
669   for(int pw=1; pw<3; pw++) {
670     if(pw>3) break;
671     Field = pw==1? Prime : Prime*Prime;
672 
673     if(pw == 2) {
674       for(wsquare=1; wsquare<Prime; wsquare++) {
675         int roots = 0;
676         for(int a=0; a<Prime; a++) if((a*a)%Prime == wsquare) roots++;
677         if(!roots) break;
678         }
679       } else wsquare = 0;
680 
681     #if MAXMDIM >= 4
682     if(WDIM == 3) {
683       if(dual == 0 && (Prime <= limitsq || pw == 1)) {
684         int s = solve3();
685         if(s) return 0;
686         }
687       continue;
688       }
689     #endif
690 
691     if(dual == 2) {
692       if(Field <= 10) {
693         vector<matrix> all_isometries = generate_isometries();
694         for(auto& X: all_isometries)
695           if(check_order(X, rotations))
696             for(auto& Y: all_isometries)
697               if(check_order(Y, 2) && check_order(mmul(X, Y), S3)) {
698                 R = X; P = Y;
699                 return 0;
700                 }
701         }
702       continue;
703       }
704 
705 #ifdef EASY
706     std::vector<int> sqrts(Prime, 0);
707     for(int k=1-Prime; k<Prime; k++) sqrts[sqr(k)] = k;
708     int fmax = Prime;
709 #else
710     std::vector<int> sqrts(Field);
711     for(int k=0; k<Field; k++) sqrts[sqr(k)] = k;
712     int fmax = Field;
713 #endif
714 
715     R = P = X = Id;
716     X[1][1] = 0; X[2][2] = 0;
717     X[1][2] = 1; X[2][1] = Prime-1;
718 
719     for(cs=0; cs<fmax; cs++) {
720       int sb = sub(1, sqr(cs));
721 
722       sn = sqrts[sb];
723 
724       R[0][0] = cs; R[1][1] = cs;
725       R[0][1] = sn; R[1][0] = sub(0, sn);
726 
727       if(!check_order(R, dual ? S3 : rotations)) continue;
728 
729       if(R[0][0] == 1) continue;
730 
731       for(ch=2; ch<fmax; ch++) {
732         int chx = sub(mul(ch,ch), 1);
733 
734         sh = sqrts[chx];
735         P[0][0] = sub(0, ch);
736         P[0][WDIM] = sub(0, sh);
737         P[1][1] = Prime-1;
738         P[WDIM][0] = sh;
739         P[WDIM][WDIM] = ch;
740 
741         if(!check_order(mmul(P, R), dual ? rotations : S3)) continue;
742 
743         if(dual) R = mmul(P, R);
744 
745         return 0;
746         }
747       }
748     }
749     }
750 
751   return 2;
752   }
753 
order(const matrix & M)754 int fpattern::order(const matrix& M) {
755   int cnt = 1;
756   matrix Po = M;
757   while(Po != Id) Po = mmul(Po, M), cnt++;
758   return cnt;
759   }
760 
761 EX int triplet_id = 0;
762 
find_triplets()763 vector<triplet_info> fpattern::find_triplets() {
764   int N = isize(matrices);
765   auto compute_transcript = [&] (int i, int j) {
766 
767     vector<int> indices(N, -1);
768     vector<int> transcript;
769     vector<int> q;
770 
771     int qty = 0;
772 
773     auto visit = [&] (int id) {
774       transcript.push_back(indices[id]);
775       if(indices[id] == -1) {
776         indices[id] = isize(q);
777         q.push_back(id);
778         qty++;
779         }
780       };
781 
782     visit(0);
783     for(int x=0; x<isize(q); x++) {
784       int at = q[x];
785       visit(gmul(at, i));
786       visit(gmul(at, j));
787       }
788 
789     transcript.push_back(qty);
790 
791     return transcript;
792     };
793 
794   DEBB(DF_FIELD, ("looking for alternate solutions"));
795   auto orig_transcript = compute_transcript(1, S7);
796 
797   set<vector<int>> transcripts_seen;
798   transcripts_seen.insert(orig_transcript);
799 
800   set<int> conjugacy_classes;
801   vector<int> cc;
802 
803   for(int i=0; i<N; i++) conjugacy_classes.insert(i);
804   for(int i=0; i<N; i++) {
805     if(!conjugacy_classes.count(i)) continue;
806     vector<int> removals;
807     for(int j=0; j<N; j++) {
808       int c = gmul(inverses[j], gmul(i, j));
809       if(c > i) removals.push_back(c);
810       }
811     for(auto r: removals) conjugacy_classes.erase(r);
812     cc.push_back(i);
813     }
814 
815   DEBB(DF_FIELD, ("conjugacy_classes = ", cc));
816 
817   vector<triplet_info> tinf;
818   triplet_info ti;
819   ti.i = 1; ti.j = S7; ti.size = orig_transcript.back();
820   tinf.push_back(ti);
821 
822   for(int i: conjugacy_classes) if(gorder(i) == S7) {
823     DEBB(DF_FIELD, ("checking i=", i));
824     for(int j=1; j<N; j++) if(gorder(j) == 2 && gorder(gmul(i, j)) == S3) {
825       auto t = compute_transcript(i, j);
826       if(!transcripts_seen.count(t)) {
827         transcripts_seen.insert(t);
828         triplet_info ti;
829         ti.i = i; ti.j = j; ti.size = t.back();
830         tinf.push_back(ti);
831         }
832       }
833     }
834 
835   DEBB(DF_FIELD, ("solutions found = ", isize(transcripts_seen)));
836   return tinf;
837   }
838 
build()839 void fpattern::build() {
840 
841   if(WDIM == 3) return;
842 
843   for(int i=0; i<isize(qpaths); i++) {
844     matrix M = strtomatrix(qpaths[i]);
845     qcoords.push_back(M);
846     printf("Solved %s as matrix of order %d\n", qpaths[i].c_str(), order(M));
847     }
848 
849   matcode.clear(); matrices.clear();
850   add(Id);
851   if(isize(matrices) != local_group) { printf("Error: rotation crash #1 (%d)\n", isize(matrices)); exit(1); }
852 
853   connections.clear();
854 
855   for(int i=0; i<(int)matrices.size(); i++) {
856 
857     matrix M = matrices[i];
858 
859     matrix PM = mmul(P, M);
860 
861     add(PM);
862 
863     if(isize(matrices) % local_group) { printf("Error: rotation crash (%d)\n", isize(matrices)); exit(1); }
864 
865     if(!matcode.count(PM)) { printf("Error: not marked\n"); exit(1); }
866 
867     connections.push_back(matcode[PM]);
868     }
869 
870   DEBB(DF_FIELD, ("Computing inverses...\n"));
871   int N = isize(matrices);
872 
873   DEBB(DF_FIELD, ("Number of heptagons: %d\n", N));
874 
875   if(WDIM == 3) return;
876 
877   rrf.resize(N); rrf[0] = S7-1;
878   for(int i=0; i<N; i++)
879     rrf[btspin(i,1)] = btspin(rrf[i], 1),
880     rrf[connections[i]] = connections[rrf[i]];
881 
882   rpf.resize(N); rpf[0] = S7;
883   for(int i=0; i<N; i++)
884     rpf[btspin(i,1)] = btspin(rpf[i], 1),
885     rpf[connections[i]] = connections[rpf[i]];
886 
887   inverses.resize(N);
888   inverses[0] = 0;
889   for(int i=0; i<N; i++) // inverses[i] = gpow(i, N-1);
890     inverses[btspin(i,1)] = rrf[inverses[i]], // btspin(inverses[i],6),
891     inverses[connections[i]] = rpf[inverses[i]];
892 
893   int errs = 0;
894   for(int i=0; i<N; i++) if(gmul(i, inverses[i])) errs++;
895   if(errs) printf("errs = %d\n", errs);
896 
897   if(0) for(int i=0; i<isize(matrices); i++) {
898     printf("%5d/%4d", connections[i], inverses[i]);
899     if(i%S7 == S7-1) printf("\n");
900     }
901 
902   DEBB(DF_FIELD, ("triplet_id = ", triplet_id, " N = ", N));
903   if(triplet_id) {
904     auto triplets = find_triplets();
905     if(triplet_id >= 0 && triplet_id < isize(triplets)) {
906       auto ti = triplets[triplet_id];
907       R = matrices[ti.i];
908       P = matrices[ti.j];
909       dynamicval<int> t(triplet_id, 0);
910       build();
911       DEBB(DF_FIELD, ("triplet built successfully"));
912       return;
913       }
914     }
915 
916   DEBB(DF_FIELD, ("Built.\n"));
917   }
918 
getdist(pair<int,bool> a,vector<char> & dists)919 int fpattern::getdist(pair<int,bool> a, vector<char>& dists) {
920   if(!a.second) return dists[a.first];
921   int m = MAXDIST;
922   int ma = dists[a.first];
923   int mb = dists[connections[btspin(a.first, 3)]];
924   int mc = dists[connections[btspin(a.first, 4)]];
925   m = min(m, 1 + ma);
926   m = min(m, 1 + mb);
927   m = min(m, 1 + mc);
928   if(m <= 2 && ma+mb+mc <= m*3-2) return m-1; // special case
929   m = min(m, 2 + dists[connections[btspin(a.first, 2)]]);
930   m = min(m, 2 + dists[connections[btspin(a.first, 5)]]);
931   m = min(m, 2 + dists[connections[btspin(connections[btspin(a.first, 3)], 5)]]);
932   return m;
933   }
934 
getdist(pair<int,bool> a,pair<int,bool> b)935 int fpattern::getdist(pair<int,bool> a, pair<int,bool> b) {
936   if(a.first == b.first) return a.second == b.second ? 0 : 1;
937   if(b.first) a.first = gmul(a.first, inverses[b.first]), b.first = 0;
938   return getdist(a, b.second ? disthex : disthep);
939   }
940 
dijkstra(vector<char> & dists,vector<int> indist[MAXDIST])941 int fpattern::dijkstra(vector<char>& dists, vector<int> indist[MAXDIST]) {
942   int N = isize(matrices);
943   dists.resize(N);
944   for(int i=0; i<N; i++) dists[i] = MAXDIST-1;
945   int maxd = 0;
946   for(int i=0; i<MAXDIST; i++) while(!indist[i].empty()) {
947     int at = indist[i].back();
948     indist[i].pop_back();
949     if(dists[at] <= i) continue;
950     maxd = i;
951     dists[at] = i;
952     int lg = MWDIM == 4 ? local_group : S7;
953     for(int q=0; q<lg; q++) {
954       dists[at] = i;
955       if(WDIM == 3)
956         indist[i+1].push_back(gmul(at, local_group));
957       else if(PURE) // todo-variation: PURE here?
958         indist[i+1].push_back(connections[at]);
959       else {
960         indist[i+2].push_back(connections[at]);
961         indist[i+3].push_back(connections[btspin(connections[at], 2)]);
962         }
963       at = groupspin(at, 1, lg);
964       }
965     }
966   return maxd;
967   }
968 
analyze()969 void fpattern::analyze() {
970 
971   if(MWDIM == 4) {
972     /* we need to compute inverses */
973     int N = isize(matrices);
974     inverses.resize(N);
975     for(int i=0; i<N; i++) {
976       matrix M = matrices[i];
977       matrix M2 = mpow(M, N-1);
978       inverses[i] = matcode[M2];
979       }
980     }
981 
982   DEBB(DF_FIELD, ("variation = %d\n", int(variation)));
983   int N = isize(connections);
984 
985   markers.resize(N);
986 
987   vector<int> indist[MAXDIST];
988 
989   indist[0].push_back(0);
990   int md0 = dijkstra(disthep, indist);
991 
992   if(MWDIM == 4) return;
993 
994   indist[1].push_back(0);
995   indist[1].push_back(connections[3]);
996   indist[1].push_back(connections[4]);
997   indist[2].push_back(connections[btspin(connections[3], 5)]);
998   indist[2].push_back(connections[2]);
999   indist[2].push_back(connections[5]);
1000   int md1 = dijkstra(disthex, indist);
1001 
1002   maxdist = max(md0, md1);
1003 
1004   otherpole = 0;
1005 
1006   for(int i=0; i<N; i+=S7) {
1007     int mp = 0;
1008     for(int q=0; q<S7; q++) if(disthep[connections[i+q]] < disthep[i]) mp++;
1009     if(mp == S7) {
1010       bool eq = true;
1011       for(int q=0; q<S7; q++) if(disthep[connections[i+q]] != disthep[connections[i]]) eq = false;
1012       if(eq) {
1013         // for(int q=0; q<S7; q++) printf("%3d", disthep[connections[i+q]]);
1014         // printf(" (%2d) at %d\n", disthep[i], i);
1015         if(disthep[i] > disthep[otherpole]) otherpole = i;
1016         // for(int r=0; r<S7; r++) {
1017         // printf("Matrix: "); for(int a=0; a<3; a++) for(int b=0; b<3; b++)
1018         //    printf("%4d", matrices[i+r][a][b]); printf("\n");
1019         //  }
1020         }
1021       }
1022     }
1023 
1024   circrad = 99;
1025 
1026   for(int i=0; i<N; i++) for(int u=2; u<4; u++) if(disthep[i] < circrad)
1027     if(disthep[connections[i]] < disthep[i] && disthep[connections[btspin(i,u)]] < disthep[i])
1028       circrad = disthep[i];
1029 
1030   DEBB(DF_FIELD, ("maxdist = %d otherpole = %d circrad = %d\n", maxdist, otherpole, circrad));
1031 
1032   matrix PRRR = strtomatrix("PRRR");
1033   matrix PRRPRRRRR = strtomatrix("PRRPRRRRR");
1034   matrix PRRRP = strtomatrix("PRRRP");
1035   matrix PRP = strtomatrix("PRP");
1036   matrix PR = strtomatrix("PR");
1037   matrix Wall = strtomatrix("RRRPRRRRRPRRRP");
1038 
1039   wallorder = order(Wall);
1040   wallid = matcode[Wall];
1041 
1042   DEBB(DF_FIELD, ("wall order = %d\n", wallorder));
1043 
1044 #define SETDIST(X, d, it) {int c = matcode[X]; indist[d].push_back(c); if(it == itNone) ; else if(markers[c] && markers[c] != it) markers[c] = itBuggy; else markers[c] = it; }
1045 
1046   matrix W = Id;
1047   for(int i=0; i<wallorder; i++) {
1048     SETDIST(W, 0, itAmethyst)
1049     W = mmul(W, Wall);
1050     }
1051   W = P;
1052   for(int i=0; i<wallorder; i++) {
1053     SETDIST(W, 0, itEmerald)
1054     W = mmul(W, Wall);
1055     }
1056 
1057   int walldist = dijkstra(distwall, indist);
1058   DEBB(DF_FIELD, ("wall dist = %d\n", walldist));
1059 
1060 
1061   W = strtomatrix("RRRRPR");
1062   for(int j=0; j<wallorder; j++) {
1063     W = mmul(W, Wall);
1064     for(int i=0; i<wallorder; i++) {
1065       SETDIST(W, 0, itNone)
1066       SETDIST(mmul(PRRR, W), 1, itNone)
1067       W = mmul(Wall, W);
1068       }
1069     }
1070   dijkstra(distwall2, indist);
1071 
1072   int rpushid = matcode[PRRPRRRRR];
1073   riverid = 0;
1074 
1075   for(int i=0; i<N; i++) {
1076     int j = i;
1077     int ipush = gmul(rpushid, i);
1078     for(int k=0; k<wallorder; k++) {
1079       if(ipush == j) {
1080         DEBB(DF_FIELD, ("River found at %d:%d\n", i, k));
1081         riverid = i;
1082         goto riveridfound;
1083         }
1084       j = gmul(j, wallid);
1085       }
1086     }
1087 
1088   riveridfound: ;
1089 
1090   W = strtomatrix("RRRRPR");
1091   for(int j=0; j<wallorder; j++) {
1092     W = mmul(W, Wall);
1093     for(int i=0; i<wallorder; i++) {
1094       if(i == 7) SETDIST(W, 0, itCoast)
1095       if(i == 3) SETDIST(mmul(PRRRP, W), 0, itWhirlpool)
1096       W = mmul(Wall, W);
1097       }
1098     }
1099   dijkstra(PURE ? distriver : distflower, indist);
1100 
1101   W = matrices[riverid];
1102   for(int i=0; i<wallorder; i++) {
1103     SETDIST(W, 0, itStatue)
1104     W = mmul(W, Wall);
1105     }
1106   W = mmul(P, W);
1107   for(int i=0; i<wallorder; i++) {
1108     SETDIST(W, 0, itSapphire)
1109     W = mmul(W, Wall);
1110     }
1111   W = mmul(PRP, matrices[riverid]);
1112   for(int i=0; i<wallorder; i++) {
1113     SETDIST(W, 1, itShard)
1114     W = mmul(W, Wall);
1115     }
1116   W = mmul(PR, matrices[riverid]);
1117   for(int i=0; i<wallorder; i++) {
1118     SETDIST(W, 1, itGold)
1119     W = mmul(W, Wall);
1120     }
1121   int riverdist = dijkstra(PURE ? distflower : distriver, indist);
1122   DEBB(DF_FIELD, ("river dist = %d\n", riverdist));
1123 
1124   for(int i=0; i<isize(currfp.matrices); i++)
1125     if(currfp.distflower[i] == 0) {
1126       distflower0 = currfp.inverses[i]+1;
1127       break;
1128       }
1129 
1130   if(!PURE) {
1131     W = matrices[riverid];
1132     for(int i=0; i<wallorder; i++) {
1133       SETDIST(W, 0, itStatue)
1134       W = mmul(W, Wall);
1135       }
1136     W = mmul(PR, matrices[riverid]);
1137     for(int i=0; i<wallorder; i++) {
1138       SETDIST(W, 0, itGold)
1139       W = mmul(W, Wall);
1140       }
1141     W = mmul(P, matrices[riverid]);
1142     for(int i=0; i<wallorder; i++) {
1143       SETDIST(W, 1, itSapphire)
1144       W = mmul(W, Wall);
1145       }
1146     dijkstra(distriverleft, indist);
1147     W = mmul(PRP, matrices[riverid]);
1148     for(int i=0; i<wallorder; i++) {
1149       SETDIST(W, 0, itShard)
1150       W = mmul(W, Wall);
1151       }
1152     W = mmul(P, matrices[riverid]);
1153     for(int i=0; i<wallorder; i++) {
1154       SETDIST(W, 0, itSapphire)
1155       W = mmul(W, Wall);
1156       }
1157     W = matrices[riverid];
1158     for(int i=0; i<wallorder; i++) {
1159       SETDIST(W, 1, itStatue)
1160       W = mmul(W, Wall);
1161       }
1162     dijkstra(distriverright, indist);
1163     }
1164   else {
1165     W = strtomatrix("RRRRPR");
1166     for(int j=0; j<wallorder; j++) {
1167       W = mmul(W, Wall);
1168       for(int i=0; i<wallorder; i++) {
1169         if(i == 7) SETDIST(W, 0, itCoast)
1170         W = mmul(Wall, W);
1171         }
1172       }
1173     dijkstra(distriverleft, indist);
1174     W = strtomatrix("RRRRPR");
1175     for(int j=0; j<wallorder; j++) {
1176       W = mmul(W, Wall);
1177       for(int i=0; i<wallorder; i++) {
1178         if(i == 3) SETDIST(mmul(PRRRP, W), 0, itWhirlpool)
1179         W = mmul(Wall, W);
1180         }
1181       }
1182     dijkstra(distriverright, indist);
1183     }
1184 
1185   DEBB(DF_FIELD, ("wall-river distance = %d\n", distwall[riverid]));
1186   DEBB(DF_FIELD, ("river-wall distance = %d\n", distriver[0]));
1187   }
1188 
orderstats()1189 int fpattern::orderstats() {
1190   int N = isize(matrices);
1191 
1192   #define MAXORD 10000
1193   int ordcount[MAXORD];
1194   int ordsample[MAXORD];
1195 
1196   for(int i=0; i<MAXORD; i++) ordcount[i] = 0;
1197 
1198   for(int i=0; i<N; i++) {
1199     int cnt = order(matrices[i]);
1200 
1201     if(cnt < MAXORD) {
1202       if(!ordcount[cnt]) ordsample[cnt] = i;
1203       ordcount[cnt]++;
1204       }
1205     }
1206 
1207   printf("Listing:\n");
1208   for(int i=0; i<MAXORD; i++) if(ordcount[i])
1209     printf("Found %4d matrices of order %3d: %s\n", ordcount[i], i, decodepath(ordsample[i]).c_str());
1210 
1211   return ordsample[Prime];
1212   }
1213 
findsubpath()1214 void fpattern::findsubpath() {
1215   int N = isize(matrices);
1216   for(int i=1; i<N; i++)
1217     if(gpow(i, Prime) == 0) {
1218       subpathid = i;
1219       subpathorder = Prime;
1220       DEBB(DF_FIELD, ("Subpath found: %s\n", decodepath(i).c_str()));
1221       return;
1222       }
1223   }
1224 
1225 fpattern *fp43;
1226 
info()1227 EX void info() {
1228   fpattern fp(0);
1229   int cases = 0, hard = 0;
1230   for(int p=0; p<500; p++) {
1231     fp.Prime = p;
1232     if(fp.solve() == 0) {
1233       printf("%4d: wsquare=%d cs=%d sn=%d ch=%d sh=%d dual=%d\n",
1234         p, fp.wsquare, fp.cs, fp.sn, fp.ch, fp.sh, fp.dual);
1235       cases++;
1236       if(!fp.easy(fp.cs) || !fp.easy(fp.sn) || !fp.easy(fp.ch) || !fp.easy(fp.sn))
1237         hard++;
1238       #ifndef EASY
1239       neasy = 0;
1240       #endif
1241       if(WDIM == 3) continue;
1242       fp.build();
1243       #ifndef EASY
1244       printf("Not easy: %d\n", neasy);
1245       #endif
1246       int N = isize(fp.matrices);
1247       int left = N / fp.Prime;
1248       printf("Prime decomposition: %d = %d", N, fp.Prime);
1249       for(int p=2; p<=left; p++) while(left%p == 0) printf("*%d", p), left /= p;
1250       printf("\n");
1251       printf("Order of RRP is: %d\n", fp.order(fp.strtomatrix("RRP")));
1252       printf("Order of RRRP is: %d\n", fp.order(fp.strtomatrix("RRRP")));
1253       printf("Order of RRRPRRRRRPRRRP is: %d\n", fp.order(fp.strtomatrix("RRRPRRRRRPRRRP")));
1254       }
1255     }
1256   printf("cases found = %d (%d hard)\n", cases, hard);
1257   }
1258 
1259 EX fpattern current_quotient_field = fpattern(0);
1260 EX fpattern fp_invalid = fpattern(0);
1261 EX bool quotient_field_changed;
1262 
1263 // these strings contain \x00
1264 #define STR(x) string(x, sizeof(x))
1265 
getcurrfp()1266 EX struct fpattern& getcurrfp() {
1267   if(fake::in()) return *FPIU(&getcurrfp());
1268   if(geometry == gFieldQuotient && quotient_field_changed)
1269     return current_quotient_field;
1270   if(geometry == gSpace535) {
1271     // 120 cells, hash = 9EF7A9C4
1272     static fpattern fp(5);
1273     return fp;
1274     }
1275   if(geometry == gSpace534) {
1276     // 260 cells, hash = 72414D0C (not 0C62E214)
1277     static fpattern fp(0);
1278     if(fp.Prime) return fp;
1279     // fp.Prime = 5; fp.force_hash = 0x72414D0C; fp.solve();
1280 
1281     if(use_rule_fp) {
1282       fp.Prime = 11; fp.force_hash = 0x5FC4CFF0; fp.solve();
1283       }
1284     else {
1285       shstream ins(STR("\x05\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\xfc\xff\xff\xff\x01\x00\x00\x00\x04\x00\x00\x00\xfc\xff\xff\xff\x04\x00\x00\x00\xfe\xff\xff\xff\x00\x00\x00\x00\x01\x00\x00\x00\xfe\xff\xff\xff\x04\x00\x00\x00\x02\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x03\x00\x00\x00\x04\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x02\x00\x00\x00\xff\xff\xff\xff\x00\x00\x00\x00\x00\x00\x00\x00\xfc\xff\xff\xff\x02\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\xfc\xff\xff\xff\x02\x00\x00\x00\x00\x00\x00\x00\xfc\xff\xff\xff\x01\x00\x00\x00\xfd\xff\xff\xff\x00\x00\x00\x00\x02\x00\x00\x00\xfd\xff\xff\xff\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00"));
1286       hread_fpattern(ins, fp);
1287       }
1288     return fp;
1289     }
1290   if(geometry == gSpace435) {
1291     // 650 cells, hash = EB201050
1292     static fpattern fp(0);
1293     if(fp.Prime) return fp;
1294     // fp.Prime = 5; fp.force_hash = 0xEB201050; fp.solve();
1295     // what is 0x72414D0C??
1296 
1297     if(use_rule_fp) {
1298       fp.Prime = 11; fp.force_hash = 0x65CE0C00; fp.solve();
1299       }
1300     else {
1301       shstream ins(STR("\x05\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\xfc\xff\xff\xff\xff\xff\xff\xff\xfe\xff\xff\xff\xfc\xff\xff\xff\x04\x00\x00\x00\x02\x00\x00\x00\x04\x00\x00\x00\xff\xff\xff\xff\x02\x00\x00\x00\x03\x00\x00\x00\x03\x00\x00\x00\xfd\xff\xff\xff\x01\x00\x00\x00\x02\x00\x00\x00\x03\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\xfc\xff\xff\xff\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xff\xff\xff\xff\x00\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\xfd\xff\xff\xff\xfd\xff\xff\xff\x00\x00\x00\x00\xfd\xff\xff\xff\x02\x00\x00\x00\x03\x00\x00\x00\x00\x00\x00\x00\xfd\xff\xff\xff\x03\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00"));
1302       hread_fpattern(ins, fp);
1303       }
1304     return fp;
1305     }
1306   if(geometry == gSpace436) {
1307     static fpattern fp(0);
1308     // FF82A214
1309     shstream ins(STR("\x05\x00\x00\x00\x02\x00\x00\x00\x01\x00\x00\x00\xfd\xff\xff\xff\x00\x00\x00\x00\xfe\xff\xff\xff\xfd\xff\xff\xff\x01\x00\x00\x00\x01\x00\x00\x00\x03\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x01\x00\x00\x00\x04\x00\x00\x00\xfd\xff\xff\xff\x02\x00\x00\x00\x01\x00\x00\x00\x02\x00\x00\x00\x02\x00\x00\x00\xfc\xff\xff\xff\x00\x00\x00\x00\x00\x00\x00\x00\xff\xff\xff\xff\x02\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\xfd\xff\xff\xff\xfd\xff\xff\xff\x00\x00\x00\x00\xfd\xff\xff\xff\x02\x00\x00\x00\x03\x00\x00\x00\x00\x00\x00\x00\xfd\xff\xff\xff\x03\x00\x00\x00\x02\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00"));
1310     hread_fpattern(ins, fp);
1311     return fp;
1312     }
1313   if(geometry == gSpace336) {
1314     // 672 cells in E3F6B7BC
1315     // 672 cells in 885F1184
1316     // 9408 cells in C4089F34
1317     static fpattern fp(0);
1318     if(fp.Prime) return fp;
1319     // fp.Prime = 7; fp.force_hash = 0xE3F6B7BCu; fp.solve();
1320     shstream ins(STR("\x07\x00\x00\x00\x03\x00\x00\x00\xfa\xff\xff\xff\x02\x00\x00\x00\x03\x00\x00\x00\x06\x00\x00\x00\x02\x00\x00\x00\xfe\xff\xff\xff\xfb\xff\xff\xff\xfc\xff\xff\xff\x03\x00\x00\x00\xfb\xff\xff\xff\xfd\xff\xff\xff\xfb\xff\xff\xff\x01\x00\x00\x00\xfd\xff\xff\xff\xfe\xff\xff\xff\xfd\xff\xff\xff\x03\x00\x00\x00\x00\x00\x00\x00\xfd\xff\xff\xff\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\xfc\xff\xff\xff\x00\x00\x00\x00\x03\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00\x00\x00\x00\x00\xfa\xff\xff\xff\xfb\xff\xff\xff\x00\x00\x00\x00\xfa\xff\xff\xff\x02\x00\x00\x00\x06\x00\x00\x00\x00\x00\x00\x00\xfb\xff\xff\xff\x06\x00\x00\x00\x04\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x00\x01\x00\x00\x00"));
1321     hread_fpattern(ins, fp);
1322     return fp;
1323     }
1324   if(geometry == gSpace344) {
1325     // 600 cells in 558C8ED0
1326     // 2400 cells in AF042EA8
1327     // 2600 cells in D26948E0
1328     // 2600 cells in EC29DCEC
1329     static fpattern fp(0);
1330     if(fp.Prime) return fp;
1331     fp.Prime = 5; fp.force_hash = 0x558C8ED0u; fp.solve();
1332     return fp;
1333     // 4900 cells in CDCC7860 (7)
1334     }
1335   if(geometry == gSpace536) {
1336     static fpattern fp(0);
1337     if(fp.Prime) return fp;
1338     // 130 cells in 3BA5C5A4
1339     // 260 cells in 9FDE7B38
1340     fp.Prime = 5; fp.force_hash = 0x9FDE7B38u; fp.solve();
1341     return fp;
1342     }
1343   if(geometry == gSpace345) {
1344     static fpattern fp(0);
1345     if(fp.Prime) return fp;
1346     // 30 cells in 02ADCAA4 (3^2)
1347     // 650 cells in 7EFE8D98 (5^2)
1348     // 55 cells in F447F75C (11)
1349     fp.Prime = 11; fp.force_hash = 0xF447F75Cu; fp.solve();
1350     return fp;
1351     }
1352   if(geometry == gSpace353) {
1353     static fpattern fp(0);
1354     if(fp.Prime) return fp;
1355     // 130 cells in 1566EBAC (5^2)
1356     // 11 cells in 5A2E2B88 (11)
1357     fp.Prime = 11; fp.force_hash = 0x5A2E2B88u; fp.solve();
1358     return fp;
1359     }
1360   if(geometry == gSpace354) {
1361     static fpattern fp(0);
1362     if(fp.Prime) return fp;
1363     fp.Prime = 11; fp.force_hash = 0x363D8DA4u; fp.solve();
1364     return fp;
1365     }
1366   if(!hyperbolic) return fp_invalid;
1367   if(WDIM == 3 && !quotient && !hybri && !bt::in()) {
1368     static fpattern fp(0);
1369     if(fp.Prime) return fp;
1370     for(int p=2; p<20; p++) { fp.Prime = p; if(!fp.solve()) break; }
1371     DEBB(DF_FIELD, ("set prime = ", fp.Prime));
1372     return fp;
1373     }
1374   if(S7 == 8 && S3 == 3 && !bt::in()) {
1375     static fpattern fp(17);
1376     return fp;
1377     }
1378   if(S7 == 5 && S3 == 4 && !bt::in()) {
1379     static fpattern fp(11);
1380     return fp;
1381     }
1382   if(S7 == 6 && S3 == 4 && !bt::in()) {
1383     static fpattern fp(13);
1384     return fp;
1385     }
1386   if(S7 == 7 && S3 == 4 && !bt::in()) {
1387     static fpattern fp(13);
1388     return fp;
1389     }
1390   if(sphere || euclid) return fp_invalid;
1391   if(S7 == 7 && S3 == 3 && !bt::in()) {
1392     if(!fp43) fp43 = new fpattern(43);
1393     return *fp43;
1394     }
1395   return fp_invalid;
1396   }
1397 
1398 #undef STR
1399 
1400 // todo undefined behavior
1401 EX int subpathid = -1;
1402 EX int subpathorder = -1;
1403 
1404 // extra information for field quotient extra configuration
1405 
1406 EX vector<fgeomextra> fgeomextras = {
1407   fgeomextra(gNormal, 4),
1408   fgeomextra(gOctagon, 1),
1409   fgeomextra(g45, 1),
1410   fgeomextra(g46, 5),
1411   fgeomextra(g47, 1),
1412   fgeomextra(gSchmutzM3, 0),
1413 /*  fgeomextra(gSphere, 0),
1414   fgeomextra(gSmallSphere, 0), -> does not find the prime
1415   fgeomextra(gEuclid, 0),
1416   fgeomextra(gEuclidSquare, 0),
1417   fgeomextra(gTinySphere, 0) */
1418   };
1419 
1420 EX int current_extra = 0;
1421 
nextPrime(fgeomextra & ex)1422 EX void nextPrime(fgeomextra& ex) {
1423   dynamicval<eGeometry> g(geometry, ex.base);
1424   dynamicval<int> t(triplet_id, 0);
1425   int nextprime;
1426   if(isize(ex.primes))
1427     nextprime = ex.primes.back().p + 1;
1428   else
1429     nextprime = 2;
1430   while(true) {
1431     fieldpattern::fpattern fp(0);
1432     fp.Prime = nextprime;
1433     if(fp.solve() == 0) {
1434       fp.build();
1435       int cells = isize(fp.matrices) / S7;
1436       ex.primes.emplace_back(primeinfo{nextprime, cells, (bool) fp.wsquare});
1437       ex.dualval.emplace_back(fp.dual);
1438       break;
1439       }
1440     nextprime++;
1441     }
1442   }
1443 
nextPrimes(fgeomextra & ex)1444 EX void nextPrimes(fgeomextra& ex) {
1445   while(isize(ex.primes) < 6)
1446     nextPrime(ex);
1447   }
1448 
enableFieldChange()1449 EX void enableFieldChange() {
1450   fgeomextra& gxcur = fgeomextras[current_extra];
1451   fieldpattern::quotient_field_changed = true;
1452   nextPrimes(gxcur);
1453   dynamicval<eGeometry> g(geometry, gFieldQuotient);
1454   ginf[geometry].g = ginf[gxcur.base].g;
1455   ginf[geometry].sides = ginf[gxcur.base].sides;
1456   ginf[geometry].vertex = ginf[gxcur.base].vertex;
1457   ginf[geometry].distlimit = ginf[gxcur.base].distlimit;
1458   ginf[geometry].tiling_name = ginf[gxcur.base].tiling_name;
1459   ginf[geometry].default_variation = ginf[gxcur.base].default_variation;
1460   ginf[geometry].flags = qFIELD | qANYQ | qBOUNDED;
1461   fieldpattern::current_quotient_field.init(gxcur.primes[gxcur.current_prime_id].p);
1462   }
1463 
field_from_current()1464 EX void field_from_current() {
1465   auto& go = ginf[geometry];
1466   dynamicval<eGeometry> g(geometry, gFieldQuotient);
1467   auto& gg = ginf[geometry];
1468   gg.sides = go.sides;
1469   gg.vertex = go.vertex;
1470   gg.distlimit = go.distlimit;
1471   gg.tiling_name = go.tiling_name;
1472   gg.flags = go.flags | qANYQ | qFIELD | qBOUNDED;
1473   gg.g = go.g;
1474   gg.default_variation = go.default_variation;
1475   fieldpattern::quotient_field_changed = true;
1476   }
1477 
1478 #if CAP_THREAD && MAXMDIM >= 4
1479 EX map<string, discovery> discoveries;
1480 
activate()1481 void discovery::activate() {
1482   if(!discoverer) {
1483     discoverer = std::make_shared<std::thread> ( [this] {
1484       for(int p=2; p<100; p++) {
1485         experiment.Prime = p;
1486         experiment.solve();
1487         if(stop_it) break;
1488         }
1489       });
1490     }
1491   if(is_suspended) {
1492     if(1) {
1493       std::unique_lock<std::mutex> lk(lock);
1494       is_suspended = false;
1495       }
1496     cv.notify_one();
1497     }
1498   }
1499 
discovered()1500 void discovery::discovered() {
1501   std::unique_lock<std::mutex> lk(lock);
1502   auto& e = experiment;
1503   hashes_found[e.hashv] = make_tuple(e.Prime, e.wsquare, e.R, e.P, e.X, isize(e.matrices) / e.local_group);
1504   }
1505 
suspend()1506 void discovery::suspend() { is_suspended = true; }
1507 
check_suspend()1508 void discovery::check_suspend() {
1509   std::unique_lock<std::mutex> lk(lock);
1510   if(is_suspended) cv.wait(lk, [this] { return !is_suspended; });
1511   }
1512 
schedule_destruction()1513 void discovery::schedule_destruction() { stop_it = true; }
~discovery()1514 discovery::~discovery() { schedule_destruction(); if(discoverer) discoverer->join(); }
1515 #endif
1516 
1517 int hk =
1518 #if CAP_THREAD
1519 #if MAXMDIM >= 4
__anon555eb95f0902null1520   + addHook(hooks_on_geometry_change, 100, [] { for(auto& d:discoveries) if(!d.second.is_suspended) d.second.suspend(); })
__anon555eb95f0a02null1521   + addHook(hooks_final_cleanup, 100, [] {
1522       for(auto& d:discoveries) { d.second.schedule_destruction(); if(d.second.is_suspended) d.second.activate(); }
1523       discoveries.clear();
1524       })
1525 #endif
1526 #endif
1527 #if CAP_COMMANDLINE
__anon555eb95f0b02null1528   + addHook(hooks_args, 0, [] {
1529       using namespace arg;
1530       if(0) ;
1531       else if(argis("-q3-limitsq")) { shift(); limitsq = argi(); }
1532       else if(argis("-q3-limitp")) { shift(); limitp = argi(); }
1533       else if(argis("-q3-limitv")) { shift(); limitv = argi(); }
1534       else return 1;
1535       return 0;
1536       })
1537 #endif
1538   + 0;
1539 
1540 EX purehookset hooks_on_geometry_change;
1541 
field_celldistance(cell * c1,cell * c2)1542 EX int field_celldistance(cell *c1, cell *c2) {
1543   if(geometry != gFieldQuotient) return DISTANCE_UNKNOWN;
1544   if(GOLDBERG) return DISTANCE_UNKNOWN;
1545   auto v1 =fieldpattern::fieldval(c1);
1546   auto v2 =fieldpattern::fieldval(c2);
1547   int d = currfp.getdist(v1, v2);
1548   return d;
1549   }
1550 
1551 EX }
1552 
1553 #define currfp fieldpattern::getcurrfp()
1554 
currfp_gmul(int a,int b)1555 EX int currfp_gmul(int a, int b) { return currfp.gmul(a,b); }
currfp_inverses(int i)1556 EX int currfp_inverses(int i) { return currfp.inverses[i]; }
currfp_distwall(int i)1557 EX int currfp_distwall(int i) { return currfp.distwall[i]; }
currfp_n()1558 EX int currfp_n() { return isize(currfp.matrices); }
currfp_get_R()1559 EX int currfp_get_R() { return currfp.matcode[currfp.R]; }
currfp_get_P()1560 EX int currfp_get_P() { return currfp.matcode[currfp.P]; }
currfp_get_X()1561 EX int currfp_get_X() { return currfp.matcode[currfp.X]; }
1562 
hread_fpattern(hstream & hs,fieldpattern::fpattern & fp)1563 EX void hread_fpattern(hstream& hs, fieldpattern::fpattern& fp) {
1564   hread(hs, fp.Prime);
1565   hread(hs, fp.wsquare);
1566   hread(hs, fp.P);
1567   hread(hs, fp.R);
1568   hread(hs, fp.X);
1569   fp.set_field(fp.Prime, fp.wsquare);
1570   #if MAXMDIM >= 4
1571   fp.generate_all3();
1572   #endif
1573   }
1574 
hwrite_fpattern(hstream & hs,fieldpattern::fpattern & fp)1575 EX void hwrite_fpattern(hstream& hs, fieldpattern::fpattern& fp) {
1576   hwrite(hs, fp.Prime);
1577   hwrite(hs, fp.wsquare);
1578   hwrite(hs, fp.P);
1579   hwrite(hs, fp.R);
1580   hwrite(hs, fp.X);
1581   }
1582 
1583 }
1584 #endif
1585