load "Subgroups/functions.txt";

function LargeModulesSelfDual(Irr,n)

Irrsmall:=[]; 
for i in Irr do if(Dimension(i) le n) then Append(~Irrsmall,i); end if; end for;

for i in [1..#Irrsmall] do
  if(not(IsSelfDual(Irrsmall[i]))) then
    Irrsmall[i]:=DirectSum(Irrsmall[i],Dual(Irrsmall[i]));
  end if;
end for;
Irrsmall:=StripDuplicates(Irrsmall);
IrrDual:=[];
for i in Irrsmall do if(Dimension(i) le n) then Append(~IrrDual,i); end if; end for;
return OrderByDimension(IrrDual);

end function;

function ExpressElementInGenerators(G,x)

d:=#Generators(G);
Subs1:=[[1..i-1] cat [i+1..d]:i in [1..d]];
Norms:=[]; funcs:=[];
for i in Subs1 do H,f:=quo<G|[G.j:j in i]>; Append(~Norms,H); Append(~funcs,f); end for;
El:=[];
for i in [1..#Norms] do
  for j in [0..Order(G.i)-1] do if(funcs[i](G.i^j) eq funcs[i](x)) then Append(~El,j); continue i; end if; end for;
end for;

return El;

end function;

/* This function takes a sequence X of sequences and produces the Cartesian product of these sequences,
   namely all sequences whose ith term comes from the ith term of X. */

function AllSequences(X)

if(#X eq 1) then return [[i]:i in X[1]]; end if;
A:=AllSequences(Prune(X));
return [Append(i,j):i in A,j in X[#X]];
end function;

/* This function, for a homocyclic group of exponent n, takes an element El and constructs all elements whose
   mth power is El. */

function ConstructPreimages(El,n,m)
Powers:=[];
for i in El do Append(~Powers,[j:j in [0..n-1] | (m*j mod n) eq i]); end for;
return AllSequences(Powers);
end function;

function ProduceGroupElements(G,Seqs)

return [&*[G.i^j[i]:i in [1..#Generators(G)]]:j in Seqs];
end function;

function ProduceGroupElement(G,Seq)

return &*[G.i^Seq[i]:i in [1..#Generators(G)]];
end function;

function EigenvaluesOfElement(M,x)
X:=sub<Group(M)|x>;
return Eigenvalues(ActionGenerators(Restriction(M,X))[1]);
end function;

function FindOrb(Orbs,i)
for Orb in Orbs do if(i in Orb) then return Orb; end if; end for;
return {};
end function;


function OrbitReps(X,frob)

Orbs:=[];
SpanOfOrbs:=[];
for i in X do
  if(not(Sort(i) in SpanOfOrbs)) then
    Append(~Orbs,i);
    SpanOfOrbs cat:=[Sort(i^(frob^j)):j in [1..Order(frob)-1]];
  end if;
end for;
return Orbs;
end function;

q:=2^8;
n:=1;
F<w>:=GF(q);
G:=GroupOfLieType("E8",q);
V:=VectorSpace(GF(q),8);
rho1:=StandardRepresentation(G);
Over1:=GL(248,q);
g1:=elt<G|V![w^n,1,1,1,1,1,1,1]>;
g2:=elt<G|V![1,w^n,1,1,1,1,1,1]>;
g3:=elt<G|V![1,1,w^n,1,1,1,1,1]>;
g4:=elt<G|V![1,1,1,w^n,1,1,1,1]>;
g5:=elt<G|V![1,1,1,1,w^n,1,1,1]>;
g6:=elt<G|V![1,1,1,1,1,w^n,1,1]>;
g7:=elt<G|V![1,1,1,1,1,1,w^n,1]>;
g8:=elt<G|V![1,1,1,1,1,1,1,w^n]>;
W:=Reflections(G);
Mats1:=[rho1(i):i in [g1,g2,g3,g4,g5,g6,g7,g8]];
MatsE:=[rho1(i):i in W];
NGT:=sub<Over1|Mats1 cat MatsE>;
T:=sub<NGT|Mats1>;
E:=sub<NGT|MatsE>;

T17:=sub<T|[T.i^15:i in [1..8]]>;
T51:=sub<T|[T.i^5:i in [1..8]]>;
T85:=sub<T|[T.i^3:i in [1..8]]>;

R:=ResidueClassRing(255);
X:=[[R!Log((T.j)[i,i]):i in [1..248]]:j in [1..8]];


  G2:=SU(3,16);
  CC:=ConjugacyClasses(G2);


  S1:=TrivialModule(G2,GF(256));
  S31:=GModule(G2); S31d:=Dual(S31);
  S32:=GModule(G2,FrobeniusImage(S31,1)); S32d:=Dual(S32);
  S33:=GModule(G2,FrobeniusImage(S31,2)); S33d:=Dual(S33);
  S34:=GModule(G2,FrobeniusImage(S31,3)); S34d:=Dual(S34);
  S81:=TensorProduct(S31,S31d); S81:=S81/Fix(S81);
  S82:=TensorProduct(S32,S32d); S82:=S82/Fix(S82);
  S83:=TensorProduct(S33,S33d); S83:=S83/Fix(S83);
  S84:=TensorProduct(S34,S34d); S84:=S84/Fix(S84);
  S912:=TensorProduct(S31,S32); S912d:=Dual(S912);
  S921:=TensorProduct(S31,S32d); S921d:=Dual(S921);
  S913:=TensorProduct(S31,S33); S913d:=Dual(S913);
  S931:=TensorProduct(S31,S33d); S931d:=Dual(S931);
  S914:=TensorProduct(S31,S34); S914d:=Dual(S914);
  S941:=TensorProduct(S31,S34d); S941d:=Dual(S941);
  S923:=TensorProduct(S32,S33); S923d:=Dual(S923);
  S932:=TensorProduct(S32,S33d); S932d:=Dual(S932);
  S924:=TensorProduct(S32,S34); S924d:=Dual(S924);
  S942:=TensorProduct(S32,S34d); S942d:=Dual(S942);
  S934:=TensorProduct(S33,S34); S934d:=Dual(S934);
  S943:=TensorProduct(S33,S34d); S943d:=Dual(S943);
  S2412:=TensorProduct(S31,S82); S2412d:=Dual(S2412);
  S2421:=TensorProduct(S32,S81); S2421d:=Dual(S2421);
  S2413:=TensorProduct(S31,S83); S2413d:=Dual(S2413);
  S2431:=TensorProduct(S33,S81); S2431d:=Dual(S2431);
  S2414:=TensorProduct(S31,S84); S2414d:=Dual(S2414);
  S2441:=TensorProduct(S34,S81); S2441d:=Dual(S2441);
  S2423:=TensorProduct(S32,S83); S2423d:=Dual(S2423);
  S2432:=TensorProduct(S33,S82); S2432d:=Dual(S2432);
  S2424:=TensorProduct(S32,S84); S2424d:=Dual(S2424);
  S2442:=TensorProduct(S34,S82); S2442d:=Dual(S2442);
  S2434:=TensorProduct(S33,S84); S2434d:=Dual(S2434);
  S2443:=TensorProduct(S34,S83); S2443d:=Dual(S2443);
  S6412:=TensorProduct(S81,S82);
  S6413:=TensorProduct(S81,S83);
  S6414:=TensorProduct(S81,S84);
  S6423:=TensorProduct(S82,S83);
  S6424:=TensorProduct(S82,S84);
  S6434:=TensorProduct(S83,S84);

  Irr2:=[S1,S31,S31d,S32,S32d,S33,S33d,S34,S34d,S81,S82,S83,S84,S912,S912d,S921,S921d,S913,S913d,S931,
S931d,S914,S914d,S941,S941d,S923,S923d,S924,S924d,S942,S942d,S932,S932d,S934,S934d,S943,S943d,
S2412,S2412d,S2421,S2421d,S2413,S2413d,S2431,S2431d,S2414,S2414d,S2441,S2441d,S2423,S2423d,
S2432,S2432d,S2424,S2424d,S2442,S2442d,S2434,S2434d,S2443,S2443d,S6412,S6413,S6414,S6423,S6424,S6434];

  for i in Irr2 do if(Dimension(i) eq 3) then
     for j in Irr2 do if(Dimension(j) eq 9) then
       Irr2:=AppendWithoutDuplicates(Irr2,CompositionFactors(TensorProduct(i,j)));
     end if; end for;
  end if; end for;
  for i in Irr2 do if(Dimension(i) eq 8) then
     for j in Irr2 do if(Dimension(j) eq 9) then
       Irr2:=AppendWithoutDuplicates(Irr2,CompositionFactors(TensorProduct(i,j)));
     end if; end for;
  end if; end for;
  for i in Irr2 do if(Dimension(i) eq 9) then
     for j in Irr2 do if(Dimension(j) eq 9) then
       Irr2:=AppendWithoutDuplicates(Irr2,CompositionFactors(TensorProduct(i,j)));
     end if; end for;
  end if; end for;
  Irr2:=OrderByDimension(Irr2);


  Irr:=LargeModulesSelfDual(Irr2,248);
  for i in [211..274] do if(EigenvaluesOfElement(S31,CC[i,3]) eq { <w^15, 1>, <w^239, 1>, <w, 1> }) then x:=CC[i,3]; end if; end for;

AnsTot3:=[
    [ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 6, 6, 6, 6, 6, 6, 6, 6, 9, 15, 15, 15, 15, 15, 15, 15, 15, 17 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 6, 6, 6, 6, 15, 15, 15, 15, 15, 15, 17, 17, 52 ],
    [ 1, 1, 1, 1, 1, 1, 1, 1, 6, 7, 8, 10, 10, 10, 13, 13, 13, 16, 16, 16, 41 ],
    [ 1, 1, 1, 1, 1, 1, 6, 6, 10, 15, 15, 15, 17, 17, 39, 52 ],
    [ 1, 1, 1, 1, 10, 10, 15, 17, 17, 17, 21, 39, 54 ],
    [ 1, 1, 2, 2, 2, 2, 2, 3, 5, 5, 5, 5, 6, 6, 6, 9, 9, 9, 14, 14, 14, 15, 15, 26 ],
    [ 1, 1, 2, 3, 6, 6, 7, 7, 7, 7, 10, 10, 11, 11, 12, 22, 23 ],
    [ 1, 1, 2, 6, 7, 7, 8, 8, 8, 16, 16, 22, 34, 35 ],
    [ 1, 1, 2, 6, 8, 8, 9, 9, 9, 20, 20, 24, 42, 43 ],
    [ 2, 2, 2, 2, 3, 3, 4, 6, 6, 6, 8, 11, 12, 12, 13, 13, 16, 19, 25 ],
    [ 2, 2, 2, 2, 3, 3, 5, 6, 6, 6, 9, 11, 14, 14, 15, 15, 17, 18, 27 ],
    [ 2, 2, 2, 3, 5, 6, 6, 9, 9, 14, 14, 15, 15, 17, 26, 27 ],
    [ 2, 2, 3, 3, 6, 7, 9, 9, 14, 14, 30, 36, 39 ],
    [ 2, 2, 3, 4, 6, 8, 9, 9, 14, 14, 32, 42, 45 ],
    [ 2, 4, 6, 6, 8, 8, 12, 12, 13, 13, 14, 19, 24, 25 ],
    [ 2, 4, 6, 6, 8, 8, 12, 12, 13, 13, 15, 16, 24, 25 ],
    [ 2, 6, 7, 7, 9, 17, 17, 21, 22, 36, 37 ],
    [ 2, 6, 7, 7, 9, 18, 18, 20, 22, 38, 39 ],
    [ 2, 6, 7, 8, 8, 16, 16, 18, 24, 34, 35 ],
    [ 2, 6, 7, 8, 8, 17, 19, 19, 24, 40, 41 ],
    [ 2, 6, 7, 9, 9, 10, 17, 17, 26, 36, 37 ],
    [ 2, 6, 7, 9, 9, 11, 18, 18, 26, 38, 39 ],
    [ 2, 6, 8, 9, 9, 12, 20, 20, 26, 42, 43 ],
    [ 2, 6, 8, 9, 9, 13, 21, 21, 26, 44, 45 ],
    [ 6, 7, 8, 9, 34, 37, 45, 47 ],
    [ 6, 7, 8, 9, 34, 39, 43, 49 ],
    [ 6, 7, 8, 9, 35, 36, 44, 47 ],
    [ 6, 7, 8, 9, 37, 40, 42, 48 ]
];

Mods:=[DirectSum([Irr[i]:i in j]):j in AnsTot3];
Eigs1:=[EigenvaluesOfElement(i,x):i in Mods];
Eigs2:=[EigenvaluesOfElement(i,x^3):i in Mods];
Eigs4:=[EigenvaluesOfElement(i,x^15):i in Mods];



El1:=[
    [ 4, 10, 3, 7, 0, 0, 0, 0 ],
    [ 10, 0, 3, 14, 0, 0, 0, 0 ],
    [ 3, 0, 14, 3, 0, 0, 0, 0 ],
    [ 6, 0, 4, 10, 0, 0, 0, 0 ],
    [ 16, 0, 10, 14, 0, 0, 0, 0 ],
    [ 10, 14, 9, 1, 0, 7, 1, 10 ],
    [ 14, 0, 6, 11, 3, 4, 11, 0 ],
    [ 0, 4, 0, 1, 6, 1, 4, 10 ],
    [ 5, 7, 5, 12, 15, 13, 1, 16 ],
    [ 7, 7, 9, 0, 16, 3, 11, 8 ],
    [ 8, 0, 16, 16, 9, 2, 10, 13 ],
    [ 14, 2, 1, 15, 2, 8, 1, 14 ],
    [ 11, 16, 0, 15, 13, 15, 13, 2 ],
    [ 9, 6, 8, 13, 8, 11, 3, 2 ],
    [ 14, 13, 2, 4, 1, 4, 7, 3 ],
    [ 11, 14, 6, 6, 8, 4, 2, 15 ],
    [ 8, 10, 4, 3, 10, 1, 9, 5 ],
    [ 7, 8, 0, 4, 1, 11, 16, 16 ],
    [ 6, 4, 4, 14, 10, 12, 4, 2 ],
    [ 14, 15, 4, 14, 4, 9, 5, 2 ],
    [ 0, 1, 9, 15, 13, 8, 3, 11 ],
    [ 2, 14, 16, 11, 11, 10, 0, 0 ],
    [ 14, 9, 4, 7, 1, 11, 10, 0 ],
    [ 10, 10, 1, 5, 2, 12, 0, 5 ],
    [ 10, 6, 14, 4, 8, 12, 12, 10 ],
    [ 1, 2, 16, 5, 0, 13, 3, 0 ],
    [ 3, 9, 1, 3, 16, 5, 3, 2 ],
    [ 14, 6, 9, 7, 9, 5, 1, 8 ]
];

El2:=[
    [ 216, 81, 162, 21, 51, 0, 0, 0 ],
    [ 183, 0, 213, 93, 51, 51, 0, 0 ],
    [ 111, 51, 195, 111, 51, 51, 0, 0 ],
    [ 120, 51, 216, 81, 51, 51, 0, 0 ],
    [ 99, 153, 81, 93, 51, 51, 51, 0 ],
    [ 132, 93, 231, 105, 102, 21, 3, 81 ],
    [ 144, 51, 69, 84, 9, 114, 84, 102 ],
    [ 0, 63, 0, 105, 18, 3, 12, 30 ],
    [ 168, 174, 66, 87, 96, 39, 3, 48 ],
    [ 72, 21, 231, 204, 48, 60, 84, 75 ],
    [ 126, 153, 150, 48, 78, 57, 30, 39 ],
    [ 42, 57, 207, 96, 57, 24, 3, 42 ],
    [ 84, 150, 102, 147, 192, 147, 39, 6 ],
    [ 78, 120, 24, 39, 75, 33, 9, 6 ],
    [ 246, 192, 108, 12, 3, 12, 123, 9 ],
    [ 186, 195, 171, 18, 75, 12, 6, 45 ],
    [ 24, 30, 12, 111, 30, 3, 231, 66 ],
    [ 123, 24, 153, 63, 105, 33, 99, 48 ],
    [ 120, 12, 12, 42, 132, 36, 63, 6 ],
    [ 93, 198, 165, 144, 63, 27, 66, 6 ],
    [ 102, 3, 78, 45, 141, 126, 9, 33 ],
    [ 57, 144, 201, 84, 135, 81, 51, 51 ],
    [ 246, 27, 114, 21, 3, 33, 81, 51 ],
    [ 81, 132, 54, 66, 108, 36, 51, 15 ],
    [ 183, 171, 42, 63, 24, 240, 36, 30 ],
    [ 207, 108, 150, 66, 51, 39, 9, 51 ],
    [ 60, 129, 156, 111, 201, 219, 9, 57 ],
    [ 144, 69, 27, 174, 129, 219, 54, 24 ]
];

El3:=[
    [ 157, 112, 224, 7, 17, 0, 0, 0 ],
    [ 231, 0, 241, 31, 102, 17, 0, 0 ],
    [ 37, 17, 235, 207, 17, 17, 0, 0 ],
    [ 210, 17, 157, 112, 17, 102, 0, 0 ],
    [ 203, 136, 112, 31, 102, 17, 17, 85 ],
    [ 129, 31, 247, 120, 119, 7, 1, 112 ],
    [ 48, 17, 193, 28, 3, 208, 28, 34 ],
    [ 0, 106, 0, 205, 6, 1, 174, 95 ],
    [ 56, 228, 192, 199, 202, 183, 1, 16 ],
    [ 194, 177, 77, 238, 16, 20, 28, 195 ],
    [ 127, 136, 135, 16, 111, 104, 10, 13 ],
    [ 14, 104, 239, 32, 104, 8, 1, 14 ],
    [ 113, 135, 119, 134, 234, 134, 183, 87 ],
    [ 111, 125, 8, 183, 195, 11, 173, 87 ],
    [ 252, 234, 206, 4, 1, 4, 41, 3 ],
    [ 62, 235, 227, 91, 195, 4, 172, 15 ],
    [ 8, 10, 4, 122, 95, 1, 77, 192 ],
    [ 126, 178, 136, 21, 35, 96, 118, 16 ],
    [ 125, 89, 174, 184, 44, 12, 21, 172 ],
    [ 31, 66, 225, 218, 21, 179, 107, 172 ],
    [ 34, 1, 111, 15, 217, 127, 173, 96 ],
    [ 104, 218, 152, 28, 130, 112, 102, 17 ],
    [ 252, 9, 123, 7, 1, 96, 27, 187 ],
    [ 112, 44, 188, 22, 121, 12, 187, 175 ],
    [ 61, 142, 99, 21, 178, 250, 182, 10 ],
    [ 69, 206, 220, 192, 17, 183, 3, 102 ],
    [ 190, 43, 222, 122, 152, 158, 88, 189 ],
    [ 218, 108, 94, 228, 43, 73, 18, 8 ]
];

if(OneToCheck eq 10) then CE:=[E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.5*E.4*E.2*E.3*E.4*E.5,
E.2*E.4*E.3*E.1*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2,
E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7,
E.4,
E.3*E.1*E.4*E.2*E.3*E.4*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3];

elif(OneToCheck eq 15) then CE:=[E.1*E.2*E.3*E.4*E.5*E.4*E.2*E.3*E.1,
E.2*E.3*E.4*E.5*E.6*E.7*E.8*E.7*E.6*E.5*E.4*E.2*E.3,
E.2*E.4*E.2,
E.3*E.1*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3];

elif(OneToCheck eq 19) then CE:=[E.1*E.3*E.4*E.5*E.4*E.3*E.1,
E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.5*E.4*E.2*E.3*E.4*E.5,
E.1*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4,
E.2*E.4*E.3*E.1*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2,
E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.8*E.7*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3*E.5*E.4*E.2*E.6*E.5*E.4*E.3*E.1*E.7*E.6*E.5*E.4*E.2*E.3*E.4*E.5*E.6*E.7*E.8,
E.3*E.1*E.4*E.2*E.3*E.4*E.5*E.6*E.5*E.4*E.2*E.3*E.1*E.4*E.3];

else CE:=E;

end if;

Els1:=ConstructPreimages([15*i:i in El1[OneToCheck]],255,5);
E11:={*i[1]^^i[2]:i in Eigs2[OneToCheck]*};
E1:={*R!(Log(i)):i in E11*};
NextEls1:=[];
for ii in [1..#Els1] do
  hh:={*&+[Els1[ii,i]*X[i,j]:i in [1..8]]:j in [1..120]*};
  if E1 eq hh join {*-i:i in hh*} join {*0^^8*} then Append(~NextEls1,Els1[ii]); end if;
  if(ii mod 10000 eq 0) then "At stage",ii,"of",#Els1,"there are",#NextEls1,"elements with the correct eigenvalues"; end if;
end for;
#NextEls1;

gg1:=ProduceGroupElements(T,NextEls1);
g:=gg1[1]^5;
Orbs:={{i}:i in [1..#NextEls1]};
EltsToMerge:=[]; C:=sub<E|EltsToMerge>;
for e in CE do if(e*g eq g*e and not(e in C)) then
  "Found an element to merge orbits with";
  Append(~EltsToMerge,e); C:=sub<E|EltsToMerge>;
  for i in [1..#NextEls1] do
    j:=Position(gg1,gg1[i]^e);
    Orb1:=FindOrb(Orbs,i);
    if(not(j in Orb1)) then
      Orb2:=FindOrb(Orbs,j);
      Exclude(~Orbs,Orb1);
      Exclude(~Orbs,Orb2);
      Include(~Orbs,Orb1 join Orb2); //printf "Merged %o and %o.\n",Orb1,Orb2;
      delete Orb1; delete Orb2; if(#Orbs eq 1) then "All merged"; break e; end if;
    end if;
  end for;
  "So far we have merged elements into",#Orbs,"orbits";
end if; end for;

Els2:=ConstructPreimages(NextEls1[1],255,3);
E22:={*i[1]^^i[2]:i in Eigs1[OneToCheck]*};
E2:={*R!(Log(i)):i in E22*};
NextEls2:=[];
for ii in [1..#Els2] do
  hh:={*&+[Els2[ii,i]*X[i,j]:i in [1..8]]:j in [1..120]*};
  if E2 eq hh join {*-i:i in hh*} join {*0^^8*} then Append(~NextEls2,Els2[ii]); end if;
  if(ii mod 1000 eq 0) then "At stage",ii,"of",#Els2,"there are",#NextEls2,"elements with the correct eigenvalues"; end if;
end for;
#NextEls2;

if(#NextEls2 eq 1) then "Nothing to do here";
else

gg2:=ProduceGroupElements(T,NextEls2);
g:=gg2[1]^3;
Orbs:={{i}:i in [1..#NextEls2]};
EltsToMerge:=[]; C2:=sub<E|EltsToMerge>;
for e in E do if(e*g eq g*e and not(e in C2)) then
  "Found an element to merge orbits with";
  Append(~EltsToMerge,e); C2:=sub<E|EltsToMerge>;
  for i in [1..#NextEls2] do
    j:=Position(gg2,gg2[i]^e);
    Orb1:=FindOrb(Orbs,i);
    if(not(j in Orb1)) then
      Orb2:=FindOrb(Orbs,j);
      Exclude(~Orbs,Orb1);
      Exclude(~Orbs,Orb2);
      Include(~Orbs,Orb1 join Orb2); //printf "Merged %o and %o.\n",Orb1,Orb2;
      delete Orb1; delete Orb2; if(#Orbs eq 1) then "All merged"; break e; end if;
    end if;
  end for;
  "So far we have merged elements into",#Orbs,"orbits";
end if; end for;

end if;

El2[OneToCheck] eq NextEls1[1] and El3[OneToCheck] eq NextEls2[1];
