/* cat.gp by Jon McCammond and Murray Elder */
/* started on Oct 17, 2001                  */

\p 50; /*reset the precision */

establishGlobals()=
{
 /* generic global variables */
 global(true,false,infinity);
 global(gMax,verbose,testC,testA,testM,testV);
 global(p,pMods,pValues,nf,lrMatch,rotateMatch);
 global(specialRoot);
 /* global t variables */
 global(tMax,tList,tFaces,tBase,tLeft,tRight);
 global(tSqCoord,tSqShort,tSqDihedVec);
 global(tCoord,tShort,tDihedVec);
 global(tSqCosine,tSqSine,tCosine,tSine);
 global(spinTable,sMax,sMin);
 /* global g variables */
 global(testNum,done,aGalleries,mGalleries,vGalleries);
 global(alength,mlength,vlength,aRegions,mRegions,vAngles);
 global(aHalf,vHalf,g,gEdges,gCoords,glength);
 global(reg,regSides,regCorners,regCurrent);
 global(moreThanHalf,halfWayPt,halfG);
}
if(true!=1,establishGlobals;);

initGlobalVar()=
{
 true=1;
 false=0;
 infinity=10^(10^4);
 gMax=1000;
 verbose=true;
 testC=true;
 testA=true;
 testM=true;
 testV=true;
 p=listcreate(gMax);
 pMods=listcreate(gMax);
 pValues=listcreate(gMax);
 lrMatch=listcreate(gMax);
 rotateMatch=listcreate(gMax);
 tMax=0; 
}
initGlobalVar();


clearTVar()= 
{
 gMax=500;
 p=listcreate(20);
 pMods=listcreate(20);
 pValues=listcreate(20);
 tMax=0;
 tList=listcreate(gMax);
 tFaces=listcreate(gMax);
 tBase=listcreate(gMax);
 tLeft=listcreate(gMax);
 tRight=listcreate(gMax);
 tSqCoord=listcreate(gMax);
 tSqShort=listcreate(gMax);
 tSqDihedVec=listcreate(gMax);
 tSqCosine=listcreate(gMax);
 tSqSine=listcreate(gMax);
 tCoord=listcreate(gMax);
 tShort=listcreate(gMax);
 tDihedVec=listcreate(gMax); 
 tCosine=listcreate(gMax);
 tSine=listcreate(gMax);
 spinTable=listcreate(gMax);
 sMax=infinity;
 sMin=0;
 specialRoot=0;
}
clearTVar();

clearGVar()=
{
 gpMax=65533;
 gMax=500;
 testNum=0;
 done=false;
 aGalleries=listcreate(gpMax);
 mGalleries=listcreate(gpMax);
 vGalleries=listcreate(gpMax);
 aRegions=listcreate(gpMax);
 mRegions=listcreate(gpMax);
 vAngles=listcreate(gpMax);
 alength=listcreate(gpMax);
 mlength=listcreate(gpMax);
 vlength=listcreate(gpMax);
 aHalf=listcreate(gpMax);
 vHalf=listcreate(gpMax);
 g=vector(gMax);
 gEdges=List(vector(gMax));
 g[1]=-1*tMax;
 listput(gEdges,[1,2],1);    
 gCoords=List(vector(gMax+2));
 glength=1;
 reg=List(vector(gMax,i,listcreate(50)));
 regSides=List(vector(gMax,i,vector(50)));
 regCorners=List(vector(gMax,i,vector(50)));
 regCurrent=List(vector(gMax,x,false));
 moreThanHalf=false;
 halfWayPt=0;
 halfG=listcreate(gMax);
}
clearGVar();

/* catbits.gp*/

/* signNum replaced with sign(x) */
/* iFactors replaced with factor(x) */
/* sqrtPrime replaced with sqrtRat(x) */

sqrtRat(n)=
{
  local(a,b,c,d);
  a=numerator(n);
  b=denominator(n);
  c=sqrtint(a);
  d=sqrtint(b);
    if(c^2==a && d^2==b,
      return(c/d),
      return(sqrt(a/b))
    );
}

findSpecialRoot()=
{
  local(ans,i,j);
  specialRoot=0;
  ans=0;
  for(i=1,length(nf.roots),
    ans=true;
    for(j=1,length(p),
      ans=ans && sign(subst(lift(pMods[j]),x,nf.roots[i]))==1;
    );
    if(ans==1,specialRoot=i;break;);
  );
  return(specialRoot);
}

createNF(pr)=
{
  listput(pMods,x);
  listput(pValues,sqrt(pr));
  nf=nfinit(x^2-pr);
  return();  
}

extendNF(pr)=
{
 local(mat1,pol1,i,pol2);
 nf=changevar(nf,[y]);
 mat1=nffactor(nf,x^2-pr);
 pol1=rnfequation(nf,mat1[1,1],1);
 for(i=1,length(p),pMods[i]=subst(lift(pMods[i]),x,pol1[2]));
 listput(pMods,x-pol1[3]*pol1[2]);
 nf=nfinit(pol1[1]);
 pol2=polredabs(pol1[1],1);
 for(i=1,length(p)+1,pMods[i]=subst(lift(pMods[i]),x,pol2[2]));
 nf=nfinit(pol2[1]);
 return();
}

addPrimeNum(n)=
{
 if(abs(n)!=1 && !setsearch(Set(p),n),
   if(length(p)==0,
     createNF(n);
     print("");
     print("The number of distinct primes which appear under");
     print("square-roots when your tetrahedra are represented");
     print("in 3-space will greatly affect the speed of the program.\n");
     print("At this point the number of such primes has already reached:");
   ,
     extendNF(n);
   );
   findSpecialRoot();
   listput(p,n);
   print(" ",length(p));
 );
}

addPrimes(L)=
{ 
 for(i=1,length(L),
  addPrimeNum(L[i]);
 );
}

checkPrimesIn(n)=
{
  if(n!=0,
    a=factor(n);
    for(j=1,matsize(a)[1],if(Mod(a[j,2],2)==Mod(1,2),addPrimeNum(a[j,1])))
  );
}

evaluatePrimes(i)=
{
  return(subst(lift(pMods[i]),x,nf.roots[specialRoot]));
}

algToNum(u)=return(subst(lift(u),x,nf.roots[specialRoot]));

squNumToAlg(n)=
{
  local(c,d,ans,f,i,j,loc);
  num=numerator(n);
  denom=denominator(n);
  c=sqrtint(num);
  d=sqrtint(denom);
  if(c^2==num && d^2==denom,
    return(c/d);
  ,
    ans=1;
    f=factor(n);
    for(i=1,matsize(f)[1],
      ans=ans*f[i,1]^(divrem(f[i,2],2)[1]);
      if(Mod(f[i,2],2)==Mod(1,2),
        loc=0;
        for(j=1,length(p),
          if(p[j]==f[i,1],loc=j;break;);
        );
        ans=ans*pMods[loc];
      );
    );
  );
  return(ans);
} 


vecEqual(u,v)=
{
  if(length(u)!=length(v),return(false););
  for(i=1,length(u),
    if(algToNum(u[i])!=algToNum(v[i]),return(false););
  );
  return(true);
}

vecToNums(u)=return(vector(length(u),i,algToNum(u[i])));
vecAdd(u,v)=return(vector(length(u),i,u[i]+v[i]));
vecSub(u,v)=return(vector(length(u),i,u[i]-v[i]));
vecSMult(x,u)=return(vector(length(u),i,x*u[i]));

xyProj(u)=return([u[1],u[2],0]);
complexMult(u,v)=return([u[1]*v[1]-u[2]*v[2], u[1]*v[2]+u[2]*v[1]]);
crossProd(u,v)=return([u[2]*v[3]-u[3]*v[2],u[3]*v[1]-u[1]*v[3],u[1]*v[2]-u[2]*v[1]]);
dotProd(u,v)=return(u[1]*v[1]+u[2]*v[2]+u[3]*v[3]);

inEdge(t,u,v)=
{
  local(i,j);
  i=crossProd(u,v);
  j=crossProd(i,t);
  if(sign(algToNum(dotProd(i,t)))!=0,return(false));
  if(sign(algToNum(dotProd(t,u)))<=0,return(false));
  return(sign(algToNum(dotProd(j,u)*dotProd(j,v)))==-1);
}

sameSide(s,t,u,v)=return(sign(algToNum(dotProd(s,crossProd(u,v))*dotProd(t,crossProd(u,v))))==1);
inTriangle(s,t,u,v)=return(sameSide(s,t,u,v) && sameSide(s,u,v,t) && sameSide(s,v,t,u));

/* *********** */
/* cattetra.gp */
/* ************ */
in(u,v)=
{
 local(ans,i);
 ans=false;
 for(i=1,length(v),
  if(u==v[i],ans=true;break;);
 );
 return(ans);
}
 
properT(l)=
{
  if(length(l)!=6,error("Edge list must have 6 entries"););
/*for(i=1,6,
    if(type(l[i])!=t_Frac) then Error("Squared edge lengths must be rational");fi;
      if (l[i]<=0) then Error("Squared edge lengths must be positive");fi;
  od;*/
  tris=[[1,2,4],[1,3,5],[2,3,6],[4,5,6]];
  for(i=1,4,
    j=[l[tris[i][1]],l[tris[i][2]],l[tris[i][3]]];
    vecsort(j);
    if(sqrt(j[1])+sqrt(j[2])-sqrt(j[3])<=0,error("Edge list fails the triangle inequality"););
  );
}

tSqCoordinates(l)=
{
  local(a,b,c,d,e,f,i,j,k,co_a,p1,q,r,s,t);
  properT(l);
  a=l[1];b=l[2];c=l[3];d=l[4];e=l[5];f=l[6];
  i=(a+b-d)/2; j=(a+c-e)/2; k=(b+c-f)/2;
  co_a=(i*j-a*k)^2/((a*b-i^2)*(a*c-j^2));
  p1=i^2/a; q=b-i^2/a; r=j^2/a;
  s=(c-j^2/a)*co_a; t=(c-j^2/a)*(1-co_a);
  return([[a,0,0],[p1,q,0],[r,s,t]]);
}

tSqShortcut(l)=
{
  local(a,i,j,k,p1,q,r,s,t);
  a=l[1][1]; p1=l[2][1]; q=l[2][2];
  r=l[3][1]; s=l[3][2]; t=l[3][3];
  j=s/q; k=t/(a*q);
  i=(r/a)+((s*p1)/(q*a))-2*sqrtRat((r*s*p1)/(q*a^2));
  return([i,j,k]);
}

faceNum(l)=
{
  local(i);
  if(!(in(l,tFaces)),
    listput(tFaces,l);
  );
  for(i=1,length(tFaces),
    if(tFaces[i]==l,return(i););
  );
  error("Something went wrong");
}

addOrientedTetra(el)=
{
  if(!(in(el,tList)),
    c=tSqCoordinates(el);
    s=tSqShortcut(c);
    tMax=length(tList)+1;
    listput(tList,el);
    listput(tBase,faceNum([el[1],el[2],el[4]]));
    listput(tLeft,faceNum([el[1],el[3],el[5]]));
    listput(tRight,faceNum([el[3],el[2],el[6]]));
    listput(tSqCoord,c);
    listput(tSqShort,s);
    tMax=length(tList);
 );
}

baseFace(i)=return([tList[i][1],tList[i][2],tList[i][4]]);
leftFace(i)=return([tList[i][1],tList[i][3],tList[i][5]]);
rightFace(i)=return([tList[i][3],tList[i][2],tList[i][6]]);

switchLR()=
{
  local(i,j,a,b,c);
  lrMatch=vector(tMax);
  for(i=1,tMax,
    for(j=i,tMax,
      a=tList[i]; b=tList[j];
      c=[b[2],b[1],b[3],b[4],b[6],b[5]];
      if(a==c,
        lrMatch[i]=j; lrMatch[j]=i;
      );
    );
  );
  return(lrMatch);
}

calcFaceAngles()=
{
  local(a,b,c,i,ans);
  tSqCosine=vector(length(tFaces));
  tSqSine=vector(length(tFaces));
  for(i=1,length(tFaces),
    a=tFaces[i][1];
    b=tFaces[i][2];
    c=tFaces[i][3];
    tSqCosine[i]=(a+b-c)^2/(4*a*b);
    tSqSine[i]=1-tSqCosine[i];
  );
  return;
}

calcDihedralVectors()=
{
 local(i,a,b,c,d,e,f,j,k,l,m);
 tSqDihedVec=vector(tMax);
 for(i=1,tMax,
  a=tList[i][1];
  b=tList[i][2];
  c=tList[i][3];
  d=tList[i][4];
  e=tList[i][5];
  f=tList[i][6];
  j=(a+b-d)/2;
  k=(a+c-e)/2;
  l=(b+c-f)/2;
  m=(j*k-a*l)^2/((a*b-j^2)*(a*c-k^2));
  tSqDihedVec[i]=[m,1-m];
  );
}


checkPrimes()=
{
  local(i,j,k,l);
  for(i=1,length(tFaces),
    checkPrimesIn(tSqCosine[i]);
    checkPrimesIn(tSqSine[i]);
  );
  for(i=1,length(tSqCoord),
    for(j=1,3,
      for(k=1,j,
        checkPrimesIn(tSqCoord[i][j][k]);
      );
    );
  );
  return();
}

calcAlgCoord()=
{
  local(a,b,c,d,e,f,z,i);
  for(i=1,tMax,
    z=squNumToAlg(0);
    a=squNumToAlg(tSqCoord[i][1][1]);
    b=squNumToAlg(tSqCoord[i][2][1]);
    c=squNumToAlg(tSqCoord[i][2][2]);
    e=[a,b,c];
    d=squNumToAlg(tSqCoord[i][3][1]);
    e=squNumToAlg(tSqCoord[i][3][2]);
    f=squNumToAlg(tSqCoord[i][3][3]);
    listput(tCoord,[[a,z,z],[b,c,z],[d,e,f]],i);
    a=squNumToAlg(tSqShort[i][1]);
    b=squNumToAlg(tSqShort[i][2]);
    c=squNumToAlg(tSqShort[i][3]);
    listput(tShort,[a,b,c],i);
    a=squNumToAlg(tSqDihedVec[i][1]);
    b=squNumToAlg(tSqDihedVec[i][2]);
    listput(tDihedVec,[a,b],i);
  );
  for(i=1,length(tFaces),
    listput(tCosine,squNumToAlg(tSqCosine[i]),i);
    listput(tSine,squNumToAlg(tSqSine[i]),i);
  );
}

getTDihedralVec(n)=if(n>0,return(tDihedVec[n]);,return(tDihedVec[lrMatch[-1*n]]););


calcSpinBounds()=
{
  local(i,j,u,v);
  sMin=infinity;
  sMax=0;
  for(i=1,tMax,
    v=getTDihedralVec(i);
    u=v;
    j=1;
    while(algToNum(u[2])>0,
      j=j+1;
      u=complexMult(u,v);
    );
    sMin=min(sMin,j-1);
    sMax=max(sMax,j);
  );
  return([sMin,sMax]);
}


calcSpin(z)=
{
  local(i,pilength,new,spin);
  if(z[1]<0,
    for(i=1,length(z),z[i]=lrMatch[-1*z[i]];);
  );
  pilength=-1;
  i=glength;
  spin=[1,0];
  for(i=1,length(z),
    new=complexMult(spin,getTDihedralVec(z[i]));
    if(algToNum(new[2]*spin[2])<= 0,pilength=pilength+1;);
    spin=new;
  );
  return(pilength);
}

testSpin(z)=
{
  local(a);
  if(length(z) <= sMin,return(true););
  if(sMax <= length(z),return(false););
  a=calcSpin(z);
  return(a==0);
}


incPosGallery()=
{
  if(g[1]==0,
    done=true;
  ,
    if(g[glength]<tMax,
      g[glength]=g[glength]+1;
    ,
      g[glength]=0;
      glength=glength-1;
      incPosGallery();
    );
  );
}

calcSpinTable()=
{
  local(a,b,c,d,i,j,k,l,testNum1,bounds);
  testNum1=0;
  calcSpinBounds();
  glength=2;
  g=vector(sMax-1);
  g[1]=1;
  incPosGallery();
  done=false;
  spinTable=listcreate((2*tMax)^sMax);
  while(!done,
    if(testSpin(vector(glength,i,g[i])),
      if(sMin<glength,
        testNum1=testNum1+1;
        k=length(spinTable);
        listput(spinTable,vector(glength,i,g[i]));
        if(glength+1<sMax,glength=glength+1;);
      ,
        if(glength+1<sMax,glength=glength+1;);
      );
    );
    incPosGallery();
  );
  g=vector(sMax-1);
  return(spinTable);
}



addTetra(el)=
{
  local(i,l,m);
  properT(el);
  m=[[0,el[1],el[2],el[3]],\
    [el[1],0,el[4],el[5]],\
    [el[2],el[4],0,el[6]],\
    [el[3],el[5],el[6],0]];
  l=[[1,2,3,4],[1,2,4,3],[1,3,2,4],[1,3,4,2],[1,4,2,3],[1,4,3,2],\
    [2,1,3,4],[2,1,4,3],[2,3,1,4],[2,3,4,1],[2,4,1,3],[2,4,3,1],\
    [3,1,2,4],[3,1,4,2],[3,2,1,4],[3,2,4,1],[3,4,1,2],[3,4,2,1],\
    [4,1,2,3],[4,1,3,2],[4,2,1,3],[4,2,3,1],[4,3,1,2],[4,3,2,1]];
  for(i=1,24,\
    L=[m[l[i][1]][l[i][2]],m[l[i][1]][l[i][3]],m[l[i][1]][l[i][4]],\
       m[l[i][2]][l[i][3]],m[l[i][2]][l[i][4]],m[l[i][3]][l[i][4]]];
    addOrientedTetra(L);
  );
  switchLR;
  calcFaceAngles;
  calcDihedralVectors;
  checkPrimes;
  findSpecialRoot;
  calcAlgCoord;
  calcSpinBounds;
  calcSpinTable;
  return;
}

createTetra(l)=
{
 local(i);
 clearTVar();
 for(i=1,length(l),addTetra(l[i]););
 print("");
 print(" Number of Tetrahedra created: ",length(l));
 print(" Number of Oriented Tetrahedra: ",tMax);
 print(" Number of Faces up to isometry: ",length(tFaces));
 print(" Primes whose roots are needed = ",vecsort(vector(length(p),i,p[i])));
 print("");
 return;
}

regular()=createTetra([[1,1,1,1,1,1]]);
coxeter()=createTetra([[4,3,3,3,3,4]]);
thurston()=createTetra([[3,3,3,3,3,3],[4,3,3,3,3,3],[4,3,3,3,3,4]]);

/*  catgall.gp */

doTest(s)=
{
  if(s=="c",testC=true;);
  if(s=="a" && testC,testA=true;);
  if(s=="m" && testC,testM=true;);
  if(s=="v" && testC,testV=true;);
  print("Test for     coordinates = ",testC);
  print("Test for annular gallery = ",testA);
  print("Test for moebius gallery = ",testM);
  print("Test for  vertex gallery = ",testV);
  print("(note:when coordinates are not found, the");
  print(" other tests cannot be run)");
}

dontTest(s)=
{
  if(s=="c",testC=false;testA=false;testM=false;testV=false;);
  if(s=="a",testA=false;);
  if(s=="m",testM=false;);
  if(s=="v",testV=false;);
  print("Test for annular gallery = ",testA);
  print("Test for moebius gallery = ",testM);
  print("Test for  vertex gallery = ",testV);
}

whatToTest(n)=
{
  local(i,j,k);
  i= lift(Mod(n,2));
  j= lift(Mod((n-i)/2,2));
  k= lift(Mod((n-i-2*j)/4,2));
  testA=false;
  testM=false;
  testV=false;
  if(n>0, testC=true;);
  if(i==1, testA=true;);
  if(j==1, testM=true;);
  if(k==1, testV=true;);
  print("Test for annular gallery = ",testA,"");
  print("Test for moebius gallery = ",testM,"");
  print("Test for  vertex gallery = ",testV,"");
  return;
}

calcCoords(b)=
{
  testC=b;
  if(!b,print("Calculate coordinates = ",testC);whatToTest(0););
}

separates(u,v,r)=return(algToNum(dotProd(u,r)*dotProd(v,r))<0);
determinant(u,v)=return(algToNum(u[1]*v[2]-u[2]*v[1]));
flip(a)=return(-1*sign(a)*lrMatch[abs(a)]);
flipList(u)=local(i);return(vector(length(u),i,-1*sign(u[i])*lrMatch[abs(u[i])]));

lessThan(u,v)=
{
  local(i,l);
  l=min(length(u),length(v));
  i=1;
  while(u[i]==v[i],
    i=i+1;
    if(i==l,return(length(u)<length(v)););
  );
  return(u[i]<v[i]);
}


getFinalSyl()=
{
  local(i,j);
  j=glength;
  while(g[j]*g[glength]>0,
    j=j-1;
    if(j==0, break;);
  );
  return(vector(glength-j,i,g[i+j]));
}

spinOk()=
{
  local(a,i,fin);
  fin=getFinalSyl();
  if((length(fin) <= sMin),return(true););
  if((sMax <= length(fin)),return(false););
  if(fin[1]<0,fin=vector(length(fin),i,lrMatch[abs(fin[i])]));
  if(in(fin,spinTable),return(true);,return(false););
}



faceCheck(i,j)=
{
  if(i>0,return(tLeft[i]==tBase[abs(j)]););
  if(i<0,return(tRight[-1*i]==tBase[abs(j)]););
  error("something went wrong");
}

faceCheckRev(i,j)=
{
  if(i>0,return(tLeft[i]==tBase[lrMatch[abs(j)]]););
  if(i<0,return(tRight[-1*i]==tBase[lrMatch[abs(j)]]););
  error("something went wrong");
}

checkFaces(l)=
{
  local(i);
  if(length(l)==1,return(true););
  for(i=2,length(l),if(!faceCheck(l[i-1],l[i]),return(false);););
  return(true);
}
    
newPt(u,v,i)=
{
  local(a,v1,v2,v3,w,z);
  a=tShort[i];
  z=crossProd(u,v);
  v1=vecSMult(a[1],u);
  v2=vecSMult(a[2],v);
  v3=vecSMult(a[3],z);
  w=vecAdd(vecAdd(v1,v2),v3);
  return(w);
}

getNewGCoord()=
{
  local(i,j,k);
  if(glength==1,
    i=abs(g[1]);
    gCoords[1]=tCoord[i][1];
    gCoords[2]=tCoord[i][2];
    gCoords[3]=tCoord[i][3];
  ,
    i=gCoords[gEdges[glength][1]];
    j=gCoords[gEdges[glength][2]];
    k=abs(g[glength]);
    gCoords[glength+2]=newPt(i,j,k);
  );
}

getGEdges()=
{
  local(i);
  gEdges[1]=[1,2];
  for(i=1,glength,
    if(g[i]>0,gEdges[i+1]=[gEdges[i][1],i+2];);
    if(g[i]<0,gEdges[i+1]=[i+2,gEdges[i][2]];);
    if(g[i]==0,error("something went wrong"););
  );
  return(gEdges);
}

getGCoords()=
{
  local(i,l);
  l=glength;
  for(i=1,l,glength=i;getNewGCoord(););
  return(vector(glength,i,gCoords[i]));
}

minCycle(c)=
{
  local(i,j,l,xx,m);
  l=length(c);
  xx=vector(2*l);
  for(i=1,l,xx[i]=c[i];xx[i+l]=c[i];);
  m=0;
  for(i=1,l,
    j=1;
    while(xx[j+m]==xx[j+i],
      j=j+1;
      if(j>l,return(vector(l,i,xx[m+i])););
    );
    if(xx[j+m]>xx[j+i],m=i;);
  );
  return(vector(l,i,xx[m+i])); 
}

aCycleRep()=
{
  local(i,l,m1,m2);
  l=vector(glength,i,g[i]);
  m1=minCycle(l);
  if(l!=m1,return(false););
  m2=minCycle(flipList(l));
  return(!lessThan(m2,m1));
}

mCycleRep()=
{
  local(i,l,l1,l2,m);
  l=vector(2*glength);
  l2=flipList(vector(glength,i,g[i]));
  for(i=1,glength,
    l[i]=g[i];
    l[glength+i]=l2[i];);
  m=minCycle(l);
  return(l==m);
}

getRegCorners(k)=
{
  local(i);
  for(i=1,length(reg[k])-1,
    regCorners[k][i]=crossProd(regSides[k][i],regSides[k][i+1]);
  );
  regCorners[k][length(reg[k])]=regCorners[k][1];
}
    
getRegSides(k)=
{
  local(i);
  for(i=1,length(reg[k])-1,
    if(reg[k][i]>0,regSides[k][i]=gCoords[reg[k][i]];);
    if(reg[k][i]<0,regSides[k][i]=vecSMult(-1,gCoords[abs(reg[k][i])]););
    if(reg[k][i]==0,error("Something's wrong"););
  );
  regSides[k][length(reg[k])]=regSides[k][1];
}

getNewAReg(k)=
{
  local(i,newSide,newNum,t);
  if(regCurrent[k],return(true););
  if(k==1, 
    if(g[1]>0,reg[1]=[-1,2,3,-1];,reg[1]=[-1,2,-3,-1];);
  ,
    if(!getNewAReg(k-1),return(false););
    if(g[k]>0,newNum=k+2;newSide=gCoords[k+2];);
    if(g[k]<0,newNum=-k-2;newSide=vecSMult(-1,gCoords[k+2]););
    if(g[k]==0,error("Something's wrong"););
    reg[k]=listcreate(50);
    for(i=1,length(reg[k-1]),
      if(algToNum(dotProd(regCorners[k-1][i],newSide))>0,
        if(i==length(reg[k-1])+1,
          listput(reg[k],reg[k-1][1]);
        ,
          listput(reg[k],reg[k-1][i]);
        );
     ,
       if(length(reg[k])==0,
         listput(reg[k],newNum);
       ,
         if(reg[k][length(reg[k])]<>newNum,
           listput(reg[k],reg[k-1][i]);
           listput(reg[k],newNum);
         );
       );
     );
   );
   if(length(reg[k])<3,return(false););
 );
 getRegSides(k);
 getRegCorners(k);
 regCurrent[k]=true;           
 return(true);
}


getAReg(l)=
{
  local(i);
  glength=length(l);
  for(i=1,glength,
    g[i]=l[i];
    regCurrent[i]=false;
  );
  getGEdges();
  getGCoords();
  return(getNewAReg(glength));
}

/* cattest.gp */

testAGall()=
{
  local(a1,b1,c1,a2,b2,c2,i,j,k,l,ap,bp,as,ans); 
  if(glength==1,return;);
  if(!faceCheck(g[glength],g[1]),return;);
  if(!aCycleRep,return;);
  a1=gCoords[1];
  b1=gCoords[2];
  c1=gCoords[3];
  a2=gCoords[gEdges[glength+1][1]];
  b2=gCoords[gEdges[glength+1][2]];
  c2=newPt(a2,b2,abs(g[1]));
  i=vecSub(a1,a2);
  j=vecSub(b1,b2);
  k=vecSub(c1,c2);
  if(vecEqual(i,[0,0,0]) || vecEqual(j,[0,0,0]) || vecEqual(k,[0,0,0]),return;);
  ap=crossProd(i,j);
  bp=crossProd(i,k);
  if(!vecEqual(crossProd(ap,bp),[0,0,0]),return;);
  if(vecEqual(ap,[0,0,0]),ap=bp;);
  ans=true;  
  for(i=1,glength+1,
    j=gCoords[gEdges[i][1]];
    k=gCoords[gEdges[i][2]];
    ans=ans && separates(j,k,ap);
  );
  if(ans, 
    if(getNewAReg(glength),
      as=crossProd([0,0,1],ap);
      if(!inEdge(as,a1,b1),as=vecSMult(-1,as););
      for(i=1,glength,
        j=gCoords[gEdges[i][1]];
        k=gCoords[gEdges[i][2]];
        if(g[i]<0,
          l=gCoords[gEdges[i+1][1]];
          if(inEdge(as,l,k),return;);
        ,
          l=gCoords[gEdges[i+1][2]];
          if(inEdge(as,j,l),return;);
        );
        if(inTriangle(as,j,k,l),return;);
      );
      listput(aGalleries,vector(glength,i,g[i]));
      listput(aRegions,reg[glength]);
      return(true);
    ,
      error("Something's wrong");
    );
  );
 return;
}



testMMaybe()=
{
  local(i,l);
  if(mCycleRep,        
    l=flipList(vector(glength,i,g[i]));
    for(i=1,glength,
      g[glength+i]=l[i];
    );
    glength=2*glength;
    getGEdges();
    getGCoords();
    if(getNewAReg(glength),
      listput(mGalleries,vector(glength/2,i,g[i]));
      listput(mRegions,reg[glength]);
    );
    glength=glength/2;
    for(i=1,glength,
      g[glength+i]=0;
      regCurrent[glength+i]=false;
    );
  );
}

testMGall()=
{
  local(a1,b1,c1,a2,b2,c2,i,j,k,l,ap,bp,as,ans); 
  if(!faceCheck(g[glength],flip(g[1])),return("a"););
  if(!mCycleRep(),return("b"););
  a1=gCoords[1];
  b1=gCoords[2];
  c1=gCoords[3];
  b2=gCoords[gEdges[glength+1][1]];
  a2=gCoords[gEdges[glength+1][2]];
  c2=newPt(b2,a2,lrMatch[abs(g[1])]);
  /* add since it's moebius */
  i=vecAdd(a1,a2); 
  j=vecAdd(b1,b2);
  k=vecAdd(c1,c2);
  /* temp 
  print("i=",vecToNums(i));
  print("j=",vecToNums(j));
  print("k=",vecToNums(k));
  print(vecToNums(crossProd(i,j))); */
  /* 3-dim case */
  if(dotProd(crossProd(i,j),k)!=0,return("c"););
  /* 0-dim case */
  if(vecEqual(i,[0,0,0]) && vecEqual(j,[0,0,0]) && vecEqual(k,[0,0,0]),
    testMMaybe();return("d");
  );
  ap=crossProd(i,j);
  bp=crossProd(i,k);
  if(!vecEqual(crossProd(ap,bp),[0,0,0]),return("e"););
  if(vecEqual(ap,[0,0,0]),ap=bp;);
  ans=true;  
  for(i=1,glength+1,
    j=gCoords[gEdges[i][1]];
    k=gCoords[gEdges[i][2]];
    ans=ans && separates(j,k,ap);
  );
  if(ans, 
    /* calculate length */
    as=crossProd([0,0,1],ap);
    if(!inEdge(as,a1,b1),as=vecSMult(-1,as););
    for(i=1,glength,
      j=gCoords[gEdges[i][1]];
      k=gCoords[gEdges[i][2]];
      if(g[i]<0,
        l=gCoords[gEdges[i+1][1]];
        if(inEdge(as,l,k), return("f"););
      ,
        l=gCoords[gEdges[i+1][2]];
        if(inEdge(as,j,l), return("g"););
      );
      if(inTriangle(as,j,k,l), return("h"););
    );
    /* add mRegion! */
    listput(mGalleries,vector(glength,i,g[i]));
  );
  return("i");
}
/* the return values aren't used except as diagnostics*/

testVMaybe()=
{
  local(i,u,minAngle,maxAngle,as,j,k,l);
  minAngle=[gCoords[2][2],gCoords[2][3]];
  maxAngle=[gCoords[3][2],gCoords[3][3]];
  for(i=2,glength,
    u=[gCoords[i+2][2],gCoords[i+2][3]];
    if(g[i]<0,if(determinant(u,maxAngle)>0,maxAngle=u;););
    if(g[i]>0,if(determinant(minAngle,u)>0, minAngle=u;););
    if(g[i]==0,error("something's wrong"););
    if(determinant(minAngle,maxAngle)<=0, return;);
  );
  listput(vGalleries,vector(glength,i,g[i]));
  listput(vAngles,[minAngle,maxAngle]);
  return;
}

testVTypical(pole)=
{
  local(i,j,k,l,as);
  for(i=2,glength,
    j=gCoords[gEdges[i][1]];
    k=gCoords[gEdges[i][2]];
    if(!separates(j,k,pole),return("a"););
  );
  /* calculate length */
  if(algToNum(gCoords[glength][2])==0 &&\
     algToNum(gCoords[glength][3])==0,return("b"););
  as=[1,0,0];
  for(i=1,glength,
    j=gCoords[gEdges[i][1]];
    k=gCoords[gEdges[i][2]];
    l=gCoords[i+2];
    if(g[i]<0,
      if(inEdge(as,l,k),return("c"););
    ,
      if(inEdge(as,j,l),return("d"););
    );
    if(inTriangle(as,j,k,l),return("e"););
  );
  listput(vGalleries,vector(glength,i,g[i]));
  listput(vAngles,[[gCoords[glength+2][2],gCoords[glength+2][3]]]);
  return("f");
}

testVGall()=
{
  local(i,j,k,l,vp,ans); 
  if(glength==1,return(false););
  if(g[2]<0,return(false);); 
  if(g[glength]<0, return(false););
  vp=crossProd(gCoords[1],gCoords[glength+2]);
  if(vecEqual(vp,[0,0,0]) && algToNum(gCoords[glength+2][1])<0,
    testVMaybe();
  ,
    testVTypical(vp);
  );
}

/* Iterative tests */
incGallery()=
{
  if(g[1]==0,done=true;return;);
  if(moreThanHalf && glength==halfWayPt,moreThanHalf=false;);
  if(g[glength]==tMax,
    g[glength]=-1*tMax;
    regCurrent[glength]=false;
  ,
    g[glength]=g[glength]+1;
    regCurrent[glength]=false;
  );
 if(g[glength]==0,
   glength=glength-1;
   incGallery();
  );
  return;
}

feasible()=
{
  if(glength==1,return(true););
  if(faceCheck(g[glength-1],g[glength]),return(spinOk()););
  return(false);        
}

sameSigns(u,v)=
{
  local(ans,i);
  ans=true;
  for(i=1,length(u),
    if(signNum(algMult(u[i],v[i]))<>1,ans=false;);
  );
  return(ans);
}

printResults(t)=
{
  print("Results ");
  print("  tested  = ",testNum," galleries");
  print("   time   = ",t," milliseconds");
  print("");

  if(testA, print("  annular = ",length(aGalleries)););
  if(testM, print("  moebius = ",length(mGalleries)););
  if(testV, print("  vertex  = ",length(vGalleries)););
  print("");
    
  print("To see the actual lists type:");
  if(testA, print("  aGalleries;"););
  if(testM, print("  mGalleries;"););
  if(testV, print("  vGalleries;"););
  print("");
  return;
}


/* Tests */
testRange(minL,maxL)=
{
  local(s);
  s=gettime();
  clearGVar();
  while(!done,
    if(g[glength]>0,
      gEdges[glength+1]=[gEdges[glength][1],glength+2];
    ,
      gEdges[glength+1]=[glength+2,gEdges[glength][2]];
    );
   if(feasible, 
     if(testC,getNewGCoord(););
     if(glength>=minL,
       testNum=testNum+1;
       if(testA,testAGall(););
       if(testM,testMGall(););
       if(testV,testVGall(););
     );
     if((glength+1<=maxL),glength=glength+1;);
   );
   incGallery();
   if(verbose && glength==2, print("[",g[1],",",g[2],"]"););
 );
 printResults(gettime()-s);
}


testSignGalleries(minL,maxL,signList)=
{
  local(i,s);
  s=gettime();
  clearGVar();
  while(!done,
    if(g[glength]>0,
      gEdges[glength+1]=[gEdges[glength][1],glength+2];
    ,
      gEdges[glength+1]=[glength+2,gEdges[glength][2]];
    );
    if(feasible, 
      if(sameSigns(vector(glength,i,g[i]),signList), 
        if(testC,getNewGCoord(););
        if(glength>=minL,
          testNum=testNum+1;
          if(testA,testAGall(););
          if(testM,testMGall(););
          if(testV,testVGall(););
        );
        if((glength+1<=maxL),glength=glength+1;);
      );
    );
    incGallery();
  );
  printResults(gettime()-s);
  return;
}

crossingPoint(u,v)=
{
  local cp;
  /* this routine finds the point at which the arc connecting
     u and v crosses the equator.
     We only call this when the answer makes sense. */
  cp=crossProd([0,0,1],crossProd(u,v));
  /* check to make sure cp is between u and v. */
  if(vecEqual(crossProd(cp,u),[0,0,0]),return(u););
  if(vecEqual(crossProd(cp,v),[0,0,0]),return(v););
  if(!inEdge(cp,u,v),
    cp=vecSMult(-1,cp);
    if(!inEdge(cp,u,v),Error("look at crossingPoint"););
  );
  return(cp);
}

halfWay()=
{
  local(b,t,bot,top,v,w);
  top=gCoords[gEdges[glength+1][1]];
  bot=gCoords[gEdges[glength+1][2]];
  t=algToNum(top[3]); b=algToNum(bot[3]);
  if(0<=b && 0<=t,return(false););
  if(b<0 && t<0,return(true););
  v=crossingPoint(top,bot);
  /* treat mixed cases (t<0<=b and b<0<=t) */
  if(b<0, /* spinning counter clockwise */
    w=crossProd(v,gCoords[2]);
    if(algToNum(w[3])==0,return(algToNum(v[2]*gCoords[2][2])<0););
    return(0<algToNum(w[3]));
  , /* spinning clockwise */
    w=crossProd(v,gCoords[1]);
    if(algToNum(w[3])==0,return(algToNum(v[1])<0););
    return(0>algToNum(w[3]));
  );
}

tooLong()=
{
  local(b,t,bot,top,v,w);
  top=gCoords[gEdges[glength+1][1]];
  bot=gCoords[gEdges[glength+1][2]];
  t=algToNum(top[3]); b=algToNum(bot[3]);
  if(0<=b && 0<=t,return(true););
  if(b<0 && t<0,return(false););
  v=crossingPoint(top,bot);
  /* treat mixed cases (t<0<=b && b<0<=t) */
  if(b<0, /* spinning counter clockwise */
    w=crossProd(v,gCoords[2]);
    if(algToNum(w[3])==0,return(algToNum(v[2]*gCoords[2][2])>0););
    return(0>algToNum(w[3]));
  , /* spinning clockwise */
    w=crossProd(v,gCoords[1]);
    if(algToNum(w[3])==0,return(algToNum(v[1])>0););
    return(0<algToNum(w[3]));
  );
}

test()=
{
  local(a,b,c,d,e,f,i,l,n,newTime,oldTime,s);
  s=gettime();
  totalTime=0;
  clearGVar();
  testC=true;
  while(!done,
    if(g[glength]>0,
      listput(gEdges,[gEdges[glength][1],glength+2],glength+1);
    ,
      listput(gEdges,[glength+2,gEdges[glength][2]],glength+1);
    );
    if(feasible(), 
      getNewGCoord();
      if(getNewAReg(glength), 
        testNum=testNum+1;
        if(testA, testAGall();); 
        if(testM, testMGall();); 
        if(testV, testVGall();); 
        /* update this */
        if(!moreThanHalf,
          if(halfWay, 
            moreThanHalf=true;
            halfWayPt=glength;
            if(verbose,
              newTime=gettime();
              print(" split = ",newTime," half = ",vector(glength,i,g[i]));
              totalTime=totalTime+newTime;
            );
          );
          glength=glength+1;
        ,
          if(!tooLong(), glength=glength+1;);
        );
      );
    );
    incGallery();
  );
  printResults(totalTime+gettime());
  return;
}

testInit(m)=
{
  local(a,b,c,d,e,f,i,l,n,newTime,oldTime,s);
  s=gettime();
  totalTime=0;
  clearGVar();
  testC=true;
  g[1]=-1*m;
  glength=1;
  while(g[1]==-1*m,
    if(g[glength]>0,
      listput(gEdges,[gEdges[glength][1],glength+2],glength+1);
    ,
      listput(gEdges,[glength+2,gEdges[glength][2]],glength+1);
    );
    if(feasible(), 
      getNewGCoord();
      if(getNewAReg(glength), 
        testNum=testNum+1;
        if(testA, testAGall();); 
        if(testM, testMGall();); 
        if(testV, testVGall();); 
        /* update this */
        if(!moreThanHalf,
          if(halfWay, 
            moreThanHalf=true;
            halfWayPt=glength;
            if(verbose,
              newTime=gettime();
              print(" split = ",newTime," half = ",vector(glength,i,g[i]));
              totalTime=totalTime+newTime;
            );
          );
          glength=glength+1;
        ,
          if(!tooLong(), glength=glength+1;);
        );
      );
    );
    incGallery();
  );
  printResults(totalTime+gettime());
  return;
}

/*
testHalfGalleries() 
    local a,b,c,d,e,f,i,l,n,s;
    s=gettime();
    clearGVar();
    testC=true;
    while !done do
        if(g[glength]>0, gEdges[glength+1]=[gEdges[glength][1],glength+2];);
        if(g[glength]<0, gEdges[glength+1]=[glength+2,gEdges[glength][2]];);
        if(feasible(), 
            getNewGCoord();
            if(getNewAReg(glength), 
                testNum=testNum+1;
                if(testA, testAGall(););
                if(testM, testMGall(););
                if(testV, testVGall(););
                # update this
                if(halfWay(), 
                    halfG[length(halfG)+1]=List([1..glength],i->g[i]);    
                else
                    glength=glength+1;
                );
            );
        );
        incGallery();
    od;
    printResults(gettime()-s);
    return;
}

singleTest(l)
    local ansA,ansM,ansV;
    g=l;
    glength=length(l);
    if(!checkFaces(g),
        print("Faces don't match");
    else
        print(g,"");#
        getGEdges();
        getGCoords();
        if(testA, testAGall(););
        if(testM, testMGall(););
        if(testV, testVGall(););
    );
}

## Summarize signs
showVAngles()
    local i,j;
    for i in [1..length(vAngles)] do
        if(length(vAngles[i])=1,
            print(vecToText(vAngles[i][1]),"");
        else
            print("[",vecToText(vAngles[i][1]),",",vecToText(vAngles[i][2]),"]");
        );
    od;
}

gallerySigns(l)=
{
  local a,i,j,s,t;
  s=[];
  t=[];
  for i in l do
        a=List([1..length(i)],j->signNum(i[j]));
        if(!a in s,
            s[length(s)+1]=a;
            t[length(t)+1]=1;
        else
            t[Position(s,a)]=t[Position(s,a)]+1;
        );
    od;
    return(List([1..length(s)],i->[s[i],t[i]]));
}

gallSigns(l)
    local a,i,j;
    for a in l do
        for i in [1..length(a)] do
            a[i]=signNum(a[i]);
        od;
    od;
    return(Collected(l));
}

gallerylengths(l)
    local a,b,i,j,s,t;
    a=List([1..length(l)],i->length(l[i]));
    b=Collected(a);
    return(b);
}

getFullGalleries()
    local a,b,i,k,l,m,aL,bL,match,offset,new,num,s;
    s=gettime();
    num=0;#
    for a in halfG do
        aL=length(a);
        print(a,"");#
        for b in halfG do
            bL=length(b);
            if(aL<=bL, 
                m=aL;
            else
                m=bL;
            );
            for k in [1..m-1] do
                match=true;
                offset=aL-k;
                for l in [1..k] do
                    if(a[l+offset]<>b[l], match=false;);
                od;
                if(match,
                    num=num+1;#
                    if(num=1000, return;);#
                    new=List([1..aL],i->a[i]);
                    for i in [k+1..bL] do
                        new[offset+i]=b[i];
                    od;
                    singleTest(new);
                );
            od;
        od;
    od;
    printResults(gettime()-s);
}
*/


/* misc */
/*
timeMult(n)
    local a,b,i,s;
    s=gettime();
    pNumber=4;
    p=[2,3,5,7];
    a=[1,2,3,4,5,6,7,8,1,2,3,4,5,6,7,8];
    b=[9,8,7,6,5,4,3,2,1,2,3,4,5,6,7,8];
    for i in [1..n] do
        algMult(a,b);
    od;
    print("time = ",gettime()-s," milliseconds");
}

testTShort(i)
    local c;
    c=tCoord[i];
    return(c[3]=newPt(c[1],c[2],i));
}

testTCoord(L)
  local a,b,c,d;
  a=[0,0,0];b=L[1];c=L[2];d=L[3];
#  return([SquEucDist(a,b),SquEucDist(a,c),SquEucDist(a,d),
#          SquEucDist(b,c),SquEucDist(b,d),SquEucDist(c,d)]);
}

testLR()
  local a,b,c,d,i,j,n,L,MD;
  n=length(tList);
  L=[];
  MD=[];
  for i in [1..n] do
    for j in [1..n] do
      if(baseFace(j)=leftFace(i),
        L=[op(L),[i,j]];
        a=tCoord[i][1];
        b=tCoord[i][2];
        c=tCoord[i][3];
        d=newPt(a,c,tShort[j]);
#        MD=[op(MD),MinDist(a,b,c,d)];
      );
    od;
  od;
#  return(min(op(MD)));
}

testThurstonDi() 
    local i,j,k,l;
    for i in [1..length(tDihedVec)] do
        j=tCoord[i][1][1];
        l=algToNum(algMult(j,j));
        k=algToNum(tDihedVec[i][1]);
        if((l=3), 
            print([l,(k<=1/2)],"");
        elif((l=4), 
            print([l,(k<=((sqrtPrime(5)-1)/4))],"");  
        else
            error("something's wrong");
        );
    od;
}

#################
### test data ###
#################
*/


consider(L)=
{
  g=L;
  glength=length(g);
  getGEdges();
}

deeplyConsider(L)=
{
  consider(L);
  getGCoords();
}

/*
sortBySize(L)
    local i,j,l,s;
    s=gettime();
    l=List([1..gMax],i->[]);
    for i in [1..length(L)] do 
        Add(l[length(L[i])],L[i]);
    od; 
    print(List([1..gMax],i->length(l[i])),"");
    return(l);
}

sortBylength(L)
    local i,j,k,l,s;
    s=gettime();
    l=List([1..7],i->[]);
    for i in [1..length(L)] do
        deeplyConsider(L[i]);
        j=gCoords[glength+2];
        if(signNum(j[1])=1 && signNum(j[2])=1,
            Add(l[1],g);
        elif(signNum(j[1])=0 && signNum(j[2])=1,
            Add(l[2],g);
        elif(signNum(j[1])=-1 && signNum(j[2])=1,
            Add(l[3],g);
        elif(signNum(j[1])=-1 && signNum(j[2])=0,
            Add(l[4],g);
        elif(signNum(j[1])=-1 && signNum(j[2])=-1,
            Add(l[5],g);
        elif(signNum(j[1])=0 && signNum(j[2])=-1,
            Add(l[6],g);
        elif(signNum(j[1])=1 && signNum(j[2])=-1,
            Add(l[7],g);
        else
            error("something's wrong");
        );
    od;
    print("time = ",gettime()-s," milliseconds");
    print(List([1..7],i->length(l[i])));
    return(l);
}

recordlength(L)
    local i,j,k,l,s;
    s=gettime();
    l=[];
    for i in [1..length(L)] do
        deeplyConsider(L[i]);
        j=gCoords[glength+2];
        Add(l,j);
    od;
    print("time = ",gettime()-s," milliseconds");
    return(l);
}

endQuad(L)
    local i,j,k,l,s;
    s=gettime();
    for i in [1..length(L)] do
        consider(L[i]);
        j=gCoords[glength+2];
        print([signNum(j[1]),signNum(j[2])],"");
    od;
    print("time = ",gettime()-s," milliseconds");
    return;
}

sortByElength(L)
    local i,j,k,l,s;
    s=gettime();
    l=[[[],[]],[[],[]]];
    for i in [1..length(L)] do
        consider(L[i]);
        j=tList[abs(g[1])][1]-2;
        k=tList[g[glength]][3]-2;
        Add(l[j][k],g);
    od;
    print([length(l[1][1]),length(l[1][2]),length(l[2][1]),length(l[2][2])],"");
    return(l);
}

*/

/*
findSqlength(L)
    local i,j,k,l,s;
    s=gettime();
    l=List([1..4],i->[]);
    for i in [1..length(L)] do
        deeplyConsider(L[i]);
        j=dotProd(gCoords[1],gCoords[glength+2]);
        k=algToNum(algSMult(1/(tList[abs(g[1])][1]*tList[abs(g[glength])][3]),algMult(j,j)));
        if(k<1/9, 
            Add(l[1],g);
        elif(k<1/4,
            Add(l[2],g);
        elif(k<1/3,
            Add(l[3],g);
        else
            Add(l[4],g);
        );
    od;
    print("time = ",gettime()-s," milliseconds");
    print(List([1..4],i->length(l[i])),"");
    return(l);
}

selectSigns(S,L)=
{
  local(i,j,k,a);
  for(i=1,length(S),
    for j in [1..length(L)] do
            if(length(S[i])=length(L[j]),
                a=true;
                for k in [1..length(S[i])] do
                    if(S[i][k]*L[j][k]<=0, a=false;);
                od;
                if(a=true, print(i," ",S[i]," ",j," ",L[j],""););
            );
        od;
    od;
}

                    
*/



/*        */
/*  main  */
/*        */


clv()=
{
 true=false;
 read("cat.gp");
}

cl()=
{
  system("clear");
  print("----------------------------------------------------");
  print("cat.gp version 0.1 by Jon McCammond and Murray Elder");
  print("----------------------------------------------------");
  print("");
}
cl();