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 then              
   min:=(i-1)/2;
   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, escalate(mlist) takes a list of dimension k
// forms, and returns a list of all isomorphism classes of possible 
// escalations of these forms.  The list of escalations is sorted by 
// discriminant.  For k<4, an escalation of a dimension k form always 
// turns out to have dimension k+1.

function escalate(mlist)
 elist:=[ [] : i in [1..10000] ];
 for i in [1..#mlist] do
  m:=mlist[i];
  t:=truant(m,290);
  n:=DiagonalJoin(m,M1![t]);
  nc:=NumberOfColumns(m);
  if nc eq 0 then Append(~elist[Z!(2*Determinant(n))],n); end if;
  if nc eq 1 then 
   for x2 in [0..2*Floor(Sqrt(t*n[1,1]))] do x:=x2/2;
    n[1,2]:=x; n[2,1]:=x;
    if IsPositiveSemiDefinite(n) then 
     d:=Z!(4*Determinant(n)); Append(~elist[d],LLLGram(n));  
     if isrepeat(elist[d]) then Remove(~elist[d],#elist[d]);
    end if;  end if;
   end for;
  end if;
  if nc eq 2 then 
   for x2 in [-2*Floor(Sqrt(t*n[1,1]))..2*Floor(Sqrt(t*n[1,1]))] do x:=x2/2;
   for y2 in [0..2*Floor(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 
     d:=Z!(8*Determinant(n)); Append(~elist[d],LLLGram(n));  
     if isrepeat(elist[d]) then Remove(~elist[d],#elist[d]);
    end if;  end if;
   end for; end for;
  end if;
  if nc eq 3 then 
   for x2 in [-2*Floor(Sqrt(t*n[1,1]))..2*Floor(Sqrt(t*n[1,1]))] do x:=x2/2;
   for y2 in [-2*Floor(Sqrt(t*n[2,2]))..2*Floor(Sqrt(t*n[2,2]))] do y:=y2/2;
   for z2 in [0..2*Floor(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 
     d:=Z!(16*Determinant(n)); Append(~elist[d],n);  
     if isrepeat(elist[d]) then Remove(~elist[d],#elist[d]);
    end if;  end if;
   end for; end for; end for;
  end if;
 end for;
 return &cat elist;  
 // return Sort(&cat elist,func<x,y|Determinant(x)-Determinant(y)>);
end function;

// The escalators of dimension 1, 2, 3, 4 are computed and stored in
// e1, e2, e3, e4 respectively.  

e1:=[M1![1]];
print "done e1";
e2:=escalate(e1);
print "done e2";
e3:=escalate(e2);
print "done e3";
e4:=escalate(e3);
print "done e4";