// This MAGMA script is to deal with the loose ends -
// namely forms among the 24888 whose smallest
// odd exception is not known to be required
// It computes all escalations of these and
// checks to see if any of these fail to represent
// numbers on checknumlist.

R<q> := PowerSeriesRing(RationalField());
load "readable_quaternaries.txt";
//load "nonlocuniv.txt";

//printf "Enter lattice number to check:";
//readi ind;
ind := 8819;
checknumlist := [133];

M := quatver[ind];

thetaprec := 2*checknumlist[#checknumlist]+2;
/*
done := false;
while done eq false do
  printf "Enter number to check (0 to quit):";
  readi numinp;
  if (numinp ne 0) then
    Append(~checknumlist,numinp);
  else
    done := true;
  end if;
end while;
*/
//M := nonlocuniv[3];

// U(d) operator

function U(f,d)
  start := Valuation(f);
  en := AbsolutePrecision(f)-1;
  ret := 0;
  for n in [start..en] do
    if n mod d eq 0 then
      ret := ret + Coefficient(f,n)*q^(Floor(n/d));
    end if;
  end for;
  ret := ret + BigO(q^Floor((AbsolutePrecision(f)+d-1)/d));
  return ret;
end function;

// V(d) operator

function V(f,d)
  start := Valuation(f);
  en := AbsolutePrecision(f)-1;
  ret := 0;  
  for n in [start..en] do
    ret := ret + Coefficient(f,n)*q^(n*d);
  end for;
  ret := ret + BigO(q^(AbsolutePrecision(f)*d));
  return ret;
end function;

// IsNorm returns yes if there is a vector v of norm n
// in the lattice L

function IsNorm(L,n)
  ret := false;
  P := ShortVectorsProcess(L,n,n);
  done := false;
  found := false;
  while done eq false do
    if IsEmpty(P) then
      done := true;
    else
      v, norm := NextVector(P);
      if norm eq n then
        done := true;
        found := true;
      end if;
    end if;
  end while;
  if found eq true then
    ret := true;
  end if;
  return ret;
end function;

L := LatticeWithGram(2*M);
T := U(R!ThetaSeries(L,5001),2);

RRR<x,y,z,w> := PolynomialRing(RationalField(),4);
poly := M[1][1]*x^2 + 2*M[1][2]*x*y + 2*M[1][3]*x*z + 2*M[1][4]*x*w;
poly := poly + M[2][2]*y^2 + 2*M[2][3]*y*z + 2*M[2][4]*y*w;
poly := poly + M[3][3]*z^2 + 2*M[3][4]*z*w + M[4][4]*w^2;
printf "The form is %o.\n",poly;

exceplist := [];
for n in [1..2500] do
  if (Coefficient(T,n) eq 0) and (n mod 3 ne 0) then
    Append(~exceplist,n);
    printf "The form fails to represent %o.\n",n;
  end if;
end for;

t := exceplist[1];

printf "We will escalate by %o.\n",t;

// Option 1 - Escalate the form by t

MM1 := RMatrixSpace(RationalField(),1,1);

M2 := DiagonalJoin(M,MM1![t]);
numnewlats := (&*[ (Floor(2*Sqrt(t*M2[i,i]))+1) : i in [1..4]]);
printf "There are %o lattices that we must consider.\n",numnewlats;

count := 0;
actual := 0;
numexcep := 0;
foundtrus := {};
exceplats := [];
for x2 in [0..Floor(2*Sqrt(t*M2[1,1]))] do
  x := x2/2;
  M2[1,5] := x;
  M2[5,1] := x;
  for y2 in [0..Floor(2*Sqrt(t*M2[2,2]))] do
    y := y2/2;
    M2[2,5] := y;
    M2[5,2] := y;
    for z2 in [0..Floor(2*Sqrt(t*M2[3,3]))] do
      z := z2/2;
      M2[3,5] := z;
      M2[5,3] := z;
      for w2 in [0..Floor(2*Sqrt(t*M2[4,4]))] do
        w := w2/2; 
        M2[4,5] := w;
        M2[5,4] := w; 
        count := count + 1;
        if count mod 100000 eq 0 then
          printf "Checking lattice %o.\n",count;
        end if;
        if IsPositiveSemiDefinite(M2) then
          if Determinant(M2) eq 0 then
            n2, temp, rank := LLLGram(M2);
            n2 := Submatrix(n2,1,1,rank,rank);
            n2 := LLLGram(n2);
          else
            n2 := M2;
          end if;
          actual := actual + 1;
          fails := {};
          L := LatticeWithGram(2*n2);
          for mm in [1..#checknumlist] do
            if not IsNorm(L,2*checknumlist[mm]) then
              Include(~fails,checknumlist[mm]);
            end if;
          end for;
          if #fails gt 0 then
            numexcep := numexcep + 1;
            printf "Lattice %o has gram matrix:\n",count;
            printf "M = %o.\n",M2;
            printf "Exceptions are %o.\n",fails;
            Append(~exceplats,<M2,fails>);
          end if;
        end if;
      end for;
    end for;
  end for;
end for;
printf "There were %o total options for escalators.\n",numnewlats;
printf "Of these %o were positive definite.\n",actual;
printf "Exceptions were found for %o of these.\n",numexcep;
printf "They are %o.\n",exceplats;
quit;