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;

// makeunilist() makes a list of all possible integer-coefficient universal
// quaternary quadratic forms, by use of the 290-Theorem.  A running count
// of universal forms is printed as they are found.  The total number of
// universal forms is discovered to be 6436.

makeunilist:= function()
// unicount is a running count of the universal forms found so far;
// unilist[d] will contain all universal forms of discriminant d
 unicount:=0;  
 unilist:= [ [] : i in [1..5000] ];   
 for n11 in [1..1] do
 for n22 in [n11..2] do   
 for n33 in [n22..5] do
 for n44 in [n33..31] do
 for m12 in [0..n11] do           n12:=m12/2; 
 for m13 in [n11..-n11 by -1] do  n13:=m13/2; 
 for m14 in [n11..-n11 by -1] do  n14:=m14/2;
 for m23 in [0..n22] do           n23:=m23/2; 
 for m24 in [n22..-n22 by -1] do  n24:=m24/2;
 for m34 in [0..n33] do           n34:=m34/2;  
  n:=M4![n11,n12,n13,n14,n12,n22,n23,n24,n13,n23,n33,n34,n14,n24,n34,n44];
  if IsPositiveDefinite(n) and (truant(n,290) eq 291) then 
   d:=Z!(16*Determinant(n)); Append(~unilist[d],n);
   if isrepeat(unilist[d]) then Remove(~unilist[d],#unilist[d]); 
    else unicount:=unicount+1; print unicount; end if;
  end if;
 end for; end for; end for; end for; 
 end for; end for; end for;
 end for; end for;
 end for; 
 return &cat unilist;
end function;

// The list of universal forms is computed below and assigned to unilist.

unilist:=makeunilist();