makeUnitPair(i) = [unitlist[i],th(unitlist[i])]

list_of_middle=[];

one_is_coprime(a,b)=gcd(a,q)==1||gcd(b,q)==1

/*for(i=0,q-1,for(j=0,q-1,for(k=0,length(ceez)-1,if(one_is_coprime(i,j),list_of_middle=concat(list_of_middle,[[myc_at_2(k),[i;j]]])))))
*/

myXX=sqrt(XX+O(p^10)) % q;

/*act_on_left(i, j) = 
mUP=makeUnitPair(i);[myc_at_2(testquat(quat_mult(mUP[1],list_of_middle[j][1]))),subst(mUP[2]*list_of_middle[j][2],x0,myXX)%q]

act_on_leftB(muP,j) 
=[myc_at_2(testquat(quat_mult(mUP[1],list_of_middle[j][1]))),subst(mUP[2]*list_of_middle[j][2],x0,myXX)%q]

getOrbit(jj)=aOLp=[];for(i=1,24,aOLp=concat(aOLp,[act_on_left(i,jj)]));aOLp

findCoset(ii)=aOLi=getOrbit(ii);coseti=[];for(i=1,length(list_of_middle),if(Set(aOLi)==Set(getOrbit(i)),coseti=concat(coseti,[i])));coseti
*/
lom=[];

for(i=0,q-1,for(j=0,q-1,if(gcd(i,q)==1||gcd(j,q)==1,lom=concat(lom,[[i;j]]))))

unitlist_H=[];

/*forstep(a=-2, 2, 1, forstep(b=-2, 2, 1, \
forstep(cc=-2, 2, 1, forstep(d=-2, 2, 1, \
if(a^2+b^2+cc^2+d^2==1, \
unitlist_H = concat(unitlist_H, [[a,b,cc,d]]))))))*/

/* experimental switching code */

if(flag,forstep(a=-2, 2, 1/2, forstep(b=-2, 2, 1/2, \
forstep(cc=-2, 2, 1/2, forstep(d=-2, 2, 1/2, \
if(a^2+b^2+cc^2+d^2==1, \
unitlist_H = concat(unitlist_H, [[a,b,cc,d]])))))),forstep(a=-2, 2, 1, \
forstep(b=-2, 2, 1/2, \
forstep(cc=-2, 2, 1/2, forstep(d=-2, 2, 1/2, \
if(a^2+b^2+cc^2+d^2==1, \
unitlist_H = concat(unitlist_H, [[a,b,cc,d]])))))))

unitlist_H=[[1,0,0,0],[-1,0,0,0]]

matunitlist_H = vector(length(unitlist_H),i,th(unitlist_H[i]));

getOrbit2(jj)=aOLp=[];for(i=1,length(matunitlist_H),aOLp=concat(aOLp,[matunitlist_H[i]*lom[jj]]));subst(aOLp,x0,myXX)%q

findSquares()=temp=[];for(i=1,q-1,if(issquare(i+O(p^2)),temp=concat(temp,[i])));temp

getOrbit_new(jj)=gO2=getOrbit2(jj);fS=findSquares();temp=[];for(i=1,length(gO2),for(j=1,length(fS),temp=concat(temp,[gO2[i]*fS[j] % q])));temp

findCoset2(ii)=aOLi=getOrbit2(ii);coseti=[];\
for(i=1,length(lom),if(Set(aOLi)==Set(getOrbit2(i)),coseti=concat(coseti,[i])));coseti

findCoset_new(ii)=aOLi=getOrbit_new(ii);coseti=[];\
for(i=1,length(lom),if(Set(aOLi)==Set(getOrbit_new(i)),coseti=concat(coseti,[i])));coseti

findRepresentatives()=S=[];T=Set(vector(length(lom),ii,ii));\
while(length(T)!=0,t1=eval(T[1]);S=concat(S,[t1]);T=setminus(T,Set(findCoset2(t1))));\
S

findRepresentatives_new()=S=[];T=Set(vector(length(lom),ii,ii));\
while(length(T)!=0,t1=eval(T[1]);S=concat(S,[t1]);T=setminus(T,Set(findCoset_new(t1))));\
S

getInverse(a)=tempnum=0;for(i=0,q-1,if((a*i)%q == 1,tempnum=i));tempnum

makeLift(ii)=ll=lom[ii];if(gcd(ll[1,1],q)==1,cc=0;dd=getInverse(ll[1,1]),cc=(-getInverse(ll[2,1]))%q;dd=0);[ll[1,1],cc;ll[2,1],dd]

findCosetMatrices()=S2=[];fR=findRepresentatives();for(i=1,length(fR),S2=concat(S2,[makeLift(fR[i])]));S2

findCosetMatrices_new()=S2=[];fR=findRepresentatives_new();for(i=1,length(fR),S2=concat(S2,[makeLift(fR[i])]));S2