// This MAGMA script is for classifying
// quaternaries that represent all odd numbers.

load "isembeddable.txt";
load "univ3.txt";

M := ZeroMatrix(Rationals(),4,4);
M[1][1] := 1;
chklist := [1,3,5,7,11,13,15,17,19,21,23,29,31,33,35,37,39,41,47,51,53,57,59,77,83,85,87,89,91,93,105,119,123,133,137,143,145,187,195,203,205,209,231,319,385,451];

// Generate matrices in Minkowski-reduced form
// This means that a22 >= a11, a33 >= a22, a44 >= a33,
// a12 >= 0, a23 >= 0, a34 >= 0 as well as other inequalities
// of the form f(m_1,m_2,m_3,m_4) >= a_ii
// where m_i = 1 and all other m_js are 1, 0 or -1

ternlatlist := [ LatticeWithGram(2*univ3[m]) : m in [1..#univ3]];

a11 := 1;
matcount1 := 0; 
matcount2 := 0;
matcount3 := 0;  
matcount4 := 0;  
matcount5 := 0;
latlist := [ [] : j in [1..50000]];
matlist := [];
for a22 in [a11..3] do
  M[2][2] := a22;
  for a12 in [0..a11] do
    M[1][2] := a12/2;
    M[2][1] := a12/2;
    M22 := Submatrix(M,1,1,2,2);
    if IsPositiveDefinite(M22) then
      for a33 in [a22..7] do
        M[3][3] := a33;
        for a13 in [-a11..a11] do
          M[1][3] := a13/2;
          M[3][1] := a13/2;
          m23min := Max([0,-a22,-a11-a22-a12-a13,-a11-a22+a12+a13]);
          m23max := Min([a22,a11+a22-a12+a13,a11+a22+a12-a13]);
          for a23 in [m23min..m23max] do
            M[2][3] := a23/2;
            M[3][2] := a23/2;
            M33 := Submatrix(M,1,1,3,3);
            if IsPositiveDefinite(M33) then
              for a44 in [a33..77] do
                M[4][4] := a44;
                for a14 in [-a11..a11] do
                  M[1][4] := a14/2;
                  M[4][1] := a14/2;
                  a24lowbound := Max([-a22,-a11-a22-a12-a14,-a11-a22+a12+a14]);
                  a24uppbound := Min([a22,a11+a22-a12+a14,a11+a22+a12-a14]);
                  for a24 in [a24lowbound..a24uppbound] do
                    M[2][4] := a24/2;
                    M[4][2] := a24/2;
                    a34lowbound := Max([0,-a33,-a11-a33-a13-a14,-a22-a33-a23-a24,-a11-a33+a13+a14,-a22-a33+a23+a24,-a11-a22-a33-a12-a13-a14-a23-a24,-a11-a22-a33+a12+a13+a14-a23-a24,-a11-a22-a33+a12-a13-a14+a23+a24,-a11-a22-a33-a12+a13+a14+a23+a24]);
                    a34uppbound := Min([a33,a11+a33-a13+a14,a11+a33+a13-a14,a22+a33-a23+a24,a22+a33+a23-a24,a11+a22+a33+a12-a13+a14-a23+a24,a11+a22+a33-a12+a13-a14-a23+a24,a11+a22+a33-a12-a13+a14+a23-a24,a11+a22+a33+a12+a13-a14+a23-a24]);
                    for a34 in [a34lowbound..a34uppbound] do
                      M[3][4] := a34/2;
                      M[4][3] := a34/2;
                      matcount1 := matcount1 + 1;
                      if (matcount1 mod 10^3) eq 0 then
                        printf "Checked %o matrices : %o positive definite, %o universal, %o not ternary.\n",matcount1,matcount2,matcount4,matcount5;
                        printf "Current gram is %o.\n",M;
                      end if;
                      if IsPositiveDefinite(M) then
                        matcount2 := matcount2 + 1;
                        // Check odd universality
                        L := LatticeWithGram(2*M);
                        T := ThetaSeries(L,903);  
                        isodduniv := true;
                        for m in [1..#chklist] do
                          if Coefficient(T,2*chklist[m]) eq 0 then
                            isodduniv := false;
                          end if;
                        end for;
                        if isodduniv eq true then
                          matcount4 := matcount4 + 1;
                          found := false;
                          done := false;
                          it := 1;
                          while done eq false do
                            if IsEmbeddable(ternlatlist[it],L) then
                              found := true;
                              done := true;
                            else
                              it := it + 1;
                              if (it gt 23) then
                                done := true;
                              end if;
                            end if;
                          end while;
                          if (found eq false) then
// Check to see if it is a repeat
                            D := Integers()!Determinant(2*M);
                            isiso := false;
                            done2 := false;
                            it := 0;
                            while done2 eq false do
                              it := it + 1;
                              if (it gt #latlist[D]) then
                                done2 := true;
                              else
                                if IsIsometric(L,latlist[D][it]) then
                                  isiso := true;
                                  done2 := true;
                                end if;
                              end if;
                            end while;
                            if isiso eq false then
                              matcount5 := matcount5 + 1;
                              Append(~latlist[D],L);
                              Append(~matlist,M);
                            end if;
                          end if;
                        end if;
                      end if;
                    end for;
                  end for;
                end for;
              end for;
            end if;
          end for;
        end for;
      end for;
    end if;
  end for;
end for;
printf "Total number of matrices is %o.\n",matcount1;
printf "Total number of positive definite is %o.\n",matcount2;
printf "Total number of universal is %o.\n",matcount4;
printf "Total number of non-ternary is %o.\n",matcount5;
PrintFile("univquats3.txt","univquats := ");
PrintFileMagma("univquats3.txt",matlist);
PrintFile("univquats3.txt",";");