1 /* schreier.c - procedures for manipulating a permutation group using
2 * the random schreier algorithm. There is a separate file schreier.txt
3 * which describes the usage.
4 *
5 * Written for nauty and traces, Brendan McKay 2010-2013.
6 */
7
8 #include "schreier.h"
9
10 TLS_ATTR long long multcount = 0;
11 TLS_ATTR long long filtercount = 0;
12
13 static permnode id_permnode;
14 /* represents identity, no actual content, doesn't need TLS_ATTR */
15 #define ID_PERMNODE (&id_permnode)
16
17 #if !MAXN
18 DYNALLSTAT(int,workperm,workperm_sz);
19 DYNALLSTAT(int,workperm2,workperm2_sz);
20 DYNALLSTAT(int,workpermA,workpermA_sz);
21 DYNALLSTAT(int,workpermB,workpermB_sz);
22 DYNALLSTAT(set,workset,workset_sz);
23 DYNALLSTAT(set,workset2,workset2_sz);
24 #else
25 static TLS_ATTR int workperm[MAXN];
26 static TLS_ATTR int workperm2[MAXN];
27 static TLS_ATTR int workpermA[MAXN];
28 static TLS_ATTR int workpermB[MAXN];
29 static TLS_ATTR set workset[MAXM];
30 static TLS_ATTR set workset2[MAXM];
31 #endif
32
33 static TLS_ATTR schreier *schreier_freelist = NULL;
34 /* Freelist of scheier structures connected by next field.
35 * vec, pwr and orbits fields are assumed allocated. */
36 static TLS_ATTR permnode *permnode_freelist = NULL;
37 /* Freelist of permnode structures connected by next field.
38 * p[] is assumed extended. */
39
40 static TLS_ATTR int schreierfails = SCHREIERFAILS;
41
42 #define TMP
43
44 static boolean filterschreier(schreier*,int*,permnode**,boolean,int,int);
45 #define PNCODE(x) ((int)(((size_t)(x)>>3)&0xFFFUL))
46
47 /* #define TESTP(id,p,n) testispermutation(id,p,n) */
48 #define TESTP(id,p,n)
49
50 /************************************************************************/
51
52 static void
testispermutation(int id,int * p,int n)53 testispermutation(int id, int *p, int n)
54 /* For debugging purposes, crash with a message if p[0..n-1] is
55 not a permutation. */
56 {
57 int i,m;
58 DYNALLSTAT(set,seen,seen_sz);
59
60 for (i = 0; i < n; ++i)
61 if (p[i] < 0 || p[i] > n) break;
62
63 if (i < n)
64 {
65 fprintf(stderr,">E Bad permutation (id=%d): n=%d p[%d]=%d\n",
66 id,n,i,p[i]);
67 exit(1);
68 }
69
70 m = SETWORDSNEEDED(n);
71 DYNALLOC1(set,seen,seen_sz,m,"malloc seen");
72 EMPTYSET(seen,m);
73
74 for (i = 0; i < n; ++i)
75 {
76 if (ISELEMENT(seen,p[i]))
77 {
78 fprintf(stderr,
79 ">E Bad permutation (id=%d): n=%d p[%d]=%d is a repeat\n",
80 id,n,i,p[i]);
81 exit(1);
82 }
83 ADDELEMENT(seen,p[i]);
84 }
85 }
86
87 /************************************************************************/
88
89 int
schreier_fails(int nfails)90 schreier_fails(int nfails)
91 /* Set the number of consecutive failures for filtering;
92 * A value of <= 0 defaults to SCHREIERFAILS.
93 * The function value is the previous setting. */
94 {
95 int prev;
96
97 prev = schreierfails;
98
99 if (nfails <= 0) schreierfails = SCHREIERFAILS;
100 else schreierfails = nfails;
101
102 return prev;
103 }
104
105 /************************************************************************/
106
107 static void
clearfreelists(void)108 clearfreelists(void)
109 /* Clear the schreier and permnode freelists */
110 {
111 schreier *sh,*nextsh;
112 permnode *p,*nextp;
113
114 nextsh = schreier_freelist;
115 while (nextsh)
116 {
117 sh = nextsh;
118 nextsh = sh->next;
119 free(sh->vec);
120 free(sh->pwr);
121 free(sh->orbits);
122 free(sh);
123 }
124 schreier_freelist = NULL;
125
126 nextp = permnode_freelist;
127 while (nextp)
128 {
129 p = nextp;
130 nextp = p->next;
131 free(p);
132 }
133 permnode_freelist = NULL;
134 }
135
136 /************************************************************************/
137
138 static permnode
newpermnode(int n)139 *newpermnode(int n)
140 /* Allocate a new permode structure, with initialized next/prev fields */
141 {
142 permnode *p;
143
144 while (permnode_freelist)
145 {
146 p = permnode_freelist;
147 permnode_freelist = p->next;
148 if (p->nalloc >= n && p->nalloc <= n+100)
149 {
150 p->next = p->prev = NULL;
151 p->mark = 0;
152 return p;
153 }
154 else
155 free(p);
156 }
157
158 p = (permnode*) malloc(sizeof(permnode)+(n-2)*sizeof(int));
159
160 if (p == NULL)
161 {
162 fprintf(ERRFILE,">E malloc failed in newpermnode()\n");
163 exit(1);
164 }
165
166 p->next = p->prev = NULL;
167 p->nalloc = n;
168
169 return p;
170 }
171
172 /************************************************************************/
173
174 static schreier
newschreier(int n)175 *newschreier(int n)
176 /* Allocate a new schreier structure, with initialised next field */
177 {
178 schreier *sh;
179
180 while (schreier_freelist)
181 {
182 sh = schreier_freelist;
183 schreier_freelist = sh->next;
184 if (sh->nalloc >= n && sh->nalloc <= n+100)
185 {
186 sh->next = NULL;
187 return sh;
188 }
189 else
190 {
191 free(sh->vec);
192 free(sh->pwr);
193 free(sh->orbits);
194 free(sh);
195 }
196 }
197
198 sh = (schreier*) malloc(sizeof(schreier));
199
200 if (sh == NULL)
201 {
202 fprintf(ERRFILE,">E malloc failed in newschreier()\n");
203 exit(1);
204 }
205
206 sh->vec = (permnode**) malloc(sizeof(permnode*)*n);
207 sh->pwr = (int*) malloc(sizeof(int)*n);
208 sh->orbits = (int*) malloc(sizeof(int)*n);
209
210 if (sh->vec == NULL || sh->pwr == NULL || sh->orbits == NULL)
211 {
212 fprintf(ERRFILE,">E malloc failed in newschreier()\n");
213 exit(1);
214 }
215
216 sh->next = NULL;
217 sh->nalloc = n;
218
219 return sh;
220 }
221
222 /************************************************************************/
223
224 void
freeschreier(schreier ** gp,permnode ** gens)225 freeschreier(schreier **gp, permnode **gens)
226 /* Free schreier structure and permutation ring. Assume this is everything. */
227 /* Use NULL for arguments which don't need freeing. */
228 {
229 schreier *sh,*nextsh;
230 permnode *p,*nextp;
231
232 if (gp && *gp)
233 {
234 nextsh = *gp;
235 while (nextsh)
236 {
237 sh = nextsh;
238 nextsh = sh->next;
239 sh->next = schreier_freelist;
240 schreier_freelist = sh;
241 }
242 *gp = NULL;
243 }
244
245 if (gens && *gens)
246 {
247 p = *gens;
248 do
249 {
250 nextp = p->next;
251 p->next = permnode_freelist;
252 permnode_freelist = p;
253 p = nextp;
254 } while (p != *gens);
255 *gens = NULL;
256 }
257 }
258
259 /************************************************************************/
260
261 permnode*
findpermutation(permnode * pn,int * p,int n)262 findpermutation(permnode *pn, int *p, int n)
263 /* Return a pointer to permutation p in the circular list,
264 * or NULL if it isn't present. */
265 {
266 permnode *rn;
267 int i;
268
269 if (!pn) return NULL;
270
271 rn = pn;
272 do
273 {
274 for (i = 0; i < n; ++i)
275 if (rn->p[i] != p[i]) break;
276 if (i == n) return rn;
277 rn = rn->next;
278 } while (rn != pn);
279
280 return NULL;
281 }
282
283 /************************************************************************/
284
285 void
addpermutation(permnode ** ring,int * p,int n)286 addpermutation(permnode **ring, int *p, int n)
287 /* Add new permutation to circular list, marked.
288 * and return pointer to it in *ring. */
289 {
290 permnode *pn,*rn;
291
292 pn = newpermnode(n);
293 rn = *ring;
294
295 memcpy(pn->p,p,n*sizeof(int));
296
297 if (!rn)
298 pn->next = pn->prev = pn;
299 else
300 {
301 pn->next = rn->next;
302 pn->prev = rn;
303 rn->next = pn->next->prev = pn;
304 }
305
306 pn->refcount = 0;
307 pn->mark = 1;
308 *ring = pn;
309 }
310
311 /************************************************************************/
312
313 static void
addpermutationunmarked(permnode ** ring,int * p,int n)314 addpermutationunmarked(permnode **ring, int *p, int n)
315 /* Add new permutation to circular list, not marked.
316 * and return pointer to it in *ring. */
317 {
318 TESTP(3,p,n);
319 addpermutation(ring,p,n);
320 (*ring)->mark = 0;
321 }
322
323 /************************************************************************/
324
325 boolean
addgenerator(schreier ** gp,permnode ** ring,int * p,int n)326 addgenerator(schreier **gp, permnode **ring, int *p, int n)
327 /* Add new permutation to group, unless it is discovered to be
328 * already in the group. It is is possible to be in the group
329 * and yet this fact is not discovered.
330 * Return TRUE if the generator (or an equivalent) is added or the
331 * group knowledge with the current partial base is improved. */
332 {
333 TESTP(2,p,n);
334 return filterschreier(*gp,p,ring,FALSE,-1,n);
335 }
336
337 /************************************************************************/
338
339 boolean
condaddgenerator(schreier ** gp,permnode ** ring,int * p,int n)340 condaddgenerator(schreier **gp, permnode **ring, int *p, int n)
341 /* Add new permutation to group, unless it is discovered to be
342 * already in the group. It is is possible to be in the group
343 * and yet this fact is not discovered, but this version will
344 * always notice if this permutation precisely is present.
345 * Return TRUE if the generator (or an equivalent) is added or the
346 * group knowledge with the current partial base is improved. */
347 {
348 TESTP(4,p,n);
349 if (findpermutation(*ring,p,n))
350 return FALSE;
351 else
352 return filterschreier(*gp,p,ring,FALSE,-1,n);
353 }
354
355 /************************************************************************/
356
357 static void
delpermnode(permnode ** ring)358 delpermnode(permnode **ring)
359 /* Delete permnode at head of circular list, making the next node head. */
360 {
361 permnode *newring;
362
363 if (!*ring) return;
364
365 if ((*ring)->next == *ring)
366 newring = NULL;
367 else
368 {
369 newring = (*ring)->next;
370 newring->prev = (*ring)->prev;
371 (*ring)->prev->next = newring;
372 }
373
374 (*ring)->next = permnode_freelist;
375 permnode_freelist = *ring;
376
377 *ring = newring;
378 }
379
380 /************************************************************************/
381
382 void
deleteunmarked(permnode ** ring)383 deleteunmarked(permnode **ring)
384 /* Delete all permutations in the ring that are not marked */
385 {
386 permnode *pn,*firstmarked;
387
388 pn = *ring;
389 firstmarked = NULL;
390
391 while (pn != NULL && pn != firstmarked)
392 {
393 if (pn->mark)
394 {
395 if (!firstmarked) firstmarked = pn;
396 pn = pn->next;
397 }
398 else
399 delpermnode(&pn);
400 }
401
402 *ring = pn;
403 }
404
405 /************************************************************************/
406
407 static void
clearvector(permnode ** vec,permnode ** ring,int n)408 clearvector(permnode **vec, permnode **ring, int n)
409 /* clear vec[0..n-1], freeing permnodes that have no other references
410 * and are not marked */
411 {
412 int i;
413
414 for (i = 0; i < n; ++i)
415 if (vec[i])
416 {
417 if (vec[i] != ID_PERMNODE)
418 {
419 --(vec[i]->refcount);
420 if (vec[i]->refcount == 0 && !vec[i]->mark)
421 {
422 *ring = vec[i];
423 delpermnode(ring);
424 }
425 }
426 vec[i] = NULL;
427 }
428 }
429
430 /************************************************************************/
431
432 static void
initschreier(schreier * sh,int n)433 initschreier(schreier *sh, int n)
434 /* Initialise schreier structure to trivial orbits and empty vector */
435 {
436 int i;
437
438 sh->fixed = -1;
439 for (i = 0; i < n; ++i)
440 {
441 sh->vec[i] = NULL;
442 sh->orbits[i] = i;
443 }
444 }
445
446 /************************************************************************/
447
448 void
newgroup(schreier ** sh,permnode ** ring,int n)449 newgroup(schreier **sh, permnode **ring, int n)
450 /* Make the trivial group, allow for ring to be set elsewhere */
451 {
452 *sh = newschreier(n);
453 initschreier(*sh,n);
454 if (ring) *ring = NULL;
455 }
456
457 /************************************************************************/
458
459 static void
applyperm(int * wp,int * p,int k,int n)460 applyperm(int *wp, int *p, int k, int n)
461 /* Apply the permutation p, k times to each element of wp */
462 {
463 int i,j,cyclen,kk,m;
464
465 TESTP(1,p,n);
466
467 if (k <= 5)
468 {
469 if (k == 0)
470 return;
471 else if (k == 1)
472 for (i = 0; i < n; ++i) wp[i] = p[wp[i]];
473 else if (k == 2)
474 for (i = 0; i < n; ++i) wp[i] = p[p[wp[i]]];
475 else if (k == 3)
476 for (i = 0; i < n; ++i) wp[i] = p[p[p[wp[i]]]];
477 else if (k == 4)
478 for (i = 0; i < n; ++i) wp[i] = p[p[p[p[wp[i]]]]];
479 else if (k == 5)
480 for (i = 0; i < n; ++i) wp[i] = p[p[p[p[p[wp[i]]]]]];
481 }
482 else if (k <= 19)
483 {
484 #if !MAXN
485 DYNALLOC1(int,workpermA,workpermA_sz,n,"applyperm");
486 #endif
487 for (i = 0; i < n; ++i) workpermA[i] = p[p[p[i]]];
488 for (; k >= 6; k -= 6)
489 for (i = 0; i < n; ++i) wp[i] = workpermA[workpermA[wp[i]]];
490 if (k == 1)
491 for (i = 0; i < n; ++i) wp[i] = p[wp[i]];
492 else if (k == 2)
493 for (i = 0; i < n; ++i) wp[i] = p[p[wp[i]]];
494 else if (k == 3)
495 for (i = 0; i < n; ++i) wp[i] = workpermA[wp[i]];
496 else if (k == 4)
497 for (i = 0; i < n; ++i) wp[i] = p[workpermA[wp[i]]];
498 else if (k == 5)
499 for (i = 0; i < n; ++i) wp[i] = p[p[workpermA[wp[i]]]];
500 }
501 else
502 {
503 m = SETWORDSNEEDED(n);
504 #if !MAXN
505 DYNALLOC1(int,workpermA,workpermA_sz,n,"applyperm");
506 DYNALLOC1(int,workpermB,workpermB_sz,n,"applyperm");
507 DYNALLOC1(set,workset2,workset2_sz,m,"applyperm");
508 #endif
509
510 EMPTYSET(workset2,m);
511
512 /* We will construct p^k in workpermB one cycle at a time. */
513
514 for (i = 0; i < n; ++i)
515 {
516 if (ISELEMENT(workset2,i)) continue;
517 if (p[i] == i)
518 workpermB[i] = i;
519 else
520 {
521 cyclen = 1;
522 workpermA[0] = i;
523 for (j = p[i]; j != i; j = p[j])
524 {
525 workpermA[cyclen++] = j;
526 ADDELEMENT(workset2,j);
527 }
528 kk = k % cyclen;
529 for (j = 0; j < cyclen; ++j)
530 {
531 workpermB[workpermA[j]] = workpermA[kk];
532 if (++kk == cyclen) kk = 0;
533 }
534 }
535 }
536 for (i = 0; i < n; ++i) wp[i] = workpermB[wp[i]];
537 }
538 }
539
540 /************************************************************************/
541
542 static boolean
filterschreier(schreier * gp,int * p,permnode ** ring,boolean ingroup,int maxlevel,int n)543 filterschreier(schreier *gp, int *p, permnode **ring,
544 boolean ingroup, int maxlevel, int n)
545 /* Filter permutation p up to level maxlevel of gp.
546 * Use ingroup=TRUE if p is known to be in the group, otherwise
547 * at least one equivalent generator is added unless it is proved
548 * (nondeterministically) that it is in the group already.
549 * maxlevel < 0 means no limit, maxlevel=0 means top level only, etc.
550 * Return TRUE iff some change is made. */
551 {
552 int i,j,j1,j2,lev;
553 int ipwr;
554 schreier *sh;
555 int *orbits,*pwr;
556 permnode **vec,*curr;
557 boolean changed,lchanged,ident;
558 #if !MAXN
559 DYNALLOC1(int,workperm,workperm_sz,n,"filterschreier");
560 #endif
561
562 ++filtercount;
563
564 memcpy(workperm,p,n*sizeof(int));
565
566 if (*ring && p == (*ring)->p)
567 {
568 ingroup = TRUE;
569 curr = *ring;
570 }
571 else
572 curr = NULL;
573
574 /* curr is the location of workperm in ring, if anywhere */
575
576 sh = gp;
577 changed = FALSE;
578 if (maxlevel < 0) maxlevel = n+1;
579
580 for (lev = 0; lev <= maxlevel; ++lev)
581 {
582 for (i = 0; i < n; ++i) if (workperm[i] != i) break;
583 ident = (i == n);
584 if (ident) break;
585
586 lchanged = FALSE;
587 orbits = sh->orbits;
588 vec = sh->vec;
589 pwr = sh->pwr;
590 for (i = 0; i < n; ++i)
591 {
592 j1 = orbits[i];
593 while (orbits[j1] != j1) j1 = orbits[j1];
594 j2 = orbits[workperm[i]];
595 while (orbits[j2] != j2) j2 = orbits[j2];
596
597 if (j1 != j2)
598 {
599 lchanged = TRUE;
600 if (j1 < j2) orbits[j2] = j1;
601 else orbits[j1] = j2;
602 }
603 }
604 if (lchanged)
605 for (i = 0; i < n; ++i) orbits[i] = orbits[orbits[i]];
606
607 if (lchanged) changed = TRUE;
608
609 if (sh->fixed >= 0)
610 {
611 for (i = 0; i < n; ++i)
612 if (vec[i] && !vec[workperm[i]])
613 {
614 changed = TRUE;
615 ipwr = 0;
616 for (j = workperm[i]; !vec[j] ; j = workperm[j]) ++ipwr;
617
618 for (j = workperm[i]; !vec[j] ; j = workperm[j])
619 {
620 if (!curr)
621 {
622 if (!ingroup) addpermutation(ring,workperm,n);
623 else addpermutationunmarked(ring,workperm,n);
624 ingroup = TRUE;
625 curr = *ring;
626 }
627 vec[j] = curr;
628 pwr[j] = ipwr--;
629 ++curr->refcount;
630 }
631 }
632
633 j = workperm[sh->fixed];
634
635 while (j != sh->fixed)
636 {
637 applyperm(workperm,vec[j]->p,pwr[j],n);
638 ++multcount;
639 curr = NULL;
640 j = workperm[sh->fixed];
641 }
642 sh = sh->next;
643 }
644 else
645 break;
646 }
647
648 if (!ident && !ingroup)
649 {
650 changed = TRUE;
651 addpermutation(ring,p,n);
652 }
653
654 return changed;
655 }
656
657 /************************************************************************/
658
659 boolean
expandschreier(schreier * gp,permnode ** ring,int n)660 expandschreier(schreier *gp, permnode **ring, int n)
661 /* filter random elements until schreierfails failures.
662 * Return true if it ever expanded. */
663 {
664 int i,j,nfails,wordlen,skips;
665 boolean changed;
666 permnode *pn;
667 #if !MAXN
668 DYNALLOC1(int,workperm2,workperm2_sz,n,"expandschreier");
669 #endif
670
671 pn = *ring;
672 if (pn == NULL) return FALSE;
673
674 nfails = 0;
675 changed = FALSE;
676
677 for (skips = KRAN(17); --skips >= 0; ) pn = pn->next;
678
679 memcpy(workperm2,pn->p,n*sizeof(int));
680
681 while (nfails < schreierfails)
682 {
683 wordlen = 1 + KRAN(3);
684 for (j = 0; j < wordlen; ++j)
685 {
686 for (skips = KRAN(17); --skips >= 0; ) pn = pn->next;
687 for (i = 0; i < n; ++i) workperm2[i] = pn->p[workperm2[i]];
688 }
689 if (filterschreier(gp,workperm2,ring,TRUE,-1,n))
690 {
691 changed = TRUE;
692 nfails = 0;
693 }
694 else
695 ++nfails;
696 }
697
698 return changed;
699 }
700
701 /************************************************************************/
702
703 int*
getorbits(int * fix,int nfix,schreier * gp,permnode ** ring,int n)704 getorbits(int *fix, int nfix, schreier *gp, permnode **ring, int n)
705 /* Get a pointer to the orbits for this partial base. The pointer
706 * remains valid until pruneset(), getorbits(), getorbitsmin()
707 * or grouporder() is called with an incompatible base (neither a
708 * prefix nor an extension). The contents of the array pointed to
709 * MUST NOT BE MODIFIED by the calling program.
710 */
711 {
712 int k;
713 schreier *sh,*sha;
714
715 sh = gp;
716 for (k = 0; k < nfix; ++k)
717 {
718 if (sh->fixed != fix[k]) break;
719 sh = sh->next;
720 }
721
722 if (k == nfix) return sh->orbits;
723
724 sh->fixed = fix[k];
725 clearvector(sh->vec,ring,n);
726 sh->vec[fix[k]] = ID_PERMNODE;
727
728 for (sha = sh->next; sha ; sha = sha->next) clearvector(sha->vec,ring,n);
729
730 for (++k; k <= nfix; ++k)
731 {
732 if (!sh->next) sh->next = newschreier(n);
733 sh = sh->next;
734 initschreier(sh,n);
735 if (k < nfix)
736 {
737 sh->fixed = fix[k];
738 sh->vec[fix[k]] = ID_PERMNODE;
739 }
740 else
741 sh->fixed = -1;
742 }
743
744 if (*ring) expandschreier(gp,ring,n);
745 return sh->orbits;
746 }
747
748 /************************************************************************/
749
750 int
getorbitsmin(int * fix,int nfix,schreier * gp,permnode ** ring,int ** orbits,int * cell,int ncell,int n,boolean changed)751 getorbitsmin(int *fix, int nfix, schreier *gp, permnode **ring,
752 int **orbits, int *cell, int ncell, int n, boolean changed)
753 /* If the basis elements fix[0..nfix-1] are minimal in their orbits,
754 * as far as we know, return value nfix and set *orbits to point
755 * to orbits fixing fix[0..nfix-1]. If fix[i] is seen to be not
756 * minimal for some i <= nfix-1, return i and set *orbits to point
757 * to orbits fixing fix[0..i-1]. If the partial base is already
758 * known, or fix[0..nfix-1] can already be seen to be non-minimal,
759 * do this work without more filtering. This shortcut is turned
760 * off if changed==TRUE. Otherwise, filter until schreierfails
761 * failures.
762 * The pointer returned remains valid until pruneset(), getorbits(),
763 * getorbitsmin() or grouporder() is called with an incompatible base
764 * (neither a prefix nor an extension). The contents of the array
765 * pointed to MUST NOT BE MODIFIED by the calling program.
766 * If cell != NULL, return early if possible when cell[0..ncell-1]
767 * are all in the same orbit fixing fix[0..nfix-1]. Otherwise
768 * cell,ncell play no part in the computation.
769 */
770 {
771 schreier *sh,*sha;
772 int *fixorbs;
773 int i,j,k,icell,nfails,wordlen,skips;
774 permnode *pn;
775 #if !MAXN
776 DYNALLOC1(int,workperm2,workperm2_sz,n,"expandschreier");
777 #endif
778
779 sh = gp;
780 k = 0;
781 if (!changed)
782 for (k = 0; k < nfix; ++k)
783 {
784 if (sh->orbits[fix[k]] != fix[k])
785 {
786 *orbits = sh->orbits;
787 return k;
788 }
789 if (sh->fixed != fix[k]) break;
790 sh = sh->next;
791 }
792
793 if (k == nfix)
794 {
795 *orbits = sh->orbits;
796 return nfix;
797 }
798
799 sh->fixed = fix[k];
800 clearvector(sh->vec,ring,n);
801 sh->vec[fix[k]] = ID_PERMNODE;
802
803 for (sha = sh->next; sha ; sha = sha->next) clearvector(sha->vec,ring,n);
804
805 for (++k; k <= nfix; ++k)
806 {
807 if (!sh->next) sh->next = newschreier(n);
808 sh = sh->next;
809 initschreier(sh,n);
810 if (k < nfix)
811 {
812 sh->fixed = fix[k];
813 sh->vec[fix[k]] = ID_PERMNODE;
814 }
815 else
816 sh->fixed = -1;
817 }
818 *orbits = fixorbs = sh->orbits;
819
820 if (cell)
821 {
822 for (icell = 1; icell < ncell; ++icell)
823 if (fixorbs[cell[icell]] != fixorbs[cell[0]]) break;
824
825 if (icell >= ncell) return nfix;
826 }
827
828 if (*ring)
829 {
830 pn = *ring;
831
832 nfails = 0;
833
834 for (skips = KRAN(17); --skips >= 0; ) pn = pn->next;
835
836 memcpy(workperm2,pn->p,n*sizeof(int));
837
838 while (nfails < schreierfails)
839 {
840 wordlen = 1 + KRAN(3);
841 for (j = 0; j < wordlen; ++j)
842 {
843 for (skips = KRAN(17); --skips >= 0; ) pn = pn->next;
844 for (i = 0; i < n; ++i) workperm2[i] = pn->p[workperm2[i]];
845 }
846 if (filterschreier(gp,workperm2,ring,TRUE,-1,n))
847 {
848 nfails = 0;
849 sh = gp;
850 for (k = 0; k < nfix; ++k)
851 {
852 if (sh->orbits[fix[k]] != fix[k])
853 {
854 *orbits = sh->orbits;
855 return k;
856 }
857 sh = sh->next;
858 }
859 if (cell)
860 {
861 for ( ; icell < ncell; ++icell)
862 if (fixorbs[cell[icell]] != fixorbs[cell[0]]) break;
863
864 if (icell >= ncell) return nfix;
865 }
866 }
867 else
868 ++nfails;
869 }
870 }
871
872 return nfix;
873 }
874
875 /************************************************************************/
876
877 void
pruneset(set * fixset,schreier * gp,permnode ** ring,set * x,int m,int n)878 pruneset(set *fixset, schreier *gp, permnode **ring, set *x, int m, int n)
879 /* Remove from x any point not minimal for the orbits for this base.
880 * If the base is already known, just provide the orbits without
881 * more filtering. Otherwise, filter until schreierfails failures.
882 */
883 {
884 int i,k;
885 schreier *sh,*sha;
886 int *orbits;
887
888 #if !MAXN
889 DYNALLOC1(set,workset,workset_sz,m,"pruneset");
890 #endif
891 for (i = 0; i < m; ++i) workset[i] = fixset[i];
892
893 sh = gp;
894 while (sh->fixed >= 0 && ISELEMENT(workset,sh->fixed))
895 {
896 DELELEMENT(workset,sh->fixed);
897 sh = sh->next;
898 }
899
900 k = nextelement(workset,m,-1);
901 if (k < 0)
902 orbits = sh->orbits;
903 else
904 {
905 sh->fixed = k;
906 clearvector(sh->vec,ring,n);
907 sh->vec[k] = ID_PERMNODE;
908
909 for (sha = sh->next; sha ; sha = sha->next)
910 clearvector(sha->vec,ring,n);
911
912 while ((k = nextelement(workset,m,k)) >= 0)
913 {
914 if (!sh->next) sh->next = newschreier(n);
915 sh = sh->next;
916 initschreier(sh,n);
917 sh->fixed = k;
918 sh->vec[k] = ID_PERMNODE;
919 }
920 if (!sh->next) sh->next = newschreier(n);
921 sh = sh->next;
922 initschreier(sh,n);
923 sh->fixed = -1;
924
925 if (*ring) expandschreier(gp,ring,n);
926 orbits = sh->orbits;
927 }
928
929 for (k = -1; (k = nextelement(x,m,k)) >= 0; )
930 if (orbits[k] != k) DELELEMENT(x,k);
931 }
932
933 /************************************************************************/
934
935 int
schreier_gens(permnode * ring)936 schreier_gens(permnode *ring)
937 /* Returns the number of generators in the ring */
938 {
939 int j;
940 permnode *pn;
941
942 if (!ring) j = 0;
943 else for (j = 1, pn = ring->next; pn != ring; pn = pn->next) ++j;
944
945 return j;
946 }
947
948 /************************************************************************/
949
950 void
dumpschreier(FILE * f,schreier * gp,permnode * ring,int n)951 dumpschreier(FILE *f, schreier *gp, permnode *ring, int n)
952 /* Dump the whole schreier structure to file f. */
953 {
954 schreier *sh;
955 permnode *pn;
956 int i,j,jj,k;
957
958
959 fprintf(f,"Schreier structure n=%d; ",n);
960
961 jj = -1;
962 for (j = 0, sh = gp; sh; sh = sh->next)
963 {
964 ++j;
965 if (sh->fixed < 0 && jj < 0) jj = j;
966 }
967 fprintf(f," levels=%d (%d used); ",j,jj);
968
969 if (!ring) j = 0;
970 else for (j = 1, pn = ring->next; pn != ring; pn = pn->next) ++j;
971 fprintf(f,"gens=%d; ",j);
972
973 for (j = 0, sh = schreier_freelist; sh; sh = sh->next) ++j;
974 for (k = 0, pn = permnode_freelist; pn; pn = pn->next) ++k;
975 fprintf(f,"freelists: %d,%d\n",j,k);
976
977 if (ring)
978 {
979 fprintf(f,"Generators:\n");
980 pn = ring;
981 do
982 {
983 fprintf(f," %03x ref=%lu mk=%d alloc=%d p=",PNCODE(pn),
984 pn->refcount,pn->mark,pn->nalloc);
985 for (i = 0; i < n; ++i) fprintf(f," %d",pn->p[i]);
986 fprintf(f,"\n");
987 pn = pn->next;
988 } while (pn != ring);
989 }
990
991 if (gp)
992 {
993 fprintf(f,"Levels:\n");
994 for (sh = gp; sh; sh = sh->next)
995 {
996 fprintf(f,"fixed=%2d alloc=%d vec=",sh->fixed,sh->nalloc);
997 for (i = 0; i < n; ++i)
998 {
999 if (sh->vec[i] == ID_PERMNODE) fprintf(f," %d=e",i);
1000 else if (sh->vec[i])
1001 {
1002 k = sh->pwr[i];
1003 j = (sh->vec[i])->p[i];
1004 fprintf(f," %03x",PNCODE(sh->vec[i]));
1005 if (k == 1)
1006 fprintf(f,"(%d,%d)",i,j);
1007 else
1008 {
1009 fprintf(f,"^%d",k);
1010 while (--k >= 1) j = (sh->vec[i])->p[j];
1011 fprintf(f,"(%d,%d)",i,j);
1012 }
1013 }
1014 }
1015 fprintf(f,"\n Orb=");
1016 j = 0;
1017 for (i = 0; i < n; ++i)
1018 {
1019 fprintf(f," %d",sh->orbits[i]);
1020 if (sh->orbits[i] == i) ++j;
1021 }
1022 fprintf(f," [%d]\n",j);
1023 if (sh->fixed < 0) break;
1024 }
1025 }
1026 }
1027
1028 /************************************************************************/
1029
1030 void
grouporder(int * fix,int nfix,schreier * gp,permnode ** ring,double * grpsize1,int * grpsize2,int n)1031 grouporder(int *fix, int nfix, schreier *gp, permnode **ring,
1032 double *grpsize1, int *grpsize2, int n)
1033 /* process the base like in getorbits(), then return the product of the
1034 * orbits along the base, using the largest orbit at the end if the
1035 * base is not complete.
1036 */
1037 {
1038 schreier *sh;
1039 int i,j,k,fx;
1040 int *orb;
1041
1042 #if !MAXN
1043 DYNALLOC1(int,workperm,workperm_sz,n,"grouporder");
1044 #endif
1045
1046 getorbits(fix,nfix,gp,ring,n);
1047 expandschreier(gp,ring,n);
1048 expandschreier(gp,ring,n);
1049 *grpsize1 = 1.0; *grpsize2 = 0;
1050
1051 for (i = 0, sh = gp; i < nfix; ++i, sh = sh->next)
1052 {
1053 orb = sh->orbits;
1054 fx = orb[sh->fixed];
1055 k = 0;
1056 for (j = fx; j < n; ++j) if (orb[j] == fx) ++k;
1057 MULTIPLY(*grpsize1,*grpsize2,k);
1058 }
1059
1060 orb = sh->orbits;
1061 k = 1;
1062 for (i = 0; i < n; ++i)
1063 if (orb[i] == i)
1064 workperm[i] = 1;
1065 else
1066 {
1067 ++workperm[orb[i]];
1068 if (workperm[orb[i]] > k) k = workperm[orb[i]];
1069 }
1070 MULTIPLY(*grpsize1,*grpsize2,k);
1071 }
1072
1073 /*****************************************************************************
1074 * *
1075 * schreier_check() checks that this file is compiled compatibly with the *
1076 * given parameters. If not, call exit(1). *
1077 * *
1078 *****************************************************************************/
1079
1080 void
schreier_check(int wordsize,int m,int n,int version)1081 schreier_check(int wordsize, int m, int n, int version)
1082 {
1083 if (wordsize != WORDSIZE)
1084 {
1085 fprintf(ERRFILE,"Error: WORDSIZE mismatch in schreier.c\n");
1086 exit(1);
1087 }
1088
1089 #if MAXN
1090 if (m > MAXM)
1091 {
1092 fprintf(ERRFILE,"Error: MAXM inadequate in schreier.c\n");
1093 exit(1);
1094 }
1095
1096 if (n > MAXN)
1097 {
1098 fprintf(ERRFILE,"Error: MAXN inadequate in schreier.c\n");
1099 exit(1);
1100 }
1101 #endif
1102
1103 if (version < NAUTYREQUIRED)
1104 {
1105 fprintf(ERRFILE,"Error: schreier.c version mismatch\n");
1106 exit(1);
1107 }
1108 }
1109
1110 /************************************************************************/
1111
1112 void
schreier_freedyn(void)1113 schreier_freedyn(void)
1114 {
1115 #if !MAXN
1116 DYNFREE(workperm,workperm_sz);
1117 DYNFREE(workperm2,workperm2_sz);
1118 DYNFREE(workpermA,workpermA_sz);
1119 DYNFREE(workpermB,workpermB_sz);
1120 DYNFREE(workset,workset_sz);
1121 DYNFREE(workset2,workset2_sz);
1122 #endif
1123 clearfreelists();
1124 }
1125