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