//Minimal classes
//Oliver Braun

//Algorithm to compute minimal classes and maximal finite subgroups 

clear;

//n less or equal 3, K imaginary quadratic number field, fractional ideal arbitrary, to be chosen below
//via the variable steinitz and the numbering of the ideal classes provided by magma

n:=2;

d:=-21;
steinitz:=1;

print "************************************************************";
print "Voronoi algorithm and minimal classes";
print "For imaginary quadratic field " cat IntegerToString(d);
print "Lattice " cat IntegerToString(steinitz);
print "************************************************************";


load initialize;
load functions;

//Find a first perfect form

Pini:=MatrixRing(K,n)!1;
Pini[n][n]:=1/Norm(p2);
Pinie:=spurform(Pini);
Pinie2:=spurform(w*Pini);

Lini:=LatticeWithGram(Pinie);
Sini:=minvecs(Pini);
Rini:=prank(Pini);

count:=1;

while Rini lt n^2 and count lt 100 do
 count:=count+1;
 dir:=findperp(Sini);

 tsup:=1000;
 tinf:=0;
 t:=(tsup+tinf)/2;

 bool:=false;
 count2:=1;

 while not bool and count2 lt 100 do

  count2:=count2+1;
  M:=1;
  Pt:=Pini+t*dir;
  while M eq 1 do
   if IsPositiveDefinite(spurform(Pt)) then
    Lt:=LatticeWithGram(spurform(Pt));
    M:=hermitianmin(Pt);
    if M eq 1 then
     tinf:=t;
     t:=(tinf+tsup)/2;
     Pt:=Pini+t*dir;
    end if;
   else
    tsup:=t;
    t:=(tinf+tsup)/2;
    Pt:=Pini+t*dir;
   end if;
  end while;

  St:=minvecs(Pt);

  tt:=Rationals()!
      Min([(idealnorm(v)-K!((v*Pini*HermitianTranspose(v))[1,1]))/(K!((v*dir*HermitianTranspose(v))[1,1])) : v in St]);
  bool:=false;
  if tt lt t and tt gt 0 then
   Pc:=Pini+tt*dir;
   Pce:=spurform(Pc);
   Lc:=LatticeWithGram(Pce);
   M:=hermitianmin(Pc);
   if M eq 1 then
    bool:=true;
   else
    tsup:=tt;
    t:=(tinf+tsup)/2;
    Pt:=Pini+t*dir;
   end if;
  else
   tsup:=t;
   t:=(tsup+tinf)/2;
   Pt:=Pini+t*dir;
  end if;
 end while;

 Pini:=Pc;
 Pinie:=spurform(Pini);

 Lini:=LatticeWithGram(Pinie);
 Sini:=minvecs(Pini);
 Rini:=prank(Pini);
end while;

if Rini ne n^2 then
 error "In FirstPerfect: the form Rini is not perfect.";
end if;

//Enumerate perfect neighbours in order to obtain a set of representatives
//of perfect Hermitian forms

perfectlist:=[Pini];         //List of representatives of perfect forms
vectlist:=[**];              //List of shortest vectors of perfect forms
facelist:=[**];              //List of facets of V-domains of p. forms; given by shortest vectors
facevectList:=[**];          //Perpendicular form to shortest vectors defining the respective facet
Dim2facevectList:=[**];      //
FaceFormList:=[**];          //List of forms defined by those shortest vectors, which define the respective facet
AutList:=[**];               //List of Aut-Groups of the inverse FaceForms
Dim2FormList:=[**];          //
Dim2FaceList:=[**];          //
Dim2AutList:=[**];           //

numberoffaces:=[];           //List of number of faces of V-domains of p. forms
E:={**};                     //multiset encoding the Voronoi graph of perfect forms
Todo:=[Pini];                //List of perfect forms to be treated with Voronoi

PerfectNeighbourList:=[**];  //List of perfect neighbours of all (mod GL) perfect forms

CriticalValueList:=[**];     //List of critical rho values (from Voronoi's algorithm)
FacetVectorList:=[**];       //List of facet vectors (from Voronoi's algorithm)

while(#Todo gt 0) do
 P:=Todo[1];
 Pe:=spurform(P);
 L:=LatticeWithGram(Pe);
 m:=hermitianmin(P);
 Sk:=minvecs(P);
 Sproj:=[projzeileNorm(v) : v in Sk];
 Append(~vectlist,Sk);

 Exclude(~Todo,Todo[1]);

 if perfrank(Sk) ne n^2 then
  error "In enumerating perfect neighbours: p-rank of potential neighbour is too small."; 
 end if;

 G:=aut(P);
 G:=ChangeRing(G,Rationals());

 DonQhull:=Open("DonneePourQhull","w");
 Puts(DonQhull, IntegerToString(n^2) cat " RBOX c");
 Puts(DonQhull, IntegerToString(#Sproj+1) );
 Puts(DonQhull, &cat ["0 " : i in [1..n^2] ] );
 for st in [ &cat [RealToString(Injec(n)) cat " " : n in Eltseq(X)] : X in Sproj] do
  Puts(DonQhull,st);
 end for;
 delete DonQhull;

 //INSERT DIRECTORY OF QHULL HERE
 System("/home3/castor/tmp/kirschme/qhull/bin/qhull -Fv <DonneePourQhull >SommetsParFace");

 Faces:=[];
 SomFac:=Open("SommetsParFace","r");
 nbface := StringToInteger(Gets(SomFac));
 for i in [1..nbface] do
  Faces:=Append(Faces, Remove([StringToInteger(n) : n in Split(Gets(SomFac)," ")],1));
 end for;
 delete SomFac;
 Faces:=[ Exclude(F,0) : F in Faces | 0 in F];
 Faces:=[ {n : n in F} : F in Faces];
 Append(~numberoffaces,#Faces);
 Append(~facelist,Faces);
 FaceForms:=[]; 
 AutFF:=[]; 
 facevect:=[]; 
 for F in Faces do 
  FF:=Parent(Pini) ! 0;
  for k in F do 
   Fk:= HermitianTranspose(Sk[k])*Sk[k]; 
   FF := FF+ Fk;
  end for;
  gL:=findperp2([Sk[k] : k in F]);
  if #gL eq 1 then 
   Append(~facevect,gL);
  end if;
 end for;
 Append(~facevectList,facevect);
 
 count:=0;

 //print Faces;

 PerfectNeighbours:=[**];    //List of perfect neighbours of P being treated
 CriticalValues:=[**];       //List of critical rho-values of P
 while(#Faces gt 0) do
  count:=count+1;
  F1:=findperp1([Sk[n] : n in Faces[1]]);
  Append(~FacetVectorList,F1);  //[??]
  Exclude(~Faces,Faces[1]);

  sgn:=Sign(&+ [Rationals()!auswerten(F1,x) : x in Sk]);
  F1:=sgn*F1;

  tsup:=10000;
  tinf:=0;
  t:=(tinf+tsup)/2;
  minimcont:=0;
  while minimcont ne 1 do
   coherent:=false;
   while not coherent do
    Pt:=P+t*F1;
    M:=1;
    while M eq 1 do
     if IsPositiveDefinite(spurform(Pt)) then
      M:=hermitianmin(Pt);
      if M eq 1 then
       tinf:=t;
       t:=(tinf+tsup)/2;
       Pt:=P+t*F1;
      end if;
     else
      tsup:=t;
      t:=(tinf+tsup)/2;
      Pt:=P+t*F1;
     end if;
    end while;
    St:=minvecs(Pt);
    SFace:=[ s : s in Sk | auswerten(F1,s) eq 0];

    Cond:=[projzeileNorm(s) : s in SFace] cat [projzeileNorm(s) : s in St];
    Uns:=Vector( #Cond , [ F!(idealnorm(v)) : v in SFace ] cat [F!(idealnorm(v)) : v in St] );
    Cond:=Transpose(Matrix(Cond));

    coherent:=IsConsistent(Cond,Uns);
    if not coherent then
     tsup:=t;
     t:=(tinf+tsup)/2;
     Pt:=P+t*F1;
    end if;
   end while;
   Pcont:=ListToSmallMatrixNorm(Solution(Cond,Uns));
   Pconte:=spurform(Pcont);
   Lcont:=LatticeWithGram(Pconte);
   //Scont:=ShortVectors(Lcont,hermitianmin(Pcont)/mmax,hermitianmin(Pcont)*mmax);
   Scontk:=minvecs(Pcont);

   minimcont:=hermitianmin(Pcont);

   tsup:=t;
   t:=(tinf+tsup)/2;
   Pt:=P+t*F1;
  end while;

  Append(~PerfectNeighbours,Pcont);

  //Determine critical value rho:
  C:=Pcont-P;
  I:=0; J:=0;
  for i in [1..n] do
   for j in [1..n] do
    if C[i][j] ne 0 then
     I:=i; J:=j;
     break i;
    end if;
   end for;
  end for;
  Append(~CriticalValues , sgn*(C[I][J])/(F1[I][J]) );
  

  iso:=false;
  for i in [1..#perfectlist] do
   if isisom(Pcont,perfectlist[i]) then
    iso:=true;
    Include(~E,<Position(perfectlist,P),i>);
   end if;
  end for;
  if not iso then
   Append(~perfectlist,Pcont);
   Append(~Todo,Pcont);
   Include(~E,<Position(perfectlist,P),Position(perfectlist,Pcont)>);
  end if;
 end while;
 Append(~PerfectNeighbourList,PerfectNeighbours);
 Append(~CriticalValueList,CriticalValues);
end while;

//Calculation of perfect forms ends here

autlist:=[#aut(p) : p in perfectlist];
//print autlist;
//Write lists in order to identify automorphism groups with GAP
if Maximum(autlist) lt 2000 then
 L1:=Open("Liste1","w");
 Puts(L1, "Liste1:=[");
 for i in [1..#perfectlist] do
  if i ne #perfectlist then
   Puts(L1, IntegerToString(IdentifyGroup(aut(perfectlist[i]))[1]) cat " , ");
  else
   Puts(L1, IntegerToString(IdentifyGroup(aut(perfectlist[i]))[1]));
  end if;
 end for;
 Puts(L1, " ];");
 delete L1;

 L2:=Open("Liste2","w");
 Puts(L2, "Liste2:=[");
 for i in [1..#perfectlist] do
  if i ne #perfectlist then
   Puts(L2, IntegerToString(IdentifyGroup(aut(perfectlist[i]))[2]) cat " , ");
  else
   Puts(L2, IntegerToString(IdentifyGroup(aut(perfectlist[i]))[2]));
  end if;
 end for;
 Puts(L2, " ];");
 delete L2;

end if;

FaceTupleList:=[**];                           //This will contain the tuples of faces of codim>1
Representatives:=[* perfectlist *];            //Representatives of all minimal classes
FaceTupleListOfRepresentatives:=[* [] *];      //We need this to check the dimension of the intersection

print "Starting the computation of minimal classes. Please be patient.";

//n=2

if n eq 2 then

 //Generate the tuples
 for i in [1..#perfectlist] do
  FaceTuples:=[**];
  S:=minvecs(perfectlist[i]);
  for j in [1..#facelist[i]] do
   for k in [j+1..#facelist[i]] do
    Intersection:=facelist[i][j] meet facelist[i][k];
    if #Intersection gt 1 then
     if #findperp2([S[l] : l in Intersection]) eq 2 then
      Append(~FaceTuples,[j,k]);
      if k gt #facelist[i] then print "Error."; end if;
     end if;
    end if;
   end for;
  end for;
  Append(~FaceTupleList,FaceTuples);
 end for;
 FaceTupleList:= [* FaceTupleList *];       //This is a bit clumsy; now FTL[1] is the data for codim 2

 //Compute corank 1 classes:

 TempList:=[**];
 for i in [1..#perfectlist] do
  for j in [1..#CriticalValueList[i]] do
   Append(~TempList , perfectlist[i]+(CriticalValueList[i][j]/2)*facevectList[i][j][1] );
  end for;
 end for;

 MinClassReps:=[TempList[1]];
 for x in TempList do
  isi:=false;
  for y in MinClassReps do
   if AreEquivalentMinimalClasses(x,y) then
    isi:=true;
   end if;
  end for;
  if not isi then
   Append(~MinClassReps,x);
  end if;
 end for;

 Append(~Representatives,MinClassReps);
 print "Corank 1 done.";

 //Compute classes of corank >= 2

 codim:=2;

 while n^2-codim ge n do
  TempList:=[**];
  for i in [1..#perfectlist] do
   for j in [1..#FaceTupleList[codim-1][i]] do
    T:=MatrixRing(K,n)!perfectlist[i];
    for k in FaceTupleList[codim-1][i][j] do
     T:=T+(CriticalValueList[i][k]/(2*codim))*facevectList[i][k][1];
    end for;
    Append(~TempList,T);
   end for;
  end for;

  MinClassReps:=[TempList[1]];
  for x in TempList do
   isi:=false;
   for y in MinClassReps do
    if AreEquivalentMinimalClasses(x,y) then
     isi:=true;
    end if;
   end for;
   if not isi then
    Append(~MinClassReps,x);
   end if;
  end for;

  Append(~Representatives,MinClassReps);
  codim:=codim+1;
 end while;
end if;


//n=3

if n eq 3 then
 //Generate the tuples
 print "n=3";
 print "Starting to assemble the tuples";

 for i in [1..#perfectlist] do
  FaceTuples:=[**];
  S:=minvecs(perfectlist[i]);
  for j in [1..#facelist[i]] do
   for k in [j+1..#facelist[i]] do
    Intersection:=facelist[i][j] meet facelist[i][k];
    if #Intersection gt 1 then
     if #findperp2([S[l] : l in Intersection]) eq 2 then
      Append(~FaceTuples,[j,k]);
      if k gt #facelist[i] then print "Error."; end if;
     end if;
    end if;
   end for;
  end for;
  Append(~FaceTupleList,FaceTuples);
 end for;
 FaceTupleList:= [* FaceTupleList *];       //This is a bit clumsy; now FTL[1] is the data for codim 2

 codim:=3;
 print "Now I've set codim to 3. FaceTupleList has " cat IntegerToString(#FaceTupleList) cat " entries at present.";
 print "Representatives has " cat IntegerToString(#Representatives) cat " entries at present.";
 while n^2-codim ge limit do
  print "Now doing it for Codim " cat IntegerToString(codim);
  FaceTuplesInCodim:=[**];
  for i in [1..#perfectlist] do
   FaceTuples:=[**];
   S:=minvecs(perfectlist[i]);
   for j in [1..#FaceTupleList[codim-2][i]] do
    Intersection:=facelist[i][FaceTupleList[codim-2][i][j][1]];
    for l in FaceTupleList[codim-2][i][j] do
     Intersection:=facelist[i][l] meet Intersection;
    end for;
    for k in [1..#facelist[i]] do
     IntersectionTemp:=facelist[i][k] meet Intersection;
     if #IntersectionTemp gt 1 then
      if #findperp2([S[l] : l in IntersectionTemp]) eq codim then
       L:=FaceTupleList[codim-2][i][j];
       Append(~L,k);
       Append(~FaceTuples, L  );
      end if;
     end if;
    end for;
   end for;
   Append(~FaceTuplesInCodim,FaceTuples);
  end for;
  Append(~FaceTupleList,FaceTuplesInCodim);
  codim:=codim+1;
 end while; 
 print "Tuples done. FaceTupleList has " cat IntegerToString(#FaceTupleList) cat " entries.";
 //Compute corank 1 classes:
 print "Starting the computation of codim 1 classes.";
 TempList:=[**];
 for i in [1..#perfectlist] do
  for j in [1..#CriticalValueList[i]] do
   Append(~TempList , perfectlist[i]+(CriticalValueList[i][j]/2)*facevectList[i][j][1] );
  end for;
 end for;

 print "TempList done. Starting equivalence testing for codim 1.";

 MinClassReps:=[TempList[1]];
 for x in TempList do
  isi:=false;
  for y in MinClassReps do
   if AreEquivalentMinimalClasses(x,y) then
    isi:=true;
   end if;
  end for;
  if not isi then
   Append(~MinClassReps,x);
  end if;
 end for;

 Append(~Representatives,MinClassReps);
 print "Corank 1 done.";
 //Compute classes of corank >= 2

 codim:=2;

 print "Starting codim" cat IntegerToString(codim) cat " computations.";

 while n^2-codim ge limit do
  TempList:=[**];
  for i in [1..#perfectlist] do
   for j in [1..#FaceTupleList[codim-1][i]] do
    T:=MatrixRing(K,n)!perfectlist[i];
    for k in FaceTupleList[codim-1][i][j] do
     T:=T+(CriticalValueList[i][k]/(2*codim))*facevectList[i][k][1];
    end for;
    if not IsPositiveDefinite(spurform(T)) then print "Error, not pos.def.";end if;
    if pcorank(T) ne codim then print "Wrong pcorank at " cat IntegerToString(j); break i; codim:=n^3; end if;
    Append(~TempList,T);
   end for;
  end for;

  print "TempList for Codim " cat IntegerToString(codim) cat " done. It has " cat IntegerToString(#TempList) cat " entries. Starting equivalence testing.";

  MinClassReps:=[TempList[1]];
  counter:=1;
  for x in TempList do
   counter+:=1;
   isi:=false;
   for y in MinClassReps do
    if AreEquivalentMinimalClasses(x,y) then
     isi:=true;
    end if;
   end for;
   if not isi then
    Append(~MinClassReps,x);
   end if;
   if counter mod Floor(#TempList/5) eq 0 then
    print "Ca. #TempList/5 forms checked. Please remain patient.";
   end if;
  end for;

  Append(~Representatives,MinClassReps);
  print "Codim " cat IntegerToString(codim) cat " done.";
  codim:=codim+1;
 end while;
end if;

//Create a generating set of GL(L) as Z-matrices
X:=[];
for p in perfectlist do
 X:=X cat [MatrixRing(Integers(),2*n)!x : x in Generators(aut(p))];
end for;
for L in PerfectNeighbourList do
 for A in L do
  for p in perfectlist do
   a,b:=isisom(A,p);
   if a then
    Append(~X,MatrixRing(Integers(),2*n)!b);
   end if;
  end for;
 end for;
end for;
//"Clean up" X:
while MatrixRing(Integers(),2*n)!1 in X do
 Remove(~X,Position(X,MatrixRing(Integers(),2*n)!1));
end while;
ZGENS:=X;             //Z-generating system
OKGENS:=matbas(X);    //OK-generating system

load mydatafunctions;
print "Data assembled.";
load testconjugacy;