1 /****************************************************************************
2 **
3 ** This file is part of GAP, a system for computational discrete algebra.
4 **
5 ** Copyright of GAP belongs to its developers, whose names are too numerous
6 ** to list here. Please refer to the COPYRIGHT file for details.
7 **
8 ** SPDX-License-Identifier: GPL-2.0-or-later
9 **
10 ** This file contains the functions for permutations (small and large).
11 **
12 ** Mathematically a permutation is a bijective mapping of a finite set onto
13 ** itself. In \GAP\ this subset must always be of the form [ 1, 2, .., N ],
14 ** where N is at most $2^32$.
15 **
16 ** Internally a permutation <p> is viewed as a mapping of [ 0, 1, .., N-1 ],
17 ** because in C indexing of arrays is done with the origin 0 instead of 1.
18 ** A permutation is represented by a bag of type 'T_PERM2' or 'T_PERM4' of
19 ** the form
20 **
21 ** (CONST_)ADDR_OBJ(p)
22 ** |
23 ** v
24 ** +------+-------+-------+-------+-------+- - - -+-------+-------+
25 ** | inv. | image | image | image | image | | image | image |
26 ** | perm | of 0 | of 1 | of 2 | of 3 | | of N-2| of N-1|
27 ** +------+-------+-------+-------+-------+- - - -+-------+-------+
28 ** ^
29 ** |
30 ** (CONST_)ADDR_PERM2(p) resp. (CONST_)ADDR_PERM4(p)
31 **
32 ** The first entry of the bag <p> is either zero, or a reference to another
33 ** permutation representing the inverse of <p>. The remaining entries of the
34 ** bag form an array describing the permutation. For bags of type 'T_PERM2',
35 ** the entries are of type 'UInt2' (defined in 'system.h' as an 16 bit wide
36 ** unsigned integer type), for type 'T_PERM4' the entries are of type
37 ** 'UInt4' (defined as a 32bit wide unsigned integer type). The first of
38 ** these entries is the image of 0, the second is the image of 1, and so on.
39 ** Thus, the entry at C index <i> is the image of <i>, if we view the
40 ** permutation as mapping of [ 0, 1, 2, .., N-1 ] as described above.
41 **
42 ** Permutations are never shortened. For example, if the product of two
43 ** permutations of degree 100 is the identity, it is nevertheless stored as
44 ** array of length 100, in which the <i>-th entry is of course simply <i>.
45 ** Testing whether a product has trailing fixpoints would be pretty costly,
46 ** and permutations of equal degree can be handled by the functions faster.
47 */
48
49 extern "C" {
50
51 #include "permutat.h"
52
53 #include "ariths.h"
54 #include "bool.h"
55 #include "error.h"
56 #include "gapstate.h"
57 #include "integer.h"
58 #include "io.h"
59 #include "listfunc.h"
60 #include "lists.h"
61 #include "modules.h"
62 #include "opers.h"
63 #include "plist.h"
64 #include "precord.h"
65 #include "range.h"
66 #include "records.h"
67 #include "saveload.h"
68 #include "sysfiles.h"
69 #include "trans.h"
70
71 } // extern "C"
72
73 #include "permutat_intern.hh"
74
75
76 /****************************************************************************
77 **
78 *F IMAGE(<i>,<pt>,<dg>) . . . . . . image of <i> under <pt> of degree <dg>
79 **
80 ** 'IMAGE' returns the image of the point <i> under the permutation of
81 ** degree <dg> pointed to by <pt>. If the point <i> is greater than or
82 ** equal to <dg> the image is <i> itself.
83 **
84 ** 'IMAGE' is implemented as a macro so do not use it with arguments that
85 ** have side effects.
86 */
87 #define IMAGE(i,pt,dg) (((i) < (dg)) ? (pt)[(i)] : (i))
88
89
90 /****************************************************************************
91 **
92 *V IdentityPerm . . . . . . . . . . . . . . . . . . . identity permutation
93 **
94 ** 'IdentityPerm' is an identity permutation.
95 */
96 Obj IdentityPerm;
97
98
99
100 static ModuleStateOffset PermutatStateOffset = -1;
101
102 typedef struct {
103
104 /****************************************************************************
105 **
106 *V TmpPerm . . . . . . . handle of the buffer bag of the permutation package
107 **
108 ** 'TmpPerm' is the handle of a bag of type 'T_PERM4', which is created at
109 ** initialization time of this package. Functions in this package can use
110 ** this bag for whatever purpose they want. They have to make sure of
111 ** course that it is large enough.
112 ** The buffer is *not* guaranteed to have any particular value, routines
113 ** that require a zero-initialization need to do this at the start.
114 ** This buffer is only constructed once it is needed, to avoid startup
115 ** costs (particularly when starting new threads).
116 ** Use the UseTmpPerm(<size>) utility function to ensure it is constructed!
117 */
118 Obj TmpPerm;
119
120 } PermutatModuleState;
121
122 #define TmpPerm MODULE_STATE(Permutat).TmpPerm
123
124
UseTmpPerm(UInt size)125 static UInt1 * UseTmpPerm(UInt size)
126 {
127 if (TmpPerm == (Obj)0)
128 TmpPerm = NewBag(T_PERM4, size);
129 else if (SIZE_BAG(TmpPerm) < size)
130 ResizeBag(TmpPerm, size);
131 return (UInt1 *)(ADDR_OBJ(TmpPerm) + 1);
132 }
133
134 template <typename T>
ADDR_TMP_PERM()135 static inline T * ADDR_TMP_PERM()
136 {
137 // no GAP_ASSERT here on purpose
138 return (T *)(ADDR_OBJ(TmpPerm) + 1);
139 }
140
141
142 /****************************************************************************
143 **
144 *F TypePerm( <perm> ) . . . . . . . . . . . . . . . . type of a permutation
145 **
146 ** 'TypePerm' returns the type of permutations.
147 **
148 ** 'TypePerm' is the function in 'TypeObjFuncs' for permutations.
149 */
150 static Obj TYPE_PERM2;
151
152 static Obj TYPE_PERM4;
153
TypePerm2(Obj perm)154 static Obj TypePerm2(Obj perm)
155 {
156 return TYPE_PERM2;
157 }
158
TypePerm4(Obj perm)159 static Obj TypePerm4(Obj perm)
160 {
161 return TYPE_PERM4;
162 }
163
164 // Forward declarations
165 template <typename T>
166 static inline UInt LargestMovedPointPerm_(Obj perm);
167
168 template <typename T>
169 static Obj InvPerm(Obj perm);
170
171
172 /****************************************************************************
173 **
174 *F PrintPerm( <perm> ) . . . . . . . . . . . . . . . . . print a permutation
175 **
176 ** 'PrintPerm' prints the permutation <perm> in the usual cycle notation. It
177 ** uses the degree to print all points with same width, which looks nicer.
178 ** Linebreaks are prefered most after cycles and next most after commas.
179 **
180 ** It does not remember which points have already been printed. To avoid
181 ** printing a cycle twice each is printed with the smallest element first.
182 ** This may in the worst case, for (1,2,..,n), take n^2/2 steps, but is fast
183 ** enough to keep a terminal at 9600 baud busy for all but the extrem cases.
184 */
185 template <typename T>
PrintPerm(Obj perm)186 static void PrintPerm(Obj perm)
187 {
188 UInt degPerm; /* degree of the permutation */
189 const T * ptPerm; /* pointer to the permutation */
190 UInt p, q; /* loop variables */
191 UInt isId; /* permutation is the identity? */
192 const char * fmt1; /* common formats to print points */
193 const char * fmt2; /* common formats to print points */
194
195 /* set up the formats used, so all points are printed with equal width */
196 // Use LargestMovedPoint so printing is consistent
197 degPerm = LargestMovedPointPerm_<T>(perm);
198 if ( degPerm < 10 ) { fmt1 = "%>(%>%1d%<"; fmt2 = ",%>%1d%<"; }
199 else if ( degPerm < 100 ) { fmt1 = "%>(%>%2d%<"; fmt2 = ",%>%2d%<"; }
200 else if ( degPerm < 1000 ) { fmt1 = "%>(%>%3d%<"; fmt2 = ",%>%3d%<"; }
201 else if ( degPerm < 10000 ) { fmt1 = "%>(%>%4d%<"; fmt2 = ",%>%4d%<"; }
202 else { fmt1 = "%>(%>%5d%<"; fmt2 = ",%>%5d%<"; }
203
204 /* run through all points */
205 isId = 1;
206 ptPerm = CONST_ADDR_PERM<T>(perm);
207 for ( p = 0; p < degPerm; p++ ) {
208
209 /* find the smallest element in this cycle */
210 q = ptPerm[p];
211 while ( p < q ) q = ptPerm[q];
212
213 /* if the smallest is the one we started with lets print the cycle */
214 if ( p == q && ptPerm[p] != p ) {
215 isId = 0;
216 Pr(fmt1,(Int)(p+1),0L);
217 ptPerm = CONST_ADDR_PERM<T>(perm);
218 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
219 Pr(fmt2,(Int)(q+1),0L);
220 ptPerm = CONST_ADDR_PERM<T>(perm);
221 }
222 Pr("%<)",0L,0L);
223 /* restore pointer, in case Pr caused a garbage collection */
224 ptPerm = CONST_ADDR_PERM<T>(perm);
225 }
226
227 }
228
229 /* special case for the identity */
230 if ( isId ) Pr("()",0L,0L);
231 }
232
233
234 /****************************************************************************
235 **
236 *F EqPerm( <opL>, <opR> ) . . . . . . . test if two permutations are equal
237 **
238 ** 'EqPerm' returns 'true' if the two permutations <opL> and <opR> are equal
239 ** and 'false' otherwise.
240 **
241 ** Two permutations may be equal, even if the two sequences do not have the
242 ** same length, if the larger permutation fixes the exceeding points.
243 **
244 */
245 template <typename TL, typename TR>
EqPerm(Obj opL,Obj opR)246 static Int EqPerm(Obj opL, Obj opR)
247 {
248 UInt degL; /* degree of the left operand */
249 const TL * ptL; /* pointer to the left operand */
250 UInt degR; /* degree of the right operand */
251 const TR * ptR; /* pointer to the right operand */
252 UInt p; /* loop variable */
253
254 /* get the degrees */
255 degL = DEG_PERM<TL>(opL);
256 degR = DEG_PERM<TR>(opR);
257
258 /* set up the pointers */
259 ptL = CONST_ADDR_PERM<TL>(opL);
260 ptR = CONST_ADDR_PERM<TR>(opR);
261
262 /* search for a difference and return False if you find one */
263 if ( degL <= degR ) {
264 for ( p = 0; p < degL; p++ )
265 if ( *(ptL++) != *(ptR++) )
266 return 0L;
267 for ( p = degL; p < degR; p++ )
268 if ( p != *(ptR++) )
269 return 0L;
270 }
271 else {
272 for ( p = 0; p < degR; p++ )
273 if ( *(ptL++) != *(ptR++) )
274 return 0L;
275 for ( p = degR; p < degL; p++ )
276 if ( *(ptL++) != p )
277 return 0L;
278 }
279
280 /* otherwise they must be equal */
281 return 1L;
282 }
283
284
285 /****************************************************************************
286 **
287 *F LtPerm( <opL>, <opR> ) . test if one permutation is smaller than another
288 **
289 ** 'LtPerm' returns 'true' if the permutation <opL> is strictly less than
290 ** the permutation <opR>. Permutations are ordered lexicographically with
291 ** respect to the images of 1,2,.., etc.
292 */
293 template <typename TL, typename TR>
LtPerm(Obj opL,Obj opR)294 static Int LtPerm(Obj opL, Obj opR)
295 {
296 UInt degL; /* degree of the left operand */
297 const TL * ptL; /* pointer to the left operand */
298 UInt degR; /* degree of the right operand */
299 const TR * ptR; /* pointer to the right operand */
300 UInt p; /* loop variable */
301
302 /* get the degrees of the permutations */
303 degL = DEG_PERM<TL>(opL);
304 degR = DEG_PERM<TR>(opR);
305
306 /* set up the pointers */
307 ptL = CONST_ADDR_PERM<TL>(opL);
308 ptR = CONST_ADDR_PERM<TR>(opR);
309
310 /* search for a difference and return if you find one */
311 if ( degL <= degR ) {
312
313 for (p = 0; p < degL; p++, ptL++, ptR++)
314 if (*ptL != *ptR) {
315 return *ptL < *ptR;
316 }
317 for (p = degL; p < degR; p++, ptR++)
318 if (p != *ptR) {
319 return p < *ptR;
320 }
321 }
322 else {
323 for (p = 0; p < degR; p++, ptL++, ptR++)
324 if (*ptL != *ptR) {
325 return *ptL < *ptR;
326 }
327 for (p = degR; p < degL; p++, ptL++)
328 if (*ptL != p) {
329 return *ptL < p;
330 }
331 }
332
333 /* otherwise they must be equal */
334 return 0;
335 }
336
337
338 /****************************************************************************
339 **
340 *F ProdPerm( <opL>, <opR> ) . . . . . . . . . . . . product of permutations
341 **
342 ** 'ProdPerm' returns the product of the two permutations <opL> and <opR>.
343 **
344 ** This is a little bit tuned but should be sufficiently easy to understand.
345 */
346 template <typename TL, typename TR>
ProdPerm(Obj opL,Obj opR)347 static Obj ProdPerm(Obj opL, Obj opR)
348 {
349 typedef typename ResultType<TL,TR>::type Res;
350
351 Obj prd; /* handle of the product (result) */
352 UInt degP; /* degree of the product */
353 Res * ptP; /* pointer to the product */
354 UInt degL; /* degree of the left operand */
355 const TL * ptL; /* pointer to the left operand */
356 UInt degR; /* degree of the right operand */
357 const TR * ptR; /* pointer to the right operand */
358 UInt p; /* loop variable */
359
360 degL = DEG_PERM<TL>(opL);
361 degR = DEG_PERM<TR>(opR);
362
363 // handle trivial cases
364 if (degL == 0) {
365 return opR;
366 }
367 if (degR == 0) {
368 return opL;
369 }
370
371 // compute the size of the result and allocate a bag
372 degP = degL < degR ? degR : degL;
373 prd = NEW_PERM<Res>( degP );
374
375 /* set up the pointers */
376 ptL = CONST_ADDR_PERM<TL>(opL);
377 ptR = CONST_ADDR_PERM<TR>(opR);
378 ptP = ADDR_PERM<Res>(prd);
379
380 /* if the left (inner) permutation has smaller degree, it is very easy */
381 if ( degL <= degR ) {
382 for ( p = 0; p < degL; p++ )
383 *(ptP++) = ptR[ *(ptL++) ];
384 for ( p = degL; p < degR; p++ )
385 *(ptP++) = ptR[ p ];
386 }
387
388 /* otherwise we have to use the macro 'IMAGE' */
389 else {
390 for ( p = 0; p < degL; p++ )
391 *(ptP++) = IMAGE( ptL[ p ], ptR, degR );
392 }
393
394 /* return the result */
395 return prd;
396 }
397
398
399 /****************************************************************************
400 **
401 *F QuoPerm( <opL>, <opR> ) . . . . . . . . . . . . quotient of permutations
402 **
403 ** 'QuoPerm' returns the quotient of the permutations <opL> and <opR>, i.e.,
404 ** the product '<opL>\*<opR>\^-1'.
405 **
406 ** Unfortunatly this can not be done in <degree> steps, we need 2 * <degree>
407 ** steps.
408 */
QuoPerm(Obj opL,Obj opR)409 static Obj QuoPerm(Obj opL, Obj opR)
410 {
411 return PROD(opL, INV(opR));
412 }
413
414
415 /****************************************************************************
416 **
417 *F LQuoPerm( <opL>, <opR> ) . . . . . . . . . left quotient of permutations
418 **
419 ** 'LQuoPerm' returns the left quotient of the two permutations <opL> and
420 ** <opR>, i.e., the value of '<opL>\^-1*<opR>', which sometimes comes handy.
421 **
422 ** This can be done as fast as a single multiplication or inversion.
423 */
424 template <typename TL, typename TR>
LQuoPerm(Obj opL,Obj opR)425 static Obj LQuoPerm(Obj opL, Obj opR)
426 {
427 typedef typename ResultType<TL,TR>::type Res;
428
429 Obj mod; /* handle of the quotient (result) */
430 UInt degM; /* degree of the quotient */
431 Res * ptM; /* pointer to the quotient */
432 UInt degL; /* degree of the left operand */
433 const TL * ptL; /* pointer to the left operand */
434 UInt degR; /* degree of the right operand */
435 const TR * ptR; /* pointer to the right operand */
436 UInt p; /* loop variable */
437
438 degL = DEG_PERM<TL>(opL);
439 degR = DEG_PERM<TR>(opR);
440
441 // handle trivial cases
442 if (degL == 0) {
443 return opR;
444 }
445 if (degR == 0) {
446 return InvPerm<TL>(opL);
447 }
448
449 // compute the size of the result and allocate a bag
450 degM = degL < degR ? degR : degL;
451 mod = NEW_PERM<Res>( degM );
452
453 /* set up the pointers */
454 ptL = CONST_ADDR_PERM<TL>(opL);
455 ptR = CONST_ADDR_PERM<TR>(opR);
456 ptM = ADDR_PERM<Res>(mod);
457
458 /* it is one thing if the left (inner) permutation is smaller */
459 if ( degL <= degR ) {
460 for ( p = 0; p < degL; p++ )
461 ptM[ *(ptL++) ] = *(ptR++);
462 for ( p = degL; p < degR; p++ )
463 ptM[ p ] = *(ptR++);
464 }
465
466 /* and another if the right (outer) permutation is smaller */
467 else {
468 for ( p = 0; p < degR; p++ )
469 ptM[ *(ptL++) ] = *(ptR++);
470 for ( p = degR; p < degL; p++ )
471 ptM[ *(ptL++) ] = p;
472 }
473
474 /* return the result */
475 return mod;
476 }
477
478
479 /****************************************************************************
480 **
481 *F InvPerm( <perm> ) . . . . . . . . . . . . . . . inverse of a permutation
482 */
483 template <typename T>
InvPerm(Obj perm)484 static Obj InvPerm(Obj perm)
485 {
486 Obj inv; /* handle of the inverse (result) */
487 T * ptInv; /* pointer to the inverse */
488 const T * ptPerm; /* pointer to the permutation */
489 UInt deg; /* degree of the permutation */
490 UInt p; /* loop variables */
491
492 inv = STOREDINV_PERM(perm);
493 if (inv != 0)
494 return inv;
495
496 deg = DEG_PERM<T>(perm);
497 inv = NEW_PERM<T>(deg);
498
499 // get pointer to the permutation and the inverse
500 ptPerm = CONST_ADDR_PERM<T>(perm);
501 ptInv = ADDR_PERM<T>(inv);
502
503 // invert the permutation
504 for ( p = 0; p < deg; p++ )
505 ptInv[ *ptPerm++ ] = p;
506
507 // store and return the inverse
508 SET_STOREDINV_PERM(perm, inv);
509 return inv;
510 }
511
512
513 /****************************************************************************
514 **
515 *F PowPermInt( <opL>, <opR> ) . . . . . . . integer power of a permutation
516 **
517 ** 'PowPermInt' returns the <opR>-th power of the permutation <opL>. <opR>
518 ** must be a small integer.
519 **
520 ** This repeatedly applies the permutation <opR> to all points which seems
521 ** to be faster than binary powering, and does not need temporary storage.
522 */
523 template <typename T>
PowPermInt(Obj opL,Obj opR)524 static Obj PowPermInt(Obj opL, Obj opR)
525 {
526 Obj pow; /* handle of the power (result) */
527 T * ptP; /* pointer to the power */
528 const T * ptL; /* pointer to the permutation */
529 UInt1 * ptKnown; /* pointer to temporary bag */
530 UInt deg; /* degree of the permutation */
531 Int exp, e; /* exponent (right operand) */
532 UInt len; /* length of cycle (result) */
533 UInt p, q, r; /* loop variables */
534
535
536 /* handle zeroth and first powers and stored inverses separately */
537 if ( opR == INTOBJ_INT(0))
538 return IdentityPerm;
539 if ( opR == INTOBJ_INT(1))
540 return opL;
541 if (opR == INTOBJ_INT(-1))
542 return InvPerm<T>(opL);
543
544 deg = DEG_PERM<T>(opL);
545
546 // handle trivial permutations
547 if (deg == 0) {
548 return IdentityPerm;
549 }
550
551 // allocate a result bag
552 pow = NEW_PERM<T>(deg);
553
554 /* compute the power by repeated mapping for small positive exponents */
555 if ( IS_INTOBJ(opR)
556 && 2 <= INT_INTOBJ(opR) && INT_INTOBJ(opR) < 8 ) {
557
558 /* get pointer to the permutation and the power */
559 exp = INT_INTOBJ(opR);
560 ptL = CONST_ADDR_PERM<T>(opL);
561 ptP = ADDR_PERM<T>(pow);
562
563 /* loop over the points of the permutation */
564 for ( p = 0; p < deg; p++ ) {
565 q = p;
566 for ( e = 0; e < exp; e++ )
567 q = ptL[q];
568 ptP[p] = q;
569 }
570
571 }
572
573 /* compute the power by raising the cycles individually for large exps */
574 else if ( IS_INTOBJ(opR) && 8 <= INT_INTOBJ(opR) ) {
575
576 /* make sure that the buffer bag is large enough */
577 UseTmpPerm(SIZE_OBJ(opL));
578 ptKnown = ADDR_TMP_PERM<UInt1>();
579
580 /* clear the buffer bag */
581 memset(ptKnown, 0, DEG_PERM<T>(opL));
582
583 /* get pointer to the permutation and the power */
584 exp = INT_INTOBJ(opR);
585 ptL = CONST_ADDR_PERM<T>(opL);
586 ptP = ADDR_PERM<T>(pow);
587
588 /* loop over all cycles */
589 for ( p = 0; p < deg; p++ ) {
590
591 /* if we haven't looked at this cycle so far */
592 if ( ptKnown[p] == 0 ) {
593
594 /* find the length of this cycle */
595 len = 1;
596 for ( q = ptL[p]; q != p; q = ptL[q] ) {
597 len++; ptKnown[q] = 1;
598 }
599
600 /* raise this cycle to the power <exp> mod <len> */
601 r = p;
602 for ( e = 0; e < exp % len; e++ )
603 r = ptL[r];
604 ptP[p] = r;
605 r = ptL[r];
606 for ( q = ptL[p]; q != p; q = ptL[q] ) {
607 ptP[q] = r;
608 r = ptL[r];
609 }
610
611 }
612
613 }
614
615
616 }
617
618 /* compute the power by raising the cycles individually for large exps */
619 else if ( TNUM_OBJ(opR) == T_INTPOS ) {
620
621 /* make sure that the buffer bag is large enough */
622 UseTmpPerm(SIZE_OBJ(opL));
623 ptKnown = ADDR_TMP_PERM<UInt1>();
624
625 /* clear the buffer bag */
626 memset(ptKnown, 0, DEG_PERM<T>(opL));
627
628 /* get pointer to the permutation and the power */
629 ptL = CONST_ADDR_PERM<T>(opL);
630 ptP = ADDR_PERM<T>(pow);
631
632 /* loop over all cycles */
633 for ( p = 0; p < deg; p++ ) {
634
635 /* if we haven't looked at this cycle so far */
636 if ( ptKnown[p] == 0 ) {
637
638 /* find the length of this cycle */
639 len = 1;
640 for ( q = ptL[p]; q != p; q = ptL[q] ) {
641 len++; ptKnown[q] = 1;
642 }
643
644 /* raise this cycle to the power <exp> mod <len> */
645 r = p;
646 exp = INT_INTOBJ( ModInt( opR, INTOBJ_INT(len) ) );
647 for ( e = 0; e < exp; e++ )
648 r = ptL[r];
649 ptP[p] = r;
650 r = ptL[r];
651 for ( q = ptL[p]; q != p; q = ptL[q] ) {
652 ptP[q] = r;
653 r = ptL[r];
654 }
655
656 }
657
658 }
659
660 }
661
662 /* compute the power by repeated mapping for small negative exponents */
663 else if ( IS_INTOBJ(opR)
664 && -8 < INT_INTOBJ(opR) && INT_INTOBJ(opR) < 0 ) {
665
666 /* get pointer to the permutation and the power */
667 exp = -INT_INTOBJ(opR);
668 ptL = CONST_ADDR_PERM<T>(opL);
669 ptP = ADDR_PERM<T>(pow);
670
671 /* loop over the points */
672 for ( p = 0; p < deg; p++ ) {
673 q = p;
674 for ( e = 0; e < exp; e++ )
675 q = ptL[q];
676 ptP[q] = p;
677 }
678
679 }
680
681 /* compute the power by raising the cycles individually for large exps */
682 else if ( IS_INTOBJ(opR) && INT_INTOBJ(opR) <= -8 ) {
683
684 /* make sure that the buffer bag is large enough */
685 UseTmpPerm(SIZE_OBJ(opL));
686 ptKnown = ADDR_TMP_PERM<UInt1>();
687
688 /* clear the buffer bag */
689 memset(ptKnown, 0, DEG_PERM<T>(opL));
690
691 /* get pointer to the permutation and the power */
692 exp = -INT_INTOBJ(opR);
693 ptL = CONST_ADDR_PERM<T>(opL);
694 ptP = ADDR_PERM<T>(pow);
695
696 /* loop over all cycles */
697 for ( p = 0; p < deg; p++ ) {
698
699 /* if we haven't looked at this cycle so far */
700 if ( ptKnown[p] == 0 ) {
701
702 /* find the length of this cycle */
703 len = 1;
704 for ( q = ptL[p]; q != p; q = ptL[q] ) {
705 len++; ptKnown[q] = 1;
706 }
707
708 /* raise this cycle to the power <exp> mod <len> */
709 r = p;
710 for ( e = 0; e < exp % len; e++ )
711 r = ptL[r];
712 ptP[r] = p;
713 r = ptL[r];
714 for ( q = ptL[p]; q != p; q = ptL[q] ) {
715 ptP[r] = q;
716 r = ptL[r];
717 }
718
719 }
720
721 }
722
723 }
724
725 /* compute the power by raising the cycles individually for large exps */
726 else if ( TNUM_OBJ(opR) == T_INTNEG ) {
727 /* do negation first as it can cause a garbage collection */
728 opR = AInvInt(opR);
729
730 /* make sure that the buffer bag is large enough */
731 UseTmpPerm(SIZE_OBJ(opL));
732 ptKnown = ADDR_TMP_PERM<UInt1>();
733
734 /* clear the buffer bag */
735 memset(ptKnown, 0, DEG_PERM<T>(opL));
736
737 /* get pointer to the permutation and the power */
738 ptL = CONST_ADDR_PERM<T>(opL);
739 ptP = ADDR_PERM<T>(pow);
740
741 /* loop over all cycles */
742 for ( p = 0; p < deg; p++ ) {
743
744 /* if we haven't looked at this cycle so far */
745 if ( ptKnown[p] == 0 ) {
746
747 /* find the length of this cycle */
748 len = 1;
749 for ( q = ptL[p]; q != p; q = ptL[q] ) {
750 len++; ptKnown[q] = 1;
751 }
752
753 /* raise this cycle to the power <exp> mod <len> */
754 r = p;
755 exp = INT_INTOBJ( ModInt( opR, INTOBJ_INT(len) ) );
756 for ( e = 0; e < exp % len; e++ )
757 r = ptL[r];
758 ptP[r] = p;
759 r = ptL[r];
760 for ( q = ptL[p]; q != p; q = ptL[q] ) {
761 ptP[r] = q;
762 r = ptL[r];
763 }
764
765 }
766
767 }
768
769 }
770
771 /* return the result */
772 return pow;
773 }
774
775
776 /****************************************************************************
777 **
778 *F PowIntPerm( <opL>, <opR> ) . . . image of an integer under a permutation
779 **
780 ** 'PowIntPerm' returns the image of the positive integer <opL> under the
781 ** permutation <opR>. If <opL> is larger than the degree of <opR> it is a
782 ** fixpoint of the permutation and thus simply returned.
783 */
784 template <typename T>
PowIntPerm(Obj opL,Obj opR)785 static Obj PowIntPerm(Obj opL, Obj opR)
786 {
787 Int img; /* image (result) */
788
789 GAP_ASSERT(TNUM_OBJ(opL) == T_INTPOS || TNUM_OBJ(opL) == T_INT);
790
791 /* large positive integers (> 2^28-1) are fixed by any permutation */
792 if ( TNUM_OBJ(opL) == T_INTPOS )
793 return opL;
794
795 img = INT_INTOBJ(opL);
796 RequireArgumentConditionEx("PowIntPerm", opL, "<point>", img > 0,
797 "must be a positive integer");
798
799 /* compute the image */
800 if ( img <= DEG_PERM<T>(opR) ) {
801 img = (CONST_ADDR_PERM<T>(opR))[img-1] + 1;
802 }
803
804 /* return it */
805 return INTOBJ_INT(img);
806 }
807
808
809 /****************************************************************************
810 **
811 *F QuoIntPerm( <opL>, <opR> ) . preimage of an integer under a permutation
812 **
813 ** 'QuoIntPerm' returns the preimage of the preimage integer <opL> under the
814 ** permutation <opR>. If <opL> is larger than the degree of <opR> is is a
815 ** fixpoint, and thus simply returned.
816 **
817 ** There are basically two ways to find the preimage. One is to run through
818 ** <opR> and look for <opL>. The index where it's found is the preimage.
819 ** The other is to find the image of <opL> under <opR>, the image of that
820 ** point and so on, until we come back to <opL>. The last point is the
821 ** preimage of <opL>. This is faster because the cycles are usually short.
822 */
823 static Obj PERM_INVERSE_THRESHOLD;
824
825 template <typename T>
QuoIntPerm(Obj opL,Obj opR)826 static Obj QuoIntPerm(Obj opL, Obj opR)
827 {
828 T pre; /* preimage (result) */
829 Int img; /* image (left operand) */
830 const T * ptR; /* pointer to the permutation */
831
832 GAP_ASSERT(TNUM_OBJ(opL) == T_INTPOS || TNUM_OBJ(opL) == T_INT);
833
834 /* large positive integers (> 2^28-1) are fixed by any permutation */
835 if ( TNUM_OBJ(opL) == T_INTPOS )
836 return opL;
837
838 img = INT_INTOBJ(opL);
839 RequireArgumentConditionEx("QuoIntPerm", opL, "<point>", img > 0,
840 "must be a positive integer");
841
842 Obj inv = STOREDINV_PERM(opR);
843
844 if (inv == 0 && PERM_INVERSE_THRESHOLD != 0 &&
845 IS_INTOBJ(PERM_INVERSE_THRESHOLD) &&
846 DEG_PERM<T>(opR) <= INT_INTOBJ(PERM_INVERSE_THRESHOLD))
847 inv = InvPerm<T>(opR);
848
849 if (inv != 0)
850 return INTOBJ_INT(
851 IMAGE(img - 1, CONST_ADDR_PERM<T>(inv), DEG_PERM<T>(inv)) + 1);
852
853 /* compute the preimage */
854 if ( img <= DEG_PERM<T>(opR) ) {
855 pre = T(img - 1);
856 ptR = CONST_ADDR_PERM<T>(opR);
857 while (ptR[pre] != T(img - 1))
858 pre = ptR[pre];
859 /* return it */
860 return INTOBJ_INT(pre + 1);
861 }
862 else
863 return INTOBJ_INT(img);
864 }
865
866
867 /****************************************************************************
868 **
869 *F PowPerm( <opL>, <opR> ) . . . . . . . . . . . conjugation of permutations
870 **
871 ** 'PowPerm' returns the conjugation of the two permutations <opL> and
872 ** <opR>, that s defined as the following product '<opR>\^-1 \*\ <opL> \*\
873 ** <opR>'.
874 */
875 template <typename TL, typename TR>
PowPerm(Obj opL,Obj opR)876 static Obj PowPerm(Obj opL, Obj opR)
877 {
878 typedef typename ResultType<TL,TR>::type Res;
879
880 Obj cnj; /* handle of the conjugation (res) */
881 UInt degC; /* degree of the conjugation */
882 Res * ptC; /* pointer to the conjugation */
883 UInt degL; /* degree of the left operand */
884 const TL * ptL; /* pointer to the left operand */
885 UInt degR; /* degree of the right operand */
886 const TR * ptR; /* pointer to the right operand */
887 UInt p; /* loop variable */
888
889 degL = DEG_PERM<TL>(opL);
890 degR = DEG_PERM<TR>(opR);
891
892 // handle trivial cases
893 if (degL == 0) {
894 return IdentityPerm;
895 }
896
897 if (degR == 0) {
898 return opL;
899 }
900
901 // compute the size of the result and allocate a bag
902 degC = degL < degR ? degR : degL;
903 cnj = NEW_PERM<Res>( degC );
904
905 /* set up the pointers */
906 ptL = CONST_ADDR_PERM<TL>(opL);
907 ptR = CONST_ADDR_PERM<TR>(opR);
908 ptC = ADDR_PERM<Res>(cnj);
909
910 /* it is faster if the both permutations have the same size */
911 if ( degL == degR ) {
912 for ( p = 0; p < degC; p++ )
913 ptC[ ptR[p] ] = ptR[ ptL[p] ];
914 }
915
916 /* otherwise we have to use the macro 'IMAGE' three times */
917 else {
918 for ( p = 0; p < degC; p++ )
919 ptC[ IMAGE(p,ptR,degR) ] = IMAGE( IMAGE(p,ptL,degL), ptR, degR );
920 }
921
922 /* return the result */
923 return cnj;
924 }
925
926
927 /****************************************************************************
928 **
929 *F CommPerm( <opL>, <opR> ) . . . . . . . . commutator of two permutations
930 **
931 ** 'CommPerm' returns the commutator of the two permutations <opL> and
932 ** <opR>, that is defined as '<hd>\^-1 \*\ <opR>\^-1 \*\ <opL> \*\ <opR>'.
933 */
934 template <typename TL, typename TR>
CommPerm(Obj opL,Obj opR)935 static Obj CommPerm(Obj opL, Obj opR)
936 {
937 typedef typename ResultType<TL,TR>::type Res;
938
939 Obj com; /* handle of the commutator (res) */
940 UInt degC; /* degree of the commutator */
941 Res * ptC; /* pointer to the commutator */
942 UInt degL; /* degree of the left operand */
943 const TL * ptL; /* pointer to the left operand */
944 UInt degR; /* degree of the right operand */
945 const TR * ptR; /* pointer to the right operand */
946 UInt p; /* loop variable */
947
948 degL = DEG_PERM<TL>(opL);
949 degR = DEG_PERM<TR>(opR);
950
951 // handle trivial cases
952 if (degL == 0 || degR == 0) {
953 return IdentityPerm;
954 }
955
956 // compute the size of the result and allocate a bag
957 degC = degL < degR ? degR : degL;
958 com = NEW_PERM<Res>( degC );
959
960 /* set up the pointers */
961 ptL = CONST_ADDR_PERM<TL>(opL);
962 ptR = CONST_ADDR_PERM<TR>(opR);
963 ptC = ADDR_PERM<Res>(com);
964
965 /* it is faster if the both permutations have the same size */
966 if ( degL == degR ) {
967 for ( p = 0; p < degC; p++ )
968 ptC[ ptL[ ptR[ p ] ] ] = ptR[ ptL[ p ] ];
969 }
970
971 /* otherwise we have to use the macro 'IMAGE' four times */
972 else {
973 for ( p = 0; p < degC; p++ )
974 ptC[ IMAGE( IMAGE(p,ptR,degR), ptL, degL ) ]
975 = IMAGE( IMAGE(p,ptL,degL), ptR, degR );
976 }
977
978 /* return the result */
979 return com;
980 }
981
982
983 /****************************************************************************
984 **
985 *F OnePerm( <perm> )
986 */
OnePerm(Obj op)987 static Obj OnePerm(Obj op)
988 {
989 return IdentityPerm;
990 }
991
992
993
994 /****************************************************************************
995 **
996 *F FiltIS_PERM( <self>, <val> ) . . . . . . test if a value is a permutation
997 **
998 ** 'FiltIS_PERM' implements the internal function 'IsPerm'.
999 **
1000 ** 'IsPerm( <val> )'
1001 **
1002 ** 'IsPerm' returns 'true' if the value <val> is a permutation and 'false'
1003 ** otherwise.
1004 */
1005 static Obj IsPermFilt;
1006
FiltIS_PERM(Obj self,Obj val)1007 static Obj FiltIS_PERM(Obj self, Obj val)
1008 {
1009 /* return 'true' if <val> is a permutation and 'false' otherwise */
1010 if ( TNUM_OBJ(val) == T_PERM2 || TNUM_OBJ(val) == T_PERM4 ) {
1011 return True;
1012 }
1013 else if ( TNUM_OBJ(val) < FIRST_EXTERNAL_TNUM ) {
1014 return False;
1015 }
1016 else {
1017 return DoFilter( self, val );
1018 }
1019 }
1020
1021
1022 /****************************************************************************
1023 **
1024 *F FuncPermList( <self>, <list> ) . . . . . convert a list to a permutation
1025 **
1026 ** 'FuncPermList' implements the internal function 'PermList'
1027 **
1028 ** 'PermList( <list> )'
1029 **
1030 ** Converts the list <list> into a permutation, which is then returned.
1031 **
1032 ** 'FuncPermList' simply copies the list pointwise into a permutation bag.
1033 ** It also does some checks to make sure that the list is a permutation.
1034 */
1035 template <typename T>
PermList(Obj list)1036 static inline Obj PermList(Obj list)
1037 {
1038 Obj perm; /* handle of the permutation */
1039 T * ptPerm; /* pointer to the permutation */
1040 UInt degPerm; /* degree of the permutation */
1041 const Obj * ptList; /* pointer to the list */
1042 T * ptTmp; /* pointer to the buffer bag */
1043 Int i, k; /* loop variables */
1044
1045 PLAIN_LIST( list );
1046
1047 degPerm = LEN_LIST( list );
1048
1049 /* make sure that the global buffer bag is large enough for checkin*/
1050 UseTmpPerm(SIZEBAG_PERM<T>(degPerm));
1051
1052 /* allocate the bag for the permutation and get pointer */
1053 perm = NEW_PERM<T>(degPerm);
1054 ptPerm = ADDR_PERM<T>(perm);
1055 ptList = CONST_ADDR_OBJ(list);
1056 ptTmp = ADDR_TMP_PERM<T>();
1057
1058 /* make the buffer bag clean */
1059 for ( i = 1; i <= degPerm; i++ )
1060 ptTmp[i-1] = 0;
1061
1062 /* run through all entries of the list */
1063 for ( i = 1; i <= degPerm; i++ ) {
1064
1065 /* get the <i>th entry of the list */
1066 if ( ptList[i] == 0 ) {
1067 return Fail;
1068 }
1069 if ( !IS_INTOBJ(ptList[i]) ) {
1070 return Fail;
1071 }
1072 k = INT_INTOBJ(ptList[i]);
1073 if ( k <= 0 || degPerm < k ) {
1074 return Fail;
1075 }
1076
1077 /* make sure we haven't seen this entry yet */
1078 if ( ptTmp[k-1] != 0 ) {
1079 return Fail;
1080 }
1081 ptTmp[k-1] = 1;
1082
1083 /* and finally copy it into the permutation */
1084 ptPerm[i-1] = k-1;
1085 }
1086
1087 /* return the permutation */
1088 return perm;
1089 }
1090
FuncPermList(Obj self,Obj list)1091 static Obj FuncPermList(Obj self, Obj list)
1092 {
1093 /* check the arguments */
1094 RequireSmallList("PermList", list);
1095
1096 UInt len = LEN_LIST( list );
1097 if ( len <= 65536 ) {
1098 return PermList<UInt2>(list);
1099 }
1100 else if (len <= MAX_DEG_PERM4) {
1101 return PermList<UInt4>(list);
1102 }
1103 else {
1104 ErrorMayQuit("PermList: list length %i exceeds maximum permutation degree\n",
1105 len, 0);
1106 }
1107 }
1108
1109 /****************************************************************************
1110 **
1111 *F LargestMovedPointPerm( <perm> ) largest point moved by perm
1112 **
1113 ** 'LargestMovedPointPerm' returns the largest positive integer that is
1114 ** moved by the permutation <perm>.
1115 **
1116 ** This is easy, except that permutations may contain trailing fixpoints.
1117 */
1118 template <typename T>
LargestMovedPointPerm_(Obj perm)1119 static inline UInt LargestMovedPointPerm_(Obj perm)
1120 {
1121 UInt sup;
1122 const T * ptPerm;
1123
1124 ptPerm = CONST_ADDR_PERM<T>(perm);
1125 for (sup = DEG_PERM<T>(perm); 1 <= sup; sup--) {
1126 if (ptPerm[sup - 1] != sup - 1)
1127 break;
1128 }
1129 return sup;
1130 }
1131
LargestMovedPointPerm(Obj perm)1132 UInt LargestMovedPointPerm(Obj perm)
1133 {
1134 GAP_ASSERT(TNUM_OBJ(perm) == T_PERM2 || TNUM_OBJ(perm) == T_PERM4);
1135
1136 if (TNUM_OBJ(perm) == T_PERM2)
1137 return LargestMovedPointPerm_<UInt2>(perm);
1138 else
1139 return LargestMovedPointPerm_<UInt4>(perm);
1140 }
1141
1142
1143 /****************************************************************************
1144 **
1145 *F FuncLARGEST_MOVED_POINT_PERM( <self>, <perm> ) largest point moved by perm
1146 **
1147 ** GAP-level wrapper for 'LargestMovedPointPerm'.
1148 */
FuncLARGEST_MOVED_POINT_PERM(Obj self,Obj perm)1149 static Obj FuncLARGEST_MOVED_POINT_PERM(Obj self, Obj perm)
1150 {
1151
1152 /* check the argument */
1153 RequirePermutation("LargestMovedPointPerm", perm);
1154
1155 return INTOBJ_INT(LargestMovedPointPerm(perm));
1156 }
1157
1158
1159 /****************************************************************************
1160 **
1161 *F FuncCYCLE_LENGTH_PERM_INT( <self>, <perm>, <point> ) . . . . . . . . . .
1162 *F . . . . . . . . . . . . . . . . . . length of a cycle under a permutation
1163 **
1164 ** 'FuncCycleLengthInt' implements the internal function
1165 ** 'CycleLengthPermInt'
1166 **
1167 ** 'CycleLengthPermInt( <perm>, <point> )'
1168 **
1169 ** 'CycleLengthPermInt' returns the length of the cycle of <point>, which
1170 ** must be a positive integer, under the permutation <perm>.
1171 **
1172 ** Note that the order of the arguments to this function has been reversed.
1173 */
1174 template <typename T>
CYCLE_LENGTH_PERM_INT(Obj perm,UInt pnt)1175 static inline Obj CYCLE_LENGTH_PERM_INT(Obj perm, UInt pnt)
1176 {
1177 const T * ptPerm; /* pointer to the permutation */
1178 UInt deg; /* degree of the permutation */
1179 UInt len; /* length of cycle (result) */
1180 UInt p; /* loop variable */
1181
1182 /* get pointer to the permutation, the degree, and the point */
1183 ptPerm = CONST_ADDR_PERM<T>(perm);
1184 deg = DEG_PERM<T>(perm);
1185
1186 /* now compute the length by looping over the cycle */
1187 len = 1;
1188 if ( pnt < deg ) {
1189 for ( p = ptPerm[pnt]; p != pnt; p = ptPerm[p] )
1190 len++;
1191 }
1192
1193 return INTOBJ_INT(len);
1194 }
1195
FuncCYCLE_LENGTH_PERM_INT(Obj self,Obj perm,Obj point)1196 static Obj FuncCYCLE_LENGTH_PERM_INT(Obj self, Obj perm, Obj point)
1197 {
1198 UInt pnt; /* value of the point */
1199
1200 /* evaluate and check the arguments */
1201 RequirePermutation("CycleLengthPermInt", perm);
1202 pnt = GetPositiveSmallInt("CycleLengthPermInt", point) - 1;
1203
1204 if ( TNUM_OBJ(perm) == T_PERM2 ) {
1205 return CYCLE_LENGTH_PERM_INT<UInt2>(perm, pnt);
1206 }
1207 else {
1208 return CYCLE_LENGTH_PERM_INT<UInt4>(perm, pnt);
1209 }
1210 }
1211
1212
1213 /****************************************************************************
1214 **
1215 *F FuncCYCLE_PERM_INT( <self>, <perm>, <point> ) . . cycle of a permutation
1216 *
1217 ** 'FuncCYCLE_PERM_INT' implements the internal function 'CyclePermInt'.
1218 **
1219 ** 'CyclePermInt( <perm>, <point> )'
1220 **
1221 ** 'CyclePermInt' returns the cycle of <point>, which must be a positive
1222 ** integer, under the permutation <perm> as a list.
1223 */
1224 template <typename T>
CYCLE_PERM_INT(Obj perm,UInt pnt)1225 static inline Obj CYCLE_PERM_INT(Obj perm, UInt pnt)
1226 {
1227 Obj list; /* handle of the list (result) */
1228 Obj * ptList; /* pointer to the list */
1229 const T * ptPerm; /* pointer to the permutation */
1230 UInt deg; /* degree of the permutation */
1231 UInt len; /* length of the cycle */
1232 UInt p; /* loop variable */
1233
1234 /* get pointer to the permutation, the degree, and the point */
1235 ptPerm = CONST_ADDR_PERM<T>(perm);
1236 deg = DEG_PERM<T>(perm);
1237
1238 /* now compute the length by looping over the cycle */
1239 len = 1;
1240 if ( pnt < deg ) {
1241 for ( p = ptPerm[pnt]; p != pnt; p = ptPerm[p] )
1242 len++;
1243 }
1244
1245 /* allocate the list */
1246 list = NEW_PLIST( T_PLIST, len );
1247 SET_LEN_PLIST( list, len );
1248 ptList = ADDR_OBJ(list);
1249 ptPerm = CONST_ADDR_PERM<T>(perm);
1250
1251 /* copy the points into the list */
1252 len = 1;
1253 ptList[len++] = INTOBJ_INT( pnt+1 );
1254 if ( pnt < deg ) {
1255 for ( p = ptPerm[pnt]; p != pnt; p = ptPerm[p] )
1256 ptList[len++] = INTOBJ_INT( p+1 );
1257 }
1258
1259 return list;
1260 }
1261
FuncCYCLE_PERM_INT(Obj self,Obj perm,Obj point)1262 static Obj FuncCYCLE_PERM_INT(Obj self, Obj perm, Obj point)
1263 {
1264 UInt pnt; /* value of the point */
1265
1266 /* evaluate and check the arguments */
1267 RequirePermutation("CyclePermInt", perm);
1268 pnt = GetPositiveSmallInt("CyclePermInt", point) - 1;
1269
1270 if ( TNUM_OBJ(perm) == T_PERM2 ) {
1271 return CYCLE_PERM_INT<UInt2>(perm, pnt);
1272 }
1273 else {
1274 return CYCLE_PERM_INT<UInt4>(perm, pnt);
1275 }
1276 }
1277
1278 /****************************************************************************
1279 **
1280 *F FuncCYCLE_STRUCT_PERM( <self>, <perm> ) . . . . . cycle structure of perm
1281 *
1282 ** 'FuncCYCLE_STRUCT_PERM' implements the internal function
1283 ** `CycleStructPerm'.
1284 **
1285 ** `CycleStructPerm( <perm> )'
1286 **
1287 ** 'CycleStructPerm' returns a list of the form as described under
1288 ** `CycleStructure'.
1289 */
1290 template <typename T>
CYCLE_STRUCT_PERM(Obj perm)1291 static inline Obj CYCLE_STRUCT_PERM(Obj perm)
1292 {
1293 Obj list; /* handle of the list (result) */
1294 Obj * ptList; /* pointer to the list */
1295 const T * ptPerm; /* pointer to the permutation */
1296 T * scratch;
1297 T * offset;
1298 UInt deg; /* degree of the permutation */
1299 UInt pnt; /* value of the point */
1300 UInt len; /* length of the cycle */
1301 UInt p; /* loop variable */
1302 UInt max; /* maximal cycle length */
1303 UInt cnt;
1304 UInt ende;
1305 UInt bytes;
1306 UInt1 * clr;
1307
1308 /* make sure that the buffer bag is large enough */
1309 UseTmpPerm(SIZE_OBJ(perm) + 8);
1310
1311 /* find the largest moved point */
1312 ptPerm = CONST_ADDR_PERM<T>(perm);
1313 for (deg = DEG_PERM<T>(perm); 1 <= deg; deg--) {
1314 if (ptPerm[deg - 1] != deg - 1)
1315 break;
1316 }
1317 if (deg == 0) {
1318 /* special treatment of identity */
1319 list = NEW_PLIST(T_PLIST, 0);
1320 return list;
1321 }
1322
1323 scratch = ADDR_TMP_PERM<T>();
1324
1325 /* the first deg bytes of TmpPerm hold a bit list of points done
1326 * so far. The remaining bytes will form the lengths of nontrivial
1327 * cycles (as numbers of type T). As every nontrivial cycle requires
1328 * at least 2 points, this is guaranteed to fit. */
1329 bytes = ((deg / sizeof(T)) + 1) * sizeof(T); // ensure alignment
1330 offset = (T *)((UInt)scratch + (bytes));
1331 clr = (UInt1 *)scratch;
1332
1333 /* clear out the bits */
1334 for (cnt = 0; cnt < bytes; cnt++) {
1335 clr[cnt] = (UInt1)0;
1336 }
1337
1338 cnt = 0;
1339 clr = (UInt1 *)scratch;
1340 max = 0;
1341 for (pnt = 0; pnt < deg; pnt++) {
1342 if (clr[pnt] == 0) {
1343 len = 0;
1344 clr[pnt] = 1;
1345 for (p = ptPerm[pnt]; p != pnt; p = ptPerm[p]) {
1346 clr[p] = 1;
1347 len++;
1348 }
1349
1350 if (len > 0) {
1351 offset[cnt] = (T)len;
1352 cnt++;
1353 if (len > max) {
1354 max = len;
1355 }
1356 }
1357 }
1358 }
1359
1360 ende = cnt;
1361
1362 /* create the list */
1363 list = NEW_PLIST(T_PLIST, max);
1364 SET_LEN_PLIST(list, max);
1365 ptList = ADDR_OBJ(list);
1366
1367 /* Recalculate after possible GC */
1368 scratch = ADDR_TMP_PERM<T>();
1369 offset = (T *)((UInt)scratch + (bytes));
1370
1371 for (cnt = 0; cnt < ende; cnt++) {
1372 pnt = (UInt)offset[cnt];
1373 if (ptList[pnt] == 0) {
1374 ptList[pnt] = INTOBJ_INT(1);
1375 }
1376 else {
1377 ptList[pnt] = INTOBJ_INT(INT_INTOBJ(ptList[pnt]) + 1);
1378 }
1379 }
1380
1381 return list;
1382 }
1383
FuncCYCLE_STRUCT_PERM(Obj self,Obj perm)1384 static Obj FuncCYCLE_STRUCT_PERM(Obj self, Obj perm)
1385 {
1386 /* evaluate and check the arguments */
1387 RequirePermutation("CycleStructPerm", perm);
1388
1389 if (TNUM_OBJ(perm) == T_PERM2) {
1390 return CYCLE_STRUCT_PERM<UInt2>(perm);
1391 }
1392 else {
1393 return CYCLE_STRUCT_PERM<UInt4>(perm);
1394 }
1395 }
1396
1397 /****************************************************************************
1398 **
1399 *F FuncORDER_PERM( <self>, <perm> ) . . . . . . . . . order of a permutation
1400 **
1401 ** 'FuncORDER_PERM' implements the internal function 'OrderPerm'.
1402 **
1403 ** 'OrderPerm( <perm> )'
1404 **
1405 ** 'OrderPerm' returns the order of the permutation <perm>, i.e., the
1406 ** smallest positive integer <n> such that '<perm>\^<n>' is the identity.
1407 **
1408 ** Since the largest element in S(65536) has order greater than 10^382 this
1409 ** computation may easily overflow. So we have to use arbitrary precision.
1410 */
1411 template <typename T>
ORDER_PERM(Obj perm)1412 static inline Obj ORDER_PERM(Obj perm)
1413 {
1414 const T * ptPerm; /* pointer to the permutation */
1415 Obj ord; /* order (result), may be huge */
1416 T * ptKnown; /* pointer to temporary bag */
1417 UInt len; /* length of one cycle */
1418 UInt p, q; /* loop variables */
1419
1420 /* make sure that the buffer bag is large enough */
1421 UseTmpPerm(SIZE_OBJ(perm));
1422
1423 /* get the pointer to the bags */
1424 ptPerm = CONST_ADDR_PERM<T>(perm);
1425 ptKnown = ADDR_TMP_PERM<T>();
1426
1427 /* clear the buffer bag */
1428 for ( p = 0; p < DEG_PERM<T>(perm); p++ )
1429 ptKnown[p] = 0;
1430
1431 /* start with order 1 */
1432 ord = INTOBJ_INT(1);
1433
1434 /* loop over all cycles */
1435 for ( p = 0; p < DEG_PERM<T>(perm); p++ ) {
1436
1437 /* if we haven't looked at this cycle so far */
1438 if ( ptKnown[p] == 0 && ptPerm[p] != p ) {
1439
1440 /* find the length of this cycle */
1441 len = 1;
1442 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
1443 len++; ptKnown[q] = 1;
1444 }
1445
1446 ord = LcmInt( ord, INTOBJ_INT( len ) );
1447
1448 // update bag pointers, in case a garbage collection happened
1449 ptPerm = CONST_ADDR_PERM<T>(perm);
1450 ptKnown = ADDR_TMP_PERM<T>();
1451
1452 }
1453
1454 }
1455
1456 /* return the order */
1457 return ord;
1458 }
1459
FuncORDER_PERM(Obj self,Obj perm)1460 static Obj FuncORDER_PERM(Obj self, Obj perm)
1461 {
1462 /* check arguments and extract permutation */
1463 RequirePermutation("OrderPerm", perm);
1464
1465 if ( TNUM_OBJ(perm) == T_PERM2 ) {
1466 return ORDER_PERM<UInt2>(perm);
1467 }
1468 else {
1469 return ORDER_PERM<UInt4>(perm);
1470 }
1471 }
1472
1473
1474 /****************************************************************************
1475 **
1476 *F FuncSIGN_PERM( <self>, <perm> ) . . . . . . . . . . sign of a permutation
1477 **
1478 ** 'FuncSIGN_PERM' implements the internal function 'SignPerm'.
1479 **
1480 ** 'SignPerm( <perm> )'
1481 **
1482 ** 'SignPerm' returns the sign of the permutation <perm>. The sign is +1 if
1483 ** <perm> is the product of an *even* number of transpositions, and -1 if
1484 ** <perm> is the product of an *odd* number of transpositions. The sign
1485 ** is a homomorphism from the symmetric group onto the multiplicative group
1486 ** $\{ +1, -1 \}$, the kernel of which is the alternating group.
1487 */
1488 template <typename T>
SIGN_PERM(Obj perm)1489 static inline Obj SIGN_PERM(Obj perm)
1490 {
1491 const T * ptPerm; /* pointer to the permutation */
1492 Int sign; /* sign (result) */
1493 T * ptKnown; /* pointer to temporary bag */
1494 UInt len; /* length of one cycle */
1495 UInt p, q; /* loop variables */
1496
1497 /* make sure that the buffer bag is large enough */
1498 UseTmpPerm(SIZE_OBJ(perm));
1499
1500 /* get the pointer to the bags */
1501 ptPerm = CONST_ADDR_PERM<T>(perm);
1502 ptKnown = ADDR_TMP_PERM<T>();
1503
1504 /* clear the buffer bag */
1505 for ( p = 0; p < DEG_PERM<T>(perm); p++ )
1506 ptKnown[p] = 0;
1507
1508 /* start with sign 1 */
1509 sign = 1;
1510
1511 /* loop over all cycles */
1512 for ( p = 0; p < DEG_PERM<T>(perm); p++ ) {
1513
1514 /* if we haven't looked at this cycle so far */
1515 if ( ptKnown[p] == 0 && ptPerm[p] != p ) {
1516
1517 /* find the length of this cycle */
1518 len = 1;
1519 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
1520 len++; ptKnown[q] = 1;
1521 }
1522
1523 /* if the length is even invert the sign */
1524 if ( len % 2 == 0 )
1525 sign = -sign;
1526
1527 }
1528
1529 }
1530
1531 /* return the sign */
1532 return INTOBJ_INT( sign );
1533 }
1534
FuncSIGN_PERM(Obj self,Obj perm)1535 static Obj FuncSIGN_PERM(Obj self, Obj perm)
1536 {
1537 /* check arguments and extract permutation */
1538 RequirePermutation("SignPerm", perm);
1539
1540 if ( TNUM_OBJ(perm) == T_PERM2 ) {
1541 return SIGN_PERM<UInt2>(perm);
1542 }
1543 else {
1544 return SIGN_PERM<UInt4>(perm);
1545 }
1546 }
1547
1548
1549 /****************************************************************************
1550 **
1551 *F FuncSMALLEST_GENERATOR_PERM( <self>, <perm> ) . . . . . . . . . . . . . .
1552 *F . . . . . . . smallest generator of cyclic group generated by permutation
1553 **
1554 ** 'FuncSMALLEST_GENERATOR_PERM' implements the internal function
1555 ** 'SmallestGeneratorPerm'.
1556 **
1557 ** 'SmallestGeneratorPerm( <perm> )'
1558 **
1559 ** 'SmallestGeneratorPerm' returns the smallest generator of the cyclic
1560 ** group generated by the permutation <perm>. That is the result is a
1561 ** permutation that generates the same cyclic group as <perm> and is with
1562 ** respect to the lexicographical order defined by '\<' the smallest such
1563 ** permutation.
1564 */
1565 template <typename T>
SMALLEST_GENERATOR_PERM(Obj perm)1566 static inline Obj SMALLEST_GENERATOR_PERM(Obj perm)
1567 {
1568 Obj small; /* handle of the smallest gen */
1569 T * ptSmall; /* pointer to the smallest gen */
1570 const T * ptPerm; /* pointer to the permutation */
1571 T * ptKnown; /* pointer to temporary bag */
1572 Obj ord; /* order, may be huge */
1573 Obj pow; /* power, may also be huge */
1574 UInt len; /* length of one cycle */
1575 UInt gcd, s, t; /* gcd( len, ord ), temporaries */
1576 UInt min; /* minimal element in a cycle */
1577 UInt p, q; /* loop variables */
1578 UInt l, n, x, gcd2; /* loop variable */
1579
1580 /* make sure that the buffer bag is large enough */
1581 UseTmpPerm(SIZE_OBJ(perm));
1582
1583 /* allocate the result bag */
1584 small = NEW_PERM<T>( DEG_PERM<T>(perm) );
1585
1586 /* get the pointer to the bags */
1587 ptPerm = CONST_ADDR_PERM<T>(perm);
1588 ptKnown = ADDR_TMP_PERM<T>();
1589 ptSmall = ADDR_PERM<T>(small);
1590
1591 /* clear the buffer bag */
1592 for ( p = 0; p < DEG_PERM<T>(perm); p++ )
1593 ptKnown[p] = 0;
1594
1595 /* we only know that we must raise <perm> to a power = 0 mod 1 */
1596 ord = INTOBJ_INT(1); pow = INTOBJ_INT(0);
1597
1598 /* loop over all cycles */
1599 for ( p = 0; p < DEG_PERM<T>(perm); p++ ) {
1600
1601 /* if we haven't looked at this cycle so far */
1602 if ( ptKnown[p] == 0 ) {
1603
1604 /* find the length of this cycle */
1605 len = 1;
1606 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
1607 len++; ptKnown[q] = 1;
1608 }
1609
1610 /* compute the gcd with the previously order ord */
1611 /* Note that since len is single precision, ord % len is to*/
1612 gcd = len; s = INT_INTOBJ( ModInt( ord, INTOBJ_INT(len) ) );
1613 while ( s != 0 ) {
1614 t = s; s = gcd % s; gcd = t;
1615 }
1616
1617 /* we must raise the cycle into a power = pow mod gcd */
1618 x = INT_INTOBJ( ModInt( pow, INTOBJ_INT( gcd ) ) );
1619
1620 /* find the smallest element in the cycle at such a positio*/
1621 min = DEG_PERM<T>(perm)-1;
1622 n = 0;
1623 for ( q = p, l = 0; l < len; l++ ) {
1624 gcd2 = len; s = l;
1625 while ( s != 0 ) { t = s; s = gcd2 % s; gcd2 = t; }
1626 if ( l % gcd == x && gcd2 == 1 && q <= min ) {
1627 min = q;
1628 n = l;
1629 }
1630 q = ptPerm[q];
1631 }
1632
1633 /* raise the cycle to that power and put it in the result */
1634 ptSmall[p] = min;
1635 for ( q = ptPerm[p]; q != p; q = ptPerm[q] ) {
1636 min = ptPerm[min]; ptSmall[q] = min;
1637 }
1638
1639 /* compute the new order and the new power */
1640 while ( INT_INTOBJ( ModInt( pow, INTOBJ_INT(len) ) ) != n )
1641 pow = SumInt( pow, ord );
1642 ord = ProdInt( ord, INTOBJ_INT( len / gcd ) );
1643
1644 }
1645
1646 }
1647
1648 /* return the smallest generator */
1649 return small;
1650 }
1651
FuncSMALLEST_GENERATOR_PERM(Obj self,Obj perm)1652 static Obj FuncSMALLEST_GENERATOR_PERM(Obj self, Obj perm)
1653 {
1654 /* check arguments and extract permutation */
1655 RequirePermutation("SmallestGeneratorPerm", perm);
1656
1657 if ( TNUM_OBJ(perm) == T_PERM2 ) {
1658 return SMALLEST_GENERATOR_PERM<UInt2>(perm);
1659 }
1660 else {
1661 return SMALLEST_GENERATOR_PERM<UInt4>(perm);
1662 }
1663 }
1664
1665 /****************************************************************************
1666 **
1667 *F FuncRESTRICTED_PERM( <self>, <perm>, <dom>, <test> ) . . RestrictedPerm
1668 **
1669 ** 'FuncRESTRICTED_PERM' implements the internal function
1670 ** 'RESTRICTED_PERM'.
1671 **
1672 ** 'RESTRICTED_PERM( <perm>, <dom>, <test> )'
1673 **
1674 ** 'RESTRICTED_PERM' returns the restriction of <perm> to <dom>. If <test>
1675 ** is set to `true' is is verified that <dom> is the union of cycles of
1676 ** <perm>.
1677 */
1678 template <typename T>
RESTRICTED_PERM(Obj perm,Obj dom,Obj test)1679 static inline Obj RESTRICTED_PERM(Obj perm, Obj dom, Obj test)
1680 {
1681 Obj rest;
1682 T * ptRest;
1683 const T * ptPerm;
1684 const Obj * ptDom;
1685 Int i,inc,len,p,deg;
1686
1687 /* make sure that the buffer bag is large enough */
1688 UseTmpPerm(SIZE_OBJ(perm));
1689
1690 /* allocate the result bag */
1691 deg = DEG_PERM<T>(perm);
1692 rest = NEW_PERM<T>(deg);
1693
1694 /* get the pointer to the bags */
1695 ptPerm = CONST_ADDR_PERM<T>(perm);
1696 ptRest = ADDR_PERM<T>(rest);
1697
1698 /* create identity everywhere */
1699 for ( p = 0; p < deg; p++ ) {
1700 ptRest[p]=(T)p;
1701 }
1702
1703 if ( ! IS_RANGE(dom) ) {
1704 if ( ! IS_PLIST( dom ) ) {
1705 return Fail;
1706 }
1707 /* domain is list */
1708 ptPerm = CONST_ADDR_PERM<T>(perm);
1709 ptRest = ADDR_PERM<T>(rest);
1710 ptDom = CONST_ADDR_OBJ(dom);
1711 len = LEN_LIST(dom);
1712 for (i = 1; i <= len; i++) {
1713 if (IS_POS_INTOBJ(ptDom[i])) {
1714 p = INT_INTOBJ(ptDom[i]);
1715 if (p <= deg) {
1716 p -= 1;
1717 ptRest[p] = ptPerm[p];
1718 }
1719 }
1720 else {
1721 return Fail;
1722 }
1723 }
1724 }
1725 else {
1726 len = GET_LEN_RANGE(dom);
1727 p = GET_LOW_RANGE(dom);
1728 inc = GET_INC_RANGE(dom);
1729 while (p<1) {
1730 p+=inc;
1731 len=-1;
1732 }
1733 i=p+(inc*len)-1;
1734 while (i>deg) {
1735 i-=inc;
1736 }
1737 p-=1;
1738 i-=1;
1739 while (p<=i) {
1740 ptRest[p]=ptPerm[p];
1741 p+=inc;
1742 }
1743 }
1744
1745 if (test==True) {
1746
1747 T * ptTmp = ADDR_TMP_PERM<T>();
1748
1749 /* cleanout */
1750 for ( p = 0; p < deg; p++ ) {
1751 ptTmp[p]=0;
1752 }
1753
1754 /* check whether the result is a permutation */
1755 for (p=0;p<deg;p++) {
1756 inc=ptRest[p];
1757 if (ptTmp[inc]==1) return Fail; /* point was known */
1758 else ptTmp[inc]=1; /* now point is known */
1759 }
1760
1761 }
1762
1763 /* return the restriction */
1764 return rest;
1765 }
1766
FuncRESTRICTED_PERM(Obj self,Obj perm,Obj dom,Obj test)1767 static Obj FuncRESTRICTED_PERM(Obj self, Obj perm, Obj dom, Obj test)
1768 {
1769 /* check arguments and extract permutation */
1770 RequirePermutation("RestrictedPerm", perm);
1771
1772 if ( TNUM_OBJ(perm) == T_PERM2 ) {
1773 return RESTRICTED_PERM<UInt2>(perm, dom, test);
1774 }
1775 else {
1776 return RESTRICTED_PERM<UInt4>(perm, dom, test);
1777 }
1778 }
1779
1780 /****************************************************************************
1781 **
1782 *F FuncTRIM_PERM( <self>, <perm>, <n> ) . . . . . . . . . trim a permutation
1783 **
1784 ** 'TRIM_PERM' trims a permutation to the first <n> points. This can be
1785 ## useful to save memory
1786 */
FuncTRIM_PERM(Obj self,Obj perm,Obj n)1787 static Obj FuncTRIM_PERM(Obj self, Obj perm, Obj n)
1788 {
1789 RequirePermutation("TRIM_PERM", perm);
1790 RequireNonnegativeSmallInt("TRIM_PERM", n);
1791 UInt newDeg = INT_INTOBJ(n);
1792 UInt oldDeg =
1793 (TNUM_OBJ(perm) == T_PERM2) ? DEG_PERM2(perm) : DEG_PERM4(perm);
1794 if (newDeg > oldDeg)
1795 newDeg = oldDeg;
1796 TrimPerm(perm, newDeg);
1797 return 0;
1798 }
1799
1800 /****************************************************************************
1801 **
1802 *F FuncSPLIT_PARTITION( <Ppoints>, <Qnum>,<j>,<g>,<l>)
1803 ** <l> is a list [<a>,<b>,<max>] -- needed because of large parameter number
1804 **
1805 ** This function is used in the partition backtrack to split a partition.
1806 ** The points <i> in the list Ppoints between a (start) and b (end) such
1807 ** that Qnum[i^g]=j will be moved at the end.
1808 ** At most <max> points will be moved.
1809 ** The function returns the start point of the new partition (or -1 if too
1810 ** many are moved).
1811 ** Ppoints and Qnum must be plain lists of small integers.
1812 */
1813 template <typename T>
1814 static inline Obj
SPLIT_PARTITION(Obj Ppoints,Obj Qnum,Obj jval,Obj g,Obj lst)1815 SPLIT_PARTITION(Obj Ppoints, Obj Qnum, Obj jval, Obj g, Obj lst)
1816 {
1817 Int a;
1818 Int b;
1819 Int cnt;
1820 Int max;
1821 Int blim;
1822 UInt deg;
1823 const T * ptPerm;
1824 Obj tmp;
1825
1826
1827 a=INT_INTOBJ(ELM_PLIST(lst,1))-1;
1828 b=INT_INTOBJ(ELM_PLIST(lst,2))+1;
1829 max=INT_INTOBJ(ELM_PLIST(lst,3));
1830 cnt=0;
1831 blim=b-max-1;
1832
1833 deg=DEG_PERM<T>(g);
1834 ptPerm=CONST_ADDR_PERM<T>(g);
1835 while ( (a<b)) {
1836 do {
1837 b--;
1838 if (b<blim) {
1839 /* too many points got moved out */
1840 return INTOBJ_INT(-1);
1841 }
1842 } while (ELM_PLIST(Qnum,
1843 IMAGE(INT_INTOBJ(ELM_PLIST(Ppoints,b))-1,ptPerm,deg)+1)==jval);
1844 do {
1845 a++;
1846 } while ((a<b)
1847 &&(!(ELM_PLIST(Qnum,
1848 IMAGE(INT_INTOBJ(ELM_PLIST(Ppoints,a))-1,ptPerm,deg)+1)==jval)));
1849 /* swap */
1850 if (a<b) {
1851 tmp=ELM_PLIST(Ppoints,a);
1852 SET_ELM_PLIST(Ppoints,a,ELM_PLIST(Ppoints,b));
1853 SET_ELM_PLIST(Ppoints,b,tmp);
1854 cnt++;
1855 }
1856 }
1857
1858 /* list is not necc. sorted wrt. \< (any longer) */
1859 RESET_FILT_LIST(Ppoints, FN_IS_SSORT);
1860 RESET_FILT_LIST(Ppoints, FN_IS_NSORT);
1861
1862 return INTOBJ_INT(b+1);
1863 }
1864
1865 static Obj
FuncSPLIT_PARTITION(Obj self,Obj Ppoints,Obj Qnum,Obj jval,Obj g,Obj lst)1866 FuncSPLIT_PARTITION(Obj self, Obj Ppoints, Obj Qnum, Obj jval, Obj g, Obj lst)
1867 {
1868 if (TNUM_OBJ(g)==T_PERM2) {
1869 return SPLIT_PARTITION<UInt2>(Ppoints, Qnum, jval, g, lst);
1870 }
1871 else {
1872 return SPLIT_PARTITION<UInt4>(Ppoints, Qnum, jval, g, lst);
1873 }
1874 }
1875
1876 /*****************************************************************************
1877 **
1878 *F FuncDISTANCE_PERMS( <perm1>, <perm2> )
1879 **
1880 ** 'DistancePerms' returns the number of points moved by <perm1>/<perm2>
1881 **
1882 */
1883 template <typename TL, typename TR>
DISTANCE_PERMS(Obj opL,Obj opR)1884 static inline Obj DISTANCE_PERMS(Obj opL, Obj opR)
1885 {
1886 UInt dist = 0;
1887 const TL * ptL = CONST_ADDR_PERM<TL>(opL);
1888 const TR * ptR = CONST_ADDR_PERM<TR>(opR);
1889 UInt degL = DEG_PERM<TL>(opL);
1890 UInt degR = DEG_PERM<TR>(opR);
1891 UInt i;
1892 if (degL < degR) {
1893 for (i = 0; i < degL; i++)
1894 if (ptL[i] != ptR[i])
1895 dist++;
1896 for (; i < degR; i++)
1897 if (ptR[i] != i)
1898 dist++;
1899 }
1900 else {
1901 for (i = 0; i < degR; i++)
1902 if (ptL[i] != ptR[i])
1903 dist++;
1904 for (; i < degL; i++)
1905 if (ptL[i] != i)
1906 dist++;
1907 }
1908
1909 return INTOBJ_INT(dist);
1910 }
1911
FuncDISTANCE_PERMS(Obj self,Obj opL,Obj opR)1912 static Obj FuncDISTANCE_PERMS(Obj self, Obj opL, Obj opR)
1913 {
1914 UInt type = (TNUM_OBJ(opL) == T_PERM2 ? 20 : 40) + (TNUM_OBJ(opR) == T_PERM2 ? 2 : 4);
1915 switch (type) {
1916 case 22:
1917 return DISTANCE_PERMS<UInt2, UInt2>(opL, opR);
1918 case 24:
1919 return DISTANCE_PERMS<UInt2, UInt4>(opL, opR);
1920 case 42:
1921 return DISTANCE_PERMS<UInt4, UInt2>(opL, opR);
1922 case 44:
1923 return DISTANCE_PERMS<UInt4, UInt4>(opL, opR);
1924 }
1925 return Fail;
1926 }
1927
1928 /****************************************************************************
1929 **
1930 *F FuncSMALLEST_IMG_TUP_PERM( <tup>, <perm> )
1931 **
1932 ** `SmallestImgTuplePerm' returns the smallest image of the tuple <tup>
1933 ** under the permutation <perm>.
1934 */
1935 template <typename T>
SMALLEST_IMG_TUP_PERM(Obj tup,Obj perm)1936 static inline Obj SMALLEST_IMG_TUP_PERM(Obj tup, Obj perm)
1937 {
1938 UInt res; /* handle of the image, result */
1939 const Obj * ptTup; /* pointer to the tuple */
1940 const T * ptPrm; /* pointer to the permutation */
1941 UInt tmp; /* temporary handle */
1942 UInt lmp; /* largest moved point */
1943 UInt i, k; /* loop variables */
1944
1945 res = MAX_DEG_PERM4; /* ``infty''. */
1946
1947 /* get the pointer */
1948 ptTup = CONST_ADDR_OBJ(tup) + LEN_LIST(tup);
1949 ptPrm = CONST_ADDR_PERM<T>(perm);
1950 lmp = DEG_PERM<T>(perm);
1951
1952 /* loop over the entries of the tuple */
1953 for ( i = LEN_LIST(tup); 1 <= i; i--, ptTup-- ) {
1954 k = INT_INTOBJ( *ptTup );
1955 if ( k <= lmp )
1956 tmp = ptPrm[k-1] + 1;
1957 else
1958 tmp = k;
1959 if (tmp<res) res = tmp;
1960 }
1961
1962 /* return the result */
1963 return INTOBJ_INT(res);
1964
1965 }
1966
FuncSMALLEST_IMG_TUP_PERM(Obj self,Obj tup,Obj perm)1967 static Obj FuncSMALLEST_IMG_TUP_PERM(Obj self, Obj tup, Obj perm)
1968 {
1969 if ( TNUM_OBJ(perm) == T_PERM2 ) {
1970 return SMALLEST_IMG_TUP_PERM<UInt2>(tup, perm);
1971 }
1972 else {
1973 return SMALLEST_IMG_TUP_PERM<UInt4>(tup, perm);
1974 }
1975 }
1976
1977 /****************************************************************************
1978 **
1979 *F OnTuplesPerm( <tup>, <perm> ) . . . . operations on tuples of points
1980 **
1981 ** 'OnTuplesPerm' returns the image of the tuple <tup> under the
1982 ** permutation <perm>. It is called from 'FuncOnTuples'.
1983 **
1984 ** The input <tup> must be a non-empty and dense plain list. This is is not
1985 ** verified.
1986 */
1987 template <typename T>
OnTuplesPerm_(Obj tup,Obj perm)1988 static inline Obj OnTuplesPerm_(Obj tup, Obj perm)
1989 {
1990 Obj res; /* handle of the image, result */
1991 Obj * ptRes; /* pointer to the result */
1992 const Obj * ptTup; /* pointer to the tuple */
1993 const T * ptPrm; /* pointer to the permutation */
1994 Obj tmp; /* temporary handle */
1995 UInt lmp; /* largest moved point */
1996 UInt i, k; /* loop variables */
1997
1998 GAP_ASSERT(IS_PLIST(tup));
1999 GAP_ASSERT(LEN_PLIST(tup) > 0);
2000
2001 const UInt len = LEN_PLIST(tup);
2002
2003 /* make a bag for the result and initialize pointers */
2004 res = NEW_PLIST_WITH_MUTABILITY(IS_PLIST_MUTABLE(tup), T_PLIST, len);
2005 SET_LEN_PLIST(res, len);
2006
2007 /* get the pointer */
2008 ptTup = CONST_ADDR_OBJ(tup) + len;
2009 ptRes = ADDR_OBJ(res) + len;
2010 ptPrm = CONST_ADDR_PERM<T>(perm);
2011 lmp = DEG_PERM<T>(perm);
2012
2013 /* loop over the entries of the tuple */
2014 for ( i = len; 1 <= i; i--, ptTup--, ptRes-- ) {
2015 if (IS_POS_INTOBJ(*ptTup)) {
2016 k = INT_INTOBJ( *ptTup );
2017 if (k <= lmp)
2018 tmp = INTOBJ_INT( ptPrm[k-1] + 1 );
2019 else
2020 tmp = *ptTup;
2021 *ptRes = tmp;
2022 }
2023 else {
2024 if (*ptTup == NULL) {
2025 ErrorQuit("OnTuples for perm: list must not contain holes",
2026 0L, 0L);
2027 }
2028 tmp = POW( *ptTup, perm );
2029 ptTup = CONST_ADDR_OBJ(tup) + i;
2030 ptRes = ADDR_OBJ(res) + i;
2031 ptPrm = CONST_ADDR_PERM<T>(perm);
2032 *ptRes = tmp;
2033 CHANGED_BAG( res );
2034 }
2035 }
2036
2037 return res;
2038 }
2039
OnTuplesPerm(Obj tup,Obj perm)2040 Obj OnTuplesPerm (
2041 Obj tup,
2042 Obj perm )
2043 {
2044 if ( TNUM_OBJ(perm) == T_PERM2 ) {
2045 return OnTuplesPerm_<UInt2>(tup, perm);
2046 }
2047 else {
2048 return OnTuplesPerm_<UInt4>(tup, perm);
2049 }
2050 }
2051
2052
2053 /****************************************************************************
2054 **
2055 *F OnSetsPerm( <set>, <perm> ) . . . . . . . . operations on sets of points
2056 **
2057 ** 'OnSetsPerm' returns the image of the tuple <set> under the permutation
2058 ** <perm>. It is called from 'FuncOnSets'.
2059 **
2060 ** The input <set> must be a non-empty set, i.e., plain, dense and strictly
2061 ** sorted. This is is not verified.
2062 */
2063 template <typename T>
OnSetsPerm_(Obj set,Obj perm)2064 static inline Obj OnSetsPerm_(Obj set, Obj perm)
2065 {
2066 Obj res; /* handle of the image, result */
2067 Obj * ptRes; /* pointer to the result */
2068 const Obj * ptTup; /* pointer to the tuple */
2069 const T * ptPrm; /* pointer to the permutation */
2070 Obj tmp; /* temporary handle */
2071 UInt lmp; /* largest moved point */
2072 UInt isint; /* <set> only holds integers */
2073 UInt i, k; /* loop variables */
2074
2075 GAP_ASSERT(IS_PLIST(set));
2076 GAP_ASSERT(LEN_PLIST(set) > 0);
2077
2078 const UInt len = LEN_PLIST(set);
2079
2080 /* make a bag for the result and initialize pointers */
2081 res = NEW_PLIST_WITH_MUTABILITY(IS_PLIST_MUTABLE(set), T_PLIST, len);
2082 SET_LEN_PLIST(res, len);
2083
2084 /* get the pointer */
2085 ptTup = CONST_ADDR_OBJ(set) + len;
2086 ptRes = ADDR_OBJ(res) + len;
2087 ptPrm = CONST_ADDR_PERM<T>(perm);
2088 lmp = DEG_PERM<T>(perm);
2089
2090 /* loop over the entries of the tuple */
2091 isint = 1;
2092 for ( i = len; 1 <= i; i--, ptTup--, ptRes-- ) {
2093 if (IS_POS_INTOBJ(*ptTup)) {
2094 k = INT_INTOBJ( *ptTup );
2095 if (k <= lmp)
2096 tmp = INTOBJ_INT( ptPrm[k-1] + 1 );
2097 else
2098 tmp = INTOBJ_INT( k );
2099 *ptRes = tmp;
2100 }
2101 else {
2102 isint = 0;
2103 tmp = POW( *ptTup, perm );
2104 ptTup = CONST_ADDR_OBJ(set) + i;
2105 ptRes = ADDR_OBJ(res) + i;
2106 ptPrm = CONST_ADDR_PERM<T>(perm);
2107 *ptRes = tmp;
2108 CHANGED_BAG( res );
2109 }
2110 }
2111
2112 // sort the result
2113 if (isint) {
2114 SortPlistByRawObj(res);
2115 RetypeBagSM(res, T_PLIST_CYC_SSORT);
2116 }
2117 else {
2118 SortDensePlist(res);
2119 }
2120
2121 /* return the result */
2122 return res;
2123 }
2124
OnSetsPerm(Obj set,Obj perm)2125 Obj OnSetsPerm (
2126 Obj set,
2127 Obj perm )
2128 {
2129 if ( TNUM_OBJ(perm) == T_PERM2 ) {
2130 return OnSetsPerm_<UInt2>(set, perm);
2131 }
2132 else {
2133 return OnSetsPerm_<UInt4>(set, perm);
2134 }
2135 }
2136
2137
2138 /****************************************************************************
2139 **
2140 *F SavePerm2( <perm2> )
2141 **
2142 */
SavePerm2(Obj perm)2143 static void SavePerm2(Obj perm)
2144 {
2145 SaveSubObj(STOREDINV_PERM(perm));
2146 UInt len = DEG_PERM2(perm);
2147 const UInt2 *ptr = CONST_ADDR_PERM2(perm);
2148 for (UInt i = 0; i < len; i++)
2149 SaveUInt2( *ptr++);
2150 }
2151
2152 /****************************************************************************
2153 **
2154 *F SavePerm4( <perm4> )
2155 **
2156 */
SavePerm4(Obj perm)2157 static void SavePerm4(Obj perm)
2158 {
2159 SaveSubObj(STOREDINV_PERM(perm));
2160 UInt len = DEG_PERM4(perm);
2161 const UInt4 *ptr = CONST_ADDR_PERM4(perm);
2162 for (UInt i = 0; i < len; i++)
2163 SaveUInt4( *ptr++);
2164 }
2165
2166 /****************************************************************************
2167 **
2168 *F LoadPerm2( <perm2> )
2169 **
2170 */
LoadPerm2(Obj perm)2171 static void LoadPerm2(Obj perm)
2172 {
2173 ADDR_OBJ(perm)[0] = LoadSubObj(); // stored inverse
2174 UInt len = DEG_PERM2(perm);
2175 UInt2 *ptr = ADDR_PERM2(perm);
2176 for (UInt i = 0; i < len; i++)
2177 *ptr++ = LoadUInt2();
2178 }
2179
2180 /****************************************************************************
2181 **
2182 *F LoadPerm4( <perm4> )
2183 **
2184 */
LoadPerm4(Obj perm)2185 static void LoadPerm4(Obj perm)
2186 {
2187 ADDR_OBJ(perm)[0] = LoadSubObj(); // stored inverse
2188 UInt len = DEG_PERM4(perm);
2189 UInt4 *ptr = ADDR_PERM4(perm);
2190 for (UInt i = 0; i < len; i++)
2191 *ptr++ = LoadUInt4( );
2192 }
2193
2194
2195 /****************************************************************************
2196 **
2197 *F Array2Perm( <array> ) . . . . . . . . . convert array of cycles into perm
2198 **
2199 ** This function may be used by C code generated by the GAP compiler.
2200 */
Array2Perm(Obj array)2201 Obj Array2Perm (
2202 Obj array )
2203 {
2204 Obj perm; /* permutation, result */
2205 UInt m; /* maximal entry in permutation */
2206 Obj cycle; /* one cycle of permutation */
2207 UInt i; /* loop variable */
2208
2209 /* special case for identity permutation */
2210 if ( LEN_LIST(array) == 0 ) {
2211 return IdentityPerm;
2212 }
2213
2214 /* allocate the new permutation */
2215 m = 0;
2216 perm = NEW_PERM4( 0 );
2217
2218 /* loop over the cycles */
2219 for ( i = 1; i <= LEN_LIST(array); i++ ) {
2220 cycle = ELM_LIST( array, i );
2221 RequireSmallList("Array2Perm", cycle);
2222
2223 m = ScanPermCycle(perm, m, cycle, LEN_LIST(cycle), ELM_LIST);
2224 }
2225
2226 /* if possible represent the permutation with short entries */
2227 TrimPerm(perm, m);
2228
2229 /* return the permutation */
2230 return perm;
2231 }
2232
TrimPerm(Obj perm,UInt m)2233 void TrimPerm(Obj perm, UInt m)
2234 {
2235 CLEAR_STOREDINV_PERM(perm);
2236 if (TNUM_OBJ(perm) == T_PERM2) {
2237 GAP_ASSERT(m <= DEG_PERM2(perm));
2238 ResizeBag(perm, SIZEBAG_PERM2(m));
2239 }
2240 else if (m <= 65536) {
2241 GAP_ASSERT(m <= DEG_PERM4(perm));
2242 UInt2 * ptr2 = ADDR_PERM2(perm);
2243 const UInt4 * ptr4 = CONST_ADDR_PERM4(perm);
2244 for (UInt k = 1; k <= m; k++) {
2245 ptr2[k - 1] = ptr4[k - 1];
2246 };
2247 RetypeBag(perm, T_PERM2);
2248 ResizeBag(perm, SIZEBAG_PERM2(m));
2249 }
2250 else {
2251 GAP_ASSERT(m <= DEG_PERM4(perm));
2252 ResizeBag(perm, SIZEBAG_PERM4(m));
2253 }
2254 }
2255
ScanPermCycle(Obj perm,UInt m,Obj cycle,UInt len,Obj (* readElm)(Obj,Int))2256 UInt ScanPermCycle(
2257 Obj perm, UInt m, Obj cycle, UInt len, Obj (*readElm)(Obj, Int))
2258 {
2259 UInt4 * ptr4; /* pointer into perm */
2260 Obj val; /* one entry as value */
2261 UInt c, p, l; /* entries in permutation */
2262 UInt j, k; /* loop variable */
2263
2264 GAP_ASSERT(len >= 1);
2265
2266 /* loop over the entries of the cycle */
2267 c = p = l = 0;
2268 for (j = len; 1 <= j; j--) {
2269
2270 /* get and check current entry for the cycle */
2271 val = readElm(cycle, j);
2272 c = GetPositiveSmallIntEx("Permutation", val, "<expr>");
2273 if (c > MAX_DEG_PERM4)
2274 ErrorMayQuit(
2275 "Permutation literal exceeds maximum permutation degree", 0,
2276 0);
2277
2278 /* if necessary resize the permutation */
2279 if (DEG_PERM4(perm) < c) {
2280 ResizeBag(perm, SIZEBAG_PERM4((c + 1023) / 1024 * 1024));
2281 ptr4 = ADDR_PERM4(perm);
2282 for (k = m + 1; k <= DEG_PERM4(perm); k++) {
2283 ptr4[k - 1] = k - 1;
2284 }
2285 }
2286 else {
2287 ptr4 = ADDR_PERM4(perm);
2288 }
2289 if (m < c) {
2290 m = c;
2291 }
2292
2293 /* check that the cycles are disjoint */
2294 if ((p != 0 && p == c) || (ptr4[c - 1] != c - 1)) {
2295 ErrorMayQuit(
2296 "Permutation: cycles must be disjoint and duplicate-free", 0,
2297 0);
2298 }
2299
2300 /* enter the previous entry at current location */
2301 if (p != 0) {
2302 ptr4[c - 1] = p - 1;
2303 }
2304 else {
2305 l = c;
2306 }
2307
2308 /* remember current entry for next round */
2309 p = c;
2310 }
2311
2312 /* enter first (last popped) entry at last (first popped) location */
2313 ptr4 = ADDR_PERM4(perm);
2314 if (ptr4[l - 1] != l - 1) {
2315 ErrorMayQuit(
2316 "Permutation: cycles must be disjoint and duplicate-free", 0, 0);
2317 }
2318 ptr4[l - 1] = p - 1;
2319
2320 return m;
2321 }
2322
2323
myquo(Obj pt,Obj perm)2324 static inline Int myquo(Obj pt, Obj perm)
2325 {
2326 if (TNUM_OBJ(perm) == T_PERM2)
2327 return INT_INTOBJ(QuoIntPerm<UInt2>(pt, perm));
2328 else if (TNUM_OBJ(perm) == T_PERM4)
2329 return INT_INTOBJ(QuoIntPerm<UInt4>(pt, perm));
2330 else
2331 return INT_INTOBJ(QUO(pt, perm));
2332 }
2333
2334
2335 /* Stabilizer chain helper implements AddGeneratorsExtendSchreierTree Inner loop */
FuncAGESTC(Obj self,Obj args)2336 static Obj FuncAGESTC(Obj self, Obj args)
2337 {
2338 Int i,j;
2339 Obj pt;
2340 Obj oj, lj;
2341 Int img;
2342 Obj orbit = ELM_PLIST(args,1);
2343 Obj newlabs = ELM_PLIST(args,2);
2344 Obj cycles = ELM_PLIST(args,3);
2345 Obj labels = ELM_PLIST(args,4);
2346 Obj translabels = ELM_PLIST(args,5);
2347 Obj transversal = ELM_PLIST(args, 6);
2348 Obj genlabels = ELM_PLIST(args,7);
2349 Int len = LEN_PLIST(orbit);
2350 Int len2 = len;
2351 Int lenn = LEN_PLIST(newlabs);
2352 Int lenl = LEN_PLIST(genlabels);
2353 for (i = 1; i<= len; i++) {
2354 pt = ELM_PLIST(orbit, i);
2355 for (j = 1; j <= lenn; j++) {
2356 oj = ELM_PLIST(newlabs,j);
2357 lj = ELM_PLIST(labels, INT_INTOBJ(oj));
2358 img = myquo(pt, lj);
2359 if (img <= LEN_PLIST(translabels) && (Obj)0 != ELM_PLIST(translabels,img))
2360 ASS_LIST(cycles, i, True);
2361 else {
2362 ASS_LIST(translabels, img, oj);
2363 ASS_LIST(transversal, img, lj);
2364 ASS_LIST(orbit, ++len2, INTOBJ_INT(img));
2365 ASS_LIST(cycles, len2, False);
2366 }
2367 }
2368 }
2369 while (i <= len2) {
2370 pt = ELM_PLIST(orbit, i);
2371 for (j = 1; j <= lenl; j++) {
2372 oj = ELM_PLIST(genlabels, j);
2373 lj = ELM_PLIST(labels, INT_INTOBJ(oj));
2374 img = myquo(pt, lj);
2375 if (img <= LEN_PLIST(translabels) && (Obj)0 != ELM_PLIST(translabels,img))
2376 ASS_LIST(cycles, i, True);
2377 else {
2378 ASS_LIST(translabels, img, oj);
2379 ASS_LIST(transversal, img, lj);
2380 ASS_LIST(orbit, ++len2, INTOBJ_INT(img));
2381 ASS_LIST(cycles, len2, False);
2382 }
2383 }
2384 i++;
2385 }
2386 return (Obj) 0;
2387 }
2388
2389 /* Stabilizer chain helper implements AddGeneratorsExtendSchreierTree Inner loop */
FuncAGEST(Obj self,Obj orbit,Obj newlabs,Obj labels,Obj translabels,Obj transversal,Obj genlabels)2390 static Obj FuncAGEST(Obj self,
2391 Obj orbit,
2392 Obj newlabs,
2393 Obj labels,
2394 Obj translabels,
2395 Obj transversal,
2396 Obj genlabels)
2397 {
2398 Int i,j;
2399 Int len = LEN_PLIST(orbit);
2400 Int len2 = len;
2401 Int lenn = LEN_PLIST(newlabs);
2402 Int lenl = LEN_PLIST(genlabels);
2403 Obj pt;
2404 Obj oj,lj;
2405 Int img;
2406 for (i = 1; i<= len; i++) {
2407 pt = ELM_PLIST(orbit, i);
2408 for (j = 1; j <= lenn; j++) {
2409 oj = ELM_PLIST(newlabs,j);
2410 lj = ELM_PLIST(labels, INT_INTOBJ(oj));
2411 img = myquo(pt,lj);
2412 if (img > LEN_PLIST(translabels) || (Obj)0 == ELM_PLIST(translabels,img)) {
2413 ASS_LIST(translabels, img, oj);
2414 ASS_LIST(transversal, img, lj);
2415 ASS_LIST(orbit, ++len2, INTOBJ_INT(img));
2416 }
2417 }
2418 }
2419 while (i <= len2) {
2420 pt = ELM_PLIST(orbit, i);
2421 for (j = 1; j <= lenl; j++) {
2422 oj = ELM_PLIST(genlabels, j);
2423 lj = ELM_PLIST(labels, INT_INTOBJ(oj));
2424 img = myquo(pt, lj);
2425 if (img > LEN_PLIST(translabels) || (Obj)0 == ELM_PLIST(translabels,img)) {
2426 ASS_LIST(translabels, img, oj);
2427 ASS_LIST(transversal, img, lj);
2428 ASS_LIST(orbit, ++len2, INTOBJ_INT(img));
2429 }
2430 }
2431 i++;
2432 }
2433 return (Obj) 0;
2434 }
2435
2436
2437 /****************************************************************************
2438 **
2439 *F MappingPermListList( <src>, <dst> ) . . . . . return a perm mapping src to
2440 ** dst
2441 **
2442 */
2443
2444 #define DEGREELIMITONSTACK 512
2445
FuncMappingPermListList(Obj self,Obj src,Obj dst)2446 static Obj FuncMappingPermListList(Obj self, Obj src, Obj dst)
2447 {
2448 Int l;
2449 Int i;
2450 Int d;
2451 Int next;
2452 Obj out;
2453 Obj tabdst, tabsrc;
2454 Int x;
2455 Obj obj;
2456 Int mytabs[DEGREELIMITONSTACK+1];
2457 Int mytabd[DEGREELIMITONSTACK+1];
2458
2459 RequireDenseList("AddRowVector", src);
2460 RequireDenseList("AddRowVector", dst);
2461 RequireSameLength("AddRowVector", src, dst);
2462
2463 l = LEN_LIST(src);
2464 d = 0;
2465 for (i = 1;i <= l;i++) {
2466 obj = ELM_LIST(src, i);
2467 if (!IS_POS_INTOBJ(obj)) {
2468 ErrorMayQuit("<src> must be a dense list of positive small integers", 0L, 0L);
2469 }
2470 x = INT_INTOBJ(obj);
2471 if (x > d) d = x;
2472 }
2473 for (i = 1;i <= l;i++) {
2474 obj = ELM_LIST(dst, i);
2475 if (!IS_POS_INTOBJ(obj)) {
2476 ErrorMayQuit("<dst> must be a dense list of positive small integers", 0L, 0L);
2477 }
2478 x = INT_INTOBJ(obj);
2479 if (x > d) d = x;
2480 }
2481 if (d <= DEGREELIMITONSTACK) {
2482 /* Small case where we work on the stack: */
2483 memset(&mytabs,0,sizeof(mytabs));
2484 memset(&mytabd,0,sizeof(mytabd));
2485 for (i = 1;i <= l;i++) {
2486 Int val = INT_INTOBJ(ELM_LIST(src, i));
2487 if (mytabs[val]) {
2488 // Already read where this value maps, check it is the same
2489 if (ELM_LIST(dst, mytabs[val]) != ELM_LIST(dst, i)) {
2490 return Fail;
2491 }
2492 }
2493 mytabs[val] = i;
2494 }
2495 for (i = 1;i <= l;i++) {
2496 Int val = INT_INTOBJ(ELM_LIST(dst, i));
2497 if (mytabd[val]) {
2498 // Already read where this value is mapped from, check it is
2499 // the same
2500 if (ELM_LIST(src, mytabd[val]) != ELM_LIST(src, i)) {
2501 return Fail;
2502 }
2503 }
2504 mytabd[val] = i;
2505 }
2506
2507 out = NEW_PLIST(T_PLIST_CYC,d);
2508 SET_LEN_PLIST(out,d);
2509 /* No garbage collection from here ... */
2510 next = 1;
2511 for (i = 1;i <= d;i++) {
2512 if (mytabs[i]) { /* if i is in src */
2513 SET_ELM_PLIST(out,i, ELM_LIST(dst,mytabs[i]));
2514 } else {
2515 if (mytabd[i]) {
2516 // Skip things in dst:
2517 while (mytabd[next] ||
2518 (mytabs[next] == 0 && mytabd[next] == 0))
2519 next++;
2520 SET_ELM_PLIST(out, i, INTOBJ_INT(next));
2521 next++;
2522 }
2523 else { // map elements in neither list to themselves
2524 SET_ELM_PLIST(out, i, INTOBJ_INT(i));
2525 }
2526 }
2527 }
2528 /* ... to here! No CHANGED_BAG needed since this is a new object! */
2529 } else {
2530 /* Version with intermediate objects: */
2531 tabsrc = NEW_PLIST(T_PLIST,d);
2532 /* No garbage collection from here ... */
2533 for (i = 1;i <= l;i++) {
2534 Int val = INT_INTOBJ(ELM_LIST(src, i));
2535 if (ELM_PLIST(tabsrc, val)) {
2536 if (ELM_LIST(dst, INT_INTOBJ(ELM_PLIST(tabsrc, val))) !=
2537 ELM_LIST(dst, i)) {
2538 return Fail;
2539 }
2540 }
2541 else {
2542 SET_ELM_PLIST(tabsrc, val, INTOBJ_INT(i));
2543 }
2544 }
2545 /* ... to here! No CHANGED_BAG needed since this is a new object! */
2546 tabdst = NEW_PLIST(T_PLIST,d);
2547 /* No garbage collection from here ... */
2548 for (i = 1;i <= l;i++) {
2549 int val = INT_INTOBJ(ELM_LIST(dst, i));
2550 if (ELM_PLIST(tabdst, val)) {
2551 if (ELM_LIST(src, INT_INTOBJ(ELM_PLIST(tabdst, val))) !=
2552 ELM_LIST(src, i)) {
2553 return Fail;
2554 }
2555 }
2556 else {
2557 SET_ELM_PLIST(tabdst, val, INTOBJ_INT(i));
2558 }
2559 }
2560 /* ... to here! No CHANGED_BAG needed since this is a new object! */
2561 out = NEW_PLIST(T_PLIST_CYC,d);
2562 SET_LEN_PLIST(out,d);
2563 /* No garbage collection from here ... */
2564 next = 1;
2565 for (i = 1;i <= d;i++) {
2566 if (ELM_PLIST(tabsrc,i)) { /* if i is in src */
2567 SET_ELM_PLIST(out,i,
2568 ELM_LIST(dst,INT_INTOBJ(ELM_PLIST(tabsrc,i))));
2569 } else {
2570 if (ELM_PLIST(tabdst, i)) {
2571 // Skip things in dst:
2572 while (ELM_PLIST(tabdst, next) ||
2573 (ELM_PLIST(tabdst, next) == 0 &&
2574 ELM_PLIST(tabsrc, next) == 0)) {
2575 next++;
2576 }
2577 SET_ELM_PLIST(out, i, INTOBJ_INT(next));
2578 next++;
2579 }
2580 else {
2581 SET_ELM_PLIST(out, i, INTOBJ_INT(i));
2582 }
2583 }
2584 }
2585 /* ... to here! No CHANGED_BAG needed since this is a new object! */
2586 }
2587 return FuncPermList(self,out);
2588 }
2589
2590 /* InstallGlobalFunction( SCRSift, function ( S, g ) */
2591 /* local stb, # the stabilizer of S we currently work with */
2592 /* bpt; # first point of stb.orbit */
2593
2594 /* stb := S; */
2595 /* while IsBound( stb.stabilizer ) do */
2596 /* bpt := stb.orbit[1]; */
2597 /* if IsBound( stb.transversal[bpt^g] ) then */
2598 /* while bpt <> bpt^g do */
2599 /* g := g*stb.transversal[bpt^g]; */
2600 /* od; */
2601 /* stb := stb.stabilizer; */
2602 /* else */
2603 /* #current g witnesses that input was not in S */
2604 /* return g; */
2605 /* fi; */
2606 /* od; */
2607
2608 /* return g; */
2609 /* end ); */
2610
2611
2612 static UInt RN_stabilizer = 0;
2613 static UInt RN_orbit = 0;
2614 static UInt RN_transversal = 0;
2615
2616 template <typename TG, typename Res>
SCR_SIFT_HELPER(Obj stb,Obj g,UInt nn)2617 static Obj SCR_SIFT_HELPER(Obj stb, Obj g, UInt nn)
2618 {
2619 int i;
2620 Obj result = NEW_PERM<Res>(nn);
2621 UInt dg = DEG_PERM<TG>(g);
2622
2623 if (dg > nn) /* In this case the caller has messed up or
2624 g just ends with a lot of fixed points which we can
2625 ignore */
2626 dg = nn;
2627
2628 /* Copy g into the buffer */
2629 Res * ptR = ADDR_PERM<Res>(result);
2630 const TG * ptG = CONST_ADDR_PERM<TG>(g);
2631 if (sizeof(TG) == sizeof(Res)) {
2632 memcpy(ptR, ptG, sizeof(TG) * dg);
2633 }
2634 else {
2635 for (i = 0; i < dg; i++)
2636 ptR[i] = (Res)ptG[i];
2637 }
2638 for (i = dg; i < nn; i++)
2639 ptR[i] = (Res)i;
2640
2641
2642 while (IsbPRec(stb, RN_stabilizer)) {
2643 Obj trans = ElmPRec(stb, RN_transversal);
2644 Obj orb = ElmPRec(stb, RN_orbit);
2645 Int bpt = INT_INTOBJ(ELM_LIST(orb, 1)) - 1;
2646
2647 const Res * ptrResult = CONST_ADDR_PERM<Res>(result);
2648 UInt degResult = DEG_PERM<Res>(result);
2649 Int im = (Int)(IMAGE(bpt, ptrResult, degResult));
2650 Obj t = ELM0_LIST(trans, im + 1);
2651
2652 if (!t)
2653 break;
2654
2655 while (bpt != im) {
2656
2657 ptR = ADDR_PERM<Res>(result);
2658 if (IS_PERM2(t)) {
2659 const UInt2 * ptT = CONST_ADDR_PERM2(t);
2660 UInt dt = DEG_PERM2(t);
2661 if (dt >= nn)
2662 for (i = 0; i < nn; i++)
2663 ptR[i] = ptT[ptR[i]];
2664 else
2665 for (i = 0; i < nn; i++)
2666 ptR[i] = IMAGE(ptR[i], ptT, dt);
2667 }
2668 else {
2669 const UInt4 * ptT = CONST_ADDR_PERM4(t);
2670 UInt dt = DEG_PERM4(t);
2671 if (dt >= nn)
2672 for (i = 0; i < nn; i++)
2673 ptR[i] = ptT[ptR[i]];
2674 else
2675 for (i = 0; i < nn; i++)
2676 ptR[i] = IMAGE(ptR[i], ptT, dt);
2677 }
2678 im = (Int)ptR[bpt];
2679 t = ELM_PLIST(trans, im + 1);
2680 }
2681 stb = ElmPRec(stb, RN_stabilizer);
2682 }
2683 /* so we're done sifting, and now we just have to clean up result */
2684 return result;
2685 }
2686
FuncSCR_SIFT_HELPER(Obj self,Obj stb,Obj g,Obj n)2687 static Obj FuncSCR_SIFT_HELPER(Obj self, Obj stb, Obj g, Obj n)
2688 {
2689 if (!IS_PREC(stb))
2690 RequireArgument("SCR_SIFT_HELPER", stb, "must be a plain record");
2691 RequirePermutation("SCR_SIFT_HELPER", g);
2692 UInt nn = GetSmallInt("SCR_SIFT_HELPER", n);
2693
2694 /* Setup the result, sort out which rep we are going to work in */
2695 if (nn > 65535) {
2696 if (IS_PERM2(g))
2697 return SCR_SIFT_HELPER<UInt2, UInt4>(stb, g, nn);
2698 else
2699 return SCR_SIFT_HELPER<UInt4, UInt4>(stb, g, nn);
2700 }
2701 else {
2702 if (IS_PERM2(g))
2703 return SCR_SIFT_HELPER<UInt2, UInt2>(stb, g, nn);
2704 else
2705 return SCR_SIFT_HELPER<UInt4, UInt2>(stb, g, nn);
2706 }
2707 }
2708
2709
2710 /****************************************************************************
2711 **
2712 *F * * * * * * * * * * * * * initialize module * * * * * * * * * * * * * * *
2713 */
2714
2715
2716 /****************************************************************************
2717 **
2718 *V BagNames . . . . . . . . . . . . . . . . . . . . . . . list of bag names
2719 */
2720 static StructBagNames BagNames[] = {
2721 { T_PERM2, "permutation (small)" },
2722 { T_PERM4, "permutation (large)" },
2723 { -1, "" }
2724 };
2725
2726
2727 /****************************************************************************
2728 **
2729 *V GVarFilts . . . . . . . . . . . . . . . . . . . list of filters to export
2730 */
2731 static StructGVarFilt GVarFilts[] = {
2732
2733 GVAR_FILT(IS_PERM, "obj", &IsPermFilt),
2734 { 0, 0, 0, 0, 0 }
2735
2736 };
2737
2738
2739 /****************************************************************************
2740 **
2741 *V GVarFuncs . . . . . . . . . . . . . . . . . . list of functions to export
2742 */
2743 static StructGVarFunc GVarFuncs [] = {
2744
2745 GVAR_FUNC(PermList, 1, "list"),
2746 GVAR_FUNC(LARGEST_MOVED_POINT_PERM, 1, "perm"),
2747 GVAR_FUNC(CYCLE_LENGTH_PERM_INT, 2, "perm, point"),
2748 GVAR_FUNC(CYCLE_PERM_INT, 2, "perm, point"),
2749 GVAR_FUNC(CYCLE_STRUCT_PERM, 1, "perm"),
2750 GVAR_FUNC(ORDER_PERM, 1, "perm"),
2751 GVAR_FUNC(SIGN_PERM, 1, "perm"),
2752 GVAR_FUNC(SMALLEST_GENERATOR_PERM, 1, "perm"),
2753 GVAR_FUNC(RESTRICTED_PERM, 3, "perm,domain,test"),
2754 GVAR_FUNC(TRIM_PERM, 2, "perm, degree"),
2755 GVAR_FUNC(SPLIT_PARTITION, 5, "Ppoints,Qn,j,g,a_b_max"),
2756 GVAR_FUNC(SMALLEST_IMG_TUP_PERM, 2, "tuple, perm"),
2757 GVAR_FUNC(DISTANCE_PERMS, 2, "perm1, perm2"),
2758 GVAR_FUNC(AGEST, 6, "orbit, newlabels, labels, translabels, transversal,genblabels"),
2759 GVAR_FUNC(AGESTC, -1, "orbit, newlabels, cycles, labels, translabels, transversal, genlabels"),
2760 GVAR_FUNC(MappingPermListList, 2, "src, dst"),
2761 GVAR_FUNC(SCR_SIFT_HELPER, 3, "stabrec, perm, n"),
2762 { 0, 0, 0, 0, 0 }
2763
2764 };
2765
2766
2767 /****************************************************************************
2768 **
2769 *F InitKernel( <module> ) . . . . . . . . initialise kernel data structures
2770 */
2771
InitKernel(StructInitInfo * module)2772 static Int InitKernel (
2773 StructInitInfo * module )
2774 {
2775 // set the bag type names (for error messages and debugging)
2776 InitBagNamesFromTable( BagNames );
2777
2778 /* install the marking function */
2779 InitMarkFuncBags(T_PERM2, MarkOneSubBags);
2780 InitMarkFuncBags(T_PERM4, MarkOneSubBags);
2781
2782 #ifdef HPCGAP
2783 MakeBagTypePublic( T_PERM2);
2784 MakeBagTypePublic( T_PERM4);
2785 #endif
2786
2787 ImportGVarFromLibrary("PERM_INVERSE_THRESHOLD", &PERM_INVERSE_THRESHOLD);
2788
2789 /* install the type functions */
2790 ImportGVarFromLibrary( "TYPE_PERM2", &TYPE_PERM2 );
2791 ImportGVarFromLibrary( "TYPE_PERM4", &TYPE_PERM4 );
2792
2793
2794 TypeObjFuncs[ T_PERM2 ] = TypePerm2;
2795 TypeObjFuncs[ T_PERM4 ] = TypePerm4;
2796
2797 /* init filters and functions */
2798 InitHdlrFiltsFromTable( GVarFilts );
2799 InitHdlrFuncsFromTable( GVarFuncs );
2800
2801 /* make the buffer bag */
2802 #ifndef HPCGAP
2803 InitGlobalBag( &TmpPerm, "src/permutat.cc:TmpPerm" );
2804 #endif
2805
2806 /* make the identity permutation */
2807 InitGlobalBag( &IdentityPerm, "src/permutat.cc:IdentityPerm" );
2808
2809 /* install the saving functions */
2810 SaveObjFuncs[ T_PERM2 ] = SavePerm2;
2811 SaveObjFuncs[ T_PERM4 ] = SavePerm4;
2812 LoadObjFuncs[ T_PERM2 ] = LoadPerm2;
2813 LoadObjFuncs[ T_PERM4 ] = LoadPerm4;
2814
2815 /* install the printing functions */
2816 PrintObjFuncs[ T_PERM2 ] = PrintPerm<UInt2>;
2817 PrintObjFuncs[ T_PERM4 ] = PrintPerm<UInt4>;
2818
2819 /* install the comparison methods */
2820 EqFuncs [ T_PERM2 ][ T_PERM2 ] = EqPerm<UInt2, UInt2>;
2821 EqFuncs [ T_PERM2 ][ T_PERM4 ] = EqPerm<UInt2, UInt4>;
2822 EqFuncs [ T_PERM4 ][ T_PERM2 ] = EqPerm<UInt4, UInt2>;
2823 EqFuncs [ T_PERM4 ][ T_PERM4 ] = EqPerm<UInt4, UInt4>;
2824 LtFuncs [ T_PERM2 ][ T_PERM2 ] = LtPerm<UInt2, UInt2>;
2825 LtFuncs [ T_PERM2 ][ T_PERM4 ] = LtPerm<UInt2, UInt4>;
2826 LtFuncs [ T_PERM4 ][ T_PERM2 ] = LtPerm<UInt4, UInt2>;
2827 LtFuncs [ T_PERM4 ][ T_PERM4 ] = LtPerm<UInt4, UInt4>;
2828
2829 /* install the binary operations */
2830 ProdFuncs[ T_PERM2 ][ T_PERM2 ] = ProdPerm<UInt2, UInt2>;
2831 ProdFuncs[ T_PERM2 ][ T_PERM4 ] = ProdPerm<UInt2, UInt4>;
2832 ProdFuncs[ T_PERM4 ][ T_PERM2 ] = ProdPerm<UInt4, UInt2>;
2833 ProdFuncs[ T_PERM4 ][ T_PERM4 ] = ProdPerm<UInt4, UInt4>;
2834 QuoFuncs[T_PERM2][T_PERM2] = QuoPerm;
2835 QuoFuncs[T_PERM2][T_PERM4] = QuoPerm;
2836 QuoFuncs[T_PERM4][T_PERM2] = QuoPerm;
2837 QuoFuncs[T_PERM4][T_PERM4] = QuoPerm;
2838 LQuoFuncs[ T_PERM2 ][ T_PERM2 ] = LQuoPerm<UInt2, UInt2>;
2839 LQuoFuncs[ T_PERM2 ][ T_PERM4 ] = LQuoPerm<UInt2, UInt4>;
2840 LQuoFuncs[ T_PERM4 ][ T_PERM2 ] = LQuoPerm<UInt4, UInt2>;
2841 LQuoFuncs[ T_PERM4 ][ T_PERM4 ] = LQuoPerm<UInt4, UInt4>;
2842 PowFuncs [ T_PERM2 ][ T_INT ] = PowPermInt<UInt2>;
2843 PowFuncs [ T_PERM2 ][ T_INTPOS ] = PowPermInt<UInt2>;
2844 PowFuncs [ T_PERM2 ][ T_INTNEG ] = PowPermInt<UInt2>;
2845 PowFuncs [ T_PERM4 ][ T_INT ] = PowPermInt<UInt4>;
2846 PowFuncs [ T_PERM4 ][ T_INTPOS ] = PowPermInt<UInt4>;
2847 PowFuncs [ T_PERM4 ][ T_INTNEG ] = PowPermInt<UInt4>;
2848 PowFuncs [ T_INT ][ T_PERM2 ] = PowIntPerm<UInt2>;
2849 PowFuncs [ T_INTPOS ][ T_PERM2 ] = PowIntPerm<UInt2>;
2850 PowFuncs [ T_INT ][ T_PERM4 ] = PowIntPerm<UInt4>;
2851 PowFuncs [ T_INTPOS ][ T_PERM4 ] = PowIntPerm<UInt4>;
2852 QuoFuncs [ T_INT ][ T_PERM2 ] = QuoIntPerm<UInt2>;
2853 QuoFuncs [ T_INTPOS ][ T_PERM2 ] = QuoIntPerm<UInt2>;
2854 QuoFuncs [ T_INT ][ T_PERM4 ] = QuoIntPerm<UInt4>;
2855 QuoFuncs [ T_INTPOS ][ T_PERM4 ] = QuoIntPerm<UInt4>;
2856 PowFuncs [ T_PERM2 ][ T_PERM2 ] = PowPerm<UInt2, UInt2>;
2857 PowFuncs [ T_PERM2 ][ T_PERM4 ] = PowPerm<UInt2, UInt4>;
2858 PowFuncs [ T_PERM4 ][ T_PERM2 ] = PowPerm<UInt4, UInt2>;
2859 PowFuncs [ T_PERM4 ][ T_PERM4 ] = PowPerm<UInt4, UInt4>;
2860 CommFuncs[ T_PERM2 ][ T_PERM2 ] = CommPerm<UInt2, UInt2>;
2861 CommFuncs[ T_PERM2 ][ T_PERM4 ] = CommPerm<UInt2, UInt4>;
2862 CommFuncs[ T_PERM4 ][ T_PERM2 ] = CommPerm<UInt4, UInt2>;
2863 CommFuncs[ T_PERM4 ][ T_PERM4 ] = CommPerm<UInt4, UInt4>;
2864
2865 /* install the 'ONE' function for permutations */
2866 OneFuncs[ T_PERM2 ] = OnePerm;
2867 OneFuncs[ T_PERM4 ] = OnePerm;
2868 OneMutFuncs[ T_PERM2 ] = OnePerm;
2869 OneMutFuncs[ T_PERM4 ] = OnePerm;
2870
2871 /* install the 'INV' function for permutations */
2872 InvFuncs[ T_PERM2 ] = InvPerm<UInt2>;
2873 InvFuncs[ T_PERM4 ] = InvPerm<UInt4>;
2874 InvMutFuncs[ T_PERM2 ] = InvPerm<UInt2>;
2875 InvMutFuncs[ T_PERM4 ] = InvPerm<UInt4>;
2876
2877 /* return success */
2878 return 0;
2879 }
2880
2881
2882 /****************************************************************************
2883 **
2884 *F PostRestore( <module> ) . . . . . . . . . . . . . after restore workspace
2885 */
PostRestore(StructInitInfo * module)2886 static Int PostRestore(StructInitInfo * module)
2887 {
2888 RN_stabilizer = RNamName("stabilizer");
2889 RN_orbit = RNamName("orbit");
2890 RN_transversal = RNamName("transversal");
2891
2892 /* return success */
2893 return 0;
2894 }
2895
2896
2897 /****************************************************************************
2898 **
2899 *F InitLibrary( <module> ) . . . . . . . initialise library data structures
2900 */
InitLibrary(StructInitInfo * module)2901 static Int InitLibrary (
2902 StructInitInfo * module )
2903 {
2904 /* init filters and functions */
2905 InitGVarFiltsFromTable( GVarFilts );
2906 InitGVarFuncsFromTable( GVarFuncs );
2907
2908 /* make the identity permutation */
2909 IdentityPerm = NEW_PERM2(0);
2910
2911 /* return success */
2912 return PostRestore(module);
2913 }
2914
2915
InitModuleState(void)2916 static Int InitModuleState(void)
2917 {
2918 /* make the buffer bag */
2919 TmpPerm = 0;
2920
2921 // return success
2922 return 0;
2923 }
2924
2925
2926 /****************************************************************************
2927 **
2928 *F InitInfoPermutat() . . . . . . . . . . . . . . . table of init functions
2929 */
2930 static StructInitInfo module = {
2931 /* type = */ MODULE_BUILTIN,
2932 /* name = */ "permutat",
2933 /* revision_c = */ 0,
2934 /* revision_h = */ 0,
2935 /* version = */ 0,
2936 /* crc = */ 0,
2937 /* initKernel = */ InitKernel,
2938 /* initLibrary = */ InitLibrary,
2939 /* checkInit = */ 0,
2940 /* preSave = */ 0,
2941 /* postSave = */ 0,
2942 /* postRestore = */ PostRestore,
2943 /* moduleStateSize = */ sizeof(PermutatModuleState),
2944 /* moduleStateOffsetPtr = */ &PermutatStateOffset,
2945 /* initModuleState = */ InitModuleState,
2946 /* destroyModuleState = */ 0,
2947 };
2948
InitInfoPermutat(void)2949 StructInitInfo * InitInfoPermutat ( void )
2950 {
2951 return &module;
2952 }
2953