1 /* $Id: perm.c 7821 2006-04-04 16:50:47Z kb $
2
3 Copyright (C) 2000-2003 The PARI group.
4
5 This file is part of the PARI/GP package.
6
7 PARI/GP is free software; you can redistribute it and/or modify it under the
8 terms of the GNU General Public License as published by the Free Software
9 Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10 ANY WARRANTY WHATSOEVER.
11
12 Check the License for details. You should have received a copy of it, along
13 with the package; see the file 'COPYING'. If not, write to the Free Software
14 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15
16 #include "pari.h"
17 #include "paripriv.h"
18
19 /*************************************************************************/
20 /** **/
21 /** Routines for handling VEC/COL **/
22 /** **/
23 /*************************************************************************/
24
25 /* The functions below are shallow*/
26
27 GEN
vec_setconst(GEN v,GEN x)28 vec_setconst(GEN v, GEN x)
29 {
30 long i, l = lg(v);
31 for (i = 1; i < l; i++) gel(v,i) = x;
32 return v;
33 }
34
35 GEN
const_vec(long n,GEN x)36 const_vec(long n, GEN x)
37 {
38 GEN v = cgetg(n+1, t_VEC);
39 long i;
40 for (i = 1; i <= n; i++) gel(v,i) = x;
41 return v;
42 }
43
44 int
vec_isconst(GEN v)45 vec_isconst(GEN v)
46 {
47 long i, l=lg(v);
48 if (l==1) return 1;
49 for(i=2;i<l;i++)
50 if (!gequal(gel(v,i), gel(v,1)))
51 return 0;
52 return 1;
53 }
54
55 /* Check if all the elements of v are different.
56 * Use a quadratic algorithm.
57 * It could be done in n*log(n) by sorting.
58 */
59
60 int
vec_is1to1(GEN v)61 vec_is1to1(GEN v)
62 {
63 long i,j;
64 long l=lg(v);
65 for (i=1; i<l; i++)
66 for(j=i+1; j<l; j++)
67 if (gequal(gel(v,i), gel(v,j)))
68 return 0;
69 return 1;
70 }
71
72 GEN
const_col(long n,GEN x)73 const_col(long n, GEN x)
74 {
75 GEN v = cgetg(n+1, t_COL);
76 long i;
77 for (i = 1; i <= n; i++) gel(v,i) = x;
78 return v;
79 }
80
81 GEN
vec_shorten(GEN v,long n)82 vec_shorten(GEN v, long n)
83 {
84 long i;
85 GEN V = cgetg(n+1,t_VEC);
86 for(i=1;i<=n;i++) V[i] = v[i];
87 return V;
88 }
89
90 GEN
vec_lengthen(GEN v,long n)91 vec_lengthen(GEN v, long n)
92 {
93 long i;
94 long l=lg(v);
95 GEN V = cgetg(n+1,t_VEC);
96 for(i=1;i<l;i++) V[i] = v[i];
97 return V;
98 }
99
100 GEN
vec_to_vecsmall(GEN z)101 vec_to_vecsmall(GEN z)
102 {
103 long i, l = lg(z);
104 GEN x = cgetg(l, t_VECSMALL);
105 for (i=1; i<l; i++) x[i] = itos(gel(z,i));
106 return x;
107 }
108
109 /*************************************************************************/
110 /** **/
111 /** Routines for handling VECSMALL **/
112 /** **/
113 /*************************************************************************/
114
115 GEN
vecsmall_to_vec(GEN z)116 vecsmall_to_vec(GEN z)
117 {
118 long i, l = lg(z);
119 GEN x = cgetg(l,t_VEC);
120 for (i=1; i<l; i++) gel(x,i) = stoi(z[i]);
121 return x;
122 }
123
124 GEN
vecsmall_to_col(GEN z)125 vecsmall_to_col(GEN z)
126 {
127 long i, l = lg(z);
128 GEN x = cgetg(l,t_COL);
129 for (i=1; i<l; i++) gel(x,i) = stoi(z[i]);
130 return x;
131 }
132
133 GEN
vecsmall_copy(GEN x)134 vecsmall_copy(GEN x)
135 {
136 long i, l = lg(x);
137 GEN z = cgetg(l, t_VECSMALL);
138 if (typ(x)!=t_VECSMALL) pari_err(typeer,"vecsmall_copy");
139 for (i=1; i<l; i++) z[i] = x[i];
140 return z;
141 }
142
143 GEN
const_vecsmall(long n,long c)144 const_vecsmall(long n, long c)
145 {
146 long i;
147 GEN V = cgetg(n+1,t_VECSMALL);
148 for(i=1;i<=n;i++) V[i] = c;
149 return V;
150 }
151
152 GEN
vecsmall_shorten(GEN v,long n)153 vecsmall_shorten(GEN v, long n)
154 {
155 long i;
156 GEN V = cgetg(n+1,t_VECSMALL);
157 for(i=1;i<=n;i++) V[i] = v[i];
158 return V;
159
160 }
161
162 GEN
vecsmall_lengthen(GEN v,long n)163 vecsmall_lengthen(GEN v, long n)
164 {
165 long i;
166 long l=lg(v);
167 GEN V = cgetg(n+1,t_VECSMALL);
168 for(i=1;i<l;i++) V[i] = v[i];
169 return V;
170 }
171
172 /* Sort v[0]...v[n-1] and put result in w[0]...w[n-1].
173 * We accept v==w. w must be allocated.
174 */
175
176 static void
vecsmall_sortspec(GEN v,long n,GEN w)177 vecsmall_sortspec(GEN v, long n, GEN w)
178 {
179 pari_sp ltop=avma;
180 long nx=n>>1, ny=n-nx;
181 long m, ix, iy;
182 GEN x, y;
183 if (n<=2)
184 {
185 if (n==1)
186 w[0]=v[0];
187 else if (n==2)
188 {
189 long v0=v[0], v1=v[1];
190 if (v0<=v1) { w[0]=v0; w[1]=v1; }
191 else { w[0]=v1; w[1]=v0; }
192 }
193 return;
194 }
195 x=new_chunk(nx); y=new_chunk(ny);
196 vecsmall_sortspec(v,nx,x);
197 vecsmall_sortspec(v+nx,ny,y);
198 for (m=0, ix=0, iy=0; ix<nx && iy<ny; )
199 if (x[ix]<=y[iy])
200 w[m++]=x[ix++];
201 else
202 w[m++]=y[iy++];
203 for(;ix<nx;)
204 w[m++]=x[ix++];
205 for(;iy<ny;)
206 w[m++]=y[iy++];
207 avma=ltop;
208 }
209
210 /*in place sort.*/
211 void
vecsmall_sort(GEN V)212 vecsmall_sort(GEN V)
213 {
214 long l=lg(V)-1;
215 if (l<=1) return;
216 vecsmall_sortspec(V+1,l,V+1);
217 }
218
219 static GEN
vecsmall_indexsortspec(GEN v,long n)220 vecsmall_indexsortspec(GEN v, long n)
221 {
222 long nx=n>>1, ny=n-nx;
223 long m, ix, iy;
224 GEN x, y;
225 GEN w=cgetg(n+1,t_VECSMALL);
226 if (n<=2)
227 {
228 if (n==1)
229 w[1]=1;
230 else if (n==2)
231 {
232 if (v[1]<=v[2]) { w[1]=1; w[2]=2; }
233 else { w[1]=2; w[2]=1; }
234 }
235 return w;
236 }
237 x=vecsmall_indexsortspec(v,nx);
238 y=vecsmall_indexsortspec(v+nx,ny);
239 for (m=1, ix=1, iy=1; ix<=nx && iy<=ny; )
240 if (v[x[ix]] <= v[y[iy]+nx])
241 w[m++] = x[ix++];
242 else
243 w[m++] = y[iy++]+nx;
244 for(;ix<=nx;) w[m++]=x[ix++];
245 for(;iy<=ny;) w[m++]=y[iy++]+nx;
246 avma = (pari_sp) w;
247 return w;
248 }
249
250 /*indirect sort.*/
251 GEN
vecsmall_indexsort(GEN V)252 vecsmall_indexsort(GEN V)
253 {
254 long l=lg(V)-1;
255 if (l==0) return V;
256 return vecsmall_indexsortspec(V,l);
257 }
258
259 /* assume V sorted */
260 GEN
vecsmall_uniq(GEN V)261 vecsmall_uniq(GEN V)
262 {
263 GEN W;
264 long i,j, l = lg(V);
265 if (l == 1) return vecsmall_copy(V);
266 W = cgetg(l,t_VECSMALL);
267 W[1] = V[1];
268 for(i=j=2; i<l; i++)
269 if (V[i] != W[j-1]) W[j++] = V[i];
270 stackdummy((pari_sp)(W + l), (pari_sp)(W + j));
271 setlg(W, j); return W;
272 }
273
274 int
vecsmall_lexcmp(GEN x,GEN y)275 vecsmall_lexcmp(GEN x, GEN y)
276 {
277 long lx,ly,l,i;
278 lx = lg(x);
279 ly = lg(y); l = min(lx,ly);
280 for (i=1; i<l; i++)
281 if (x[i] != y[i]) return x[i]<y[i]? -1: 1;
282 if (lx == ly) return 0;
283 return (lx < ly)? -1 : 1;
284 }
285
286 int
vecsmall_prefixcmp(GEN x,GEN y)287 vecsmall_prefixcmp(GEN x, GEN y)
288 {
289 long i, lx = lg(x), ly = lg(y), l = min(lx,ly);
290 for (i=1; i<l; i++)
291 if (x[i] != y[i]) return x[i]<y[i]? -1: 1;
292 return 0;
293 }
294
295 /*Can be used on t_VEC, but coeffs not gcopy-ed*/
296 GEN
vecsmall_prepend(GEN V,long s)297 vecsmall_prepend(GEN V, long s)
298 {
299 long i, l2 = lg(V);
300 GEN res = cgetg(l2+1, typ(V));
301 res[1] = s;
302 for (i = 2; i <= l2; ++i) res[i] = V[i - 1];
303 return res;
304 }
305
306 /*Can be used on t_VEC, but coeffs not gcopy-ed*/
307 GEN
vecsmall_append(GEN V,long s)308 vecsmall_append(GEN V, long s)
309 {
310 long i, l2 = lg(V);
311 GEN res = cgetg(l2+1, typ(V));
312 for (i = 1; i < l2; ++i) res[i] = V[i];
313 res[l2] = s; return res;
314 }
315
316 GEN
vecsmall_concat(GEN u,GEN v)317 vecsmall_concat(GEN u, GEN v)
318 {
319 long i, l1 = lg(u)-1, l2 = lg(v)-1;
320 GEN res = cgetg(l1+l2+1, t_VECSMALL);
321 for (i = 1; i <= l1; ++i) res[i] = u[i];
322 for (i = 1; i <= l2; ++i) res[i+l1] = v[i];
323 return res;
324 }
325
326 /* return the number of indices where u and v are equal */
327 long
vecsmall_coincidence(GEN u,GEN v)328 vecsmall_coincidence(GEN u, GEN v)
329 {
330 long i, s = 0, l = min(lg(u),lg(v));
331 for(i=1; i<l; i++)
332 if(u[i] == v[i]) s++;
333 return s;
334 }
335
336 /* returns the first index i<=n such that x=v[i] if it exits, 0 otherwise */
337 long
vecsmall_isin(GEN v,long x)338 vecsmall_isin(GEN v, long x)
339 {
340 long i, l = lg(v);
341 for (i = 1; i < l; i++)
342 if (v[i] == x) return i;
343 return 0;
344 }
345
346
347 long
vecsmall_pack(GEN V,long base,long mod)348 vecsmall_pack(GEN V, long base, long mod)
349 {
350 long i, s = 0;
351 for(i=1; i<lg(V); i++) s = (base*s + V[i]) % mod;
352 return s;
353 }
354
355 /*************************************************************************/
356 /** **/
357 /** Routines for handling bit vector **/
358 /** **/
359 /*************************************************************************/
360
361 GEN
bitvec_alloc(long n)362 bitvec_alloc(long n)
363 {
364 long l = 1 + (n>>TWOPOTBITS_IN_LONG);
365 return const_vecsmall(l,0);
366 }
367
368
369 GEN
bitvec_shorten(GEN bitvec,long n)370 bitvec_shorten(GEN bitvec, long n)
371 {
372 long l = 1 + (n>>TWOPOTBITS_IN_LONG);
373 return vecsmall_shorten(bitvec,l);
374 }
375
376 long
bitvec_test(GEN bitvec,long b)377 bitvec_test(GEN bitvec, long b)
378 {
379 long q = b>>TWOPOTBITS_IN_LONG;
380 long r = b&(BITS_IN_LONG-1);
381 return (bitvec[1+q]>>r) & 1L;
382 }
383
384 void
bitvec_set(GEN bitvec,long b)385 bitvec_set(GEN bitvec, long b)
386 {
387 long q = b>>TWOPOTBITS_IN_LONG;
388 long r = b&(BITS_IN_LONG-1);
389 bitvec[1+q] |= 1L<<r;
390 }
391
392 void
bitvec_clear(GEN bitvec,long b)393 bitvec_clear(GEN bitvec, long b)
394 {
395 long q = b>>TWOPOTBITS_IN_LONG;
396 long r = b&(BITS_IN_LONG-1);
397 bitvec[1+q] &= ~(1L<<r);
398 }
399
400 /*************************************************************************/
401 /** **/
402 /** Routines for handling vectors of VECSMALL **/
403 /** **/
404 /*************************************************************************/
405
406 GEN
vecvecsmall_sort(GEN x)407 vecvecsmall_sort(GEN x)
408 {
409 return gen_sort(x, 0, vecsmall_lexcmp);
410 }
411
412 GEN
vecvecsmall_indexsort(GEN x)413 vecvecsmall_indexsort(GEN x)
414 {
415 return gen_sort(x, cmp_IND | cmp_C, vecsmall_lexcmp);
416 }
417
418 long
vecvecsmall_search(GEN x,GEN y,long flag)419 vecvecsmall_search(GEN x, GEN y, long flag)
420 {
421 return gen_search(x,y,flag,vecsmall_prefixcmp);
422 }
423
424 /*************************************************************************/
425 /** **/
426 /** Routines for handling permutations **/
427 /** **/
428 /*************************************************************************/
429
430 /* Permutations may be given by
431 * perm (VECSMALL): a bijection from 1...n to 1...n i-->perm[i]
432 * cyc (VEC of VECSMALL): a product of disjoint cycles. */
433
434 /* Identity permutation. Not a good name : l is not a perm */
435 GEN
perm_identity(long l)436 perm_identity(long l)
437 {
438 GEN perm = cgetg(l + 1, t_VECSMALL);
439 long i;
440 for (i = 1; i <= l; i++) perm[i] = i;
441 return perm;
442 }
443
444 GEN
cyclicperm(long l,long d)445 cyclicperm(long l, long d)
446 {
447 GEN perm = cgetg(l + 1, t_VECSMALL);
448 long i;
449 for (i = 1; i <= l-d; i++) perm[i] = i+d;
450 for (i = l-d+1; i <= l; i++) perm[i] = i-l+d;
451 return perm;
452 }
453
454 /* Multiply (compose) two permutations.
455 * Can be used if s is a t_VEC but no copy */
456 GEN
perm_mul(GEN s,GEN t)457 perm_mul(GEN s, GEN t)
458 {
459 GEN u;
460 long i, l = lg(t);
461 u = cgetg(l, typ(s));
462 for (i = 1; i < l; i++) u[i] = s[ t[i] ];
463 return u;
464 }
465
466 /* Multiply (compose) two permutations, putting the result in the second one. */
467 static GEN
perm_mul_inplace(GEN s,GEN t)468 perm_mul_inplace(GEN s, GEN t)
469 {
470 long i;
471 for (i = 1; i < lg(s); i++) t[i] = s[t[i]];
472 return t;
473 }
474 /* Compute the inverse (reciprocal) of a permutation. */
475 GEN
perm_inv(GEN x)476 perm_inv(GEN x)
477 {
478 long i, lx = lg(x);
479 GEN y = cgetg(lx,t_VECSMALL);
480 for (i=1; i<lx; i++) y[ x[i] ] = i;
481 return y;
482 }
483
484 /* Return s*t*s^-1 */
485 GEN
perm_conj(GEN s,GEN t)486 perm_conj(GEN s, GEN t)
487 {
488 long i, l = lg(s);
489 GEN v = cgetg(l, t_VECSMALL);
490 for (i = 1; i < l; i++) v[ s[i] ] = s[ t[i] ];
491 return v;
492 }
493
494 /* Orbits of the subgroup generated by v on {1,..,n} */
495 static GEN
vecperm_orbits_i(GEN v,long n)496 vecperm_orbits_i(GEN v, long n)
497 {
498 long mj = 1, j, k, l, m, o, p, flag;
499 GEN bit, cycle, cy;
500 cycle = cgetg(n+1, t_VEC);
501 bit = bitvec_alloc(n);
502 for (k = 1, l = 1; k <= n;)
503 {
504 for ( ; bitvec_test(bit,mj); mj++) /*empty*/;
505 cy = cgetg(n+1, t_VECSMALL);
506 m = 1;
507 k++;
508 cy[m++] = mj;
509 bitvec_set(bit, mj++);
510 do
511 {
512 flag = 0;
513 for (o = 1; o < lg(v); o++)
514 for (p = 1; p < m; p++) /* m increases! */
515 {
516 j = mael(v,o,cy[p]);
517 if (!bitvec_test(bit,j))
518 {
519 flag = 1;
520 bitvec_set(bit,j);
521 k++;
522 cy[m++] = j;
523 }
524 }
525 }
526 while (flag);
527 setlg(cy, m);
528 gel(cycle,l++) = cy;
529 }
530 setlg(cycle, l); return cycle;
531 }
532 /* memory clean version */
533 GEN
vecperm_orbits(GEN v,long n)534 vecperm_orbits(GEN v, long n)
535 {
536 pari_sp av = avma;
537 return gerepilecopy(av, vecperm_orbits_i(v, n));
538 }
539
540 /* Compute the cyclic decomposition of a permutation */
541 GEN
perm_cycles(GEN v)542 perm_cycles(GEN v)
543 {
544 pari_sp av = avma;
545 return gerepilecopy(av, vecperm_orbits_i(mkvec(v), lg(v)-1));
546 }
547
548 /* Output the order of p */
549 long
perm_order(GEN v)550 perm_order(GEN v)
551 {
552 pari_sp ltop = avma;
553 GEN c;
554 long i, d;
555 c = vecperm_orbits_i(mkvec(v), lg(v)-1);
556 for(i=1, d=1; i<lg(c); i++) d = clcm(d, lg(c[i])-1);
557 avma = ltop; return d;
558 }
559
560 GEN
cyc_pow(GEN cyc,long exp)561 cyc_pow(GEN cyc, long exp)
562 {
563 long i, j, k, l, r;
564 GEN c;
565 for (r = j = 1; j < lg(cyc); j++)
566 {
567 long n = lg(cyc[j]) - 1;
568 r += cgcd(n, exp);
569 }
570 c = cgetg(r, t_VEC);
571 for (r = j = 1; j < lg(cyc); j++)
572 {
573 GEN v = gel(cyc,j);
574 long n = lg(v) - 1, e = smodss(exp,n), g = cgcd(n, e), m = n / g;
575 for (i = 0; i < g; i++)
576 {
577 GEN p = cgetg(m+1, t_VECSMALL);
578 gel(c,r++) = p;
579 for (k = 1, l = i; k <= m; k++)
580 {
581 p[k] = v[l+1];
582 l += e; if (l >= n) l -= n;
583 }
584 }
585 }
586 return c;
587 }
588
589 /* Compute the power of a permutation given by product of cycles
590 * Ouput a perm, not a cyc */
591 GEN
cyc_pow_perm(GEN cyc,long exp)592 cyc_pow_perm(GEN cyc, long exp)
593 {
594 long e, j, k, l, n;
595 GEN p;
596 for (n = 0, j = 1; j < lg(cyc); j++) n += lg(cyc[j])-1;
597 p = cgetg(n + 1, t_VECSMALL);
598 for (j = 1; j < lg(cyc); j++)
599 {
600 GEN v = gel(cyc,j);
601 n = lg(v) - 1; e = smodss(exp, n);
602 for (k = 1, l = e; k <= n; k++)
603 {
604 p[v[k]] = v[l+1];
605 if (++l == n) l = 0;
606 }
607 }
608 return p;
609 }
610
611 /* Compute the power of a permutation.
612 * TODO: make it more clever for small exp */
613 GEN
perm_pow(GEN perm,long exp)614 perm_pow(GEN perm, long exp)
615 {
616 return cyc_pow_perm(perm_cycles(perm), exp);
617 }
618
619 GEN
perm_to_GAP(GEN p)620 perm_to_GAP(GEN p)
621 {
622 pari_sp ltop=avma;
623 GEN gap;
624 GEN x;
625 long i;
626 long nb, c=0;
627 char *s;
628 long sz;
629 long lp=lg(p)-1;
630 if (typ(p) != t_VECSMALL) pari_err(typeer, "perm_to_GAP");
631 x = perm_cycles(p);
632 sz = (long) ((bfffo(lp)+1) * L2SL10 + 1);
633 /*Dry run*/
634 for (i = 1, nb = 1; i < lg(x); ++i)
635 {
636 GEN z = gel(x,i);
637 long lz = lg(z)-1;
638 nb += 1+lz*(sz+2);
639 }
640 nb++;
641 /*Real run*/
642 gap = cgetg(nchar2nlong(nb) + 1, t_STR);
643 s = GSTR(gap);
644 for (i = 1; i < lg(x); ++i)
645 {
646 long j;
647 GEN z = gel(x,i);
648 if (lg(z) > 2)
649 {
650 s[c++] = '(';
651 for (j = 1; j < lg(z); ++j)
652 {
653 if (j > 1)
654 {
655 s[c++] = ','; s[c++] = ' ';
656 }
657 sprintf(s+c,"%ld",z[j]);
658 while(s[c++]) /* empty */;
659 c--;
660 }
661 s[c++] = ')';
662 }
663 }
664 if (!c) { s[c++]='('; s[c++]=')'; }
665 s[c++] = 0;
666 return gerepileupto(ltop,gap);
667 }
668
669 int
perm_commute(GEN p,GEN q)670 perm_commute(GEN p, GEN q)
671 {
672 pari_sp ltop = avma;
673 int r = gequal(perm_mul(p,q), perm_mul(q,p));
674 avma = ltop; return r;
675 }
676
677 /*************************************************************************/
678 /** **/
679 /** Routines for handling groups **/
680 /** **/
681 /*************************************************************************/
682
683 /* Groups are vectors [gen,orders]
684 * gen (vecvecsmall): list of generators given by permutations
685 * orders (vecsmall): relatives orders of generators. */
686
687 static GEN
trivialsubgroups(void)688 trivialsubgroups(void)
689 {
690 GEN p2 = cgetg(2, t_VEC);
691 gel(p2,1) = mkvec2(cgetg(1,t_VEC), cgetg(1,t_VECSMALL));
692 return p2;
693 }
694
695 /* Compute the orders of p modulo the group S given by a set */
696 long
perm_relorder(GEN p,GEN S)697 perm_relorder(GEN p, GEN S)
698 {
699 pari_sp ltop = avma;
700 long n = 1;
701 GEN q = p;
702 while (!vecvecsmall_search(S, q, 0)) { q = perm_mul(q, p); n++; }
703 avma = ltop; return n;
704 }
705
706 GEN
perm_generate(GEN S,GEN H,long o)707 perm_generate(GEN S, GEN H, long o)
708 {
709 long i, k, n = lg(H)-1;
710 GEN L = cgetg(1+n*o, t_VEC);
711 for(i=1; i<=n; i++)
712 gel(L,i) = vecsmall_copy(gel(H,i));
713 for(k=n+1; k <= n*o; ++k)
714 gel(L,k) = perm_mul(gel(L,k-n), S);
715 return L;
716 }
717
718 /*Return the order (cardinal) of a group */
719 long
group_order(GEN G)720 group_order(GEN G)
721 {
722 GEN ord = gel(G,2);
723 long i, card = 1;
724 for (i = 1; i < lg(ord); i++) card *= ord[i];
725 return card;
726 }
727
728 /* G being a subgroup of S_n, output n */
729 long
group_domain(GEN G)730 group_domain(GEN G)
731 {
732 if (lg(G[1]) < 2) pari_err(talker,"empty group in group_domain");
733 return lg(gmael(G,1,1)) - 1;
734 }
735
736 /*Compute the left coset of g mod G: gG*/
737 GEN
group_leftcoset(GEN G,GEN g)738 group_leftcoset(GEN G, GEN g)
739 {
740 GEN res;
741 long i, j, k;
742 GEN gen = gel(G,1);
743 GEN ord = gel(G,2);
744 res = cgetg(group_order(G)+1, t_VEC);
745 gel(res,1) = vecsmall_copy(g);
746 k = 1;
747 for (i = 1; i < lg(gen); i++)
748 {
749 long c = k * (ord[i] - 1);
750 for (j = 1; j <= c; j++)
751 gel(res,++k) = perm_mul(gel(res,j), gel(gen,i));
752 }
753 return res;
754 }
755
756 /*Compute the right coset of g mod G: Gg*/
757 GEN
group_rightcoset(GEN G,GEN g)758 group_rightcoset(GEN G, GEN g)
759 {
760 GEN res;
761 long i, j, k;
762 GEN gen = gel(G,1);
763 GEN ord = gel(G,2);
764 res = cgetg(group_order(G)+1, t_VEC);
765 gel(res,1) = vecsmall_copy(g);
766 k = 1;
767 for (i = 1; i < lg(gen); i++)
768 {
769 long c = k * (ord[i] - 1);
770 for (j = 1; j <= c; j++)
771 gel(res,++k) = perm_mul(gel(gen,i), gel(res,j));
772 }
773 return res;
774 }
775
776 /*Compute the elements of a group from the generators*/
777 /*Not stack clean!*/
778 GEN
group_elts(GEN G,long n)779 group_elts(GEN G, long n)
780 {
781 return group_leftcoset(G, perm_identity(n));
782 }
783
784 /*Cyclic group generated by g of order s*/
785 GEN
cyclicgroup(GEN g,long s)786 cyclicgroup(GEN g, long s)
787 {
788 GEN p2 = cgetg(3, t_VEC);
789 gel(p2,1) = mkvec( vecsmall_copy(g) );
790 gel(p2,2) = mkvecsmall(s); return p2;
791 }
792
793 /*Return the group generated by g1,g2 of rel orders s1,s2*/
794
795 GEN
dicyclicgroup(GEN g1,GEN g2,long s1,long s2)796 dicyclicgroup(GEN g1, GEN g2, long s1, long s2)
797 {
798 GEN H = cgetg(3, t_VEC);
799 GEN p3,p4;
800 p3 = cgetg(3, t_VEC);
801 gel(p3,1) = vecsmall_copy((GEN)g1);
802 gel(p3,2) = vecsmall_copy((GEN)g2);
803 p4 = cgetg(3,t_VECSMALL);
804 p4[1] = s1;
805 p4[2] = s2;
806 gel(H,1) = p3;
807 gel(H,2) = p4;
808 return H;
809 }
810
811 /* return the quotient map G --> G/H */
812 /*The ouput is [gen,hash]*/
813 /* gen (vecvecsmall): coset generators
814 * hash (vecvecsmall): sorted vecsmall of concat(perm,coset number) */
815 GEN
group_quotient(GEN G,GEN H)816 group_quotient(GEN G, GEN H)
817 {
818 pari_sp ltop = avma;
819 GEN p1, p2, p3;
820 long i, j, k, a = 1;
821 long n = group_domain(G), o = group_order(H);
822 GEN elt = vecvecsmall_sort(group_elts(G,n));
823 long le = lg(elt)-1;
824 GEN used = bitvec_alloc(le+1);
825 long l = le/o;
826 p2 = cgetg(l+1, t_VEC);
827 p3 = cgetg(le+1, t_VEC);
828 for (i = 1, k = 1; i <= l; ++i)
829 {
830 GEN V;
831 while(bitvec_test(used,a)) a++;
832 V = group_leftcoset(H,gel(elt,a));
833 p2[i] = V[1];
834 for(j=1;j<lg(V);j++)
835 {
836 long b=vecvecsmall_search(elt,gel(V,j),0);
837 bitvec_set(used,b);
838 }
839 for (j = 1; j <= o; j++)
840 gel(p3,k++) = vecsmall_append(gel(V,j),i);
841 }
842 p1 = cgetg(3,t_VEC);
843 gel(p1,1) = gcopy(p2);
844 gel(p1,2) = vecvecsmall_sort(p3);
845 return gerepileupto(ltop,p1);
846 }
847
848 /*Find in which coset a perm lie.*/
849 long
cosets_perm_search(GEN C,GEN p)850 cosets_perm_search(GEN C, GEN p)
851 {
852 long n = gen_search(gel(C,2),p,0,vecsmall_prefixcmp);
853 if (!n) pari_err(talker, "coset not found in cosets_perm_search");
854 return mael3(C,2,n,lg(p));
855 }
856
857 /*Compute the image of a permutation by a quotient map.*/
858 GEN
quotient_perm(GEN C,GEN p)859 quotient_perm(GEN C, GEN p)
860 {
861 long j, l2 = lg(C[1]);
862 GEN p3 = cgetg(l2, t_VECSMALL);
863 for (j = 1; j < l2; ++j)
864 p3[j] = cosets_perm_search(C, perm_mul(p, gmael(C,1,j)));
865 return p3;
866 }
867
868 /* H is a subgroup of G, C is the quotient map G --> G/H
869 *
870 * Lift a subgroup S of G/H to a subgroup of G containing H */
871 GEN
quotient_subgroup_lift(GEN C,GEN H,GEN S)872 quotient_subgroup_lift(GEN C, GEN H, GEN S)
873 {
874 long l1 = lg(H[1])-1;
875 long l2 = lg(S[1])-1, j;
876 GEN p1 = cgetg(3, t_VEC), L = cgetg(l1+l2+1, t_VEC);
877 for (j = 1; j <= l1; ++j) gel(L,j) = gmael(H,1,j);
878 for (j = 1; j <= l2; ++j) gel(L,l1+j) = gmael(C,1,mael3(S,1,j,1));
879 gel(p1,1) = L;
880 gel(p1,2) = vecsmall_concat(gel(H,2),gel(S,2));
881 return p1;
882 }
883
884 /* Let G a group and H a quotient map G --> G/H
885 * Assume H is normal, return the group G/H */
886 GEN
quotient_group(GEN C,GEN G)887 quotient_group(GEN C, GEN G)
888 {
889 pari_sp ltop = avma;
890 GEN Qgen, Qord, Qelt, Q;
891 long i,j, n = lg(C[1])-1, l = lg(G[1]);
892 Qord = cgetg(l, t_VECSMALL);
893 Qgen = cgetg(l, t_VEC);
894 Qelt = mkvec(perm_identity(n));
895 for (i = 1, j = 1; i < l; ++i)
896 {
897 gel(Qgen,j) = quotient_perm(C, gmael(G,1,i));
898 Qord[j] = perm_relorder(gel(Qgen,j), vecvecsmall_sort(Qelt));
899 if (Qord[j] != 1)
900 {
901 Qelt = perm_generate(gel(Qgen,j), Qelt, Qord[j]);
902 j++;
903 }
904 }
905 setlg(Qgen,j);
906 setlg(Qord,j); Q = mkvec2(Qgen, Qord);
907 if (group_order(Q) != n) pari_err(talker,"galoissubgroup: not a WSS group");
908 return gerepilecopy(ltop,Q);
909 }
910
911 /* Test if g normalize N*/
912 long
group_perm_normalize(GEN N,GEN g)913 group_perm_normalize(GEN N, GEN g)
914 {
915 pari_sp ltop = avma;
916 long r = gequal(vecvecsmall_sort(group_leftcoset(N, g)),
917 vecvecsmall_sort(group_rightcoset(N, g)));
918 avma = ltop; return r;
919 }
920
921 /* L is a list of subgroups, C is a coset and r a rel. order.*/
922 static GEN
liftlistsubgroups(GEN L,GEN C,long r)923 liftlistsubgroups(GEN L, GEN C, long r)
924 {
925 pari_sp ltop = avma;
926 GEN R;
927 long c = lg(C)-1;
928 long l = lg(L)-1, n = lg(C[1])-1, i, k;
929 if ( !l )
930 return cgetg(1,t_VEC);
931 R = cgetg(l*c+1, t_VEC);
932 for (i = 1, k = 1; i <= l; ++i)
933 {
934 GEN S = gel(L,i), Selt = vecvecsmall_sort(group_elts(S,n));
935 long j;
936 for (j = 1; j <= c; ++j)
937 if ((perm_relorder(gel(C,j), Selt) == r)
938 && group_perm_normalize(S, gel(C,j)))
939 {
940 GEN g = cgetg(3, t_VEC);
941 gel(g,1) = vecsmall_append(gel(S,1), C[j]);
942 gel(g,2) = vecsmall_append(gel(S,2), r);
943 gel(R,k++) = g;
944 }
945 }
946 setlg(R, k);
947 return gerepilecopy(ltop, R);
948 }
949
950 /* H is a normal subgroup, C is the quotient map G -->G/H,
951 * S is a subgroup of G/H, and G is embedded in Sym(l)
952 * Return all the subgroups K of G such that
953 * S= K mod H and K inter H={1} */
954 static GEN
liftsubgroup(GEN C,GEN H,GEN S)955 liftsubgroup(GEN C, GEN H, GEN S)
956 {
957 pari_sp ltop = avma;
958 GEN V = trivialsubgroups();
959 long n = lg(S[1]), i;
960 for (i = 1; i < n; ++i)
961 { /*loop over generators of S*/
962 GEN W = group_leftcoset(H, gmael(C, 1, mael3(S, 1, i, 1)));
963 V = liftlistsubgroups(V, W, mael(S, 2, i));
964 }
965 return gerepilecopy(ltop,V);
966 }
967
968 /* 1:A4 2:S4 0: other */
969 long
group_isA4S4(GEN G)970 group_isA4S4(GEN G)
971 {
972 GEN elt = gel(G,1);
973 GEN ord = gel(G,2);
974 long n = lg(ord);
975 if (n != 4 && n != 5) return 0;
976 if (ord[1]!=2 || ord[2]!=2 || ord[3]!=3) return 0;
977 if (perm_commute(gel(elt,1),gel(elt,3))) return 0;
978 if (n==4) return 1;
979 if (ord[4]!=2) return 0;
980 if (perm_commute(gel(elt,3),gel(elt,4))) return 0;
981 return 2;
982 }
983 /* compute all the subgroups of a group G */
984 GEN
group_subgroups(GEN G)985 group_subgroups(GEN G)
986 {
987 pari_sp ltop = avma;
988 GEN p1, H, C, Q, M, sg1, sg2, sg3, gen = gel(G,1), ord = gel(G,2);
989 long lM, i, j, n = lg(gen);
990 if (n == 1) return trivialsubgroups();
991 if (group_isA4S4(G))
992 {
993 GEN u = perm_mul(gel(gen,1),gel(gen,2));
994 H = dicyclicgroup(gel(gen,1),gel(gen,2),2,2);
995 /* sg3 is the list of subgroups intersecting only partially with H*/
996 sg3 = cgetg((n==4)?4: 10, t_VEC);
997 gel(sg3,1) = cyclicgroup(gel(gen,1), 2);
998 gel(sg3,2) = cyclicgroup(gel(gen,2), 2);
999 gel(sg3,3) = cyclicgroup(u, 2);
1000 if (n==5)
1001 {
1002 GEN s = gel(gen,1); /*s=(1,2)(3,4)*/
1003 GEN t = gel(gen,2); /*t=(1,3)(2,4)*/
1004 GEN u = gel(gen,3);
1005 GEN v = gel(gen,4), st = perm_mul(s,t), w, u2;
1006 if (gequal(perm_conj(u,s), t)) /*u=(2,3,4)*/
1007 u2 = perm_mul(u,u);
1008 else
1009 {
1010 u2 = u;
1011 u = perm_mul(u,u);
1012 }
1013 if (perm_order(v)==2)
1014 {
1015 if (!perm_commute(s,v)) /*v=(1,2)*/
1016 {
1017 v = perm_conj(u,v);
1018 if (!perm_commute(s,v)) v = perm_conj(u,v);
1019 }
1020 w = perm_mul(v,t); /*w=(1,4,2,3)*/
1021 }
1022 else
1023 {
1024 w = v;
1025 if (!gequal(perm_mul(w,w), s)) /*w=(1,4,2,3)*/
1026 {
1027 w = perm_conj(u,w);
1028 if (!gequal(perm_mul(w,w), s)) w = perm_conj(u,w);
1029 }
1030 v = perm_mul(w,t); /*v=(1,2)*/
1031 }
1032 gel(sg3,4) = dicyclicgroup(s,v,2,2);
1033 gel(sg3,5) = dicyclicgroup(t,perm_conj(u,v),2,2);
1034 gel(sg3,6) = dicyclicgroup(st,perm_conj(u2,v),2,2);
1035 gel(sg3,7) = dicyclicgroup(s,w,2,2);
1036 gel(sg3,8) = dicyclicgroup(t,perm_conj(u,w),2,2);
1037 gel(sg3,9) = dicyclicgroup(st,perm_conj(u2,w),2,2);
1038 }
1039 }
1040 else
1041 {
1042 long osig = mael(factoru(ord[1]), 1, 1);
1043 GEN sig = perm_pow(gel(gen,1), ord[1]/osig);
1044 H = cyclicgroup(sig,osig);
1045 sg3 = NULL;
1046 }
1047 C = group_quotient(G,H);
1048 Q = quotient_group(C,G);
1049 M = group_subgroups(Q);
1050 lM = lg(M);
1051 /* sg1 is the list of subgroups containing H*/
1052 sg1 = cgetg(lM, t_VEC);
1053 for (i = 1; i < lM; ++i)
1054 gel(sg1,i) = quotient_subgroup_lift(C,H,gel(M,i));
1055 /*sg2 is a list of lists of subgroups not intersecting with H*/
1056 sg2 = cgetg(lM, t_VEC);
1057 /* Loop over all subgroups of G/H */
1058 for (j = 1; j < lM; ++j)
1059 gel(sg2,j) = liftsubgroup(C, H, gel(M,j));
1060 p1 = concat(sg1, concat(sg2,NULL));
1061 if (sg3)
1062 {
1063 p1 = concat(p1, sg3);
1064 /*Fixup to avoid the D4 subgroups of S4 to be in non supersolvable format*/
1065 if (n==5)
1066 for(j = 3; j <= 5; j++)
1067 {
1068 GEN c = gmael(p1,j,1);
1069 if (!perm_commute(gel(c,1),gel(c,3)))
1070 {
1071 if (perm_commute(gel(c,2),gel(c,3))) { lswap(c[1], c[2]); }
1072 else
1073 (void)perm_mul_inplace(gel(c,2), gel(c,1));
1074 }
1075 }
1076 }
1077 return gerepileupto(ltop,p1);
1078 }
1079
1080 /*return 1 if G is abelian, else 0*/
1081 long
group_isabelian(GEN G)1082 group_isabelian(GEN G)
1083 {
1084 long i, j;
1085 GEN g = gel(G,1);
1086 for(i=2; i<lg(g); i++)
1087 for(j=1; j<i; j++)
1088 if (!perm_commute(gel(g,i), gel(g,j))) return 0;
1089 return 1;
1090 }
1091
1092 /*If G is abelian, return its HNF matrix*/
1093 GEN
group_abelianHNF(GEN G,GEN S)1094 group_abelianHNF(GEN G, GEN S)
1095 {
1096 GEN M, g = gel(G,1), o = gel(G,2);
1097 long i, j, k, n = lg(g);
1098 if (!group_isabelian(G)) return NULL;
1099 if (n==1) return cgetg(1,t_MAT);
1100 if (!S) S = group_elts(G, group_domain(G));
1101 M = cgetg(n,t_MAT);
1102 for(i=1; i<n; i++)
1103 {
1104 GEN P, C = cgetg(n,t_COL);
1105 pari_sp av = avma;
1106 gel(M,i) = C;
1107 P = perm_pow(gel(g,i), o[i]);
1108 for(j=1; j<lg(S); j++)
1109 if (gequal(P, gel(S,j))) break;
1110 avma = av;
1111 if (j==lg(S)) pari_err(talker,"wrong argument in galoisisabelian");
1112 j--;
1113 for(k=1; k<i; k++)
1114 {
1115 long q = j / o[k];
1116 gel(C,k) = stoi(j - q*o[k]);
1117 j = q;
1118 }
1119 gel(C,k) = stoi(o[i]);
1120 for (k++; k<n; k++) gel(C,k) = gen_0;
1121 }
1122 return M;
1123 }
1124
1125 /*If G is abelian, return its abstract SNF matrix*/
1126 GEN
group_abelianSNF(GEN G,GEN L)1127 group_abelianSNF(GEN G, GEN L)
1128 {
1129 pari_sp ltop = avma;
1130 GEN H = group_abelianHNF(G,L);
1131 if (!H) return NULL;
1132 return gerepileupto(ltop, smithclean( smith(H) ));
1133 }
1134
1135 GEN
abelian_group(GEN v)1136 abelian_group(GEN v)
1137 {
1138 GEN G = cgetg(3,t_VEC);
1139 long card, i, d = 1;
1140 gel(G,1) = cgetg(lg(v),t_VEC);
1141 gel(G,2) = vecsmall_copy(v);
1142 card = group_order(G);
1143 for(i=1; i<lg(v); i++)
1144 {
1145 GEN p = cgetg(card+1, t_VECSMALL);
1146 long o = v[i], u = d*(o-1), j, k, l;
1147 gmael(G,1,i) = p;
1148 /* The following loop is over-optimized. Remember that I wrote it for
1149 * testpermutation. Something has survived... BA */
1150 for(j=1;j<=card;)
1151 {
1152 for(k=1;k<o;k++)
1153 for(l=1;l<=d; l++,j++) p[j] = j+d;
1154 for (l=1; l<=d; l++,j++) p[j] = j-u;
1155 }
1156 d += u;
1157 }
1158 return G;
1159 }
1160
1161 GEN
groupelts_center(GEN S)1162 groupelts_center(GEN S)
1163 {
1164 pari_sp ltop = avma;
1165 long i, j, n = lg(S)-1, l = n;
1166 GEN V, elts = bitvec_alloc(n+1);
1167 for(i=1; i<=n; i++)
1168 {
1169 if (bitvec_test(elts,i)) { l--; continue; }
1170 for(j=1; j<=n; j++)
1171 if (!perm_commute(gel(S,i),gel(S,j)))
1172 {
1173 bitvec_set(elts,i);
1174 bitvec_set(elts,j); l--; break;
1175 }
1176 }
1177 V = cgetg(l+1,t_VEC);
1178 for (i=1, j=1; i<=n ;i++)
1179 if (!bitvec_test(elts,i)) gel(V,j++) = vecsmall_copy(gel(S,i));
1180 return gerepileupto(ltop,V);
1181 }
1182
1183 GEN
groupelts_abelian_group(GEN S)1184 groupelts_abelian_group(GEN S)
1185 {
1186 pari_sp ltop = avma;
1187 GEN Qgen, Qord, Qelt;
1188 long i, j, n = lg(S[1])-1, l = lg(S);
1189 Qord = cgetg(l, t_VECSMALL);
1190 Qgen = cgetg(l, t_VEC);
1191 Qelt = mkvec(perm_identity(n));
1192 for (i = 1, j = 1; i < l; ++i)
1193 {
1194 Qgen[j] = S[i];
1195 Qord[j] = perm_relorder(gel(Qgen,j), vecvecsmall_sort(Qelt));
1196 if (Qord[j] != 1)
1197 {
1198 Qelt = perm_generate(gel(Qgen,j), Qelt, Qord[j]);
1199 j++;
1200 }
1201 }
1202 setlg(Qgen,j);
1203 setlg(Qord,j);
1204 return gerepilecopy(ltop, mkvec2(Qgen, Qord));
1205 }
1206
1207 GEN
group_export_GAP(GEN G)1208 group_export_GAP(GEN G)
1209 {
1210 pari_sp ltop = avma;
1211 GEN s, g = gel(G,1);
1212 long i, l = lg(g);
1213 if (l == 1) return strtoGENstr("Group(())");
1214 s = strtoGENstr("Group(");
1215 for (i = 1; i < l; ++i)
1216 {
1217 if (i > 1) s = shallowconcat(s, strtoGENstr(", "));
1218 s = shallowconcat(s, perm_to_GAP(gel(g,i)));
1219 }
1220 s = concat(s, strtoGENstr(")"));
1221 return gerepileupto(ltop,s);
1222 }
1223
1224 GEN
group_export_MAGMA(GEN G)1225 group_export_MAGMA(GEN G)
1226 {
1227 pari_sp ltop=avma;
1228 GEN s, g = gel(G,1);
1229 long i, l = lg(g);
1230 if (l == 1) return strtoGENstr("PermutationGroup<1|>");
1231 s = strtoGENstr("PermutationGroup<");
1232 s = shallowconcat(s, stoi(group_domain(G)));
1233 s = shallowconcat(s, strtoGENstr("|"));
1234 for (i = 1; i < l; ++i)
1235 {
1236 if (i > 1) s = shallowconcat(s, strtoGENstr(", "));
1237 s = shallowconcat(s, vecsmall_to_vec(gel(g,i)));
1238 }
1239 s = concat(s, strtoGENstr(">"));
1240 return gerepileupto(ltop,s);
1241 }
1242
1243 GEN
group_export(GEN G,long format)1244 group_export(GEN G, long format)
1245 {
1246 switch(format)
1247 {
1248 case 0: return group_export_GAP(G);
1249 case 1: return group_export_MAGMA(G);
1250 }
1251 pari_err(flagerr,"galoisexport");
1252 return NULL; /*-Wall*/
1253 }
1254