# This code sets up the Galois group over R(z_1,z_2) of the field
# extension KT generated by x[1], x[2], y[1], y[2], r[1] and r[2].

#@ G64 
G64 := [
 1, L, LL, LLL, M, LM, LLM, LLLM, N, LN, LLN, LLLN, MN, LMN, LLMN, LLLMN,
 A1, LA1, LLA1, LLLA1, MA1, LMA1, LLMA1, LLLMA1,
 NA1, LNA1, LLNA1, LLLNA1, MNA1, LMNA1, LLMNA1, LLLMNA1,
 A2, LA2, LLA2, LLLA2, MA2, LMA2, LLMA2, LLLMA2,
 NA2, LNA2, LLNA2, LLLNA2, MNA2, LMNA2, LLMNA2, LLLMNA2,
 A12, LA12, LLA12, LLLA12, MA12, LMA12, LLMA12, LLLMA12,
 NA12, LNA12, LLNA12, LLLNA12, MNA12, LMNA12, LLMNA12, LLLMNA12
]:

protect(
 'A1', 'LA1', 'LLA1', 'LLLA1', 'MA1', 'LMA1', 'LLMA1', 'LLLMA1',
 'NA1', 'LNA1', 'LLNA1', 'LLLNA1', 'MNA1', 'LMNA1', 'LLMNA1', 'LLLMNA1',
 'A2', 'LA2', 'LLA2', 'LLLA2', 'MA2', 'LMA2', 'LLMA2', 'LLLMA2',
 'NA2', 'LNA2', 'LLNA2', 'LLLNA2', 'MNA2', 'LMNA2', 'LLMNA2', 'LLLMNA2',
 'A12', 'LA12', 'LLA12', 'LLLA12', 'MA12', 'LMA12', 'LLMA12', 'LLLMA12',
 'NA12', 'LNA12', 'LLNA12', 'LLLNA12', 'MNA12', 'LMNA12', 'LLMNA12', 'LLLMNA12'
):

act_A[A1 ] := (u) -> eval(subs(r=rx,subs(r[1]=-r[1],subs(root_rule,u)))):
act_A[A2 ] := (u) -> eval(subs(r=rx,subs(r[2]=-r[2],subs(root_rule,u)))):
act_A[A12] := (u) -> eval(subs(r=rx,subs({r[1]=-r[1],r[2]=-r[2]},subs(root_rule,u)))):

#@ act_AT_rule 
act_AT_rule := table():

#@ xyrstu_vars 
xyrstu_vars := [
 seq(x[i],i=1..4),
 seq(y[i],i=1..2),
 seq(r[i],i=1..2),
 seq(s[i],i=1..16),
 seq(t[i],i=1..4),
 seq(u[i],i=1..16)
];

act_AT_rule[1] := {seq(a = a,a in xyrstu_vars)};

act_AT_rule[L] := {
 x[1] = x[2], x[2] = -x[1], x[4] = -x[4], y[2] = -y[2], r[1] = r[2], r[2] = r[1],
 t[1] =  t[2], t[2] = -t[1], t[3] =  t[4], t[4] = -t[3],
 seq(s[2*i-1] = s[2*i],i=1..4), seq(s[2*i] = s[2*i-1],i=1..4),
 s[ 9] = s[10], s[10] =-s[ 9], s[11] = s[12], s[12] = -s[11],
 s[13] = -s[13], s[14] =-s[15], s[15] = -s[14], s[16] = -s[16],
 u[1] = u[2], u[2] = -u[1], u[3] = u[4], u[4] = -u[3], u[5] = u[6],
 u[6] = -u[5], u[7] = u[8], u[8] = -u[7], u[9] = u[10], u[10] = -u[9],
 u[11] = u[12], u[12] = -u[11], u[13] = u[14], u[14] = -u[13],
 u[15] = u[16], u[16] = -u[15]
}:

act_AT_rule[M] := {
 x[2] = -x[2], x[3] = -x[3], x[4] = -x[4], y[1] = -y[1],
 t[1] =  t[3], t[2] = -t[4], t[3] =  t[1], t[4] = -t[2],
 s[3] = 1/s[3], s[4] = 1/s[4],
 s[9] = s[11], s[10] = -s[12], s[11] = s[9], s[12] = -s[10],
 u[2] = -u[2], u[4] = -u[4], u[5] = -u[5], u[7] = -u[7], u[9] = u[11],
 u[10] = -u[12], u[11] = u[9], u[12] = -u[10], u[13] = -u[15], u[14] = u[16],
 u[15] = -u[13], u[16] = u[14]
}:

act_AT_rule[N] := {
 x[2] = -x[2],
 t[1] =  t[1], t[2] = -t[2], t[3] =  t[3], t[4] = -t[4],
 s[10] = -s[10], s[12] = -s[12],
 seq(s[i] = -s[i],i=13..16),
 u[2] = -u[2], u[4] = -u[4], u[6] = -u[6], u[8] = -u[8],
 u[10] = -u[10], u[12] = -u[12], u[14] = -u[14], u[16] = -u[16]
}:

act_AT_rule[A1] := {
 r[1] = -r[1],
 t[1] = -1/t[1], t[3] = -1/t[3], s[1] = 1/s[1], s[3] = 1/s[3], s[ 9] = -s[ 9], s[11] = -s[11],
 s[13] = -s[15], s[14] = -s[16], s[15] = -s[13], s[16] = -s[14],
 u[1] = u[3], u[3] = u[1], u[5] = u[7], u[7] = u[5], u[9] = u[11], u[11] = u[9], u[13] = u[15], u[15] = u[13]
}:

act_AT_rule[A2] := {
 r[2] = -r[2],
 t[2] = -1/t[2], t[4] = -1/t[4], s[2] = 1/s[2], s[4] = 1/s[4], s[10] = -s[10], s[12] = -s[12],
 s[13] = -s[14], s[14] = -s[13], s[15] = -s[16], s[16] = -s[15],
 u[2] = u[4], u[4] = u[2], u[6] = u[8], u[8] = u[6], u[10] = u[12], u[12] = u[10], u[14] = u[16], u[16] = u[14]
}:

act_AT[1] := (a) -> a:
for T in [L,M,N,A1,A2] do 
 act_AT[T] := subs(rule=act_AT_rule[T],eval((a) -> subs(rule,a))); 
od:

#@ make_AT_rule 
make_AT_rule := proc(T,U,TU)
 global act_AT_rule,act_AT;
 act_AT_rule[TU] := {seq(v = act_AT[T](act_AT[U](v)),v in xyrstu_vars)};
 act_AT[TU] := subs(rule=act_AT_rule[TU],eval((a) -> subs(rule,a))); 
 NULL;
end:

make_AT_rule(L,L,LL):
make_AT_rule(L,LL,LLL):
make_AT_rule(L,M,LM):
make_AT_rule(L,LM,LLM):
make_AT_rule(L,LLM,LLLM):
make_AT_rule(L,N,LN):
make_AT_rule(L,LN,LLN):
make_AT_rule(L,LLN,LLLN):
make_AT_rule(M,N,MN):
make_AT_rule(L,MN,LMN):
make_AT_rule(L,LMN,LLMN):
make_AT_rule(L,LLMN,LLLMN):
make_AT_rule(A1,A2,A12):
make_AT_rule(L,A1,LA1):
make_AT_rule(LL,A1,LLA1):
make_AT_rule(LLL,A1,LLLA1):
make_AT_rule(M,A1,MA1):
make_AT_rule(LM,A1,LMA1):
make_AT_rule(LLM,A1,LLMA1):
make_AT_rule(LLLM,A1,LLLMA1):
make_AT_rule(N,A1,NA1):
make_AT_rule(LN,A1,LNA1):
make_AT_rule(LLN,A1,LLNA1):
make_AT_rule(LLLN,A1,LLLNA1):
make_AT_rule(MN,A1,MNA1):
make_AT_rule(LMN,A1,LMNA1):
make_AT_rule(LLMN,A1,LLMNA1):
make_AT_rule(LLLMN,A1,LLLMNA1):
make_AT_rule(L,A2,LA2):
make_AT_rule(LL,A2,LLA2):
make_AT_rule(LLL,A2,LLLA2):
make_AT_rule(M,A2,MA2):
make_AT_rule(LM,A2,LMA2):
make_AT_rule(LLM,A2,LLMA2):
make_AT_rule(LLLM,A2,LLLMA2):
make_AT_rule(N,A2,NA2):
make_AT_rule(LN,A2,LNA2):
make_AT_rule(LLN,A2,LLNA2):
make_AT_rule(LLLN,A2,LLLNA2):
make_AT_rule(MN,A2,MNA2):
make_AT_rule(LMN,A2,LMNA2):
make_AT_rule(LLMN,A2,LLMNA2):
make_AT_rule(LLLMN,A2,LLLMNA2):
make_AT_rule(L,A12,LA12):
make_AT_rule(LL,A12,LLA12):
make_AT_rule(LLL,A12,LLLA12):
make_AT_rule(M,A12,MA12):
make_AT_rule(LM,A12,LMA12):
make_AT_rule(LLM,A12,LLMA12):
make_AT_rule(LLLM,A12,LLLMA12):
make_AT_rule(N,A12,NA12):
make_AT_rule(LN,A12,LNA12):
make_AT_rule(LLN,A12,LLNA12):
make_AT_rule(LLLN,A12,LLLNA12):
make_AT_rule(MN,A12,MNA12):
make_AT_rule(LMN,A12,LMNA12):
make_AT_rule(LLMN,A12,LLMNA12):
make_AT_rule(LLLMN,A12,LLLMNA12):

#@ average_AT 
average_AT := (a,H) -> add(act_AT[T](a)/nops(H),T in H);

#@ G64_mult 
G64_mult := proc() procname(args) end:

#@ G64_inv  
G64_inv  := proc() procname(args) end:

#@ G64_lookup
for T in G64 do
 G64_lookup[map(act_AT[T],[t[1],t[2],t[3],t[4]])] := T;
od:
for T1 in G64 do
 for T2 in G64 do
  T3 := G64_lookup[[seq(act_AT[T1](act_AT[T2](t[i])),i=1..4)]];
  G64_mult(T1,T2) := T3;
  if T3 = 1 then
   G64_inv(T1) := T2;
  fi;
 od:
od:


#@ G64_conj 
G64_conj := (a,b) -> G64_mult(G64_mult(a,b),G64_inv(a)):

#@ G64_classes 
G64_classes := [op(map(T -> map(U -> G64_conj(U,T),{op(G64)}),{op(G64)}))];

stab_A  := (u) -> select(T -> act_A[T](u) = u,G16);   #@ stab_A
stab_AT := (u) -> select(T -> act_AT[T](u) = u,G64);  #@ stab_AT
orb_AT  := (v) -> {op(map(T -> act_AT[T](v),G64))};   #@ orb_AT
uorb_AT := (v) -> map(strip_sign,orb_AT(v));          #@ uorb_AT

#@ G64_subgroup_generated_by
G64_subgroup_generated_by := proc(L::set)
 local H0,H1,T,U;
 H0 := {1};
 H1 := L;

 while H1 <> H0 do
  H0 := H1;
  H1 := {op(H1),seq(seq(G64_mult(T,U),T in H0),U in H0)};
 od;

 return sort([op(H1)]);
end:

# This code is not run by default, because it is slow, and not too useful.
# There are 569 subgroups altogether.

#@ find_G64_subgroups 
#@ G64_cyclic_subgroups 
#@ G64_subgroups
#@ G64_subgroup_join;

find_G64_subgroups := proc()
 global G64_cyclic_subgroups,G64_subgroups,G64_subgroup_join;
 local T,U,C,n,SG;

 G64_cyclic_subgroups := {}:
 for T in G64 do
  U := T;
  C := {1};
  while U <> 1 do 
   C := {op(C),U};
   U := G64_mult(T,U);
  od:
  G64_cyclic_subgroups := {op(G64_cyclic_subgroups),sort(C)};
 od:

 G64_subgroup_join := proc(H1,H2)
  local U,H,HH,T1,T2;
  U := H1 union H2;
  H := U;
  while true do
   HH := {seq(seq(G64_mult(T1,T2),T1 in H),T2 in U)};
   if HH = H then break; fi;
   H := HH;
  od:

  return sort(H);
 end:

 G64_subgroups := G64_cyclic_subgroups:
 n := 6;
 while n>0 do
  SG := {seq(seq(G64_subgroup_join(H1,H2),H1 in G64_subgroups),H2 in G64_cyclic_subgroups)};
  if nops(SG) = nops(G64_subgroups) then break; fi;
  G64_subgroups := SG;
  print(nops(SG));
  n := n-1;
 od:
end: