1 /*********************************************************************
2 *********************************************************************
3 *********************************************************************
4 *
5 * xperm.c
6 *
7 * (C) Jose M. Martin-Garcia 2003-2008.
8 * jmm@iem.cfmac.csic.es, IEM, CSIC, Madrid, Spain.
9 *
10 * This is free software, distributed under the GNU GPL license.
11 * See http://metric.iem.csic.es/Martin-Garcia/xAct/
12 *
13 * These are a collection of C-functions that find Strong Generating
14 * Sets, Coset representatives and Double-coset representatives.
15 * The algorithms are based on Butler and Portugal et al.
16 *
17 * xperm can be used standalone or from the Mathematica package xPerm.
18 *
19 * 20 June - 5 July 2003. Main development.
20 * 3 September 2004. Eliminated MAX variables.
21 * 9 November 2005. Corrected treatment of SGS of group D.
22 * 6 May 2006. All arrays declared dynamically to avoid stack limits.
23 * Thanks to Kasper Peeters for spotting and correcting this.
24 * 25-28 June 2007. Large extension to included multiple dummysets and
25 * repeatedsets.
26 * 2015. Eliminate nested functions. Kasper Peeters.
27 *
28 * Main ideas:
29 * - Permutations are represented using Images notation.
30 * - Generating sets with m n-permutations are stored as lists of
31 * length m*n.
32 * - #define VERBOSE_* to turn on log output. *PPC*
33 *
34 * Comments:
35 * - Permutations are assumed to have degree n>0.
36 * - Lists can have length 0 or positive.
37 *
38 * This is ISO C99, not ANSI-C. There are some gcc extensions:
39 * - ISO C forbids nested functions
40 * - ISO C89 forbids mixed declarations and code
41 * - ISO C90 does not support `long long'
42 * - ISO C90 forbids variable-size arrays
43 *
44 *********************************************************************
45 *********************************************************************
46 *********************************************************************/
47
48 /* KP */
49 #include <stdio.h>
50 #include <stdlib.h>
51 #include <string.h>
52 #include <vector>
53 #include "xperm_new.h"
54 #include <iostream>
55
56 /*********************************************************************
57 * PROTOTYPES *
58 *********************************************************************/
59
60 /* Output */
61 void print_perm(int *p, int n, int nl);
62 void print_array_perm(int *perms, int m, int n, int nl);
63 void print_list(int *list, int n, int nl);
64 void print_array(int *array, int m, int n, int nl);
65
66 /* Lists */
67 int equal_list(int *list1, int *list2, int n);
68 void copy_list(int *list1, int *list2, int n);
69 int position(int i, int *list, int n);
70 int position_list(int *matrix, int m, int *row, int n);
71 void zeros(int *list, int n);
72 void range(int *list, int n);
73 void complement(int *all, int al, int *part, int pl, int n,
74 int *com, int *cl);
75 void sort(int *list, int *slist, int l);
76 void sortB(int *list, int *slist, int l, int *B, int Bl);
77 int minim(int *list, int n);
78 int maxim(int *list, int n);
79 void intersection(int *list1, int l1, int *list2, int l2, int *list,
80 int *l);
81
82 /* Permutations */
83 int isid(int *list, int n );
84 void product(int *p1, int *p2, int *p, int n);
85 void inverse(int *p, int *ip, int n);
86 int onpoints(int point, int *p, int n);
87 void stable_points(int *p, int n, int *list, int *m);
88 int first_nonstable_point(int *p, int n);
89 void nonstable_points(int *list1, int l1, int *GS, int m, int n,
90 int *list2, int *l2);
91 void stabilizer(int *points, int k, int *GS, int m, int n, int *subGS,
92 int *mm);
93 void one_orbit(int point, int *GS, int m, int n, int *orbit, int *ol);
94 void all_orbits(int *GS, int m, int n, int *orbits);
95 void schreier_vector(int point, int *GS, int m, int n, int *nu, int *w);
96 void trace_schreier(int point, int *nu, int *w, int *perm, int n);
97
98 /* Algorithms */
99 long long int order_of_group(int *base, int bl, int *GS, int m, int n);
100 int perm_member(int *p, int *base, int bl, int *GS, int m, int n);
101 void schreier_sims_step(int *base, int bl, int *GS, int m, int n, int i,
102 int *T, int mm, int *newbase, int *nbl, int **newGS, int *nm,
103 int *num);
104 void schreier_sims(int *base, int bl, int *GS, int m, int n,
105 int *newbase, int *nbl, int **newGS, int *nm, int *num);
106 void coset_rep(int *p, int n, int *base, int bl, int *GS, int *m,
107 int *freeps, int fl, int *cr);
108 void SGSD(int *vds, int vdsl, int *dummies, int dl, int *mQ,
109 int *vrs, int vrsl, int *repes, int rl, int n,
110 int firstd, int *KD, int *KDl, int *bD, int *bDl);
111 void double_coset_rep(int *g, int n, int *base, int bl, int *GS, int m,
112 int *vds, int vdsl, int *dummies, int dl, int *mQ,
113 int *vrs, int vrsl, int *repes, int rl, int *dcr);
114 void canonical_perm(int *perm,
115 int SGSQ, int *base, int bl, int *GS, int m, int n,
116 int *freeps, int fl, int *dummyps, int dl, int ob, int metricQ,
117 int *cperm);
118 void canonical_perm_ext(int *perm, int n,
119 int SGSQ, int *base, int bl, int *GS, int m,
120 int *frees, int fl,
121 int *vds, int vdsl, int *dummies, int dl, int *mQ,
122 int *vrs, int vrsl, int *repes, int rl,
123 int *cperm);
124
125
126 /*********************************************************************
127 * PRINTING FUNCTIONS *
128 *********************************************************************/
129
130 /**********************************************************************/
131
132 /* print_perm. JMM, 22 June 2003
133 *
134 * This function prints a permutation p of degree n. If nl=1(0) adds
135 * (does not) a newline. */
136
print_perm(int * p,int n,int nl)137 void print_perm(int *p, int n, int nl)
138 {
139 int i;
140 if (isid(p,n)) printf("id");
141 else {
142 printf("(");
143 printf("%d", p[0]); /* No comma */
144 for (i=1; i<n; i++) {
145 printf(",%d", p[i]); /* Comma separated */
146 }
147 printf(")");
148 }
149 if (nl) printf("\n");
150 }
151
152 /**********************************************************************/
153
154 /* print_array_perm. JMM, 22 June 2003
155 *
156 * This function prints an array of m n-permutations.
157 * If nl=1 (0) adds (does not) a newline after each row.
158 * There are no commas between permutations. */
159
print_array_perm(int * perms,int m,int n,int nl)160 void print_array_perm(int *perms, int m, int n, int nl)
161 {
162 int j;
163 printf("{");
164 if (nl) printf("\n");
165 for(j=0; j<m; j++) {
166 printf(" ");
167 print_perm(perms+j*n, n, nl);
168 }
169 if (nl) printf("}\n");
170 else printf(" }\n");
171 }
172
173 /**********************************************************************/
174
175 /* print_list. JMM, 22 June 2003
176 *
177 * This function prints a list of length n in curly brackets.
178 * If nl=1 (0) it adds (does not) a newline. */
179
print_list(int * list,int n,int nl)180 void print_list(int *list, int n, int nl)
181 {
182 int i;
183 printf("{");
184 if (n>0) printf("%d", list[0]); /* No comma */
185 for (i=1; i<n; i++) printf(",%d", list[i]);/* Comma separated */
186 printf("}");
187 if (nl) printf("\n");
188 }
189
190 /**********************************************************************/
191
192 /* print_array. JMM, 22 June 2003
193 *
194 * This function prints an array of dimensions m x n in curly brackets.
195 * If nl=1 (0) adds (does not) a newline after each row.
196 * There are no commas between lists. */
197
print_array(int * array,int m,int n,int nl)198 void print_array(int *array, int m, int n, int nl)
199 {
200 int j;
201 printf("{");
202 if(nl) printf("\n");
203 for(j=0; j<m; j++) {
204 printf(" ");
205 print_list(array+j*n, n, nl);
206 }
207 if (!nl) printf(" ");
208 printf("}\n");
209 }
210
211 /**********************************************************************/
212
213
214 /*********************************************************************
215 * GENERIC FUNCTIONS *
216 *********************************************************************/
217
218 /**********************************************************************/
219
220 /* equal_list. JMM, 27 June 2003
221 *
222 * Checks that list1 and list2 (both with length n) are equal.
223 * If n=0 the result is 1.
224 */
225
226 /* Old code
227 int equal_list(int *list1, int *list2, int n) {
228 while(n--) { * Run from n-1 to 0 *
229 if(*(list1+n) != *(list2+n)) return(0); * different *
230 }
231 return(1); * equal *
232 }
233 */
234
235 /* KP, 7 May 2006 */
equal_list(int * list1,int * list2,int n)236 int equal_list(int *list1, int *list2, int n)
237 {
238 if (n==0) return 1;
239 if (memcmp(list1, list2, n*sizeof(int))==0) return 1;
240 else return 0;
241 }
242
243 /**********************************************************************/
244
245 /* copy_list. JMM, 27 June 2003
246 *
247 * Copies first n elements of list1 onto list2. If n=0 nothing is done.
248 */
249
250 /* Old code
251 void copy_list(int *list1, int *list2, int n) {
252 while(n--) *(list2+n) = *(list1+n);
253 }
254 */
255
256 /* KP, 7 May 2006 */
257 /* KP, 26 Nov 2007: memcpy changed to memmove */
copy_list(int * list1,int * list2,int n)258 void copy_list(int *list1, int *list2, int n)
259 {
260 if (n==0) return;
261 memmove(list2, list1, n*sizeof(int));
262 }
263
264 /**********************************************************************/
265
266 /* position. JMM, 20 June 2003
267 *
268 * This function finds the position of the integer i in list (length n),
269 * counting from 1 to n, or else gives 0 if i is not in list. Note that
270 * it starts from the end so that if i appears several times in list it
271 * will only find the last one. If n=0 the result is 0.
272 *
273 * More than 60% of time in the Schreier-Sims algorithm is spent here!
274 *
275 * Note: it is not possible to use memchr because it only works with
276 * characters (integers up to 255).
277 */
278
position(int i,int * list,int n)279 int position(int i, int *list, int n)
280 {
281 while(n--) {
282 if (*(list+n) == i) return(n+1);
283 }
284 return(0);
285 }
286
287 /**********************************************************************/
288
289 /* position_list. JMM, 28 June 2003
290 *
291 * This function gives the position of list row (length n) in the
292 * matrix of dimensions m, n. If row is not present in matrix, 0
293 * is returned.
294 */
295
position_list(int * matrix,int m,int * row,int n)296 int position_list(int *matrix, int m, int *row, int n)
297 {
298 while(m--) {
299 if(equal_list(matrix+m*n, row, n)) return(m+1);
300 }
301 return(0);
302 }
303
304 /**********************************************************************/
305
306 /* zeros, range. JMM, 22 June 2003
307 *
308 * These functions generate lists of n zeros or the range 1...n,
309 * respectively, storing the result in list.
310 */
311
312 /* Old code
313 void zeros(int *list, int n) {
314 while(n--) *(list+n) = 0;
315 }
316 */
317
318 /* KP, 7 May 2006 */
zeros(int * list,int n)319 void zeros(int *list, int n)
320 {
321 if (n==0) return;
322 memset(list, 0, n*sizeof(int));
323 }
324
range(int * list,int n)325 void range(int *list, int n)
326 {
327 while(n--) *(list+n) = n+1;
328 }
329
330 /**********************************************************************/
331
332 /* complement. JMM, 28 June 2003
333 *
334 * This function gives all sublists of length n in list all (length al)
335 * that are not members of list part (length pl), storing the result
336 * in list com (length cl).
337 * We assume that enough space (typically al n-lists) has been already
338 * allocated for com.
339 */
340
complement(int * all,int al,int * part,int pl,int n,int * com,int * cl)341 void complement(int *all, int al, int *part, int pl, int n,
342 int *com, int *cl)
343 {
344 int i;
345 *cl = 0;
346 for(i=0; i<al; i++) {
347 if (!position_list(part, pl, all+i*n, n)) {
348 copy_list(all+i*n, com+(*cl)*n, n);
349 (*cl)++;
350 }
351 }
352 }
353
354 /**********************************************************************/
355
356 /* sort. JMM, 30 June 2003
357 *
358 * Sorts list (length l) storing the result in slist (length l, too).
359 * There are l*(l-1)/2 element comparisons. This could be improved.
360 */
361
sort(int * list,int * slist,int l)362 void sort(int *list, int *slist, int l)
363 {
364
365 int i,j;
366 int tmp, mini; /* mini is a position, not an element of list */
367
368 #ifdef VERBOSE_LISTS /*PPC*/
369 printf("sort: Sorting list ");
370 print_list(list, l, 1); /*PPC*/
371 #endif /*PPC*/
372 copy_list(list, slist, l);
373 for (i=0; i<l-1; i++) {
374 /* Assume points at positions 0...i-1 are already sorted */
375 /* Compare point at position i with points to the right */
376 mini = i;
377 for(j=i+1; j<l; j++) {
378 if ( slist[j] < slist[mini] ) mini = j;
379 }
380 /* Swap points i and mini */
381 tmp = slist[i];
382 slist[i] = slist[mini];
383 slist[mini] = tmp;
384 }
385 #ifdef VERBOSE_LISTS /*PPC*/
386 printf("sort: with result ");
387 print_list(slist, l, 1); /*PPC*/
388 #endif /*PPC*/
389 }
390
391 /* sortB. JMM, 30 June 2003
392 *
393 * Sorts list (length l) according to the order given by B (length Bl),
394 * storing the result in slist (length l, too). Elements not in B
395 * are pushed to the end, and sorted using sort.
396 */
397
sortB(int * list,int * slist,int l,int * B,int Bl)398 void sortB(int *list, int *slist, int l, int *B, int Bl)
399 {
400
401 int sl;
402 int *tmp= (int*)malloc(l*sizeof(int)), tmpl;
403 int *stmp= (int*)malloc(l*sizeof(int));
404
405 #ifdef VERBOSE_LISTS /*PPC*/
406 printf("sortB: Sorting list ");
407 print_list(list, l, 1); /*PPC*/
408 printf("sortB: using base ");
409 print_list(B, Bl, 1); /*PPC*/
410 #endif /*PPC*/
411 /* Elements of list in B, keeping order of B, moved to slist */
412 intersection(B, Bl, list, l, slist, &sl);
413 /* Other elements of list are appended to slist */
414 complement(list, l, B, Bl, 1, tmp, &tmpl);
415 /* Check that tmpl+sl==l */
416 if(tmpl+sl != l) printf("Error in sortB\n");
417 /* Sort the latter integers using sort */
418 sort(tmp, stmp, tmpl);
419 copy_list(stmp, slist+sl, tmpl);
420 #ifdef VERBOSE_LISTS /*PPC*/
421 printf("sortB: with result ");
422 print_list(slist, l, 1); /*PPC*/
423 #endif /*PPC*/
424 /* Free allocated memory */
425 free(tmp);
426 free(stmp);
427
428 }
429
430 /**********************************************************************/
431
432 /* minim, maxim. JMM, 2 July 2003
433 *
434 * These functions find the minimum and maximum elements of a given
435 * list of length n, using normal integer ordering. Names min, max are
436 * reserved.
437 */
438
minim(int * list,int n)439 int minim(int *list, int n)
440 {
441
442 int m=list[n-1];
443
444 while(n--) {
445 if (list[n]<m) m=list[n];
446 }
447 return(m);
448
449 }
450
maxim(int * list,int n)451 int maxim(int *list, int n)
452 {
453
454 int m=list[n-1];
455
456 while(n--) {
457 if (list[n]>m) m=list[n];
458 }
459 return(m);
460
461 }
462
463 /**********************************************************************/
464
465 /* intersection. JMM, 3 July 2003
466 *
467 * This function returns in list (length l) the intersection of the
468 * elements of lists list1 (length l1) and list2 (length l2). Note that
469 * we return unique (that is, non-repeated) elements. Returned elements
470 * keep the original order in list1.
471 * We assume that enough space (typical l1 or l2 integers) has been
472 * already allocated for list.
473 */
474
intersection(int * list1,int l1,int * list2,int l2,int * list,int * l)475 void intersection(int *list1, int l1, int *list2, int l2, int *list,
476 int *l)
477 {
478
479 int i, j, a;
480
481 *l = 0;
482 for(i=0; i<l1; i++) { /* Range over elements of list1 */
483 a = list1[i];
484 for(j=0; j<l2; j++) { /* Range over elements of list2 */
485 /* a is in list2 and not yet in list, append */
486 if ((list2[j] == a) && (!position(a, list, *l)))
487 list[(*l)++] = a;
488 }
489 }
490
491 }
492
493 /**********************************************************************/
494
495
496 /*********************************************************************
497 * PERMUTATIONS FUNCTIONS *
498 *********************************************************************/
499
500 /**********************************************************************/
501
502 /* isid. JMM, 20 June 2003
503 *
504 * This function detects whether the n-permutation p is the identity
505 * (returning 1) or not (returning 0). */
506
isid(int * p,int n)507 int isid(int *p, int n )
508 {
509
510 while(n--) {
511 if (*(p+n)!=n+1) return(0); /* Not the identity */
512 }
513 return(1); /* Identity */
514
515 }
516
517 /**********************************************************************/
518
519 /* product. JMM, 20 June 2003
520 *
521 * Product of n-permutations p1 and p2 (from left to right) storing the
522 * result in the n-permutation p. Note that the points of the
523 * permutations go from 1 to n, while we want to add from 0 to n-1 to
524 * the pointer p2. That's why we need to add -1.
525 * This function is highly efficient, but rather criptic. */
526
product(int * p1,int * p2,int * p,int n)527 void product(int *p1, int *p2, int *p, int n)
528 {
529
530 while(n--) *(p++) = *(p2-1+*(p1++));
531
532 /* Example:
533 * Suppose we have p1=(3,2,1) and p2=(3,1,2) of degree n=3.
534 * The following operations take place, in this order:
535 * n=3, n--
536 * *(p1++)=3
537 * *(p2-1+3)=*(p2+2)=2
538 * *(p++)=2 Computed first element of p
539 * n=2, n--
540 * *(p1++)=2
541 * *(p2-1+2)=*(p2+1)=1
542 * *(p++)=1 Computed second element of p
543 * n=1, n--
544 * *(p1++)=1
545 * *(p2-1+1)=*(p2)=3
546 * *(p++)=3 Computed third element of p
547 * n=0, END
548 * Therefore the final permutation returned is p=(2,1,3).
549 */
550
551 }
552
553 /* Taken out KP */
554 /* TAB1 is an element of S; TAB2 is an element of D */
F2(int * TAB1,int * g,int * TAB2,int * sgd,int n)555 void F2(int *TAB1, int *g, int *TAB2, int *sgd, int n)
556 {
557 int *tmp= (int*)malloc(n*sizeof(int));
558
559 product(TAB1, g, tmp, n);
560 product(tmp, TAB2, sgd, n);
561
562 free(tmp);
563 } /* End of function F2 */
564
565 /**********************************************************************/
566
567 /* inverse. JMM, 20 June 2003
568 *
569 * Inverse of the n-permutation p, storing the result in the
570 * n-permutation ip. Again we need to add -1 to the pointer ip.
571 * Conversely, note that we need to add +1 to n because n ranges from
572 * 0 to n-1 in the while loop. */
573
inverse(int * p,int * ip,int n)574 void inverse(int *p, int *ip, int n)
575 {
576
577 while(n--) *(ip-1+*(p+n)) = n+1 ;
578
579 }
580
581 /**********************************************************************/
582
583 /* onpoints. JMM, 20 June 2003
584 *
585 * Image of point under the n-permutation p. If point is larger than
586 * the degree n we return point. */
587
onpoints(int point,int * p,int n)588 int onpoints(int point, int *p, int n)
589 {
590
591 if (point <= n) return(*(p-1+point));
592 else return(point);
593
594 }
595
596 /**********************************************************************/
597
598 /* stable_points. JMM, 20 June 2003
599 *
600 * This function finds the points which remain stable under the
601 * n-permutation p, storing the result in list of length m (<=n).
602 * With n=0 we get an empty list (m=0).
603 * We assume that enough space (typically n integers) has been already
604 * allocated for list.
605 */
606
stable_points(int * p,int n,int * list,int * m)607 void stable_points(int *p, int n, int *list, int *m)
608 {
609
610 int i;
611
612 *m = 0;
613 for(i=1; i<=n; i++) {
614 if ( onpoints(i, p, n) == i ) list[(*m)++] = i;
615 }
616
617 }
618
619 /**********************************************************************/
620
621 /* first_nonstable_point. JMM, 30 June 2003
622 *
623 * This function gives the smallest moved point by the permutation p.
624 * If no point is moved the result is 0. If n=0 the result is 0.
625 * This is a restricted form of stable_points.
626 */
627
first_nonstable_point(int * p,int n)628 int first_nonstable_point(int *p, int n)
629 {
630
631 int i;
632
633 for(i=1; i<=n; i++) {
634 if (onpoints(i, p, n) != i) return(i);
635 }
636 return(0);
637 }
638
639 /**********************************************************************/
640
641 /* nonstable_points. JMM, 30 June 2003
642 *
643 * This function gives a set of points (list2, length l2) such that
644 * none of the permutations in GS fixes them all. The first l1 points
645 * of list2 will be those of list1, even if they are all stable under
646 * GS. If m=0 the resulting list2 is just list1.
647 * We assume that enough space (typically n integers) has been already
648 * allocated for list2.
649 */
650
nonstable_points(int * list1,int l1,int * GS,int m,int n,int * list2,int * l2)651 void nonstable_points(int *list1, int l1, int *GS, int m, int n,
652 int *list2, int *l2)
653 {
654
655 int i, j;
656
657 copy_list(list1, list2, l1);
658 *l2 = l1; /* Initialize list2 */
659 for(j=0; j<m; j++) { /* Range over permutations of GS */
660 int stable=1;
661 for(i=0; i<*l2; i++) {
662 if (onpoints(list2[i], GS+j*n, n) != list2[i]) {
663 stable=0;
664 break;
665 }
666 }
667 /* If all points already in list2 are stable under the
668 permutation, append the smallest nonstable point */
669 if (stable) {
670 list2[*l2] = first_nonstable_point(GS+j*n, n);
671 (*l2)++;
672 }
673 }
674 }
675
676 /**********************************************************************/
677
678 /* stabilizer. JMM, 22 June 2003
679 *
680 * This function finds the permutations that stabilize the set of k
681 * points among the GeneratingSet GS with m n-permutations. With m=0
682 * we return an empty subset (mm=0). Recall that the stabilizer of a
683 * generating set of a group is not, in general, a generating set for
684 * the corresponding stabilizer of that group.
685 * We assume that enough space (typically m*n integers) has been already
686 * allocated for subGS.
687 */
688
stabilizer(int * points,int k,int * GS,int m,int n,int * subGS,int * mm)689 void stabilizer(int *points, int k, int *GS, int m, int n,
690 int *subGS, int *mm)
691 {
692
693 int i, j;
694 int stable;
695
696 *mm=0;
697 for(j=0; j<m; j++) {
698 stable=1;
699 for(i=0; i<k; i++) {
700 if (onpoints(points[i],GS+n*j,n) != points[i]) {
701 stable = 0;
702 break;
703 }
704 }
705 if (stable) {
706 copy_list(GS+n*j, subGS+(*mm)*n, n);
707 ++(*mm);
708 }
709 }
710
711 }
712
713 /**********************************************************************/
714
715 /* one_orbit. JMM, 22 June 2003
716 *
717 * Orbit of a given point under a GeneratingSet GS with m
718 * n-permutations. The result is stored in orbit (length ol).
719 * With m=0 the orbit is just {point}.
720 * We assume that enough space (typically n integers) has been already
721 * allocated for orbit.
722 */
723
one_orbit(int point,int * GS,int m,int n,int * orbit,int * ol)724 void one_orbit(int point, int *GS, int m, int n, int *orbit, int *ol)
725 {
726
727 int np=0; /* Index of current element in the orbit */
728 int gamma; /* Current element in the orbit */
729 int mp; /* Index of current permutation in GS */
730 int newgamma;
731
732 orbit[0] = point;
733 *ol = 1;
734 while(np < *ol) {
735 gamma = orbit[np];
736 for(mp=0; mp<m; mp++) {
737 newgamma = onpoints(gamma, GS+mp*n, n);
738 if (!position(newgamma, orbit, *ol))
739 orbit[(*ol)++] = newgamma;
740 }
741 np++;
742 }
743
744 }
745
746 /**********************************************************************/
747
748 /* all_orbits. JMM, 1 July 2003
749 *
750 * This function returns all orbits of the generating set GS (with
751 * m n-permutations). Note that we use a very special notation: the
752 * result is a vector orbits of length n such that all points marked
753 * with 1 belong to the same orbit; all points marked with 2 belong to
754 * another orbit, and so on. With m=0 we return the list {1,2,...,n}.
755 */
756
all_orbits(int * GS,int m,int n,int * orbits)757 void all_orbits(int *GS, int m, int n, int *orbits)
758 {
759
760 int i, j; /* Counters */
761 int *orbit= (int*)malloc(n*sizeof(int)), ol;/* Computed orbit */
762 int orbit_index=1; /* Orbit index */
763
764 /* Initialize orbits */
765 memset(orbits, 0, n*sizeof(int));
766
767 /* Compute orbits */
768 for (i=1; i<=n; i++) { /* Points */
769 if(orbits[i-1]==0) { /* New orbit condition */
770 /* Compute new orbit */
771 one_orbit(i, GS, m, n, orbit, &ol);
772 /* Mark points of new orbit with index or */
773 for (j=0; j<ol; j++) orbits[orbit[j]-1] = orbit_index;
774 /* Increment index for next orbit */
775 orbit_index++;
776 }
777 }
778
779 /* Free allocated memory */
780 free(orbit);
781
782 }
783
784 /**********************************************************************/
785
786 /* one_schreier_orbit. JMM, 22 June 2003
787 *
788 * Orbit of point under the GeneratingSet GS with m n-permutations.
789 * The result is stored in orbit and the vectors nu of permutations and
790 * w of backward points. If init=1 both nu and w are reset to 0.
791 *
792 * Note: the profiler shows that roughly half of the time in this
793 * function is spent in the subroutine `position'.
794 */
795
one_schreier_orbit(int point,int * GS,int m,int n,int * orbit,int * ol,int * nu,int * w,int init)796 void one_schreier_orbit(int point, int *GS, int m, int n,
797 int *orbit, int *ol, int *nu, int *w, int init)
798 {
799
800 int np; /* Index of current element in the orbit */
801 int gamma; /* Current element in the orbit */
802 int mp; /* Index of current permutation in GS */
803 int *perm= (int*)malloc(n*sizeof(int));
804 int newgamma;
805
806 /* Initialize schreier with zeros if required */
807 memset(orbit, 0, n*sizeof(int));
808 if (init) {
809 memset(nu, 0, n*n*sizeof(int));
810 memset(w, 0, n*sizeof(int));
811 }
812 /* First element of orbit. There is no backward pointer */
813 orbit[0] = point;
814 *ol = 1;
815 /* Other elements of orbit */
816 np = 0;
817 while(np < *ol) {
818 gamma = orbit[np];
819 for(mp=0; mp<m; mp++) {
820 copy_list(GS+mp*n, perm, n);
821 newgamma = onpoints(gamma, perm, n);
822 if (position(newgamma, orbit, *ol));
823 else {
824 /* Append to orbit */
825 orbit[(*ol)++] = newgamma;
826 /* Perm moving gamma to newgamma */
827 copy_list(perm, nu+(newgamma-1)*n, n);
828 /* Gamma backward pointer of newgamma */
829 *(w+newgamma-1) = gamma;
830 }
831 }
832 np++;
833 }
834 /* Free allocated memory */
835 free(perm);
836
837 }
838
839 /**********************************************************************/
840
841 /* schreier_vector. JMM, 22 June 2003
842 *
843 * All orbits of a given GeneratingSet GS with m n-permutations, with
844 * combined vectors of backward permutations and pointers. Note that we
845 * do not return the orbits information, though it can be reconstructed
846 * from w and nu.
847 */
848
schreier_vector(int point,int * GS,int m,int n,int * nu,int * w)849 void schreier_vector(int point, int *GS, int m, int n, int *nu, int *w)
850 {
851
852 int i; /* Point counter (from 1 to n) */
853 int *orbit= (int*)malloc(n*sizeof(int));
854 int *usedpoints= (int*)malloc(n*sizeof(int));
855 int j=0; /* Counter of used points */
856 int ol;
857
858 /* First orbit */
859 one_schreier_orbit(point, GS, m, n, orbit, &ol, nu, w, 1);
860 while(ol--) usedpoints[j++] = orbit[ol];
861 /* Other orbits */
862 for(i=1; i<=n; i++) {
863 if (!position(i, usedpoints, j)) {
864 one_schreier_orbit(i, GS, m, n, orbit, &ol, nu, w, 0);
865 while(ol--) usedpoints[j++] = orbit[ol];
866 }
867 }
868 /* Free allocated memory */
869 free(orbit);
870 free(usedpoints);
871
872 }
873
874 /**********************************************************************/
875
876 /* trace_schreier. JMM, 22 June 2003
877 *
878 * This function traces the Schreier vector (orbits, nu, w) with the
879 * point given, returning in perm a permutation wich moves the first
880 * point of the corresponding orbit to point.
881 */
882
trace_schreier(int point,int * nu,int * w,int * perm,int n)883 void trace_schreier(int point, int *nu, int *w, int *perm, int n)
884 {
885
886 int i;
887 int found=0;
888 int *newperm= (int*)malloc(n*sizeof(int));
889
890 for(i=0; i<n; i++) {
891 if (*(w+point-1) == 0) {
892 range(perm, n);
893 found=1;
894 break;
895 }
896 }
897 if (!found) {
898 trace_schreier(*(w+point-1), nu, w, newperm, n);
899 product(newperm, nu+(point-1)*n, perm, n);
900 }
901 free(newperm);
902
903 }
904
905 /**********************************************************************/
906
907 /* order_of_group. JMM, 27 June 2003
908 *
909 * This function calculates the order of the group generated by the SGS
910 * given by base (with bl points) and GS (with m n-permutations).
911 *
912 * JMM, 12 March 2006: changed type to long long int (8 bytes).
913 * This allows numbers up to ~ 10^19.
914 */
915
order_of_group(int * base,int bl,int * GS,int m,int n)916 long long int order_of_group(int *base, int bl, int *GS, int m, int n)
917 {
918
919 if (m==0) return(1);
920 else {
921 int *stab= (int*)malloc(m*n*sizeof(int)), sl;
922 int *orbit= (int*)malloc( n*sizeof(int)), ol;
923 one_orbit(base[0], GS, m, n, orbit, &ol);
924 stabilizer(base, 1, GS, m, n, stab, &sl);
925 long long int ret=ol* order_of_group(base+1,bl-1,stab,sl,n);
926 free(stab);
927 free(orbit);
928 return ret;
929 }
930 }
931
932 /**********************************************************************/
933
934 /* perm_member. JMM, 27 June 2003
935 *
936 * Checks whether permutation p belongs (returns 1) or not (returns 0)
937 * to the group generated by the SGS given by base (length bl) and GS
938 * (m n-permutations). If bl=0 or m=0 the result is isid(p, n) .
939 *
940 * JMM, 26 June 2007: modified behaviour of m==0. It gave ret=1 before.
941 */
942
perm_member(int * p,int * base,int bl,int * GS,int m,int n)943 int perm_member(int *p, int *base, int bl, int *GS, int m, int n)
944 {
945
946 if (bl==0 || m==0) return( isid(p, n) );
947 else {
948 int *pp= (int*)malloc( n*sizeof(int));
949 int *ip= (int*)malloc( n*sizeof(int));
950 int *orbit= (int*)malloc( n*sizeof(int)), ol;
951 int *w= (int*)malloc( n*sizeof(int));
952 int *nu= (int*)malloc(n*n*sizeof(int));
953 int *stab= (int*)malloc(m*n*sizeof(int)), sl;
954 int point, ret;
955
956 one_schreier_orbit(base[0], GS,m,n, orbit,&ol, nu,w, 1);
957 point = onpoints(base[0], p, n);
958 if (position(point, orbit, ol)) {
959 trace_schreier(point, nu, w, pp, n);
960 inverse(pp, ip, n);
961 product(p, ip, pp, n);
962 stabilizer(base, 1, GS, m, n, stab, &sl);
963 ret = perm_member(pp, base+1, bl-1, stab,sl,n);
964 }
965 else ret = 0;
966
967 free(pp);
968 free(ip);
969 free(orbit);
970 free(w);
971 free(nu);
972 free(stab);
973
974 return(ret);
975 }
976
977 }
978
979 /**********************************************************************/
980
981 /* schreier_sims. JMM, 28-30 June 2003
982 *
983 * This function computes a Strong Generating Set for the group of
984 * permutations generated by GS (with m n-permutations). It is
985 * possible to give the first bl (possibly 0) elements of the base,
986 * and the code computes all the others. There are possibly (and
987 * probably) redundant points in the final base. n=0 is not consistent
988 * here.
989 * We assume that enough space (at least, and typically, m*n integers)
990 * has been already allocated for newGS. We also assume that newGS can
991 * be reallocated. That's why we do not send a pointer, but a pointer
992 * to that pointer. That is, it needs to be reallocated in a different
993 * subroutine, and that cannot be done with a normal pointer!
994 */
995
schreier_sims(int * base,int bl,int * GS,int m,int n,int * newbase,int * nbl,int ** newGS,int * nm,int * num)996 void schreier_sims(int *base, int bl, int *GS, int m, int n,
997 int *newbase, int *nbl, int **newGS, int *nm, int *num)
998 {
999
1000 #ifdef VERBOSE_SCHREIER /*PPC*/
1001 printf("******** SCHREIER-SIMS ALGORITHM ********\n"); /*PPC*/
1002 #endif /*PPC*/
1003 /* Note the cycle: base -> newbase -> base2 -> newbase -> ... */
1004
1005 /* Copy base into newbase, adding more points if needed */
1006 nonstable_points(base, bl, GS, m, n, newbase, nbl);
1007 #ifdef VERBOSE_SCHREIER /*PPC*/
1008 printf("Original base:");
1009 print_list(base, bl, 1); /*PPC*/
1010 printf("New base:");
1011 print_list(newbase, *nbl, 1); /*PPC*/
1012 #endif /*PPC*/
1013 /* Initialize newGS=GS */
1014 copy_list(GS, *newGS, m*n);
1015 *nm = m;
1016 if (*nbl==0) { /* Problem. Return input sets */
1017 copy_list(base, newbase, bl);
1018 *nbl = bl;
1019 return;
1020 }
1021
1022 /* Allocate memory for intermediate SGS and stabilizer */
1023 int i; /* Base index */
1024 int *base2= (int*)malloc( n*sizeof(int)), bl2;/* Interm base */
1025 int *GS2= (int*)malloc(m*n*sizeof(int)), m2; /* Interm GS */
1026 int *stab= (int*)malloc(m*n*sizeof(int)), mm; /* Stabilizer */
1027
1028 /* Main loop */
1029 m2 = *nm; /* Initially GS2 will be just newGS */
1030 for(i=(*nbl); i>0; i--) {
1031 #ifdef VERBOSE_SCHREIER /*PPC*/
1032 printf("\nComputing SGS for H^(%d)\n", i-1); /*PPC*/
1033 #endif /*PPC*/
1034 if (*nm > m2) { /* Reallocate GS2 and stab */
1035 GS2 = (int*)realloc(GS2, (*nm)*n*sizeof(int));
1036 stab = (int*)realloc(stab, (*nm)*n*sizeof(int));
1037 }
1038 /* Copy newbase into base2 */
1039 copy_list(newbase, base2, *nbl);
1040 bl2=*nbl;
1041 copy_list(*newGS, GS2, (*nm)*n);
1042 m2=*nm;
1043 /* Compute newbase from base2 */
1044 stabilizer(base2, i-1, GS2, m2, n, stab, &mm);
1045 schreier_sims_step(base2, bl2, GS2, m2, n, i,
1046 stab, mm, newbase, nbl, newGS, nm, num);
1047 }
1048
1049 /* Free allocated memory */
1050 free(base2);
1051 free(GS2);
1052 free(stab);
1053
1054 #ifdef VERBOSE_SCHREIER /*PPC*/
1055 printf("************ END OF ALGORITHM ***********\n"); /*PPC*/
1056 #endif /*PPC*/
1057 }
1058
1059 /* Intermediate function for schreier_sims.
1060 *
1061 * Assuming that S^(i-1) is a GenSet for H^(i-1) and that [B, S^(i)]
1062 * are a StrongGenSet for H^(i), find a StrongGenSet for H^(i-1).
1063 * T is the set of additional generators in S^(i-1) since the previous
1064 * call to the procedure with the present value of i. Assume that a
1065 * StrongGenSet of < S^(i-1) - T > (the previous value of H^(i-1)) is
1066 * included in B and S. The present value of nu^(i-1) must be an
1067 * extension of the previous value.
1068 *
1069 * base (length bl): tentative base for the whole group
1070 * GS (m n-permutations): tentative StrongGenSet for the whole group
1071 * i: current order in the hierarchy of stabilizers
1072 *
1073 * Sizes on input:
1074 * i <= bl <= n Not changed
1075 * mm <= m Not changed
1076 * nbl = bl nbl will increase
1077 * nm = m nm will increase
1078 */
1079
schreier_sims_step(int * base,int bl,int * GS,int m,int n,int i,int * T,int mm,int * newbase,int * nbl,int ** newGS,int * nm,int * num)1080 void schreier_sims_step(int *base, int bl, int *GS, int m, int n,
1081 int i, int *T, int mm,
1082 int *newbase, int *nbl, int **newGS, int *nm, int *num)
1083 {
1084
1085 /* Declarations */
1086 /* Counters */
1087 int c, j=0, jj, level;
1088 /* Intermediate permutations */
1089 int *p= (int*)malloc(n*sizeof(int));
1090 int *ip= (int*)malloc(n*sizeof(int));
1091 int *pp= (int*)malloc(n*sizeof(int));
1092 int *ppp= (int*)malloc(n*sizeof(int));
1093 /* Stabilizer of base[1...i-1] */
1094 int *Si= (int*)malloc(m*n*sizeof(int)), Sil;
1095 /* Old stabilizer. Here we could use mm*n rather than m*n */
1096 int *oldSi= (int *)malloc(m*n*sizeof(int)), oldSil;
1097 /* Orbit of base[i] */
1098 int *Deltai= (int*)malloc( n*sizeof(int)), Deltail;
1099 int *w= (int*)malloc( n*sizeof(int));
1100 int *nu= (int*)malloc(n*n*sizeof(int));
1101 /* Old orbit */
1102 int *oldDeltai= (int*)malloc( n*sizeof(int)), oldDeltail;
1103 int *oldw= (int*)malloc( n*sizeof(int));
1104 int *oldnu= (int*)malloc(n*n*sizeof(int));
1105 /* Generators to check */
1106 int *genset= (int*)malloc(m*n*sizeof(int)), gensetl;
1107 /* Loops */
1108 int gamma, gn, sn;
1109 int *s= (int*)malloc(n*sizeof(int));
1110 int *g= (int*)malloc(n*sizeof(int));
1111 /* Stabilizer */
1112 int *stab = (int*)malloc(m*n*sizeof(int)), stabl;
1113 int *stabps= (int*)malloc( n*sizeof(int)), stabpsl;
1114
1115 #ifdef VERBOSE_SCHREIER /*PPC*/
1116 printf("******** schreier_sims_step ********\n"); /*PPC*/
1117 printf("base:");
1118 print_list(base, bl, 1); /*PPC*/
1119 printf("GS (%d perms of degree %d):", m, n); /*PPC*/
1120 print_array_perm(GS, m, n, 1); /*PPC*/
1121 #endif /*PPC*/
1122 /* Initialize newbase=base and newGS=GS (and their lengths) */
1123 /* They are already equal on input. Always true? */
1124 copy_list(base, newbase, bl);
1125 *nbl = bl;
1126 copy_list(GS, *newGS, m*n);
1127 *nm = m;
1128 /* Original generating sets. We get Sil<=m and oldSil<=mm */
1129 stabilizer(base, i-1, GS, m, n, Si, &Sil);
1130 #ifdef VERBOSE_SCHREIER /*PPC*/
1131 printf("Stabilizer of first %d points of base ", i-1); /*PPC*/
1132 print_list(base, i-1, 0);
1133 printf(" :\n"); /*PPC*/
1134 print_array_perm(Si, Sil, n, 1); /*PPC*/
1135 #endif /*PPC*/
1136 complement(Si, Sil, T, mm, n, oldSi, &oldSil);
1137 #ifdef VERBOSE_SCHREIER /*PPC*/
1138 printf("Previous stabilizer of %d points:\n", i-1); /*PPC*/
1139 print_array_perm(oldSi, oldSil, n, 1); /*PPC*/
1140 #endif /*PPC*/
1141 /* Basic orbits */
1142 one_schreier_orbit(base[i-1], Si, Sil, n,
1143 Deltai, &Deltail, nu, w, 1);
1144 one_schreier_orbit(base[i-1], oldSi, oldSil, n,
1145 oldDeltai, &oldDeltail, oldnu, oldw, 1);
1146 /* Check that Deltai is an extension of oldDeltai */
1147 for(c=0; c<n; c++) {
1148 if(w[c]!=oldw[c] && oldw[c]!=0) {
1149 #ifdef VERBOSE_SCHREIER /*PPC*/
1150 printf("Deltai[%d] modified to match oldDeltai\n", c); /*PPC*/
1151 #endif /*PPC*/
1152 copy_list(oldnu+c*n, nu+c*n, n);
1153 w[c] = oldw[c];
1154 }
1155 }
1156 #ifdef VERBOSE_SCHREIER /*PPC*/
1157 printf("Orbit: ");
1158 print_list(Deltai, Deltail, 1); /*PPC*/
1159 #endif /*PPC*/
1160
1161 /* Loop gn over elements gamma of basic orbit */
1162 for(gn=0; gn<Deltail; gn++) {
1163
1164 gamma = Deltai[gn];
1165
1166 #ifdef VERBOSE_SCHREIER /*PPC*/
1167 printf(" gamma=%d\n", gamma); /*PPC*/
1168 #endif /*PPC*/
1169 /* In both cases here we get gensetl<=m */
1170 if (position(gamma, oldDeltai, oldDeltail)) {
1171 copy_list(T, genset, mm*n);
1172 gensetl=mm;
1173 }
1174 else {
1175 copy_list(Si, genset, Sil*n);
1176 gensetl=Sil;
1177 }
1178
1179 /* Loop sn over generators s in genset */
1180 for (sn=0; sn<gensetl; sn++) {
1181
1182 copy_list(genset+sn*n, s, n);
1183 (*num)++;
1184
1185 /* Compute Schreier generator */
1186 trace_schreier(gamma, nu,w, p,n);
1187 trace_schreier(onpoints(gamma,s,n), nu,w, pp,n);
1188 inverse(pp, ip, n);
1189 product(p, s, ppp, n);
1190 product(ppp, ip, g, n);
1191 #ifdef VERBOSE_SCHREIER /*PPC*/
1192 printf(" g(%d)=", *num);
1193 print_perm(g, n, 1); /*PPC*/
1194 #endif /*PPC*/
1195
1196 /* Compute stabilizer. Reallocate to maximum size */
1197 stab = (int*)realloc(stab, (*nm)*n*sizeof(int));
1198 stabilizer(newbase, i, *newGS, *nm, n, stab, &stabl);
1199 /* If g is not in subgroup H^(i) */
1200 if(!isid(g, n)) {
1201 if ((stabl==0)||(!perm_member(g, newbase+i,*nbl-i, stab,stabl, n))) {
1202 #ifdef VERBOSE_SCHREIER /*PPC*/
1203 printf(" g not in H^(%d)\n", i); /*PPC*/
1204 #endif /*PPC*/
1205 /* Enlarge newGS. We need reallocation */
1206 *newGS = (int*)realloc(*newGS, ((*nm)+1)*n*sizeof(int));
1207 copy_list(g, (*newGS)+(*nm)*n, n);
1208 (*nm)++;
1209 /* Extend newbase if needed, so that no strong
1210 * generator fixes all base points. */
1211 stable_points(g, n, stabps, &stabpsl);
1212 #ifdef VERBOSE_SCHREIER /*PPC*/
1213 printf(" Stable under g: "); /*PPC*/
1214 print_list(stabps, stabpsl, 1); /*PPC*/
1215 #endif /*PPC*/
1216 /* If g moves a point of newbase then set j */
1217 for(jj=0; jj<*nbl; jj++) {
1218 if(!position(newbase[jj], stabps, stabpsl)) {
1219 j = jj+1;
1220 #ifdef VERBOSE_SCHREIER /*PPC*/
1221 printf(" g moves %d in newbase\n", newbase[jj]); /*PPC*/
1222 #endif /*PPC*/
1223 break;
1224 }
1225 /* else choose a new point */
1226 j = *nbl+1;
1227 }
1228 if (j == *nbl+1) {
1229 for(jj=1; jj<=n; jj++) {
1230 if(!position(jj, stabps, stabpsl) &&
1231 !position(jj, newbase, *nbl)) {
1232 newbase[*nbl]=jj;
1233 (*nbl)++;
1234 #ifdef VERBOSE_SCHREIER /*PPC*/
1235 printf(" Point %d added to newbase\n", jj); /*PPC*/
1236 #endif /*PPC*/
1237 break;
1238 }
1239 }
1240 }
1241 /* Ensure we still have a base and SGS for
1242 * H^(i+1) */
1243 for(level=j; level>i; level--) {
1244 #ifdef VERBOSE_SCHREIER /*PPC*/
1245 printf("\nEnsuring H^(%d) at level %d\n", i+1, level); /*PPC*/
1246 #endif /*PPC*/
1247 schreier_sims_step(newbase,*nbl,*newGS,*nm, n,
1248 level, g, 1,
1249 newbase,nbl,newGS,nm,num);
1250 }
1251 #ifdef VERBOSE_SCHREIER /*PPC*/
1252 printf("***** Finished check of H(%d) ******\n\n", i+1);/*PPC*/
1253 #endif /*PPC*/
1254 }
1255 }
1256 }
1257 }
1258
1259 /* Free allocated memory */
1260 free(p);
1261 free(ip);
1262 free(pp);
1263 free(ppp);
1264 free(Si);
1265 free(oldSi);
1266 free(Deltai);
1267 free(w);
1268 free(nu);
1269 free(oldDeltai);
1270 free(oldw);
1271 free(oldnu);
1272 free(genset);
1273 free(s);
1274 free(g);
1275 free(stab);
1276 free(stabps);
1277
1278 }
1279
1280 /**********************************************************************
1281 *
1282 * Notations and comments:
1283 * - In xTensor, a permutation g corresponding to a given tensor
1284 * configuration is understood as acting on index-numbers giving
1285 * tensor-slot numbers. That is, p = i^g (which is p= onpoints(i, g))
1286 * means that i is the index at slot p in the configuration g.
1287 * Equivalently, i = p^ig, where ig is the inverse of g.
1288 * - That convention is precisely the opposite to the one chosen by
1289 * Renato. Everything here and in xPerm follow Renato's convention.
1290 * There are two InversePerm actions in ToCanonicalOne in xTensor to
1291 * perform the change of notation. Here canonical_perm also has two
1292 * inverse actions for backwards compatibility.
1293 * - The base of a SGS for the group S represents an ordering of the
1294 * slots of the tensor. Changing the base can be understood as a
1295 * change of priorities for index canonicalization.
1296 */
1297
1298 /**********************************************************************/
1299
1300 /* coset_rep. JMM, 30 June 2003
1301 *
1302 * This algorithm is encoded from Renato Portugal et al.
1303 *
1304 * This function gives a canonical representant of a n-permutation p
1305 * in the group described by a SGS (pair base,GS) and the result is
1306 * stored in cr. In other words, it gives the canonical representant
1307 * of the right coset S.p. The free slots are different at return time.
1308 */
1309
coset_rep(int * p,int n,int * base,int bl,int * GS,int * m,int * freeps,int fl,int * cr)1310 void coset_rep(int *p, int n,
1311 int *base, int bl, int *GS, int *m,
1312 int *freeps, int fl,
1313 int *cr)
1314 {
1315
1316 #ifdef VERBOSE_COSET /*PPC*/
1317 printf("***** RIGHT-COSET-REP ALGORITHM ****\n"); /*PPC*/
1318 printf("Permutation ");
1319 print_perm(p, n, 1); /*PPC*/
1320 printf("Base: ");
1321 print_list(base, bl, 1); /*PPC*/
1322 #endif /*PPC*/
1323 /* Trivial case without symmetries */
1324 if (*m==0) {
1325 copy_list(p, cr, n);
1326 return;
1327 }
1328 /* else */
1329
1330 int i, j, k, b, pp;
1331 int *deltap= (int*)malloc( n*sizeof(int)), deltapl;
1332 int *deltapsorted= (int*)malloc( n*sizeof(int));
1333 int *om= (int*)malloc( n*sizeof(int));
1334 int *PERM= (int*)malloc( n*sizeof(int));
1335 int *perm2= (int*)malloc( n*sizeof(int));
1336 int *orbit= (int*)malloc( n*sizeof(int)), ol;
1337 int *orbit1= (int*)malloc( n*sizeof(int)), o1l;
1338 int *w= (int*)malloc( n*sizeof(int));
1339 int *nu= (int*)malloc(n*n*sizeof(int));
1340 int *genset= (int*)malloc(*m*n*sizeof(int)), gensetl;
1341 int *stab= (int*)malloc(*m*n*sizeof(int)), mm;
1342
1343 /* Copy p to PERM and GS to genset, to avoid side effects. */
1344 copy_list(p, PERM, n);
1345 copy_list(GS, genset, *m*n);
1346 gensetl=*m;
1347
1348 /* Loop over elements of base */
1349 for(i=0; i<bl; i++) {
1350 b = base[i]; /* Slot to analyze */
1351 one_schreier_orbit(b, genset, gensetl, n,
1352 orbit, &ol, nu, w, 1);
1353 #ifdef VERBOSE_COSET /*PPC*/
1354 printf("\nAnalyzing slot %d\n", b); /*PPC*/
1355 printf("orbit: ");
1356 print_list(orbit, ol, 1); /*PPC*/
1357 printf("freeps: ");
1358 print_list(freeps, fl, 1); /*PPC*/
1359 #endif /*PPC*/
1360 intersection(orbit, ol, freeps, fl, orbit1, &o1l);
1361 #ifdef VERBOSE_COSET /*PPC*/
1362 printf("Free slots that can go to that slot: "); /*PPC*/
1363 print_list(orbit1, o1l, 1); /*PPC*/
1364 #endif /*PPC*/
1365 if (o1l==0) continue; /* Slot with no symmetries */
1366 /* else */
1367 for(j=0; j<o1l; j++) {
1368 deltap[j] = onpoints(orbit1[j], PERM, n);
1369 }
1370 deltapl=o1l;
1371 #ifdef VERBOSE_COSET /*PPC*/
1372 printf("At those slots we resp. find indices deltap: ");/*PPC*/
1373 print_list(deltap, deltapl, 1); /*PPC*/
1374 #endif /*PPC*/
1375 sortB(deltap, deltapsorted, deltapl, base, bl);
1376 k = position(deltapsorted[0], deltap, deltapl);
1377 #ifdef VERBOSE_COSET /*PPC*/
1378 printf("Least index: %d, at position k: %d of deltap\n",/*PPC*/
1379 deltap[k-1], k); /*PPC*/
1380 #endif /*PPC*/
1381 pp = orbit1[k-1];
1382 #ifdef VERBOSE_COSET /*PPC*/
1383 printf("That index is at tensor slot pp: %d\n", pp); /*PPC*/
1384 #endif /*PPC*/
1385 /* Compute permutation om such that b^om = pp */
1386 trace_schreier(pp, nu, w, om, n);
1387 #ifdef VERBOSE_COSET /*PPC*/
1388 printf("We can move slot %d to slot %d with perm om:", /*PPC*/
1389 pp, b); /*PPC*/
1390 print_perm(om, n, 0);
1391 printf(" in S\n"); /*PPC*/
1392 #endif /*PPC*/
1393 product(om, PERM, perm2, n);
1394 #ifdef VERBOSE_COSET /*PPC*/
1395 printf("New list of indices: "); /*PPC*/
1396 print_perm(perm2, n, 1); /*PPC*/
1397 #endif /*PPC*/
1398 copy_list(perm2, PERM, n);
1399 /* New positions of free indices */
1400 inverse(om, perm2, n);
1401 for(j=0; j<fl; j++) {
1402 freeps[j] = onpoints(freeps[j], perm2, n);
1403 }
1404 #ifdef VERBOSE_COSET /*PPC*/
1405 printf("Removing those perms that move slot %d\n", b); /*PPC*/
1406 #endif /*PPC*/
1407 /* Note that we do not change base to have i as first
1408 member of base. This is not general, but I think
1409 here it is valid due to the nature of our problem */
1410 stabilizer(base+i, 1, genset, gensetl, n, stab, &mm);
1411 copy_list(stab, genset, mm*n);
1412 gensetl=mm;
1413 }
1414 /* Move to result */
1415 copy_list(PERM, cr, n);
1416 /* Return new SGS */
1417 copy_list(genset, GS, gensetl*n);
1418 *m=gensetl;
1419 #ifdef VERBOSE_COSET /*PPC*/
1420 printf("************ END OF ALGORITHM ***********\n"); /*PPC*/
1421 #endif /*PPC*/
1422
1423 /* Free allocated memory */
1424 free(deltap);
1425 free(deltapsorted);
1426 free(om);
1427 free(PERM);
1428 free(perm2);
1429 free(orbit);
1430 free(orbit1);
1431 free(w);
1432 free(nu);
1433 free(genset);
1434 free(stab);
1435
1436 }
1437
1438 /**********************************************************************/
1439
1440 /* SGSD. JMM, 26 June 2007
1441 *
1442 * We construct a number of functions that generate the strong
1443 * generating set for a given list of (multiple) dummysets (all pairs
1444 * of dummies of a vbundle) and repeatedsets (the positions of repeated
1445 * indices, like components or directions). The main function is SGSD,
1446 * but there are other elementary operations: SGSofdummyset,
1447 * SGSofrepeatedset, movedummyset, moverepeatedset.
1448 * We introduce another funny name here: drummy: either a dummy or a
1449 * repeated index.
1450 */
1451
1452 /* SGS for a dummyset. List dummies not modified */
SGSofdummyset(int * dummies,int dl,int sym,int n,int * KD,int * KDl,int * bD,int * bDl)1453 void SGSofdummyset(int *dummies, int dl, int sym, int n,
1454 int *KD, int *KDl, int *bD, int *bDl)
1455 {
1456
1457 if (dl==0) {
1458 *KDl=0;
1459 *bDl=0;
1460 return;
1461 } /* else */
1462
1463 /* Number of pairs of dummies: dpl */
1464 int dpl = dl/2;
1465 int *range_perm = (int*)malloc( n*sizeof(int));
1466 int *KD1 = (int*)malloc(dpl*n*sizeof(int));
1467 int *KD2 = (int*)malloc(dpl*n*sizeof(int));
1468 int i;
1469
1470 range(range_perm, n);
1471 /* KD1: exchange indices. Ex: Cycles[{1,3},{2,4}]
1472 There are always dpl-1 permutations */
1473 for (i=0; i<dpl-1; i++) {
1474 /* Copy the identity */
1475 copy_list(range_perm, KD1+i*n, n);
1476 /* Swap elements of consecutive pairs */
1477 KD1[ i*n+dummies[2*i ]-1 ] = dummies[2*i+2];
1478 KD1[ i*n+dummies[2*i+2]-1 ] = dummies[2*i ];
1479 KD1[ i*n+dummies[2*i+1]-1 ] = dummies[2*i+3];
1480 KD1[ i*n+dummies[2*i+3]-1 ] = dummies[2*i+1];
1481 }
1482 /* KD2: exchange indices of the same pair.
1483 If sym!=0 there are dpl permutations */
1484 if (sym == 1) { /* Symmetric metric */
1485 for (i=0; i<dpl; i++) {
1486 /* Copy the identity */
1487 copy_list(range_perm, KD2+i*n, n);
1488 /* Swap elements of pair */
1489 KD2[i*n+dummies[2*i]-1] = dummies[2*i+1];
1490 KD2[i*n+dummies[2*i+1]-1] = dummies[2*i];
1491 /* Keep positive sign */
1492 }
1493 *KDl = 2*dpl-1;
1494 }
1495 else if (sym == -1) {
1496 for (i=0; i<dpl; i++) {
1497 /* Copy the identity */
1498 copy_list(range_perm, KD2+i*n, n);
1499 /* Swap elements of pair */
1500 KD2[i*n+dummies[2*i]-1] = dummies[2*i+1];
1501 KD2[i*n+dummies[2*i+1]-1] = dummies[2*i];
1502 /* Set negative sign */
1503 KD2[i*n+n-2] = n;
1504 KD2[i*n+n-1] = n-1;
1505 }
1506 *KDl = 2*dpl-1;
1507 }
1508 else if (sym == 0) {
1509 /* Do nothing */
1510 *KDl = dpl-1;
1511 }
1512 else {
1513 /* Unknown value of sym */
1514 }
1515 /* KD */
1516 copy_list(KD1, KD, (dpl-1)*n);
1517 if (sym!=0) copy_list(KD2, KD+(dpl-1)*n, dpl*n);
1518 /* base of group D */
1519 for (i=0; i<dpl; i++) {
1520 bD[i] = dummies[2*i];
1521 }
1522 *bDl=dpl;
1523 #ifdef VERBOSE_DOUBLE /*PPC*/
1524 printf("KD: ");
1525 print_array_perm(KD, *KDl, n, 1); /*PPC*/
1526 printf("bD: ");
1527 print_list(bD, *bDl, 1); /*PPC*/
1528 #endif /*PPC*/
1529 free(range_perm);
1530 free(KD1);
1531 free(KD2);
1532 }
1533
1534 /* SGS for a repeatedset. List repes not modified */
SGSofrepeatedset(int * repes,int rl,int n,int * KD,int * KDl,int * bD,int * bDl)1535 void SGSofrepeatedset(int *repes, int rl, int n,
1536 int *KD, int *KDl, int *bD, int *bDl)
1537 {
1538
1539 #ifdef VERBOSE_DOUBLE /*PPC*/
1540 printf("From SGSofrepeatedset:\n"); /*PPC*/
1541 printf("repes: ");
1542 print_list(repes, rl, 1);
1543 #endif /*PPC*/
1544 if (rl==0) {
1545 *KDl=0;
1546 *bDl=0;
1547 return;
1548 } /* else */
1549
1550 int *range_perm = (int*)malloc( n*sizeof(int));
1551 int i;
1552
1553 range(range_perm, n);
1554 for (i=0; i<rl-1; i++) {
1555 /* Copy the identity */
1556 copy_list(range_perm, KD+i*n, n);
1557 /* Swap elements of pair */
1558 KD[i*n+repes[i]-1] = repes[i+1];
1559 KD[i*n+repes[i+1]-1] = repes[i];
1560 /* Keep positive sign */
1561 }
1562 *KDl = rl-1;
1563 copy_list(repes, bD, rl-1);
1564 *bDl = rl-1;
1565 #ifdef VERBOSE_DOUBLE /*PPC*/
1566 printf("KD: ");
1567 print_array_perm(KD, *KDl, n, 1); /*PPC*/
1568 printf("bD: ");
1569 print_list(bD, *bDl, 1); /*PPC*/
1570 #endif /*PPC*/
1571 free(range_perm);
1572 }
1573
1574 /* Move index in a dummyset. List dummies reordered */
movedummyset(int firstd,int * dummies,int dl,int)1575 void movedummyset(int firstd, int *dummies, int dl, int)
1576 {
1577
1578 /* Find position of dummy and relative
1579 position of its pair */
1580 #ifdef VERBOSE_DOUBLE /*PPC*/
1581 printf("Rearrange dummies for %d. dummies: ", firstd); /*PPC*/
1582 print_list(dummies, dl, 1); /*PPC*/
1583 #endif /*PPC*/
1584 int pos, j;
1585 pos = position(firstd, dummies, dl)-1;
1586 if (pos==-1) { /* firstd not in dummies */
1587 /* Do nothing */
1588 }
1589 else {
1590 int tmp;
1591 /* Swap all pairs if firstd found as second */
1592 if (pos%2==0) {
1593 /* 1st member: Do nothing */
1594 }
1595 else {
1596 pos=pos-1;
1597 /* 2nd member: Swap all pairs */
1598 for (j=0; j<dl/2; j++) {
1599 tmp= dummies[2*j];
1600 dummies[2*j]=dummies[2*j+1];
1601 dummies[2*j+1]=tmp;
1602 }
1603 }
1604 /* Exchange position of pair with first pair */
1605 if (pos==0) {
1606 /* 1st pair: Do nothing */
1607 }
1608 else {
1609 /* Exchange two pairs of dummies */
1610 tmp=dummies[0];
1611 dummies[0]=firstd;
1612 dummies[pos]=tmp;
1613 tmp=dummies[1];
1614 dummies[1]=dummies[pos+1];
1615 dummies[pos+1]=tmp;
1616 }
1617 }
1618 #ifdef VERBOSE_DOUBLE /*PPC*/
1619 printf("Now dummies: "); /*PPC*/
1620 print_list(dummies, dl, 1); /*PPC*/
1621 #endif /*PPC*/
1622 }
1623
1624 /* Move index in a repeated set. List repes reordered */
moverepeatedset(int firstd,int * repes,int rl)1625 void moverepeatedset(int firstd, int *repes, int rl)
1626 {
1627
1628 /* Find position of dummy and relative
1629 position of its pair */
1630 #ifdef VERBOSE_DOUBLE /*PPC*/
1631 printf("Rearrange dummies for %d. repes: ", firstd); /*PPC*/
1632 print_list(repes, rl, 1); /*PPC*/
1633 #endif /*PPC*/
1634 int pos;
1635 pos = position(firstd, repes, rl)-1;
1636 if (pos==-1) { /* firstd not in repes */
1637 /* Do nothing */
1638 }
1639 else {
1640 if (pos==0) {
1641 /* 1st index: Do nothing */
1642 }
1643 else {
1644 /* Exchange two positions */
1645 repes[pos] = repes[0];
1646 repes[0] = firstd;
1647 }
1648 }
1649 #ifdef VERBOSE_DOUBLE /*PPC*/
1650 printf("Now repes: "); /*PPC*/
1651 print_list(repes, rl, 1); /*PPC*/
1652 #endif /*PPC*/
1653 }
1654
1655 /* Remove a first-pair from dummies */
dropdummyset(int firstd,int * vds,int vdsl,int * dummies,int * dl)1656 void dropdummyset(int firstd,
1657 int *vds, int vdsl, int *dummies, int *dl)
1658 {
1659
1660 int i, j, itotal=0;
1661
1662 for (i=0; i<vdsl; i++) {
1663 if (dummies[itotal]==firstd && vds[i]!=0) {
1664 for (j=itotal; j<*dl-2; j++) {
1665 dummies[j] = dummies[j+2];
1666 }
1667 vds[i] = vds[i]-2;
1668 *dl = *dl - 2;
1669 return;
1670 }
1671 itotal = itotal + vds[i];
1672 }
1673 }
1674
1675 /* Remove a first-point from repes */
droprepeatedset(int firstd,int * vrs,int vrsl,int * repes,int * rl)1676 void droprepeatedset(int firstd,
1677 int *vrs, int vrsl, int *repes, int *rl)
1678 {
1679
1680 int i, j, itotal=0;
1681
1682 for (i=0; i<vrsl; i++) {
1683 if (repes[itotal]==firstd && vrs[i]!=0) {
1684 for (j=itotal; j<*rl; j++) {
1685 repes[j] = repes[j+1];
1686 }
1687 vrs[i] = vrs[i]-1;
1688 *rl = *rl - 1;
1689 return;
1690 }
1691 itotal = itotal + vrs[i];
1692 }
1693 }
1694
1695 /* SGSD for a given list of dummies and repes */
SGSD(int * vds,int vdsl,int * dummies,int dl,int * mQ,int * vrs,int vrsl,int * repes,int rl,int n,int firstd,int * KD,int * KDl,int * bD,int * bDl)1696 void SGSD(int *vds, int vdsl, int *dummies, int dl, int *mQ,
1697 int *vrs, int vrsl, int *repes, int rl, int n,
1698 int firstd, int *KD, int *KDl, int *bD, int *bDl)
1699 {
1700
1701 if (dl==0 && rl==0) {
1702 *KDl=0;
1703 *bDl=0;
1704 return;
1705 } /* else */
1706
1707 int *tmp, tmpl;
1708 int *tmpGS = (int*)malloc(n*n*sizeof(int)), tmpGSl;
1709 int *tmpbase = (int*)malloc( n*sizeof(int)), tmpbasel;
1710 int i, itotal;
1711
1712 /* Loop over all dummysets */
1713 itotal = 0;
1714 *KDl = 0;
1715 *bDl = 0;
1716 for (i=0; i<vdsl; i++) {
1717 tmpl = vds[i];
1718 tmp = dummies + itotal;
1719 movedummyset(firstd, tmp, tmpl, mQ[i]);
1720 itotal = itotal + tmpl;
1721 SGSofdummyset(tmp, tmpl, mQ[i], n,
1722 tmpGS, &tmpGSl, tmpbase, &tmpbasel);
1723 copy_list(tmpGS, KD+(*KDl)*n, tmpGSl*n);
1724 *KDl = *KDl + tmpGSl;
1725 copy_list(tmpbase, bD+ *bDl, tmpbasel);
1726 *bDl = *bDl + tmpbasel;
1727 }
1728
1729 /* Loop over all repeatedsets */
1730 itotal=0;
1731 for (i=0; i<vrsl; i++) {
1732 tmpl = vrs[i];
1733 tmp = repes + itotal;
1734 moverepeatedset(firstd, tmp, tmpl);
1735 itotal = itotal + tmpl;
1736 SGSofrepeatedset(tmp, tmpl, n,
1737 tmpGS, &tmpGSl, tmpbase, &tmpbasel);
1738 copy_list(tmpGS, KD+(*KDl)*n, tmpGSl*n);
1739 *KDl = *KDl + tmpGSl;
1740 copy_list(tmpbase, bD+ *bDl, tmpbasel);
1741 *bDl = *bDl + tmpbasel;
1742 }
1743
1744 free(tmpGS);
1745 free(tmpbase);
1746
1747 #ifdef VERBOSE_DOUBLE
1748 printf("base of D:");
1749 print_list(bD, *bDl, 1); /*PPC*/
1750 printf("GS of D (%d perms of degree %d):", *KDl, n); /*PPC*/
1751 print_array_perm(KD, *KDl, n, 1); /*PPC*/
1752 #endif /*PPC*/
1753 }
1754
1755 /**********************************************************************/
1756
1757 /* double_coset_rep. JMM, 1-5 July 2003.
1758 * JMM, 9 November 2005. Corrected treatment of SGS of group D.
1759 * JMM, 26 June 2007. Large extension to consider several dummysets and
1760 * several repeatedsets.
1761 *
1762 * This algorithm is encoded from Renato Portugal et al, with the
1763 * extensions to repeatedsets.
1764 *
1765 * This function gives a canonical representant for the n-permutation g
1766 * in the double coset S.g.D given by the groups S, described by a SGS
1767 * (pair base, GS) and the group D, described by the direct product of
1768 * the pair symmetric dpl pairs of dummies and the S_k groups of k
1769 * repeated indices. Each dummyset has a metric_symmetry sign: if mQ=1
1770 * the metric is symmetric. If mQ=-1 the metric is antisymmetric
1771 * (spinor calculus). If mQ=0 there is no metric.
1772 * Note that n = dl + dr. The result is stored in dcr.
1773 */
1774
1775 class alphastruct {
1776 public:
1777 alphastruct();
1778 alphastruct(const alphastruct&);
1779 ~alphastruct();
1780
1781 void init(int n);
1782
1783 int *L; /* We assume that elements of L cannot
1784 be repeated. I only have experimental
1785 evidence for this */
1786 int Ll;
1787 int *s;
1788 int *d;
1789 int *o;
1790
1791 int nn;
1792 };
1793
alphastruct()1794 alphastruct::alphastruct()
1795 : L(0), Ll(0), s(0), d(0), o(0), nn(0)
1796 {
1797 }
1798
alphastruct(const alphastruct & other)1799 alphastruct::alphastruct(const alphastruct& other)
1800 {
1801 init(other.nn);
1802 for(int i=0; i<nn; ++i) {
1803 L[i]=other.L[i];
1804 s[i]=other.s[i];
1805 d[i]=other.d[i];
1806 o[i]=other.o[i];
1807 }
1808 Ll=other.Ll;
1809 }
1810
1811
init(int n)1812 void alphastruct::init(int n)
1813 {
1814 L = new int[n];
1815 s = new int[n];
1816 d = new int[n];
1817 o = new int[n];
1818 nn=n;
1819 }
1820
~alphastruct()1821 alphastruct::~alphastruct()
1822 {
1823 if(L) delete [] L;
1824 if(s) delete [] s;
1825 if(d) delete [] d;
1826 if(o) delete [] o;
1827 nn=0;
1828 }
1829
1830
TAB(std::vector<alphastruct> & ALPHA,int * L,int Ll,int * s1,int * d1,int n)1831 void TAB(std::vector<alphastruct>& ALPHA, int *L, int Ll, int *s1, int *d1, int n)
1832 {
1833 int i;
1834 int l=0;
1835
1836 /* Search ALPHA for the l corresponding to L */
1837 // std::cout << ALPHA.size() << ", " << Ll << std::endl;
1838 for (i=0; i<Ll; i++) l=ALPHA[l].o[L[i]-1];
1839 /* Copy permutations of element l */
1840 copy_list(ALPHA[l].s, s1, n);
1841 copy_list(ALPHA[l].d, d1, n);
1842 }
1843
F1(std::vector<alphastruct> & ALPHA,int * L,int Ll,int * g,int * list,int * listl,int n,int Deltabl,int * Deltab,int * DeltaD)1844 void F1(std::vector<alphastruct>& ALPHA, int *L, int Ll, int *g, int *list, int *listl, int n, int Deltabl, int *Deltab, int *DeltaD)
1845 {
1846 int c, c1, c2;
1847 int *sgd= (int*)malloc(n*sizeof(int));
1848 int *TAB1= (int*)malloc(n*sizeof(int));
1849 int *TAB2= (int*)malloc(n*sizeof(int));
1850 int *tmp= (int*)malloc(n*sizeof(int));
1851
1852 TAB(ALPHA, L, Ll, TAB1, TAB2, n);
1853
1854 /* Compute s.g.d */
1855 F2(TAB1, g, TAB2, sgd, n);
1856 #ifdef VERBOSE_DOUBLE /*PPC*/
1857 printf("With L=");
1858 print_list(L, Ll, 0); /*PPC*/
1859 printf(" we get sgd: ");
1860 print_perm(sgd, n, 1); /*PPC*/
1861 #endif /*PPC*/
1862 /* Images of Deltab under sgd. Note that tmp has length
1863 Deltabl */
1864 for (c=0; c<Deltabl; c++)
1865 tmp[c] = onpoints(Deltab[c], sgd, n);
1866 #ifdef VERBOSE_DOUBLE /*PPC*/
1867 printf(" which maps slots in Deltab to indices list: ");/*PPC*/
1868 print_list(tmp, Deltabl, 1); /*PPC*/
1869 #endif /*PPC*/
1870 /* Orbits of DeltaD containing the points in tmp */
1871 int oi;
1872 for (c1=0; c1<Deltabl; c1++) {
1873 oi = DeltaD[tmp[c1]-1];
1874 if(oi) {
1875 for(c2=0; c2<n; c2++) {
1876 if ((DeltaD[c2]==oi) &&
1877 (!position(c2+1, list, *listl)))
1878 list[(*listl)++] = c2+1;
1879 }
1880 }
1881 }
1882 #ifdef VERBOSE_DOUBLE /*PPC*/
1883 printf(" whose points belong to orbits "); /*PPC*/
1884 print_list(list, *listl, 1); /*PPC*/
1885 #endif /*PPC*/
1886 free(sgd);
1887 free(TAB1);
1888 free(TAB2);
1889 free(tmp);
1890 }
1891
1892
1893 /* Consistency check */
consistency(int * array,int m,int n)1894 int consistency(int *array, int m, int n)
1895 {
1896 int *arrayp= (int*)malloc(m*n*sizeof(int)), arraypl;
1897 int *arrayn= (int*)malloc(m*n*sizeof(int)), arraynl;
1898 int i, ip, in, ret;
1899
1900 #ifdef VERBOSE_DOUBLE /*PPC*/
1901 printf("Checking consistency in m:%d, n:%d\n", m, n); /*PPC*/
1902 print_array_perm(array, m, n, 1); /*PPC*/
1903 #endif /*PPC*/
1904
1905 /* Detect sign of permutation */
1906 arraypl=0;
1907 arraynl=0;
1908 for(i=0; i<m; i++) {
1909 if (array[i*n+n-2]<array[i*n+n-1]) /* Positive */
1910 copy_list(array+i*n, arrayp+(arraypl++)*n, n);
1911 else /* Negative */
1912 copy_list(array+i*n, arrayn+(arraynl++)*n, n);
1913 }
1914 #ifdef VERBOSE_DOUBLE /*PPC*/
1915 printf("Found positive perms: %d\n", arraypl); /*PPC*/
1916 printf("Found negative perms: %d\n", arraynl); /*PPC*/
1917 #endif /*PPC*/
1918 /* Here there are arraynl*arraypl comparisons. This
1919 should be improved with a better intersection
1920 algorithm which sorts the lists in advance */
1921 ret = 1; /* True */
1922 for (in=0; in<arraynl; in++) {
1923 for (ip=0; ip<arraypl; ip++) {
1924 if (equal_list(arrayp+ip*n, arrayn+in*n, n-2)) {
1925 ret = 0; /* False */
1926 break;
1927 }
1928 }
1929 }
1930 #ifdef VERBOSE_DOUBLE /*PPC*/
1931 if (ret) printf("Found no problem in check\n"); /*PPC*/
1932 else printf("Found perm with two signs.\n"); /*PPC*/
1933 #endif /*PPC*/
1934 free(arrayp);
1935 free(arrayn);
1936 return(ret);
1937 }
1938
1939
double_coset_rep(int * g,int n,int * base,int bl,int * GS,int m,int * vds,int vdsl,int * dummies,int dl,int * mQ,int * vrs,int vrsl,int * repes,int rl,int * dcr)1940 void double_coset_rep(int *g, int n, int *base, int bl, int *GS, int m,
1941 int *vds, int vdsl, int *dummies, int dl, int *mQ,
1942 int *vrs, int vrsl, int *repes, int rl, int *dcr)
1943 {
1944
1945 int i, j, l, jj, kk, c; /* Counters */
1946 int result;
1947 /* Inverse of permutation g */
1948 int *ig= (int*)malloc( n*sizeof(int));
1949 /* All drummies, both the pair-dummies and the repes */
1950 int *drummies= (int*)malloc( n*sizeof(int)), dril;
1951 /* The initial slots of all those dummies */
1952 int *drummyslots= (int*)malloc( n*sizeof(int));
1953 /* Temporary space for sorting */
1954 int *drummytmp= (int*)malloc( n*sizeof(int)), drummytmpl;
1955 int *drummytmp2= (int*)malloc( n*sizeof(int));
1956 /* Bases for group S */
1957 int *bS= (int*)malloc( n*sizeof(int)), bSl;
1958 int *bSsort= (int*)malloc( n*sizeof(int));
1959 /* gs for group S. It cannot be larger than GS */
1960 int *KS= (int*)malloc(m*n*sizeof(int)), KSl;
1961 /* gs for group D. The number dl+dr is an upper bound */
1962 int *KD= (int*)malloc(n*n*sizeof(int)), KDl;
1963 int *bD= (int*)malloc( n*sizeof(int)), bDl;
1964 int *nu= (int*)malloc(n*n*sizeof(int));
1965 int *w= (int*)malloc( n*sizeof(int));
1966 int *Deltab= (int*)malloc( n*sizeof(int)), Deltabl;
1967 int *DeltaD= (int*)malloc( n*sizeof(int));
1968 int *IMAGES= (int*)malloc( n*sizeof(int)), IMAGESl;
1969 int *IMAGESsorted= (int*)malloc( n*sizeof(int));
1970 int *p= (int*)malloc( n*sizeof(int));
1971 int *nuD= (int*)malloc(n*n*sizeof(int));
1972 int *wD= (int*)malloc( n*sizeof(int));
1973 int *Deltap= (int*)malloc( n*sizeof(int)), Deltapl;
1974 int *NEXT= (int*)malloc( n*sizeof(int)), NEXTl;
1975 int *L= (int*)malloc( n*sizeof(int)), Ll;
1976 int *L1= (int*)malloc( n*sizeof(int)), L1l;
1977 int *s= (int*)malloc( n*sizeof(int));
1978 int *d= (int*)malloc( n*sizeof(int));
1979 int *list1= (int*)malloc( n*sizeof(int)), list1l;
1980 int *list2= (int*)malloc( n*sizeof(int)), list2l;
1981 int *perm1= (int*)malloc( n*sizeof(int));
1982 int *perm2= (int*)malloc( n*sizeof(int));
1983 int *perm3= (int*)malloc( n*sizeof(int));
1984 int *s1= (int*)malloc( n*sizeof(int));
1985 int *d1= (int*)malloc( n*sizeof(int));
1986
1987
1988 /* We use Renato's notation, with g mapping slots to indices */
1989
1990 /**************************************************************
1991 * DEFINITIONS
1992 **************************************************************/
1993
1994 /* Given a list L={i1, ..., il} the function TAB must return
1995 a pair (s1, d1) corresponding to L. The whole collection of
1996 lists L and their pairs is stored in the array ALPHA of
1997 structures alphastruct. Each L in ALPHA is actually
1998 identified by its position l in the array. The structure
1999 for a given L contains a list of values of l corresponding
2000 to the n possible extensions of L.
2001 */
2002
2003 /* Define ALPHA and TAB */
2004
2005 int ALPHAl;
2006 int *ALPHAstep= (int*)malloc(n*sizeof(int));
2007
2008 /* Initialize ALPHA to {} and TAB to {id, id} */
2009 ALPHAl= 1;
2010 std::vector<alphastruct> ALPHA(1);
2011 // alphastruct *ALPHA= new alphastruct[1]; //(struct alphastruct*)malloc(sizeof(struct alphastruct));
2012 ALPHA[0].init(n);
2013 ALPHA[0].Ll=0;
2014 Ll=0; // ?
2015 range(ALPHA[0].s, n);
2016 range(ALPHA[0].d, n);
2017 ALPHAstep[0]=0;
2018 ALPHAstep[1]=1;
2019
2020
2021 /**************************************************************
2022 * CONSTRUCTION OF BASES
2023 **************************************************************/
2024
2025 /* Join all drummies */
2026 copy_list(dummies, drummies, dl);
2027 copy_list(repes, drummies+dl, rl);
2028 dril = dl + rl;
2029
2030 /* Slots of all those drummies */
2031 inverse(g, ig, n);
2032 for (i=0; i<dril; i++) {
2033 drummyslots[i] = onpoints(drummies[i], ig, n);
2034 }
2035
2036 /* Initialize KS */
2037 copy_list(GS, KS, m*n);
2038 KSl=m;
2039
2040 /* Extend base to bS to cover all positions of drummies.
2041 * We assume that in the intersection we get bS=base really.
2042 * We sort the missing drummies, but not the given base */
2043 intersection(base, bl, drummyslots, dril, bS, &bSl);
2044 complement(drummyslots,dril, base,bl, 1, drummytmp,&drummytmpl);
2045 sort(drummytmp, drummytmp2, drummytmpl);
2046 copy_list(drummytmp2, bS+bSl, drummytmpl);
2047 bSl = bSl + drummytmpl;
2048 #ifdef VERBOSE_DOUBLE /*PPC*/
2049 printf("All drummies: drummies: "); /*PPC*/
2050 print_list(drummies, dril, 1); /*PPC*/
2051 printf("Their slots: drummyslots: "); /*PPC*/
2052 print_list(drummyslots, dril, 1); /*PPC*/
2053 printf("base: "); /*PPC*/
2054 print_list(base, bl, 1); /*PPC*/
2055 printf("Complement: drummytmp: "); /*PPC*/
2056 print_list(drummytmp, drummytmpl, 1); /*PPC*/
2057 printf("Sort them: drummytmp2: "); /*PPC*/
2058 print_list(drummytmp2, drummytmpl, 1); /*PPC*/
2059 printf("base extended to bS: "); /*PPC*/
2060 print_list(bS, bSl, 1); /*PPC*/
2061 #endif /*PPC*/
2062 /* Generate associated base for sorting names of dummies */
2063 /* We choose a particular form of bSsort for aesthetics */
2064 sort(drummies,drummytmp2,dril);
2065 sort(bS, drummytmp, bSl);
2066 for (i=0; i<bSl; i++) {
2067 bSsort[i]=drummytmp2[position(bS[i],drummytmp,bSl)-1];
2068 }
2069
2070 #ifdef VERBOSE_DOUBLE /*PPC*/
2071 printf("drummies: "); /*PPC*/
2072 print_list(drummies, dril, 1); /*PPC*/
2073 printf("Base of SGSS: bS: "); /*PPC*/
2074 print_list(bS, bSl, 1); /*PPC*/
2075 printf("Sorted slots: "); /*PPC*/
2076 print_list(drummytmp2, bSl, 1); /*PPC*/
2077 printf("Sorted base bS: "); /*PPC*/
2078 print_list(drummytmp, bSl, 1); /*PPC*/
2079 printf("Base for sorting: bSsort: "); /*PPC*/
2080 print_list(bSsort, bSl, 1); /*PPC*/
2081 #endif /*PPC*/
2082
2083
2084 /*******************/
2085 /**** Main loop ****/
2086 /*******************/
2087
2088 /* Note we use i=1..bSl instead of the usual i=0 ..(bSl-1) */
2089 for (i=1; i<=bSl; i++) {
2090 int b=bS[i-1];
2091
2092 #ifdef VERBOSE_DOUBLE /*PPC*/
2093 printf("\n************** Loop i=%d *************\n", i);/*PPC*/
2094 printf("Analyzing slot bS[%d]=%d of tensor\n", i-1, b); /*PPC*/
2095 #endif /*PPC*/
2096 /* Schreier vector of S */
2097 schreier_vector(b, KS, KSl, n, nu, w);
2098 one_orbit(b, KS, KSl, n, Deltab, &Deltabl);
2099 #ifdef VERBOSE_DOUBLE /*PPC*/
2100 printf("Under S, slot %d go to slots Deltab: ",b); /*PPC*/
2101 print_list(Deltab, Deltabl, 1); /*PPC*/
2102 #endif /*PPC*/
2103 /* Compute SGS for group D. Do not rearrange drummies */
2104 SGSD(vds, vdsl, dummies, dl, mQ,
2105 vrs, vrsl, repes, rl, n,
2106 0, KD, &KDl, bD, &bDl);
2107 /* Orbits of D */
2108 all_orbits(KD, KDl, n, DeltaD);
2109 #ifdef VERBOSE_DOUBLE /*PPC*/
2110 printf("Orbits of indices: DeltaD: "); /*PPC*/
2111 print_list(DeltaD, n, 1); /*PPC*/
2112 #endif /*PPC*/
2113 /* Images of b under elements of S.g.D.
2114 Deltab and DeltaD are used by F1 */
2115 IMAGESl=0;
2116 for (c=ALPHAstep[i-1]; c<ALPHAstep[i]; c++)
2117 F1(ALPHA, ALPHA[c].L, ALPHA[c].Ll, g, IMAGES, &IMAGESl, n, Deltabl, Deltab, DeltaD);
2118 #ifdef VERBOSE_DOUBLE /*PPC*/
2119 printf("At slot %d we can have indices IMAGES: ", b); /*PPC*/
2120 print_list(IMAGES, IMAGESl, 1); /*PPC*/
2121 printf("IMAGESl: %d\n", IMAGESl); /*PPC*/
2122 #endif /*PPC*/
2123 /* If there are no images we have finished */
2124 if(IMAGESl==0) continue;
2125
2126 /* Find minimum index */
2127 sortB(IMAGES, IMAGESsorted, IMAGESl, bSsort, bSl);
2128 p[i-1] = IMAGESsorted[0];
2129 #ifdef VERBOSE_DOUBLE /*PPC*/
2130 printf("The least of them is p[i-1]=%d\n", p[i-1]); /*PPC*/
2131 #endif /*PPC*/
2132 /* Recompute SGS of D */
2133 if (dl>0 || rl>0) {
2134 SGSD(vds, vdsl, dummies, dl, mQ,
2135 vrs, vrsl, repes, rl, n,
2136 p[i-1], KD, &KDl, bD, &bDl);
2137 }
2138 else { /* Do nothing */
2139 #ifdef VERBOSE_DOUBLE /*PPC*/
2140 printf("Rearrangement of base of D not required.\n"); /*PPC*/
2141 #endif /*PPC*/
2142 }
2143 /* Schreier vector of D */
2144 schreier_vector(p[i-1], KD, KDl, n, nuD, wD);
2145 /* Orbit of p[i-1] under D */
2146 one_orbit(p[i-1], KD, KDl, n, Deltap, &Deltapl);
2147 #ifdef VERBOSE_DOUBLE /*PPC*/
2148 printf("The orbit of index %d is Deltap: ", p[i-1]); /*PPC*/
2149 print_list(Deltap, Deltapl, 1); /*PPC*/
2150 printf("Looking for digs moving index %d to slot %d\n", /*PPC*/
2151 p[i-1], b); /*PPC*/
2152 #endif /*PPC*/
2153 /* Calculate ALPHA and TAB */
2154 ALPHAstep[i+1]=ALPHAstep[i];
2155 for(l=ALPHAstep[i-1]; l<ALPHAstep[i]; l++) {
2156 #ifdef VERBOSE_DOUBLE /*PPC*/
2157 printf("Loop with l=%d\n", l); /*PPC*/
2158 #endif /*PPC*/
2159 Ll=ALPHA[l].Ll;
2160 copy_list(ALPHA[l].L, L, Ll);
2161 #ifdef VERBOSE_DOUBLE /*PPC*/
2162 printf("L: ");
2163 print_list(L, Ll, 1); /*PPC*/
2164 #endif /*PPC*/
2165 copy_list(ALPHA[l].s, s, n);
2166 copy_list(ALPHA[l].d, d, n);
2167 #ifdef VERBOSE_DOUBLE /*PPC*/
2168 printf("TAB[L]={"); /*PPC*/
2169 print_perm(s, n, 0); /*PPC*/
2170 printf(", "); /*PPC*/
2171 print_perm(d, n, 0); /*PPC*/
2172 printf("}\n"); /*PPC*/
2173 #endif /*PPC*/
2174 list1l=Deltabl;
2175 for (c=0; c<list1l; c++)
2176 list1[c]=onpoints(Deltab[c], s, n);
2177 product(g, d, perm1, n);
2178 inverse(perm1, perm2, n);
2179 list2l=Deltapl;
2180 for (c=0; c<list2l; c++)
2181 list2[c]=onpoints(Deltap[c], perm2, n);
2182 #ifdef VERBOSE_DOUBLE /*PPC*/
2183 printf("NEXT: intersection of sets of slots "); /*PPC*/
2184 print_list(list1, list1l, 0);
2185 printf(" and "); /*PPC*/
2186 print_list(list2, list2l, 1); /*PPC*/
2187 #endif /*PPC*/
2188 intersection(list1, list1l, list2, list2l,
2189 NEXT, &NEXTl);
2190 #ifdef VERBOSE_DOUBLE /*PPC*/
2191 printf("Intermediate slots NEXT: "); /*PPC*/
2192 print_list(NEXT, NEXTl, 1); /*PPC*/
2193 #endif /*PPC*/
2194
2195 for(jj=0; jj<NEXTl; jj++) {
2196 j = NEXT[jj];
2197 inverse(s, perm1, n);
2198 trace_schreier(onpoints(j,perm1,n), nu,w, perm2,n);
2199 product(perm2, s, s1, n);
2200 #ifdef VERBOSE_DOUBLE /*PPC*/
2201 printf("From slot %d to interm. slot %d use s1: ", /*PPC*/
2202 b, j);
2203 print_perm(s1, n, 1); /*PPC*/
2204 #endif /*PPC*/
2205 product(g, d, perm2, n);
2206 trace_schreier(onpoints(j,perm2,n),nuD,wD,perm3,n);
2207 inverse(perm3, perm1, n);
2208 product(d, perm1, d1, n);
2209 copy_list(L, L1, Ll);
2210 L1l=Ll;
2211 L1[L1l++] = j;
2212 #ifdef VERBOSE_DOUBLE /*PPC*/
2213 printf("d1: ");
2214 print_perm(d1, n, 1); /*PPC*/
2215 printf("L1: ");
2216 print_list(L1, L1l, 1); /*PPC*/
2217 #endif /*PPC*/
2218
2219 kk=ALPHAstep[i+1];
2220 ALPHA[l].o[j-1] = kk;
2221 ALPHAl++;
2222 ALPHAstep[i+1]++;
2223 // std::cout << "resizing " << ALPHAl << std::endl;
2224 // ALPHA = (struct alphastruct*)realloc(ALPHA, ALPHAl* sizeof(struct alphastruct));
2225 ALPHA.resize(ALPHAl);
2226 ALPHA.back().init(n);
2227 copy_list(L1, ALPHA[kk].L, L1l);
2228 ALPHA[kk].Ll = L1l;
2229 copy_list(s1, ALPHA[kk].s, n);
2230 copy_list(d1, ALPHA[kk].d, n);
2231
2232 F2(s1, g, d1, perm1, n);
2233 #ifdef VERBOSE_DOUBLE /*PPC*/
2234 printf("This gives the new index configuration: "); /*PPC*/
2235 inverse(perm1, perm2, n); /*PPC*/
2236 print_perm(perm2, n, 1); /*PPC*/
2237 for(ii=0; ii<i; ii++) { /*PPC*/
2238 product(s1, g, perm2, n); /*PPC*/
2239 product(perm2, d1, perm3, n); /*PPC*/
2240 printf("Checking slot %d with point %d: ", /*PPC*/
2241 bS[ii], p[ii]); /*PPC*/
2242 if(onpoints(bS[ii], perm3, n)==p[ii]) /*PPC*/
2243 printf("True\n"); /*PPC*/
2244 else printf("*************FALSE************\n");/*PPC*/
2245 } /*PPC*/
2246 #endif /*PPC*/
2247 }
2248 }
2249 {
2250 /* Verify if there are 2 equal permutations
2251 of opposite sign in SgD */
2252 int *array= (int*)malloc(n*(ALPHAstep[i+1]-ALPHAstep[i])
2253 *sizeof(int));
2254 int arrayl= 0;
2255 #ifdef VERBOSE_DOUBLE /*PPC*/
2256 printf("Astep[i-1]=%d, Astep[i]=%d, Astep[i+1]=%d\n", /*PPC*/
2257 ALPHAstep[i-1], ALPHAstep[i], ALPHAstep[i+1]); /*PPC*/
2258 #endif /*PPC*/
2259 for(l=ALPHAstep[i]; l<ALPHAstep[i+1]; l++) {
2260 F2(ALPHA[l].s, g, ALPHA[l].d, perm1, n);
2261 copy_list(perm1, array+(arrayl++)*n, n);
2262 }
2263 #ifdef VERBOSE_DOUBLE /*PPC*/
2264 printf("Perform check.\n"); /*PPC*/
2265 #endif /*PPC*/
2266 result= consistency(array, arrayl, n);
2267 #ifdef VERBOSE_DOUBLE /*PPC*/
2268 printf("Result of check: %d\n", result); /*PPC*/
2269 #endif /*PPC*/
2270 free(array);
2271 if (!result) break;
2272 }
2273 /* Find the stabilizers S^(i+1) and D^(i+1) */
2274 stabilizer(&bS[i-1], 1, KS, KSl, n, KS, &KSl);
2275 dropdummyset(p[i-1], vds, vdsl, dummies, &dl);
2276 droprepeatedset(p[i-1], vrs, vrsl, repes, &rl);
2277 #ifdef VERBOSE_DOUBLE /*PPC*/
2278 printf("Remove perms of KS moving slot %d\n", b); /*PPC*/
2279 printf("Remove perms of KD moving index %d\n", p[i-1]); /*PPC*/
2280 #endif /*PPC*/
2281
2282 } /* End of main loop */
2283
2284 /* Result */
2285 if (result==0) {
2286 zeros(perm1, n);
2287 }
2288 else {
2289 l=ALPHAstep[i-1];
2290 F2(ALPHA[l].s, g, ALPHA[l].d, perm1, n);
2291 }
2292 copy_list(perm1, dcr, n);
2293
2294 /* Free allocated memory */
2295 // free(ALPHA);
2296 free(ALPHAstep);
2297 free(ig);
2298 free(drummies);
2299 free(drummyslots);
2300 free(drummytmp);
2301 free(drummytmp2);
2302 free(bS);
2303 free(bSsort);
2304 free(KS);
2305 free(KD);
2306 free(bD);
2307 free(nu);
2308 free(w);
2309 free(Deltab);
2310 free(DeltaD);
2311 free(IMAGES);
2312 free(IMAGESsorted);
2313 free(p);
2314 free(nuD);
2315 free(wD);
2316 free(Deltap);
2317 free(NEXT);
2318 free(L);
2319 free(L1);
2320 free(s);
2321 free(d);
2322 free(list1);
2323 free(list2);
2324 free(perm1);
2325 free(perm2);
2326 free(perm3);
2327 free(s1);
2328 free(d1);
2329 }
2330
2331 /**********************************************************************/
2332
2333 /* canonical_perm. JMM, 5 July 2003
2334 *
2335 * JMM, 26 June 2007. Function extended. The old function is kept
2336 * as a call to the new function for backwards compatibility.
2337 *
2338 * The "dummy" group D is given through the list dummyps of dpl pairs
2339 * of (initial) slots of dummies. The list freeps (length fl) contains
2340 * the slots of the free indices. Clearly 2*dpl+fl+2=n.
2341 * See notes for canonical_perm_ext.
2342 * Parameter ob is now useless. Kept for backwards compatibility.
2343 */
2344
canonical_perm(int * PERM,int SGSQ,int * base,int bl,int * GS,int m,int n,int * freeps,int fl,int * dummyps,int dpl,int,int metricQ,int * CPERM)2345 void canonical_perm(int *PERM,
2346 int SGSQ, int *base, int bl, int *GS, int m, int n,
2347 int *freeps, int fl, int *dummyps, int dpl, int, int metricQ,
2348 int *CPERM)
2349 {
2350
2351 int i;
2352 int vds;
2353 int mQ;
2354 int *repes= NULL;
2355 int *vrs= NULL;
2356 int *PERM1= (int*)malloc(n *sizeof(int));
2357 int *PERM2= (int*)malloc(n *sizeof(int));
2358 int *frees= (int*)malloc(fl *sizeof(int));
2359 int *dummies= (int*)malloc(2*dpl*sizeof(int));
2360
2361 /* Construct "vectors" vds and mQ */
2362 vds = 2*dpl;
2363 mQ = metricQ;
2364
2365 /* !!!!!!!! Change to Renato's notation !!!!!!!! */
2366 inverse(PERM, PERM1, n);
2367 for (i=0; i<fl; i++) {
2368 frees[i] = onpoints(freeps[i], PERM1, n);
2369 }
2370 for (i=0; i<2*dpl; i++) {
2371 dummies[i] = onpoints(dummyps[i], PERM1, n);
2372 }
2373
2374 /* Call new, extended function */
2375 canonical_perm_ext(PERM1, n, SGSQ, base, bl, GS, m,
2376 frees, fl, &vds, 1, dummies, 2*dpl, &mQ,
2377 vrs, 0, repes, 0,
2378 PERM2);
2379
2380 /* !!!!!!!! Change back to our notation !!!!!!!! */
2381 if (PERM2[0] != 0) inverse(PERM2, CPERM, n);
2382 else copy_list(PERM2, CPERM, n);
2383
2384 /* Free allocated space */
2385 free(PERM1);
2386 free(PERM2);
2387 free(frees);
2388 free(dummies);
2389 }
2390
2391 /**********************************************************************/
2392
2393 /* canonical_perm_ext. JMM, 26 June 2007
2394 *
2395 * This function finds a canonical representant of the permutation PERM
2396 * in the double_coset with group S (slot symmetries) and group D
2397 * (dummy index or repeated index symmetries).
2398 *
2399 * There are two steps: first S, then S and D. The "slot" group S is
2400 * given through the generating set GS (with m n-permutations), that
2401 * can be a SGS (then SGSQ must be 1) or not (then SGSQ must be 0).
2402 * In the former case base (length bl) contains the associated base;
2403 * in the latter case a SGS will be constructed using the Schreier_Sims
2404 * algorithm using the points in base as the first points for the base.
2405 * The algorithm first calls coset_rep with PERM and the free indices
2406 * converting PERM into PERM1. The SGS is then stabilized with respect
2407 * to those points. Then the double_coset_rep algorithm is called with
2408 * PERM1 returning PERM2, which is finally copied to CPERM.
2409 *
2410 * The "drummy" group D is given through the specification of so-called
2411 * dummysets and repeatedsets. A dummyset is a collection of pairs of
2412 * indices (first the up-index, then the down-index) and a symmetry
2413 * switch saying if there is a symmetric metric (switch 1), an
2414 * antisymmetric metric (-1) or no metric at all (0). A repeated set
2415 * is simply a list of the names of repeated indices. Both dummies and
2416 * repes are collectively called drummies. Note that dummysets and
2417 * repeatedsets contain **names** (i.e. positions in the canonical
2418 * configuration) in the initial configuration.
2419 *
2420 * The result is stored in the permutation CPERM.
2421 * Note that we input a pointer to GS, and not a pointer to a pointer
2422 * because there is no need to return a modified GS.
2423 *
2424 * This is probably the most important function of this code, and hence
2425 * it is worth to explain the meaning of the input fields:
2426 * PERM: pointer to permutation (Images notation) to canonicalize
2427 * n: degree of perm (length of list of images)
2428 * SGSQ: switch. 1 if we supply a SGS, 0 if we only give a GS
2429 * base: sorted list of points for the SGS. Start of base if SGSQ=0
2430 * bl: length of base
2431 * GS: pointer to the list of permutations forming the GS
2432 * m: number of permutations in the GS
2433 * freeps: list of initial names of free indices
2434 * fl: length of list freeps
2435 * vds: list of lengths of dummysets
2436 * vdsl: length of vds (number of dummysets)
2437 * dummies: list with pairs of names dummies
2438 * dl: length of list dummies (sum of elements of vdsl)
2439 * mQ: list (of length vdsl) with symmetry signs
2440 * vrs: list of lengths of repeatedsets
2441 * vrsl: length of rds (number of repeatedsets)
2442 * repes: list with repeated indices
2443 * rl: length of list repes (sum of elements of vrsl)
2444 * CPERM: pointer to return the canonicalized permutation
2445 */
2446
canonical_perm_ext(int * PERM,int n,int SGSQ,int * base,int bl,int * GS,int m,int * frees,int fl,int * vds,int vdsl,int * dummies,int dl,int * mQ,int * vrs,int vrsl,int * repes,int rl,int * CPERM)2447 void canonical_perm_ext(int *PERM, int n,
2448 int SGSQ, int *base, int bl, int *GS, int m,
2449 int *frees, int fl,
2450 int *vds, int vdsl, int *dummies, int dl, int *mQ,
2451 int *vrs, int vrsl, int *repes, int rl,
2452 int *CPERM)
2453 {
2454
2455 int i;
2456 int *freeps= (int*)malloc(fl*sizeof(int));
2457 int *PERM1= (int*)malloc(n *sizeof(int));
2458 int *PERM2= (int*)malloc(n *sizeof(int));
2459 int *newbase= (int*)malloc(n *sizeof(int)), newbl;
2460 int **newGS = NULL;
2461 int *pointer;
2462 int newm;
2463 int *tmpbase= (int*)malloc(n *sizeof(int)), tmpbl;
2464 int num=0;
2465
2466 pointer= (int*)malloc(m*n*sizeof(int));
2467 newGS= &pointer;
2468
2469 #ifdef VERBOSE_CANON /*PPC*/
2470 printf("can: input base: "); /*PPC*/
2471 print_list(base, bl, 1); /*PPC*/
2472 #endif /*PPC*/
2473
2474 if (!SGSQ) { /* Compute a Strong Generating Set */
2475 nonstable_points(base, bl, GS, m, n,
2476 tmpbase, &tmpbl);
2477 schreier_sims(tmpbase, tmpbl, GS, m, n,
2478 newbase, &newbl, newGS, &newm, &num);
2479 }
2480 else {
2481 copy_list(base, newbase, bl);
2482 newbl=bl;
2483 copy_list(GS, *newGS, m*n);
2484 newm=m;
2485 }
2486
2487 #ifdef VERBOSE_CANON /*PPC*/
2488 printf("can: SGS computed.\n"); /*PPC*/
2489 printf("can: newbase: "); /*PPC*/
2490 print_list(newbase, newbl, 1); /*PPC*/
2491 #endif /*PPC*/
2492
2493 /* Compute slots of free indices */
2494 inverse(PERM, PERM1, n);
2495 for (i=0; i<fl; i++) {
2496 freeps[i] = onpoints(frees[i], PERM1, n);
2497 }
2498 /* Call coset_rep algorithm. Result in PERM1 */
2499 coset_rep(PERM, n, newbase, newbl, *newGS, &newm,
2500 freeps, fl, PERM1);
2501 #ifdef VERBOSE_CANON /*PPC*/
2502 printf("can: Canonical perm after coset algorithm: "); /*PPC*/
2503 print_perm(PERM1, n, 1); /*PPC*/
2504 printf("New positions of free indices: "); /*PPC*/
2505 print_list(freeps, fl, 1); /*PPC*/
2506 #endif /*PPC*/
2507
2508 if (dl+rl==0) { /* No drummy indices */
2509 copy_list(PERM1, CPERM, n);
2510 }
2511 else {
2512
2513 complement(newbase, newbl, freeps, fl, 1, tmpbase, &tmpbl);
2514 copy_list(tmpbase, newbase, tmpbl);
2515 newbl=tmpbl;
2516 stabilizer(freeps, fl, *newGS, newm, n, *newGS, &newm);
2517 #ifdef VERBOSE_CANON /*PPC*/
2518 printf("can: newbase after fixing: "); /*PPC*/
2519 print_list(newbase, newbl, 1); /*PPC*/
2520 printf("can: newGS after fixing: "); /*PPC*/
2521 print_array_perm(*newGS, newm, n, 1); /*PPC*/
2522 #endif /*PPC*/
2523
2524 /* Apply dummy-indices algorithm. Result in PERM2 */
2525 #ifdef VERBOSE_CANON /*PPC*/
2526 printf("can: Starting double_coset algorithm.\n"); /*PPC*/
2527 #endif /*PPC*/
2528 double_coset_rep(PERM1, n, newbase, newbl, *newGS, newm,
2529 vds, vdsl, dummies, dl, mQ,
2530 vrs, vrsl, repes, rl, PERM2);
2531 #ifdef VERBOSE_CANON /*PPC*/
2532 printf("can: Finished double_coset algorithm.\n"); /*PPC*/
2533 #endif /*PPC*/
2534
2535 /* Copy to result */
2536 copy_list(PERM2, CPERM, n);
2537
2538 }
2539
2540 /* Free allocated memory */
2541 free(freeps);
2542 free(PERM1);
2543 free(PERM2);
2544 free(newbase);
2545 free(tmpbase);
2546 free(*newGS);
2547
2548 #ifdef VERBOSE_CANON /*PPC*/
2549 printf("************ END OF ALGORITHM ***********\n"); /*PPC*/
2550 #endif /*PPC*/
2551 }
2552