# Let K be the field of rational functions on EX^*, and let KT be the biquadratic
# extension generated by r[1] = sqrt(1 - y[2]/sqrt(2)) and r[2] = sqrt(1 + y[2]/sqrt(2)).
# This file defines many elements of KT and relations between them, which are
# useful for studying the Galois theory of various intermediate fields.
# The relevant Galois group is defined in the file group64.mpl

# These elements t[i] are used when studying torus quotients of EX^*

#@ tx
tx[1] := (r[1]+1/sqrt(2))*(1-r[1]*y[1])/x[1];
tx[2] := (r[2]+1/sqrt(2))*(1-r[2]*y[1])/x[2];
tx[3] := (r[1]+1/sqrt(2))*(1+r[1]*y[1])/x[1];
tx[4] := (r[2]+1/sqrt(2))*(1+r[2]*y[1])/x[2];

# These are expressions for t[i] in terms of y[1], r[1] and r[2].
# They are correct on F16 but have the wrong sign in some other places.
#@ tr
tr[1] := sqrt((r[1]+1/sqrt(2))*(1-r[1]*y[1])/((r[1]-1/sqrt(2))*(1+r[1]*y[1])));
tr[2] := sqrt((r[2]+1/sqrt(2))*(1-r[2]*y[1])/((r[2]-1/sqrt(2))*(1+r[2]*y[1])));
tr[3] := sqrt((r[1]+1/sqrt(2))*(1+r[1]*y[1])/((r[1]-1/sqrt(2))*(1-r[1]*y[1])));
tr[4] := sqrt((r[2]+1/sqrt(2))*(1+r[2]*y[1])/((r[2]-1/sqrt(2))*(1-r[2]*y[1])));

# Some further elements with nice transformation properties
#@ st
st[ 1] := t[3]*t[1]:
st[ 2] := t[4]*t[2]:
st[ 3] := t[3]/t[1]:
st[ 4] := t[4]/t[2]:
st[ 5] := st[1]/2+1/st[1]/2-2:
st[ 6] := st[2]/2+1/st[2]/2-2:
st[ 7] := st[3]/2+1/st[3]/2:
st[ 8] := st[4]/2+1/st[4]/2:
st[ 9] := (t[1]+1/t[1])/2:
st[10] := (t[2]+1/t[2])/2:
st[11] := (t[3]+1/t[3])/2:
st[12] := (t[4]+1/t[4])/2:
st[13] := t[1]*t[2] - t[3]*t[4]:
st[14] := t[1]/t[2] - t[3]/t[4]:
st[15] := t[2]/t[1] - t[4]/t[3]:
st[16] := 1/(t[1]*t[2]) - 1/(t[3]*t[4]):

# Expressions for s[i] in terms of y[1], y[2], r[1], r[2]
#@ sx
sx[ 1] := (r[1] + 1/sqrt(2))/(r[1] - 1/sqrt(2));
sx[ 2] := (r[2] + 1/sqrt(2))/(r[2] - 1/sqrt(2));
sx[ 3] := (1 + r[1]*y[1])/(1 - r[1]*y[1]);
sx[ 4] := (1 + r[2]*y[1])/(1 - r[2]*y[1]);
sx[ 5] := (1 + y[2]*sqrt(2))/(1 - y[2]*sqrt(2));
sx[ 6] := (1 - y[2]*sqrt(2))/(1 + y[2]*sqrt(2));
sx[ 7] := (1 + y[1]^2*(1 - y[2]/sqrt(2)))/(1 - y[1]^2*(1 - y[2]/sqrt(2)));
sx[ 8] := (1 + y[1]^2*(1 + y[2]/sqrt(2)))/(1 - y[1]^2*(1 + y[2]/sqrt(2)));
sx[ 9] := 2*r[1]*x[1]*(1-y[1]/sqrt(2))/((1 - y[1]^2*(1-y[2]/sqrt(2)))*(1-y[2]*sqrt(2)));
sx[10] := 2*r[2]*x[2]*(1-y[1]/sqrt(2))/((1 - y[1]^2*(1+y[2]/sqrt(2)))*(1+y[2]*sqrt(2)));
sx[11] := 2*r[1]*x[1]*(1+y[1]/sqrt(2))/((1 - y[1]^2*(1-y[2]/sqrt(2)))*(1-y[2]*sqrt(2)));
sx[12] := 2*r[2]*x[2]*(1+y[1]/sqrt(2))/((1 - y[1]^2*(1+y[2]/sqrt(2)))*(1+y[2]*sqrt(2)));
sx[13] := -2*y[1]*(r[1]+1/sqrt(2))*(r[2]+1/sqrt(2))*(r[1]+r[2])/(x[1]*x[2]);
sx[14] := -2*y[1]*(r[1]+1/sqrt(2))*(r[2]-1/sqrt(2))*(r[1]-r[2])/(x[1]*x[2]);
sx[15] :=  2*y[1]*(r[1]-1/sqrt(2))*(r[2]+1/sqrt(2))*(r[1]-r[2])/(x[1]*x[2]);
sx[16] :=  2*y[1]*(r[1]-1/sqrt(2))*(r[2]-1/sqrt(2))*(r[1]+r[2])/(x[1]*x[2]);

# Another family of elements with nice transformation properties
#@ ut
ut[ 1] := t[1]+t[3]:
ut[ 5] := t[1]-t[3]:
ut[ 9] := t[1]-1/t[3]:
ut[13] := t[1]+1/t[3]:
for i from 0 to 3 do
 ut[4*i+2] := act_AT[L](ut[4*i+1]);
 ut[4*i+3] := act_AT[A1](ut[4*i+1]);
 ut[4*i+4] := act_AT[LA1](ut[4*i+1]);
od:

for i from 1 to 16 do
 ux[i] := FNF_y0(eval(subs(t = tx,ut[i])));
od:

#@ tu
tu[ 1] := u[1]/2 + u[5]/2;
tu[ 2] := u[2]/2 + u[6]/2;
tu[ 3] := u[1]/2 - u[5]/2;
tu[ 4] := u[2]/2 - u[6]/2;

#@ su
su[ 1] := - u[ 1]/u[ 3];
su[ 2] := - u[ 2]/u[ 4];
su[ 3] := - u[15]/u[13];
su[ 4] := - u[16]/u[14];
su[ 5] :=   u[ 9]*u[11]/2 - 1;
su[ 6] :=   u[10]*u[12]/2 - 1;
su[ 7] := - u[ 1]*u[ 3]/2 - 1;
su[ 8] := - u[ 2]*u[ 4]/2 - 1;
su[ 9] :=   ( u[ 1] - u[11])/2;
su[10] :=   ( u[ 2] - u[12])/2;
su[11] :=   (-u[ 3] + u[11])/2;
su[12] :=   (-u[ 4] + u[12])/2;
su[13] :=   (u[ 1]*u[ 6]+u[ 2]*u[ 5])/2;
su[14] := - (u[ 4]*u[ 5]+u[ 1]*u[ 8])/2;
su[15] := - (u[ 3]*u[ 6]+u[ 2]*u[ 7])/2;
su[16] :=   (u[ 3]*u[ 8]+u[ 4]*u[ 7])/2;

#@ au
au[1] := factor(1/(1/u[1]+1/u[3]));
au[2] := factor(1/(1/u[6]+1/u[8]));

#@ bu
bu[1] := u[1]*u[3];
bu[2] := u[6]*u[8];

#@ bs
bs[1] := -2*(s[7]+1);
bs[2] :=  2*(s[8]-1);

# Expressions for x[i], y[i], z[i], r[i] in terms of s[j]
zs[1]     := (s[7]*s[8]-1)/((1+s[7])*(1+s[8]));            #@ zs
ys[1]     := sqrt(2)*(s[3]-1)/(s[3]+1)*(s[1]-1)/(s[1]+1);  #@ ys
ys_alt[1] := sqrt(2)*(s[4]-1)/(s[4]+1)*(s[2]-1)/(s[2]+1);  #@ ys_alt
ys[2]     := sqrt(2)*(s[8]-s[7])/(s[7]*s[8]-1);
xs[1] := 2*sqrt(2)*t[3]/((s[1]-1)*(s[3]+1));               #@ xs
xs[2] := 2*sqrt(2)*t[4]/((s[2]-1)*(s[4]+1));
xs[3] := ys[1];
xs[4] := -ys[1]*ys[2];
rs[1] := (s[1]+1)/(s[1]-1)/sqrt(2);                        #@ rs
rs[2] := (s[2]+1)/(s[2]-1)/sqrt(2);

# Expressions for x[i], y[i], z[i], r[i] in terms of t[j]
zt[1]     := eval(subs(s = st,zs[1]    )):                 #@ zt
yt[1]     := eval(subs(s = st,ys[1]    )):                 #@ yt
yt_alt[1] := eval(subs(s = st,ys_alt[1])):                 #@ yt_alt
yt[2]     := eval(subs(s = st,ys[2]    )):
xt[1]     := eval(subs(s = st,xs[1]    )):                 #@ xt
xt[2]     := eval(subs(s = st,xs[2]    )):
xt[3]     := eval(subs(s = st,xs[3]    )):
xt[4]     := eval(subs(s = st,xs[4]    )):
rt[1]     := eval(subs(s = st,rs[1]    )):                 #@ rt
rt[2]     := eval(subs(s = st,rs[2]    )):

# Relations among the elements t[i]
#@ t_rels
t_rels[ 1] := t[1]^2*t[2]^2*t[3]^2*t[4]^2-4*t[1]^2*t[2]*t[3]^2*t[4]-4*t[1]*t[2]^2*t[3]*t[4]^2+t[1]^2*t[3]^2+12*t[1]*t[2]*t[3]*t[4]+t[2]^2*t[4]^2-4*t[1]*t[3]-4*t[2]*t[4]+1;
t_rels[ 2] := t[1]^2*t[2]*t[3]*t[4]^2-t[1]*t[2]^2*t[3]^2*t[4]+t[1]^2*t[2]*t[3]-t[1]*t[2]^2*t[4]-t[1]*t[3]^2*t[4]+t[2]*t[3]*t[4]^2-t[1]*t[4]+t[2]*t[3];

#@ ts_rels
ts_rels[1] := tan_sum(t[1], t[3]) + ((t[1]+t[3])*(r[1]-1/sqrt(2))/sqrt(2));
ts_rels[2] := tan_sum(t[2], t[4]) + ((t[2]+t[4])*(r[2]-1/sqrt(2))/sqrt(2));
ts_rels[3] := tan_sum(t[1],-t[3]) - ((t[1]-t[3])*(r[1]-1/sqrt(2))/r[1]/2);
ts_rels[4] := tan_sum(t[2],-t[4]) - ((t[2]-t[4])*(r[2]-1/sqrt(2))/r[2]/2);

# Relations among the elements s[i]
#@ s_rels
s_rels[ 1] := s[5]*s[6]-1;
s_rels[ 2] := s[5]*((s[7]-2)*(s[8]+2)+3) - ((s[7]+2)*(s[8]-2)+3);
s_rels[ 3] := s[6]*((s[7]+2)*(s[8]-2)+3) - ((s[7]-2)*(s[8]+2)+3);
s_rels[ 4] := s[9]/s[11]-s[10]/s[12];
s_rels[ 5] := s[1]*s[2]*s[3]-s[1]*s[2]*s[4]+s[1]*s[3]*s[4]-s[2]*s[3]*s[4]-s[1]+s[2]-s[3]+s[4];
s_rels[ 6] := s[1]^2*s[2]^2-4*s[1]^2*s[2]-4*s[1]*s[2]^2+s[1]^2+12*s[1]*s[2]+s[2]^2-4*s[1]-4*s[2]+1;
s_rels[ 7] := s[13]*s[14]*s[15] + s[13]*s[14]*s[16] + s[13]*s[15]*s[16] + s[14]*s[15]*s[16];
s_rels[ 8] := (s[13]-s[14]-s[15]+s[16])^2 + 8*(s[13]+s[14]+s[15]+s[16])^2 - 8*(s[13]-s[16])^2 - 8*(s[14]-s[15])^2;
s_rels[ 9] := ((1-4*s[1]+s[1]^2)*(s[2]-2) - 2*s[1])^2 - (1 - s[1])^2*(3*s[1]-1)*(s[1]-3);
s_rels[10] := ((1-4*s[2]+s[2]^2)*(s[1]-2) - 2*s[2])^2 - (1 - s[2])^2*(3*s[2]-1)*(s[2]-3);
s_rels[11] := s[1]*s[2]*s[16] + s[13];
s_rels[12] := s[1]*s[15] + s[2]*s[14];

# Quadratic relations for t[i] over K
#@ qt_rels
qt_rels[1] := t[1]^2 -2*sqrt(2)*x[1]*(1-sqrt(2)*y[1]*(1-y[2]/sqrt(2)))/((1-y[2]*sqrt(2))*(1-y[1]^2*(1-y[2]/sqrt(2)))) * t[1] - 1;
qt_rels[2] := t[2]^2 -2*sqrt(2)*x[2]*(1-sqrt(2)*y[1]*(1+y[2]/sqrt(2)))/((1+y[2]*sqrt(2))*(1-y[1]^2*(1+y[2]/sqrt(2)))) * t[2] - 1;
qt_rels[3] := t[3]^2 -2*sqrt(2)*x[1]*(1+sqrt(2)*y[1]*(1-y[2]/sqrt(2)))/((1-y[2]*sqrt(2))*(1-y[1]^2*(1-y[2]/sqrt(2)))) * t[3] - 1;
qt_rels[4] := t[4]^2 -2*sqrt(2)*x[2]*(1+sqrt(2)*y[1]*(1+y[2]/sqrt(2)))/((1+y[2]*sqrt(2))*(1-y[1]^2*(1+y[2]/sqrt(2)))) * t[4] - 1;

# The relations u_lin_rels[i] and u_quad_rels[i] depend only on the expressions 
# for the u[i] in terms of the t[j], not on the relations among the t[j]

#@ u_lin_rels
u_lin_rels := {
u[ 1] + u[ 3] - u[ 9] - u[11],
u[ 1] + u[ 5] - u[ 9] - u[13],
u[ 1] + u[ 7] - u[11] - u[13],
u[ 2] + u[ 4] - u[10] - u[12],
u[ 2] + u[ 6] - u[10] - u[14],
u[ 2] + u[ 8] - u[12] - u[14],
u[ 3] + u[ 5] - u[ 9] - u[15],
u[ 3] + u[ 7] - u[11] - u[15],
u[ 4] + u[ 6] - u[10] - u[16],
u[ 4] + u[ 8] - u[12] - u[16],
u[ 5] + u[ 7] - u[13] - u[15],
u[ 6] + u[ 8] - u[14] - u[16]
}:

#@ u_quad_rels
u_quad_rels := {
 u[ 1]*u[ 3] + u[ 5]*u[ 7] + 4,
 u[ 2]*u[ 4] + u[ 6]*u[ 8] + 4,
 u[ 9]*u[11] + u[13]*u[15] + 4,
 u[10]*u[12] + u[14]*u[16] + 4,
 u[ 1]*u[ 7] + u[ 3]*u[ 5],
 u[ 2]*u[ 8] + u[ 4]*u[ 6],
 u[ 9]*u[15] + u[11]*u[13],
 u[10]*u[16] + u[12]*u[14]
}:

# These additional relations depend on relations among the t[j].
#@ u_rels
u_rels[1] := 1/(u[9]*u[11]) + 1/(u[10]*u[12]) - 1/2;
u_rels[2] := u[1]*u[4]-u[4]*u[11]+u[11]*u[2]-u[2]*u[3]+u[3]*u[12]-u[12]*u[1];
u_rels[3] := (u[9]*u[11]+u[10]*u[12]-8)*(1/2+1/(u[1]*u[3])+1/(u[2]*u[4]))+2*(u[9]*u[11]-u[10]*u[12])*(1/(u[1]*u[3])-1/(u[2]*u[4]));

#@ ab_rels
ab_rels[1] := expand((b[1]+2*a[1]^2)*(b[2]-6*a[2]^2) + 4*a[1]^2*a[2]^2);
ab_rels[2] := 8*a[2]^2*b[1]-b[1]*b[2]+2*b[1]+2*b[2]+8;
ab_rels[3] := ((1+a[2]^2)*(3+2*a[1]^2+b[1])-  (1+a[1]^2))^2 - (a[1]^2-3*a[2]^2-2)^2 - 4*a[1]^2*a[2]^2*(1+a[1]^2)*(1+a[2]^2);
ab_rels[4] := ((1+a[1]^2)*(1+6*a[2]^2-b[2])-3*(1+a[2]^2))^2 - (a[1]^2-3*a[2]^2-2)^2 - 4*a[1]^2*a[2]^2*(1+a[1]^2)*(1+a[2]^2);
ab_rels[5] := (1+a[1]^2)*b[2] + (1+a[2]^2)*b[1] + 4*(1 - a[1]^2*a[2]^2);
ab_rels[6] := a[1]^2 + (1/4)*b[1]*(b[1]*b[2]+6*b[1]+6*b[2]+24)/(b[1]*b[2]+2*b[1]+2*b[2]+8);
ab_rels[7] := a[2]^2 - (1/8)*(b[1]*b[2]-2*b[1]-2*b[2]-8)/b[1];

#@ extra_galois_rels
extra_galois_rels := [
 tan_sum(t[1], t[3]) - 1/(1/u[1] + 1/u[3]),
 tan_sum(t[2],-t[4]) - 1/(1/u[6] + 1/u[8]),
 u[6]*u[8] - 2*(s[8]-1),
 (u[6]-u[8])^2 -  8*(s[8]-1)*(s[7]*s[8]-1)/(s[7]*s[8]-2*s[7]+2*s[8]-1),
 (u[6]+u[8])^2 - 16*(s[8]-1)^2*(s[7]+1)   /(s[7]*s[8]-2*s[7]+2*s[8]-1)
];

#@ s_to_t 
s_to_t := {seq(s[i]=st[i],i=1..16),seq(u[i]=ut[i],i=1..16)}:

#@ t_to_xyz 
t_to_xyz := {seq(t[i]=tx[i],i=1..4),seq(s[i]=sx[i],i=1..16),seq(u[i]=ux[i],i=1..16),r[1]=ry[1],r[2]=ry[2]}:

#@ t_reduce 
t_reduce := (a) -> FNF_y0(factor(subs(root_rule,expand(rationalize(simplify(a))))));

#@ t_xyz_reduce 
t_xyz_reduce := (a) -> t_reduce(subs(t_to_xyz,a));

# SNB is a basis for KT over R(s[7],s[8]).  
# The function VNB converts its argument to the list of coefficients with respect to the basis SNB

#@ SNB 
SNB := [seq(seq(seq(seq(seq(s[1]^i*s[2]^j*s[3]^k*t[1]^l*t[2]^m,m=0..1),l=0..1),k=0..1),j=0..1),i=0..1)];

#@ VNB 
VNB := (u) -> [seq(seq(seq(seq(seq(
               coeff(coeff(coeff(coeff(coeff(u,
               s[1],i),s[2],j),s[3],k),t[1],l),t[2],m),
               m=0..1),l=0..1),k=0..1),j=0..1),i=0..1)];

# The following relations can be used to express things in terms of SNB

#@ sn_rels

sn_rels[1] := s[1]^2 - ((2*(3*s[7]*s[8]+2*s[7]-2*s[8]-3))/(s[7]*s[8]+2*s[7]-2*s[8]-1))*s[1] + 1:
sn_rels[2] := s[2]^2 - ((2*(3*s[7]*s[8]-2*s[7]+2*s[8]-3))/(s[7]*s[8]-2*s[7]+2*s[8]-1))*s[2] + 1:
sn_rels[3] := s[3]^2 - 2*s[7]*s[3] + 1:
sn_rels[4] := s[4] - ((1/8)*(s[1]*s[2]*s[3]*s[7]^2*s[8]^2-s[1]*s[2]*s[7]^3*s[8]^2-3*s[1]*s[3]*s[7]^2*s[8]^2+3*s[1]*s[7]^3*s[8]^2-3*s[2]*s[3]*s[7]^2*s[8]^2+3*s[2]*s[7]^3*s[8]^2-4*s[1]*s[2]*s[3]*s[7]^2+6*s[1]*s[2]*s[3]*s[7]*s[8]-4*s[1]*s[2]*s[3]*s[8]^2+4*s[1]*s[2]*s[7]^3-6*s[1]*s[2]*s[7]^2*s[8]+4*s[1]*s[2]*s[7]*s[8]^2-4*s[1]*s[3]*s[7]^2*s[8]+4*s[1]*s[3]*s[7]*s[8]^2+4*s[1]*s[7]^3*s[8]-4*s[1]*s[7]^2*s[8]^2+4*s[2]*s[3]*s[7]^2*s[8]-4*s[2]*s[3]*s[7]*s[8]^2-4*s[2]*s[7]^3*s[8]+4*s[2]*s[7]^2*s[8]^2+9*s[3]*s[7]^2*s[8]^2-s[7]^3*s[8]^2+4*s[1]*s[3]*s[7]^2-2*s[1]*s[3]*s[7]*s[8]+4*s[1]*s[3]*s[8]^2-4*s[1]*s[7]^3+2*s[1]*s[7]^2*s[8]-4*s[1]*s[7]*s[8]^2+4*s[2]*s[3]*s[7]^2-2*s[2]*s[3]*s[7]*s[8]+4*s[2]*s[3]*s[8]^2-4*s[2]*s[7]^3+2*s[2]*s[7]^2*s[8]-4*s[2]*s[7]*s[8]^2+s[1]*s[2]*s[3]-s[1]*s[2]*s[7]+4*s[1]*s[3]*s[7]-4*s[1]*s[3]*s[8]-4*s[1]*s[7]^2+4*s[1]*s[7]*s[8]-4*s[2]*s[3]*s[7]+4*s[2]*s[3]*s[8]+4*s[2]*s[7]^2-4*s[2]*s[7]*s[8]-4*s[3]*s[7]^2-10*s[3]*s[7]*s[8]-4*s[3]*s[8]^2+4*s[7]^3+2*s[7]^2*s[8]-4*s[7]*s[8]^2-3*s[1]*s[3]+3*s[1]*s[7]-3*s[2]*s[3]+3*s[2]*s[7]+9*s[3]-9*s[7]+8*s[8])/((s[7]-1)*(s[7]+1)*(s[7]*s[8]-1))):
sn_rels[5] := t[1]^2 - s[1]*(2*s[7] - s[3]):
sn_rels[6] := t[2]^2-((1/8)*(-3*s[1]*s[2]*s[3]*s[7]^2*s[8]^2+3*s[1]*s[2]*s[7]^3*s[8]^2-4*s[1]*s[2]*s[3]*s[7]^2*s[8]+4*s[1]*s[2]*s[3]*s[7]*s[8]^2+4*s[1]*s[2]*s[7]^3*s[8]-4*s[1]*s[2]*s[7]^2*s[8]^2+s[1]*s[3]*s[7]^2*s[8]^2-s[1]*s[7]^3*s[8]^2+9*s[2]*s[3]*s[7]^2*s[8]^2-s[2]*s[7]^3*s[8]^2+4*s[1]*s[2]*s[3]*s[7]^2-2*s[1]*s[2]*s[3]*s[7]*s[8]+4*s[1]*s[2]*s[3]*s[8]^2-4*s[1]*s[2]*s[7]^3+2*s[1]*s[2]*s[7]^2*s[8]-4*s[1]*s[2]*s[7]*s[8]^2-3*s[3]*s[7]^2*s[8]^2+3*s[7]^3*s[8]^2+4*s[1]*s[2]*s[3]*s[7]-4*s[1]*s[2]*s[3]*s[8]-4*s[1]*s[2]*s[7]^2+4*s[1]*s[2]*s[7]*s[8]-4*s[1]*s[3]*s[7]^2+6*s[1]*s[3]*s[7]*s[8]-4*s[1]*s[3]*s[8]^2+4*s[1]*s[7]^3-6*s[1]*s[7]^2*s[8]+4*s[1]*s[7]*s[8]^2-4*s[2]*s[3]*s[7]^2-10*s[2]*s[3]*s[7]*s[8]-4*s[2]*s[3]*s[8]^2+4*s[2]*s[7]^3+2*s[2]*s[7]^2*s[8]-4*s[2]*s[7]*s[8]^2+4*s[3]*s[7]^2*s[8]-4*s[3]*s[7]*s[8]^2-4*s[7]^3*s[8]+4*s[7]^2*s[8]^2-3*s[1]*s[2]*s[3]+3*s[1]*s[2]*s[7]+4*s[3]*s[7]^2-2*s[3]*s[7]*s[8]+4*s[3]*s[8]^2-4*s[7]^3+2*s[7]^2*s[8]-4*s[7]*s[8]^2+s[1]*s[3]-s[1]*s[7]+9*s[2]*s[3]-9*s[2]*s[7]+8*s[2]*s[8]-4*s[3]*s[7]+4*s[3]*s[8]+4*s[7]^2-4*s[7]*s[8]-3*s[3]+3*s[7])/((s[7]-1)*(s[7]+1)*(s[7]*s[8]-1))):
sn_rels[7] := t[3] - t[1]*s[3]:
sn_rels[8] := t[4] - ((1/8)*t[2]*(s[1]*s[2]*s[3]*s[7]^2*s[8]^2-s[1]*s[2]*s[7]^3*s[8]^2-3*s[1]*s[3]*s[7]^2*s[8]^2+3*s[1]*s[7]^3*s[8]^2-3*s[2]*s[3]*s[7]^2*s[8]^2+3*s[2]*s[7]^3*s[8]^2-4*s[1]*s[2]*s[3]*s[7]^2+6*s[1]*s[2]*s[3]*s[7]*s[8]-4*s[1]*s[2]*s[3]*s[8]^2+4*s[1]*s[2]*s[7]^3-6*s[1]*s[2]*s[7]^2*s[8]+4*s[1]*s[2]*s[7]*s[8]^2-4*s[1]*s[3]*s[7]^2*s[8]+4*s[1]*s[3]*s[7]*s[8]^2+4*s[1]*s[7]^3*s[8]-4*s[1]*s[7]^2*s[8]^2+4*s[2]*s[3]*s[7]^2*s[8]-4*s[2]*s[3]*s[7]*s[8]^2-4*s[2]*s[7]^3*s[8]+4*s[2]*s[7]^2*s[8]^2+9*s[3]*s[7]^2*s[8]^2-s[7]^3*s[8]^2+4*s[1]*s[3]*s[7]^2-2*s[1]*s[3]*s[7]*s[8]+4*s[1]*s[3]*s[8]^2-4*s[1]*s[7]^3+2*s[1]*s[7]^2*s[8]-4*s[1]*s[7]*s[8]^2+4*s[2]*s[3]*s[7]^2-2*s[2]*s[3]*s[7]*s[8]+4*s[2]*s[3]*s[8]^2-4*s[2]*s[7]^3+2*s[2]*s[7]^2*s[8]-4*s[2]*s[7]*s[8]^2+s[1]*s[2]*s[3]-s[1]*s[2]*s[7]+4*s[1]*s[3]*s[7]-4*s[1]*s[3]*s[8]-4*s[1]*s[7]^2+4*s[1]*s[7]*s[8]-4*s[2]*s[3]*s[7]+4*s[2]*s[3]*s[8]+4*s[2]*s[7]^2-4*s[2]*s[7]*s[8]-4*s[3]*s[7]^2-10*s[3]*s[7]*s[8]-4*s[3]*s[8]^2+4*s[7]^3+2*s[7]^2*s[8]-4*s[7]*s[8]^2-3*s[1]*s[3]+3*s[1]*s[7]-3*s[2]*s[3]+3*s[2]*s[7]+9*s[3]-9*s[7]+8*s[8])/((s[7]-1)*(s[7]+1)*(s[7]*s[8]-1))):

# This is a basis for KT over R(z[1],z[2])
#@ BXYR 
BXYR := [seq(seq(seq(seq(seq(seq(
         x[1]^i*x[2]^j*y[1]^k*y[2]^l*r[1]^m*r[2]^n,
         i=0..1),j=0..1),k=0..1),l=0..1),m=0..1),n=0..1)];

# The action of G[64] on KT preserves the above basis up to sign, with orbits as follows.
# These are enumerated so that the first 16 orbits have size one, and the actions on these
# orbits give the 16 one-dimensional characters of G[64].  The remaining 24 orbits all have 
# size two, and so give two-dimensional characters of G[64].  These are all irreducible, and
# each character occurs twice.  Orbits with the same character are adjacent in the list.

#@ BXYR_orbits 
BXYR_orbits := [
 seq(seq(seq(seq([y[2]^i*y[1]^j*(x[1]*x[2])^k*(r[1]*r[2])^l],i=0..1),j=0..1),k=0..1),l=0..1),
 seq(seq(seq(y[2]^i*y[1]^j*(r[1]*r[2])^l*~[x[1],x[2]],i=0..1),j=0..1),l=0..1),
 seq(seq(seq(y[2]^i*y[1]^j*(x[1]*x[2])^l*~[r[1],r[2]],i=0..1),j=0..1),l=0..1),
 seq(seq(y[2]^i*y[1]^j*~[r[1]*x[1],r[2]*x[2]],i=0..1),j=0..1),
 seq(seq(y[2]^i*y[1]^j*~[r[1]*x[2],r[2]*x[1]],i=0..1),j=0..1)
];

# Characters of G64
#@ orb_char 
orb_char := proc(o)
 local c,T,eqs1,eqs2,sol1,sol2,p,q;
 if nops(o) = 1 then
  return map(T -> factor(act_AT[T](o[1])/o[1]),G64);
 fi;

 c := [];
 for T in G64 do 
  eqs1 := {coeffs(act_AT[T](o[1]) - p * o[1] - q*o[2],{x[1],x[2],y[1],y[2],r[1],r[2]})};
  eqs2 := {coeffs(act_AT[T](o[2]) - p * o[1] - q*o[2],{x[1],x[2],y[1],y[2],r[1],r[2]})};
  sol1 := solve(eqs1,{p,q});
  sol2 := solve(eqs2,{p,q});
  c := [op(c),subs(sol1,p) + subs(sol2,q)];
 od;

 return c;
end:

#@ BXYR_chars 
BXYR_chars := map(orb_char,BXYR_orbits):

#@ G64_chars 
G64_chars := [
 seq(BXYR_chars[i],i=1..16),
 seq(BXYR_chars[16+2*i],i=1..12)
];

#@ iso_AT 
iso_AT := (a,i) -> add(act_AT[G64[j]](a) * G64_chars[i][j] * G64_chars[i][1]/64,j = 1..64);

galois_check_zero := proc(a,n_)
 local n,x0,y0,z0,t0,s0,r0,u0;
 n := `if`(nargs > 1, n_, 26);
 x0 := inner_quasirational_points[n];
 y0 := combine(simplify(y_proj0(x0)));
 z0 := y0 ^~ 2;
 t0 := combine(simplify(expand(E_to_TTP(x0))));
 s0 := combine(simplify(expand(rationalize(eval(subs(t = t0,[seq(st[i],i=1..16)]))))));
 u0 := combine(simplify(expand(rationalize(eval(subs(t = t0,[seq(ut[i],i=1..16)]))))));
 r0[1] := combine(simplify(expand(sqrt(1 - y0[2]/sqrt(2)))));
 r0[2] := combine(simplify(expand(sqrt(1 + y0[2]/sqrt(2)))));
 return evalb(combine(simplify(eval(subs({x=x0,y=y0,z=z0,t=t0,s=s0,r=r0,u=u0},a)))) = 0);
end: