// This script computes all escalations
// x^2+3y^2+5z^2 by 15. These are
// auxiliary escalators. The output is
// written to 1115out.txt.

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);

// 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:=(#thetalist+1)/2;
 for i in [3..#thetalist by 2] do         
  if thetalist[i] eq 0 and i mod 4 eq 3 then              
   min:=(i-1)/2;
   break;
  end if;
 end for;
 return min;
end function;

// truant2(m,maxcheck) computes a list of all numbers less than or equal
// to maxcheck that are not represented.

truant2 := function(m,maxcheck)
  thetalist:=ElementToSequence(ThetaSeries(LatticeWithGram(2*m),2*maxcheck));
  exc := [];
  for i in [3..#thetalist by 2] do         
    if thetalist[i] eq 0 and i mod 4 eq 3 then              
      Append(~exc,Floor((i-1)/2));
    end if;
  end for;
  return exc;
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;

// escalate15 escalates the lattice x^2 + 3y^2 + 5z^2 with a vector of
// length 15 in an attempt to obtain forms that locally represent all
// odds.

function escalate15(mlist)
  // Allocate space for a list of escalators ordered by discriminant
  // For k = 3, discriminants higher than 10000 are needed.
 elist:=[ [] : i in [1..40000] ];
 elist2 := [ [] : i in [1..40000] ];
 for i in [1..#mlist] do
  printf "Escalating form %o.\n",i;
  count := 0;
  count2 := 0;
  m:=mlist[i];
  t:=15;
  n:=DiagonalJoin(m,M1![t]);
  nc:=NumberOfColumns(m);
  if nc eq 3 then 
   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
      count2 := count2 + 1;
      n2, temp, rank := LLLGram(n);
      n2 := Submatrix(n2,1,1,rank,rank);
      d := Z!(8*Determinant(n2)); Append(~elist2[d],LLLGram(n2));
      if isrepeat(elist2[d]) then 
       Remove(~elist2[d],#elist2[d]);
       count2 := count2 - 1;
      end if;
     end if;
     if Determinant(n) ne 0 then
      count := count + 1;
      d:=Z!(16*Determinant(n)); Append(~elist[d],LLLGram(n));  
     end if;
     if isrepeat(elist[d]) then 
      Remove(~elist[d],#elist[d]);
      count := count - 1;
    end if;  end if;
   end for; end for; end for;
  end if;
  printf "There were %o new forms of higher dimension.\n",count;
  if count2 ne 0 then
    printf "There were %o new forms of the same dimension.\n",count2;
  end if;
 end for;
 return &cat elist, &cat elist2;  
 // return Sort(&cat elist,func<x,y|Determinant(x)-Determinant(y)>);
end function;

// We want to escalate the form with the following matrix.

M := M3 ! Matrix([[1, 0, 0], [0, 3, 0], [0, 0, 5]]);
esc := [];
Append(~esc,M);
dim4, dim3 := escalate15(esc);

PrintFile("1115out.txt","switch1 := [");
for n in [1..#dim4] do
  if n gt 1 then
    PrintFile("1115out.txt",",");
  end if;
  PrintFile("1115out.txt","Matrix(4,4,");
  PrintFile("1115out.txt",ElementToSequence(dim4[n]));
  PrintFile("1115out.txt",")");  
end for;
PrintFile("1115out.txt","];");