\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ This program provides the compatational steps to prove Theorem 4.2


default(breakloop,0)
\\ fix of the built-in lift(x), which has a small bug
LIFT(x) =
  { local(t = type(x));

    while (t == "t_INTMOD" || t == "t_POLMOD", x = lift(x); t = type(x));
    x;
  }

z; zet1;zet2; X; r2;
\\r2;	 
r1=1;

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ CONSTRUCTION OF THE EQUATION IN ROOTS OF UNITY

om = r1/(z-xx)-r1/(z-y)+r2/(z-u)-r2/(z+u)-r1/(z+xx)+r1/(z+y)
den = (z^2-xx^2)*(z^2-y^2)*(z^2 - u^2)
nu = om*den
E1= polcoeff(nu,0,z)/xx/y/u
\\polcoeff(nu,1,z)
\\=0, hence no condition
\\polcoeff(nu,3,z)
\\=0, hence no condition
E2=(polcoeff(nu,4,z)+polcoeff(nu,2,z))/2

RE = polresultant(E1,E2,xx)

Y=zet1*(1-xx)/(1+xx);

E1s = subst(subst(subst(E1,xx,(1-X)/(1+X)),y,(1-zet1*X)/(1+zet1*X)), u, (1-zet2)/(1+zet2));
E2s = subst(subst(subst(E2,xx,(1-X)/(1+X)),y,(1-zet1*X)/(1+zet1*X)), u, (1-zet2)/(1+zet2));
\\ clear denominators

E1s = (1+X)*(1+zet1*X)*(zet2+1)*E1s
E2s = (1+X)^2*(1+zet1*X)^2*(zet2+1)^2*E2s;

PP = E2s-polcoeff(E2s,4,X)/(polcoeff(E1s,2,X))*X^2*E1s;
RES =polresultant(PP,E1s,X)/2048/r2/zet1/(zet2+1);
SQRES = zet1*(zet2+1)^2*(zet2-1)*r2^2 - (zet1^2-1)*(zet2+1)*zet2*r2 + 2*zet2*(zet2-1)*(zet1-1)^2;
SQRES^2==RES

\\ the equation appearing in Theorem 4.2
SQRES_paper = r2^2*zet1*zet2^3 - r2^2*zet1 + (4-r2^2)*zet1*zet2 + (r2^2-4)*zet1*zet2^2 + (2-r2)*zet1^2*zet2^2 + (r2+2)*zet2^2 - (r2+2)*zet1^2*zet2 + (r2-2)*zet2;
\\ Check that the equation is as written in the text
print("same equation: ", SQRES == SQRES_paper) \\ = 1, is OK.

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ List of exponents of the equation
Exps = [1,3;1,0;1,1;1,2;2,2;0,2;2,1;0,1]
for(i=1,8,for(j=i+1,8,for(k=j+1,8,M = [Exps[j,1]-Exps[i,1],Exps[k,1]-Exps[i,1];Exps[j,2]-Exps[i,2],Exps[k,2]-Exps[i,2]]; dt= matdet(M); if(abs(dt)==0,print(i,j,k,"  dt=",dt,)))))
\\ List of bad triples
\\124  dt=0
\\134  dt=0
\\234  dt=0
\\378  dt=0
\\456  dt=0
\\only those are bad as claimed in the proof of Theorem 4.2


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Check the sub-relations of length two.
E12 = subst(polresultant(SQRES,zet2^3-1,zet2),zet1,x)/(x^2-1)/r2/2
\\ (r2^2 + 12)*x^4 + (12*r2^2 - 48)*x^3 + (3*r2^4 - 26*r2^2 + 72)*x^2 + (12*r2^2 - 48)*x + (r2^2 + 12)
E58 = - subst(subst(SQRES,zet2,1/zet1^2),zet1,x)*x^5/(x^2-1)
\\ r2^2*x^4 + (r2 + 2)*x^3 + (2*r2^2 - 4)*x^2 + (r2 + 2)*x + r2^2
E67 = subst(subst(SQRES,zet2,zet1^2),zet1,x)/x/(x^2-1)
subst(E67,r2,-r2) == E58 \\ true, i.e. solutions are covered by those of E58.

\\ check the REAL quadratic irrationals that annihilate coefficients
\\ none for the first equation E12
\\ for the second equation E58:
pr1 = polroots(subst(E58,r2,sqrt(2)))
tc = pr1[3]
tc0 = 1;
for(i=1,61,tc0*=tc; if(abs(tc0-1)<0.001,print(i," ",tc0),))
\\ this checks the claims in Lemma 4.5


\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\  Reduce to (8)-Relations
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
{for(i=1,8,for(j=i+1,8,for(k=j+1,8,
e = 2^(i-1)+2^(j-1)+2^(k-1); 
M1 = vecextract(Exps~,e);
M1r = matrix(2,2,s,t,M1[s,t+1]-M1[s,1]);
d1 = matdet(M1r);
M2 = vecextract(Exps~,255-e);
M2r = matrix(2,4,s,t,M2[s,t+1]-M2[s,1]);
d=d1;
for(l=1,4,for(m=l+1,4,MM = matrix(2,2); MM[,1]= M2r[,l]; MM[,2] = M2r[,m];
  d=gcd(d,matdet(MM))));
if(d==1,,print(i,j,k,"  d=",d,));
)))}
\\ prints nothing: i.e. no problematic case


{gcm(M,r)=local(g); g=0; for(l=1,r,for(m=l+1,r,MM = matrix(2,2); MM[,1]= M[,l]; MM[,2] = M[,m];
  g=gcd(g,matdet(MM)))); g;}

\\ Case (4,4)

{for(i=1,8,for(j=i+1,8,for(k=j+1,8,for(n=k+1,8,
e = 2^(i-1)+2^(j-1)+2^(k-1)+2^(n-1); 
M1 = vecextract(Exps~,e);
M1r = matrix(2,3,s,t,M1[s,t+1]-M1[s,1]);
M2 = vecextract(Exps~,255-e);
M2r = matrix(2,3,s,t,M2[s,t+1]-M2[s,1]);
d=gcd(gcm(M1r,3),gcm(M2r,3));
if(d==1,,print(i,j,k,n,"  d=",d,));
))))}
\\1234  d=2
\\1356  d=2
\\1378  d=2
\\2456  d=2
\\2478  d=2
\\5678  d=2
\\As claimed for this case

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ check the remaining (4,4)-relations

R1A = r2^2*zet1*zet2^3  + (4-r2^2)*zet1*zet2 + (2-r2)*zet1^2*zet2^2 + (r2+2)*zet2^2;
R1B = - r2^2*zet1 + (r2^2-4)*zet1*zet2^2   - (r2+2)*zet1^2*zet2 + (r2-2)*zet2;

polresultant(R1A,R1B,zet1) == (16*zet2^4 * (r2-2)*(r2+2)*(r2-1)*(r2+1)*(zet2-1)^2*(zet2+1)^2)
\\ complete factorization into rational factors or  zet2 = \pm 1, all excluded.

R2A = r2^2*zet1*zet2^3 + (4-r2^2)*zet1*zet2  - (r2+2)*zet1^2*zet2 + (r2-2)*zet2;

R2B = - r2^2*zet1 + (r2^2-4)*zet1*zet2^2 + (2-r2)*zet1^2*zet2^2 + (r2+2)*zet2^2

polresultant(R2A,R2B,zet1) == -(zet2^2*(zet2-1)^4*(zet2+1)^4*r2^4*(r2-2)*(r2+2))
\\ complete factorization into rational factors or  zet2 = \pm 1, all excluded.


R3A = r2^2*zet1*zet2^3 - r2^2*zet1 + (4-r2^2)*zet1*zet2 + (r2^2-4)*zet1*zet2^2;
R3B =  (2-r2)*zet1^2*zet2^2 + (r2+2)*zet2^2 - (r2+2)*zet1^2*zet2 + (r2-2)*zet2;
polresultant(R3A,R3B,zet1) == (zet2*(zet2-1)^2 * (r2*(zet2+1) + 2*zet2-2)*(r2^2 + (2*r2^2 - 4)*zet2 + r2^2*zet2^2)^2)

Re3z2 = polresultant(R3A,R3B,zet2)/(8*zet1^2*(zet1^2-1)*r2^3)
\\ gives a degree three relation, irreducible unless r2=\sqrt{4/3}
\\ and \zet1 = 8-th root of unity.

subst(%,r2,sqrt(2))

\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Loop through the  case (8)

qr = r2^2 + be*r2 + ga;
PR = polresultant(qr,SQRES,r2);

\\\\\\\

\\ cross-check that the PRij pick all the coefficients from PR
PR0 = polcoeff(PR,0,be);
PR1 = polcoeff(PR,1,be);
PR2 = polcoeff(PR,2,be);
PR3 = polcoeff(PR,3,be);
PR00 = polcoeff(PR0,0,ga);
PR01 = polcoeff(PR0,1,ga);
PR02 = polcoeff(PR0,2,ga);
PR10 = polcoeff(PR1,0,ga);
PR11 = polcoeff(PR1,1,ga);
PR12 = polcoeff(PR1,2,ga);
PR20 = polcoeff(PR2,0,ga);
PR21 = polcoeff(PR2,1,ga);
PR22 = polcoeff(PR2,2,ga);
{PR ==
PR00*ga^0*be^0 + 
PR01*ga^1*be^0 + 
PR02*ga^2*be^0 + 
PR10*ga^0*be^1 + 
PR11*ga^1*be^1 + 
PR12*ga^2*be^1 + 
PR20*ga^0*be^2 + 
PR21*ga^1*be^2 + 
PR22*ga^2*be^2} \\ =1
\\\\\\\\\\\\\ End cross-check

\\ remove bad factors, since these solutions correspond to rational r (instead of quadratic)
rbf(h) = if(subst(h,ga,2*be-4)==0,h =h/(2*be-ga-4),); if(subst(h,ga,-2*be-4)==0, h=h/(2*be+ga+4),); if(subst(h,ga,-be-1)==0, h=h/(be+ga+1),); if(subst(h,ga,be-1)==0,h=h/(be-ga-1),);h;
testsol(N,i,j,b,g) = rt = (-b+sqrt(b^2-4*g))/2;  norml2(subst(subst(subst(SQRES,r2,rt), zet1,exp(2*Pi*I*i/N)),zet2,exp(2*Pi*I*j/N)))



Nstart = [1,5,7,11,13,17,5*7,5*11,7*11,5*13]
Nstart = concat(Nstart,3*Nstart)
NList = vecsort(concat(Nstart,concat(2*Nstart,concat(4*Nstart,8*Nstart))))


SolList = [];  \\ Solutions
CandList = []; \\ Candidates: 
DMList = [];   \\ Also candidates:

{
for(NN=1,length(NList), N=NList[NN]; print("Now working with N=",N); P = polcyclo(N,v); Np = poldegree(P);
for(i=1,N-1,print("current i=",i); for(j=1,N/2, \\ may restrict using N-e_xy, N-eU symmetry
if(gcd(i,gcd(j,N)) == 1,
if(j == N/2,, \\ exclude \zetA_U = -1
eqlist = [];
PR00s = subst(subst(PR00,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR01s = subst(subst(PR01,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR02s = subst(subst(PR02,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR10s = subst(subst(PR10,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR11s = subst(subst(PR11,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR12s = subst(subst(PR12,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR20s = subst(subst(PR20,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR21s = subst(subst(PR21,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));
PR22s = subst(subst(PR22,zet1,Mod(v^i,P)),zet2,Mod(v^j,P));

PRs = LIFT(PR00s)*ga^0*be^0 + LIFT(PR01s)*ga^1*be^0 + LIFT(PR02s)*ga^2*be^0 + LIFT(PR10s)*ga^0*be^1 + LIFT(PR11s)*ga^1*be^1 + LIFT(PR12s)*ga^2*be^1 + LIFT(PR20s)*ga^0*be^2 + LIFT(PR21s)*ga^1*be^2 + LIFT(PR22s)*ga^2*be^2;
print("PRs: ",PRs);

for(k=0,Np, if(polcoeff(PRs,k,v) ==0,, eqlist = concat(eqlist,rbf(polcoeff(PRs,k,v)))));
if(length(eqlist) <=1, print("Cand:",i," ",j," ",PRs); CandList=concat(CandList,[[N,i,j]]),
prv = vector(length(eqlist)-1,i,polresultant(eqlist[1],eqlist[i+1],be));
nz = 0; \\ till now no non-zero resultant found
for(k=1,length(prv),if(prv[k]==0,,nz=1; FA=factor(prv[k])~; 
     k = length(prv)+1); \\ end the loop, if a non-zero resultant was found
);
if(nz==1,
   for(l=1,length(FA),
   if(poldegree(FA[1,l],ga) == 1, \\ print(i,j," l=",l," Factor: ",FA[1,l]); 
      sga = -polcoeff(FA[1,l],0,ga)/polcoeff(FA[1,l],1,ga); \\ gamma to be substituted in
      seqlist = [];
      for(m=1,length(eqlist),s = subst(eqlist[m],ga,sga); if(s==0,,seqlist=concat(seqlist,[s])));
\\      print(i," ",j," seqlist: "seqlist);
      if(poldegree(seqlist[1],be)> 0, 
         FAB = factor(seqlist[1]);
         for(m=1,length(FAB~), 
            if(poldegree(FAB[m,1],be) == 1, 
               sbe = -polcoeff(FAB[m,1],0,be)/polcoeff(FAB[m,1],1,be); \\ beta to be substituted in
     	    issol=1;
     	    for(n=2,length(seqlist), if(subst(seqlist[n],be,sbe)==0,,issol=0; n=length(seqlist)+1));
     	    if(issol==1, 
                \\print("Eqsolution: N=",N," i="i," j=",j,"  ga=",sga," be=",sbe); 
                if((sga<0) && (testsol(N,i,j,sbe,sga)<0.0001), print("TRUESOLUTION: N=",N," i="i," j=",j,"  ga=",sga," be=",sbe); SolList = concat(SolList,[[N,i,j,sbe,sga]]),);
	    )
            , \\ do nothing with non-lin factors in FAB
            );
         ); \\ end for(m=..) 
      ,) \\ do nothing if the factor is constant (as poly in beta)
   , 
   \\ do not do anything with non-linear factors
   ) \\ end if
   ) \\ end for(l=)
   , 
   print(i," ",j," All resultants are zero, check by different method");
   DMList=concat(DMList,[[N,i,j]]);
); \\ end if nz==1
))
,)
)) \\ end for i,j looping [1,N]
) \\ end for N looping NList
}


rattest(N,i,j) = ps = subst(subst(SQRES, zet1,exp(2*Pi*I*i/N)),zet2,exp(2*Pi*I*j/N)); ps = ps/polcoeff(ps,2,r2); disc=poldisc(ps,r2); rdisc = round(10^5*disc)/10^5; issquare(rdisc)

SolList 
\\ = [[12, 6, 1, 4, -8], [[12, 6, 5, -4, -8],  [24, 3, 4, 0, -4/3], [24, 15, 4, 0, -4/3]]

CandList 
\\ = [[3, 2, 1], [6, 2, 1], [6, 3, 1], [6, 3, 2], [6, 5, 1], [6, 5, 2], [10, 2, 3], [10, 4, 1], [10, 6, 1], [10, 8, 3], [4, 1, 1], [4, 2, 1], [4, 3, 1]]

for(i=1,length(CandList),print(CandList[i]," has rationality test ", rattest(CandList[i][1],CandList[i][2],CandList[i][3])))
\\[3, 2, 1] has rationality test 1
\\[6, 2, 1] has rationality test 1
\\[6, 3, 1] has rationality test 0
\\[6, 3, 2] has rationality test 0
\\[6, 5, 1] has rationality test 0
\\[6, 5, 2] has rationality test 1
\\[10, 2, 3] has rationality test 1
\\[10, 4, 1] has rationality test 1
\\[10, 6, 1] has rationality test 1
\\[10, 8, 3] has rationality test 1
\\[4, 1, 1] has rationality test 1
\\[4, 2, 1] has rationality test 1
\\[4, 3, 1] has rationality test 1

\\DMList 
\\ = [[3, 1, 1], [6, 1, 1]]

for(i=1,length(DMList),print(DMList[i]," has rationality test ", rattest(DMList[i][1],DMList[i][2],DMList[i][3])))
\\[3, 1, 1] has rationality test 1
\\[6, 1, 1] has rationality test 0


subst(subst(subst(SQRES_paper,r2,+2+2*sqrt(3)),zet1,exp(2*Pi*I*6/12)),zet2,exp(2*Pi*I*5/12))