GrIsoFacile:=function(n);
G:=Sym(n);
GEn:=Generators(G);
T:=(SetToSequence(Set(G)));
I:=[g^(-1):g in T];
L:={[l*g: g in T]:l in GEn};
ISO:=PermutationGroup<Set(G)|Include(L, I)>;
return ISO;
end function; 



GrIsoFacileOld:=function(n);
G:=Sym(n);
GEn:=Generators(G);
T:=IndexedSetToSequence(SetToIndexedSet(Set(G)));
I:=[g^(-1):g in T];
S:=Set(G);
GG:=Sym(S);
e:=GG!I;
L:={GG![l*g: g in T]:l in GEn};
ISO:=PermutationGroup<S|Include(L, e)>;
return ISO;
end function; 


DeterminIso:=function(T);
G:=Parent(Representative(GSet(Parent(T))));
I:=Id(G);
phi:=Image(T, I);
Im1:=Representative(Support(phi^-1*Image(T,G!(1,2))) meet Support(phi^-1*Image(T, G!(1,3))));
VecAlpha:=[Im1];
for i in [2..#GSet(G)] do
Append(~VecAlpha, Representative(Exclude(Support(phi^-1 *Image(T, G!(1, i))), Im1)));
end for;
alpha:=G!VecAlpha;
beta:=phi*alpha^-1;
r:=G!(1,2,3);
if  beta*r^-1*alpha eq Image(T, r) then i:=1;
else i:=0;
end if; 
return <alpha, beta,i>;
end function;


DistPerm:=function(p,q)
r:=#Support(p*q^-1);
return r;
end function;

forward EstT1;
EstT1:=procedure(d,R,~vec);
/*le vecteur $vec$ contient le nombre d'enfants de chaque sommet*/
if not(IsEmpty(R)) then 
Append(~vec,#R);
p:=Random(R);
NewR:={q: q in R |#Support(p*q^-1) ge d};
EstT1(d,NewR, ~vec);
end if;
end procedure;



EsperEstT:=function(d,R,N);
s:=0;
for k in [1..N] do
vec:=[];
EstT1(d,R,~vec);
s1:=1;
for i in [2..#vec] do
s1:=vec[#vec-i+1]*(1+s1);
end for;
s:=s+s1;
end for;
Val:=Round(s/N);
return Val;
end function;




forward ExactT;
ExactT:=procedure(d,R,~j)
if not(IsEmpty(R)) then 
for p in R do
j:=j+1;
/*la variable $j$ est incrémentée à chaque passage par un nouveau sommet de l'arbre.*/ 
NewR:={q: q in R |#Support(p*q^-1) ge d};
ExactT(d, NewR,~j);
end for;
end if;
end procedure;

/*
n:=5;
d:=5;
G:=Sym(n);
R:={p:p in G|#Support(p) ge d};
EsperEstT(d,R,10000);
j:=0;
ExactT(d,R,~j);
j;
*/


IsRbounded:=function(PA, r);
t:=true;
N:=#GSet(Parent(Representative(PA)));
for i in [1..N] do
for j in [1..N] do
t:= t and #{p: p in PA|Image(p, j) eq i} le r;
end for;
end for;
return t;
end function;

IsRbalanced:=function(PA, r);
t:=true;
N:=#GSet(Parent(Representative(PA)));
for i in [1..N] do
for j in [1..N] do
t:= t and #{p: p in PA|Image(p, j) eq i} eq r;
end for;
end for;
return t;
end function;


BaseNupleRaff:=function(G);
SETG:=Set(G);
VSET:={CycleStructure(u): u in SETG};
D:= IndexedSetToSequence(SetToIndexedSet(VSET));
return D;
end function;


GrapheColore:=function(pa, D);
paI:=SetToIndexedSet(pa);
D:=Sort(D);
s:=#pa;
V:={1..s};
SetOfEdges:=[];
for i in [1..#D] do
SetOfEdges[i]:={{k,l} : k,l in [1..s] |CycleStructure(paI[k]/paI[l]) eq D[i]};
end for;
SetOfEdges[1]:={};
return SetOfEdges;
end function;


IsIsoGrapheColore2:=function(pa1,pa2);
t:= #pa1 eq #pa2;
D:=BaseNupleRaff(Parent(Random(pa1)));
E1:=GrapheColore(pa1,D);
E2:=GrapheColore(pa2,D);
if t then
for i in [1.. #E1] do
t:=#E1[#E1+1-i] eq #E2[#E1+1-i] and t;
if (t and #E1[#E1+1-i] ne 0) then
DSS:=Graph<#pa1|E1[#E1+1-i]>;
CQP:=Graph<#pa1|E2[#E1+1-i]>;
t:=IsIsomorphic(DSS,CQP)and t;
end if;
end for;
end if;
return t;
end function;


AffichPA:=function(PA);
n:=#GSet(Parent(Representative(PA)));
s:=#PA;
M:=Matrix(Sort([[Image(p, i) : i in [1..n]] : p in PA]));
return M;
end function;


AffichPAL:=function(PA);

/* Special pour les listes*/
n:=#GSet(Parent(Representative(PA)));
s:=#PA;
M:=Matrix([[Image(p, i) : i in [1..n]] : p in PA]);
return M;
end function;


GraphGRSPListe:=function(W);
s:=#W;
n:=#GSet(Parent(Representative(W)));
V1:={1..n};/*les n lignes du tableau*/
V2:={n+1..2*n};/* les n colonnes*/
V3:={2*n+1..2*n+s};/* Les s symboles*/
V4:={s+2*n+1..s+2*n+n*n};/* Les sommets de degré 3*/
V5:={s+2*n+n*n+1..s+2*n+n*n+3};/* les 3 sommets qui lient le reste*/
E1:={{v1, s+2*n+n*n+1}: v1 in V1};/*Tous les sommets "lignes" sont adjacent au sommet s+2n+sn+1*/
E2:={{v2, s+2*n+n*n+2}: v2 in V2};/*Tous les sommets "colonnes" sont adjacent au sommet s+2n+sn+2*/
E3:={{v3, s+2*n+n*n+3}: v3 in V3};/*Tous les sommets "symboles" sont adjacent au sommet s+2n+sn+3*/
E4:={{v1, s+n+v1*n+t}: v1 in V1, t in [1..n]};
E5:={{s+n+v1*n+t, n+t}:v1 in V1, t in [1..n]};
E6:={{2*n+u, s+2*n+(t-1)*n+Image(W[u], t)}:u in [1..s], t in [1..n]};
E7:={{s+2*n+n*n+1, s+2*n+n*n+4},{s+2*n+n*n+2, s+2*n+n*n+5},{s+2*n+n*n+3, s+2*n+n*n+6},{s+2*n+n*n+7, s+2*n+n*n+6}};
/* J'ajoute, pour créer des ancres, une arête aux sommets qui lient les sommets lignes et colonnes ,  et deux arêtes au sommet liant les symboles*/
V:={1..s+2*n+n*n+7};
E:=E1 join E2 join E3 join E4 join E5 join E6 join E7;
GR:=Graph<V| E>;
return GR;
end function;



PolType:=function(PA,n);
CC:=ConjugacyClasses(Sym(n));
CS:=[CycleStructure(CC[i, 3]): i in [1..#CC]];
Pol:=[#{{p,q}: p, q in PA| CycleStructure(p*q^-1) eq CS[i]}: i in [1..#CS]];
return Pol;
end function;

PolDist:=function(PA);
n:=#GSet(Parent(Representative(PA)));
D:=[#{{p,q}: p, q in PA|#Support(p*q^-1) eq i}: i in [0..n]];
return D;
end function;

Occurrences:=function(PA,n);
OCC:={* #{p: p in PA|Image(p, i) eq j}: i, j in [1..n] *};
return OCC;
end function;


TestPAdansD:=function(PA,n,D,GG);
s:=#PA;
D1:={PB:PB in D|#PB eq s};
/* var prend la valeur true si il existe un élément de la liste D isométrique à PA*/
var:=not(IsEmpty(D1));
if not(IsEmpty(D1)) then
PType:=PolType(PA, n);
D2:={PB: PB in D1|PolType(PB,n) eq PType};
var:=not(IsEmpty(D2));
if not(IsEmpty(D2)) then
M:=Occurrences(PA,n);
D3:={PB: PB in D2|Occurrences(PB,n) eq M};
var:=not(IsEmpty(D3));
if not(IsEmpty(D3)) then
var:= exists{PB: PB in D3|exists{T: T in GG|Image(T, PB) eq PA}};
end if;
end if;
end if;
return var;
end function;

forward GenPArecord;
procedure GenPArecord(PA, Res,n,d,GG, ~D);
if TestPAdansD(PA,n,D,GG) then
else
Include(~D, PA);
#D;
for p in Res do
NewPA:=Include(PA,p);
NewRes:={q:q in Res|#Support(q*p^-1) ge d};
GenPArecord(NewPA, NewRes,n,d,GG, ~D);
end for;
end if;
end procedure;
/*
n:=5;
d:=4;
D:={* *};
GG:=GrIsoFacile(n);
PA:={Id(Sym(n))};
Res:={p: p in Sym(n)|#Support(p) ge d}; 
time GenPArecord(PA,Res,n,d,GG,~D); 
*/



ConsEGG:=function(n);
/*--Ici pour Sym(5) */
G:=Sym(n);
GG:=GrIsoFacile(n);
EGG:=<GG>;
Append(~EGG,Stabiliser(GG, Id(Sym(n))));
if n eq 5 then
Append(~EGG,<Stabiliser(EGG[2], G!(4,5)),G!(4,5) >);
Append(~EGG,<Stabiliser(EGG[2], G!(2,3)(4,5)),G!(2,3)(4,5) >);
Append(~EGG,<Stabiliser(EGG[2], G!(3,4,5)),G!(3,4,5) >);
Append(~EGG,<Stabiliser(EGG[2], G!(2,3,4,5)),G!(2,3,4,5) >);
Append(~EGG,<Stabiliser(EGG[2], G!(1,2,3,4,5)),G!(1,2,3,4,5) >);
Append(~EGG,<Stabiliser(EGG[2], G!(1,2)(3,4,5)),G!(1,2)(3,4,5) >);
else
Append(~EGG,<Stabiliser(EGG[2], G!(1,2)(3,4)(5,6)),G!(1,2)(3,4)(5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(n-1,n)),G!(n-1,n) >);
Append(~EGG,<Stabiliser(EGG[2], G!(3,4)(5,6)),G!(3,4)(5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(1,2,3)(4,5,6)),G!(1,2,3)(4,5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(4,5,6)),G!(4,5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(1,2)(3,4,5,6)),G!(1,2)(3,4,5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(3,4,5,6)),G!(3,4,5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(2,3,4,5,6)),G!(2,3,4,5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(1,2,3,4,5,6)),G!(1,2,3,4,5,6) >);
Append(~EGG,<Stabiliser(EGG[2], G!(2,3)(4,5,6)),G!(2,3)(4,5,6) >);
end if;
return EGG;
end function;



NrPenWM:=function(C);
N:=1;
for i in [1..#C] do
N:=N*C[i,1]^C[i,2]*Factorial(C[i,2]);
end for;
return N;
end function;

NrIndPenPAWM:=function(PA);
n:=#GSet(Parent(Random(PA)));
G:=Sym(n);
ECyc:=[CycleStructure(ConjugacyClasses(G)[i,3]): i in [1..#ConjugacyClasses(G)]];
EPen:=[NrPenWM(ECyc[i]): i in [1..#ConjugacyClasses(G)]];
for i in [1..#ConjugacyClasses(G)] do
if not(exists{{p,q}: p, q in PA| CycleStructure(p*q^-1) eq CycleStructure(ConjugacyClasses(G)[i,3])})
then 
EPen[i]:=Max(EPen)+1;
end if;
end for;

l:=Min([N: N in EPen| N ne 0]);
ind:=Index(EPen, l);
return ind;
end function;

ParentCanWM:=function(PA,EGG);/*joue le role de m(G) par rapport à Kaski */

n:=#GSet(Parent(Random(PA)));
G:=Sym(n);
l:=NrIndPenPAWM(PA)+1;
EnsPaires:={[p, q]:p,q in PA|CycleStructure(p*q^-1) eq CycleStructure(EGG[l, 2])};
Trip:={};
for P in EnsPaires do      
t, h:=IsConjugate(G,P[1]^-1*P[2], EGG[l,2]) ;
if t then Include(~Trip,[P[1],P[2],h]);
end if;
end for;        
EnsAcons:=[<{Image(T,P[3]^-1*P[1]^-1*g*P[3]): g in PA}, P, T>: P in Trip, T in EGG[l, 1]];
ListMat:=[AffichPA( EnsAcons[i,1]): i in [1..#EnsAcons]];
m,k:=Min(ListMat);
perm:=G![m[NumberOfRows(m),i]:i in [1..n]];
TT:=EnsAcons[k];
PAp:={TT[2,1]*TT[2,3]*Image(TT[3]^-1,g)*TT[2,3]^-1:g in Exclude(TT[1],perm)};
return PAp;
end function;

InvariantWM:=function(PA,EGG);
n:=#GSet(Parent(Representative(PA)));
G:=Sym(n);
l:=NrIndPenPAWM(PA)+1;
EnsPaires:={[p, q]:p,q in PA|CycleStructure(p*q^-1) eq CycleStructure(EGG[l, 2])};
Trip:={};
for P in EnsPaires do      
t, h:=IsConjugate(G,P[1]^-1*P[2], EGG[l,2]) ;
if t then Include(~Trip,[P[1],P[2],h]);
end if;
end for;        
EnsAcons:=[<{Image(T,P[3]^-1*P[1]^-1*g*P[3]): g in PA}, P, T>: P in Trip, T in EGG[l, 1]];
ListMat:=[AffichPA( EnsAcons[i,1]): i in [1..#EnsAcons]];
m,k:=Min(ListMat);
TT:=EnsAcons[k,1];
return TT;
end function;


forward GenOrdonneeL;
procedure GenOrdonneeL(PAdepart,EGG,n,d,~i,Res,~EnsPAMax);
/* Attention, PAdepart est une liste.*/
if IsEmpty(Res) then Include(~EnsPAMax, PAdepart);
#EnsPAMax;
else
i:=i+1;
[#Res, i];
for p in Res do
if AffichPA(InvariantWM(Append(PAdepart,p),EGG)) eq AffichPAL(Append(PAdepart,p)) then 
NewRes:={q: q in Res|DistPerm(p, q ) ge d};
GenOrdonneeL(Append(PAdepart,p),EGG,  n, d,~i,NewRes, ~EnsPAMax);
end if;
end for;
end if;
end procedure;

/*
n:=5;
d:=4;
G:=Sym(n);
PAdepart:={Id(Sym(n)), G!(2,3,4,5)};
j:=0;
r:=2;
Res:={p: p in Sym(n)|DistPerm(p, Id(Sym(n))) ge d and DistPerm(p, G!(2,3,4,5)) 
ge d};
GG:=GrIsoFacile(n);
EGG:=ConsEGG(5);
EnsPAMax:={};
time GenOrdonnee(PAdepart,EGG,n,d,~j,Res,~EnsPAMax);*/





ParentCan:=function(PA);
/*Sur base d'un PA, renvoie le parent canonique.*/
n:=#GSet(Parent(Random(PA)));
W:=SetToIndexedSet(PA);
GRSP:=GraphGRSPListe(W);
A,B:=IsIsomorphic( GRSP,CanonicalGraph(GRSP));
e,f:=Min(B([2*n+1..2*n+#PA]));
p:=W[f];
ParentPA:=Exclude(PA, p);
return ParentPA; 
end function;





forward GenPAWeak;
procedure GenPAWeak(PA, Res,n,d,GG, ~D);
/*Programme de génération utilisation une augmentation canonique faible.*/
if IsEmpty(Res) then 
Include(~D, PA);
else
Z:={};
for p in Res do
PAp:=ParentCan(Include(PA,p));
if Occurrences(PA,n) eq Occurrences(PAp,n) then
if PolType(PA,n) eq PolType(PAp,n) then 
if exists{T: T in GG|Image(T,PA) eq PAp}
then Include(~Z, Include(PA,p));
end if;
end if;
end if;
end for;
U:={};
for PB in Z do
if not(TestPAdansD(PB,n,U, GG)) then
Include(~U, PB);
end if;
end for;
for PB in U do
p:=Representative({q: q in PB|q notin PA});
NewRes:={q:q in Res|#Support(q*p^-1) ge d};
GenPAWeak(PB, NewRes,n,d,GG, ~D);
end for;
end if;
end procedure;




forward CanAugWM;
procedure CanAugWM(PAdepart,EGG,n,d,~i,Res,~EnsPAMax);
/* PAdepart contient le tableau, EGG contient les différents sous-groupes d'$Iso(Sym(n))$ dont la procédure a besoin, Res est l'ensembles des permutations que l'on peut ajouter au PA, EnsPAMAx est l'ensembles des PA maximaux déjà enregistrés.
*/
if IsEmpty(Res) then Include(~EnsPAMax, PAdepart);
#EnsPAMax;
else
GPAdep:=Stabiliser(EGG[1], PAdepart);
CPA:={Include(PAdepart, p): p in Res};
EcycDep:={* CycleStructure(p*q^-1): p, q in PAdepart*};
OccDep:={*#{p: p in PAdepart|Image(p, i) eq j} : i, j in [1..n]*};
for ZZ in {CPA meet {Image(T, Z): T in GPAdep}: Z in CPA} do
 Z:= Representative(ZZ);
 permAjoutee:=Representative(Z diff PAdepart);
 P4Z:=ParentCanWM(Z,EGG);
EcycP4Z:={* CycleStructure(p*q^-1): p, q in P4Z*};
OccP4Z:={*#{p: p in P4Z|Image(p, i) eq j} : i, j in [1..n]*};
 if ((EcycP4Z eq EcycDep) and (OccDep eq OccP4Z))  then
i:=i+1;
 i;
 StabZ:=Stabiliser(EGG[1],Z);
 if exists{T: T in StabZ|Image(T,P4Z) eq PAdepart} then 
 NewRes:={q: q in Res|DistPerm(permAjoutee, q ) ge d};
 CanAugWM(Z,EGG,  n, d,~i,NewRes, ~EnsPAMax);
 end if;
 end if;
end for;
end if;
end procedure;


forward CanAug;
procedure CanAug(PAdepart,GG,n,d,~i,Res,~EnsPAMax);
if IsEmpty(Res) then Include(~EnsPAMax, PAdepart);
#EnsPAMax;
else
GPAdep:=Stabiliser(GG, PAdepart);
CPA:={Include(PAdepart, p): p in Res};
EcycDep:={* CycleStructure(p*q^-1): p, q in PAdepart*};
OccDep:={*#{p: p in PAdepart|Image(p, i) eq j} : i, j in [1..n]*};
/* Calculs des invariants (coefficients du polynôme énumérateur des types et entrées de la matrice des occurrences)*/ 
for ZZ in {CPA meet {Image(T, Z): T in GPAdep}: Z in CPA} do
Z:= Representative(ZZ);
permAjoutee:=Representative(Z diff PAdepart);
P4Z:=ParentCan(Z);
EcycP4Z:={* CycleStructure(p*q^-1): p, q in P4Z*};
OccP4Z:={*#{p: p in P4Z|Image(p, i) eq j} : i, j in [1..n]*};

if ((EcycP4Z eq EcycDep) and (OccDep eq OccP4Z))  then
i:=i+1;
/*Le compteur i compte le nombre total d'augmentations effectuées*/ 
StabZ:=Stabiliser(GG,Z);
if exists{T: T in StabZ|Image(T,P4Z) eq PAdepart} then 
NewRes:={q: q in Res|DistPerm(permAjoutee, q ) ge d};
CanAug(Z,GG,  n, d,~i,NewRes, ~EnsPAMax);
end if;
end if;
end for;
end if;
end procedure;





forward CanAugBal;
procedure CanAugBal(PAdepart,GG,n,d,r,~i,Res,~EnsPAMax);

if IsRbalanced(PAdepart, r) then Include(~EnsPAMax, PAdepart);
#EnsPAMax;
else
if not(IsEmpty(Res)) then
GPAdep:=Stabiliser(GG, PAdepart);
CPA:={Include(PAdepart, p): p in Res|IsRbounded(Include(PAdepart, p), r)};
for ZZ in {CPA meet {Image(T, Z): T in GPAdep}: Z in CPA} do
 Z:= Representative(ZZ);
 permAjoutee:=Representative(Z diff PAdepart);
 i:=i+1;
 NewRes:={q: q in Res|DistPerm(permAjoutee, q ) ge d};

 StabZ:=Stabiliser(GG,Z); 
 if exists{T: T in StabZ|Image(T,ParentCan(Z)) eq PAdepart} then 
CanAugBal(Z,GG,  n, d,r,~i,NewRes, ~EnsPAMax);
end if;
end for; 
end if;
end if;
end procedure;