// This MAGMA script is a modification of
// Bhargava's and Hanke's script. It
// finds integer-valued escalations of
// the zero dimensional lattice. It finds
// unary, binary, ternary and quaternary
// escalator lattices, and writes these to
// files. It also sorts the ternaries
// into universal and non-universal.

Q:=RationalField();
Z:=Integers();

V0:=RMatrixSpace(Q,1,0);
V1:=RMatrixSpace(Q,1,1);
V2:=RMatrixSpace(Q,1,2);
V3:=RMatrixSpace(Q,1,3);
V4:=RMatrixSpace(Q,1,4);
V5:=RMatrixSpace(Q,1,5);

// All quadratic forms are stored as symmetric square matrices in M0-M5.

M0:=RMatrixSpace(Q,0,0);
M1:=RMatrixSpace(Q,1,1);
M2:=RMatrixSpace(Q,2,2);
M3:=RMatrixSpace(Q,3,3);
M4:=RMatrixSpace(Q,4,4);
M5:=RMatrixSpace(Q,5,5);

// Constant specifiying how high to check

check := 319;

// P(n) is a predicate that determines whether the integer n
// is in the desired set of integers. The predicate below
// checks odd integers.

P := function(n)
  if n mod 2 eq 1 then
    ret := true;
  else
    ret := false;
  end if;
  return ret;
end function;

// truant(m,maxcheck) computes the smallest number less than or equal to
// maxcheck that is not represented by m.  If all numbers less than or
// equal to maxcheck are represented by m, then maxcheck+1 is returned.

truant := function(m,maxcheck)
  thetalist:=ElementToSequence(ThetaSeries(LatticeWithGram(2*m),2*maxcheck));
  min := maxcheck+1;
  for i in [3..#thetalist by 2] do
    chk := (i-1)/2;
    chk := Z ! chk;         
    if thetalist[i] eq 0 and P(chk) then
      min := chk;              
      break;
    end if;
  end for;
  return min;
end function;

// isrepeat(mlist) takes a list of forms and determines whether the last
// form in the list is isomorphic to any previous element in the list.
// If it is isomorphic to some previous element, then true is returned,
// otherwise false is returned.

function isrepeat(mlist)
  rep:=false;
  for i in [1..#mlist-1] do
    if IsIsometric(LatticeWithGram(2*mlist[i]),LatticeWithGram(2*mlist[#mlist])) 
    then rep:=true;  break;
    end if;
  end for;
  return rep;
end function;

// For k = 0, 1, 2, or 3, escalatek(mlist) takes a list of dimension k
// forms, and write the list of new higher dimensional escalators
// to e(k+1).txt and the list of same dimensional escalators to
// eknew.txt.

function escalate1(mlist)
  // Allocate space for a list of escalators ordered by discriminant
  elist := [ [] : i in [1..20000] ];
  elist2 := [ [] : i in [1..20000] ];
  PrintFile("e2.txt","e2:=[");
  PrintFile("e1new.txt","e1new:=[");
  count := 0;
  count2 := 0;
  for i in [1..#mlist] do
    printf "Escalating form %o.\n",i;
    m := mlist[i];
    t := truant(m,check);
    n := DiagonalJoin(m,M1![t]);
    for x2 in [0..Floor(2*Sqrt(t*n[1,1]))] do x:=x2/2;
      n[1,2]:=x; n[2,1]:=x;
      if IsPositiveSemiDefinite(n) then 
        if Determinant(n) eq 0 then
          n2, temp, rank := LLLGram(n);
          n2 := Submatrix(n2,1,1,rank,rank);
          d := Z!(2*Determinant(n2)); 
          n2 := LLLGram(n2);
          Append(~elist2[d],n2);
          if isrepeat(elist2[d]) then 
            Remove(~elist2[d],#elist2[d]);
          else
            count2 := count2 + 1;
            if count2 gt 1 then
              PrintFile("e1new.txt",",");
            end if;
            PrintFile("e1new.txt","Matrix(1,1,");
            PrintFile("e1new.txt",ElementToSequence(n2));
            PrintFile("e1new.txt",")");    
          end if;
        end if;
        if Determinant(n) ne 0 then
          d:=Z!(4*Determinant(n)); 
          n2 := LLLGram(n);
          Append(~elist[d],n2);
          if isrepeat(elist[d]) then 
            Remove(~elist[d],#elist[d]);
          else
            count := count + 1;
            if count gt 1 then
              PrintFile("e2.txt",",");
            end if;
            PrintFile("e2.txt","Matrix(2,2,");
            PrintFile("e2.txt",ElementToSequence(n2));
            PrintFile("e2.txt",")");
          end if;  
        end if;
      end if;
    end for;
  end for;
  PrintFile("e2.txt","];");
  PrintFile("e1new.txt","];");
  printf "%o new forms of dimension 2 found.\n",count;
  printf "%o new forms of dimension 1 found.\n",count2;
  return count;
end function;  

function escalate2(mlist)
  // Allocate space for a list of escalators ordered by discriminant
  elist := [ [] : i in [1..20000] ];
  elist2 := [ [] : i in [1..20000] ];
  PrintFile("e3.txt","e3:=[");
  PrintFile("e2new.txt","e2new:=[");
  count := 0;
  count2 := 0;
  for i in [1..#mlist] do
    printf "Escalating form %o.\n",i;
    m := mlist[i];
    t := truant(m,check);
    n := DiagonalJoin(m,M1![t]);
    for x2 in [-Floor(2*Sqrt(t*n[1,1]))..Floor(2*Sqrt(t*n[1,1]))] do x:=x2/2;
    for y2 in [0..Floor(2*Sqrt(t*n[2,2]))] do y:=y2/2;
      n[1,3]:=x; n[3,1]:=x; n[2,3]:=y; n[3,2]:=y;
      if IsPositiveSemiDefinite(n) then 
        if Determinant(n) eq 0 then
          n2, temp, rank := LLLGram(n);
          n2 := Submatrix(n2,1,1,rank,rank);
          d := Z!(4*Determinant(n2)); 
          n2 := LLLGram(n2);
          Append(~elist2[d],n2);
          if isrepeat(elist2[d]) then 
            Remove(~elist2[d],#elist2[d]);
          else
            count2 := count2 + 1;
            if count2 gt 1 then
              PrintFile("e2new.txt",",");
            end if;
            PrintFile("e2new.txt","Matrix(2,2,");
            PrintFile("e2new.txt",ElementToSequence(n2));
            PrintFile("e2new.txt",")");    
          end if;
        end if;
        if Determinant(n) ne 0 then
          d:=Z!(8*Determinant(n)); 
          n2 := LLLGram(n);
          Append(~elist[d],n2);
          if isrepeat(elist[d]) then 
            Remove(~elist[d],#elist[d]);
          else
            count := count + 1;
            if count gt 1 then
              PrintFile("e3.txt",",");
            end if;
            PrintFile("e3.txt","Matrix(3,3,");
            PrintFile("e3.txt",ElementToSequence(n2));
            PrintFile("e3.txt",")");
          end if;  
        end if;
      end if;
    end for;
    end for;
  end for;
  PrintFile("e3.txt","];");
  PrintFile("e2new.txt","];");
  printf "%o new forms of dimension 3 found.\n",count;
  printf "%o new forms of dimension 2 found.\n",count2;
  return count;
end function;  

function escalate3(mlist)
  // Allocate space for a list of escalators ordered by discriminant
  elist := [ [] : i in [1..20000] ];
  elist2 := [ [] : i in [1..20000] ];
  PrintFile("e4.txt","e4:=[");
  PrintFile("e3new.txt","e3new:=[");
  count := 0;
  count2 := 0;
  for i in [1..#mlist] do
    printf "Escalating form %o.\n",i;
    m := mlist[i];
    t := truant(m,check);
    n := DiagonalJoin(m,M1![t]);
    for x2 in [-Floor(2*Sqrt(t*n[1,1]))..Floor(2*Sqrt(t*n[1,1]))] do x:=x2/2;
    for y2 in [-Floor(2*Sqrt(t*n[2,2]))..Floor(2*Sqrt(t*n[2,2]))] do y:=y2/2;
    for z2 in [0..Floor(2*Sqrt(t*n[3,3]))] do z:=z2/2;
      n[1,4]:=x; n[4,1]:=x; n[2,4]:=y; n[4,2]:=y; n[3,4]:=z; n[4,3]:=z;
      if IsPositiveSemiDefinite(n) then 
        if Determinant(n) eq 0 then
          n2, temp, rank := LLLGram(n);
          n2 := Submatrix(n2,1,1,rank,rank);
          d := Z!(8*Determinant(n2)); 
          n2 := LLLGram(n2);
          Append(~elist2[d],n2);
          if isrepeat(elist2[d]) then 
            Remove(~elist2[d],#elist2[d]);
          else
            count2 := count2 + 1;
            if count2 gt 1 then
              PrintFile("e3new.txt",",");
            end if;
            PrintFile("e3new.txt","Matrix(3,3,");
            PrintFile("e3new.txt",ElementToSequence(n2));
            PrintFile("e3new.txt",")");    
          end if;
        end if;
        if Determinant(n) ne 0 then
          d:=Z!(16*Determinant(n)); 
          n2 := LLLGram(n);
          Append(~elist[d],n2);
          if isrepeat(elist[d]) then 
            Remove(~elist[d],#elist[d]);
          else
            count := count + 1;
            if count gt 1 then
              PrintFile("e4.txt",",");
            end if;
            PrintFile("e4.txt","Matrix(4,4,");
            PrintFile("e4.txt",ElementToSequence(n2));
            PrintFile("e4.txt",")");
          end if;  
        end if;
      end if;
    end for;
    end for;
    end for;
  end for;
  PrintFile("e4.txt","];");
  PrintFile("e3new.txt","];");
  printf "%o new forms of dimension 4 found.\n",count;
  printf "%o new forms of dimension 3 found.\n",count2;
  return count;
end function;  
  
// The escalators of dimension 1, 2, 3, 4 are computed and stored in
// e1, e2, e3, e4 respectively.  

for n in [1..check] do
  if P(n) then
    e1 := [M1![n]];
    break;
  end if;
end for;

printf "There is one non-universal escalator of step 1, dim 1.\n";
printf "done e1\n";

escalate1(e1);
printf "Done with step 1 escalation\n";
load "e2.txt";

escalate2(e2);
printf "Done with step 2 escalation\n";
load "e3.txt";
printf "Testing universality.\n";

univ3 := [];
nonuniv3 := [];
univ_count := 0;
nonuniv_count := 0;
PrintFile("univ3.txt","univ3 := [");
PrintFile("nonuniv3.txt","nonuniv3 := [");
for n in [1..#e3] do
  if truant(e3[n],check) gt check then
    univ_count := univ_count + 1;
    if univ_count gt 1 then
      PrintFile("univ3.txt",",");
    end if;
    PrintFileMagma("univ3.txt",e3[n]);
  else
    nonuniv_count := nonuniv_count + 1;
    if nonuniv_count gt 1 then
      PrintFile("nonuniv3.txt",",");
    end if;
    PrintFileMagma("nonuniv3.txt",e3[n]);
  end if;
end for;
PrintFile("univ3.txt","];");
PrintFile("nonuniv3.txt","];");
printf "%o universal forms found.\n",univ_count;
printf "%o non-universal forms found.\n",nonuniv_count;

load "nonuniv3.txt";
escalate3(nonuniv3);
printf "Done with step 3 escalation\n";

load "e3new.txt";

printf "Checking if the 3-dimensional lattices found are new.\n";
new := false;
for n in [1..#e3] do
  Append(~e3,e3new[n]);
  if isrepeat(e3) then
    Remove(~e3,#e3);
  else
    printf "Lattice %o is new.\n",n;
    new := true;
  end if;
end for;
if new eq false then
  printf "No new 3-dimensional lattices found.\n";
end if;