Z := Integers();

// This function coughs up a quadratic non-residue modulo
// an odd prime p.

function NonRes(p)
  done := false;
  q := 1;
  while done eq false do
    q := NextPrime(q);
    if LegendreSymbol(q,p) eq -1 then
      done := true;
    end if;
  end while;
  return q;
end function;

// This function takes a rational number a in Q_p and
// returns an integer that is in the same squareclass as a.

function MakeInt(a,p)
  pval := Valuation(a,p);
  a2 := a*p^(-pval);
  b := 1;
  if p gt 2 then
    a3 := FiniteField(p) ! a2;
    a4 := Z ! a3;
    if LegendreSymbol(a4,p) eq 1 then
      b := 1;
    else
      b := NonRes(p);
    end if;
  end if;
  if p eq 2 then
    a3 := Integers(8) ! a2;
    a4 := Z ! a3;
    b := a4 mod 8;
  end if;
  ret := b*p^(pval);
  return ret;
end function;

// This function takes as input a gram matrix M of the quadratic form
// Q = x^T M x, and outputs a matrix M2 whose entries are rational
// so that x^T M2 x is "split" and x^T M x and x^T M2 x are equivalent
// over Z_p, the ring of p-adic integers.

function GetLocal(M,p)
  n := NumberOfRows(M);
  done := false;
  while done eq false do
    // Find least p-divisible entry not already taken care of.
    // If there are none, then done = true
    minval := Valuation(0,p);
    loc := [];
    found := false;
    for a in [1..n] do
      for b in [1..n] do
        takencareof := false;
        // Check if it is taken care of
        if M[a][b] ne 0 then
          if a eq b then
            S := { M[a][c] : c in ([1..b-1] cat [b+1..n]) };
            S := S join { M[c][b] : c in ([1..a-1] cat [a+1..n]) };
            if S eq {0} or #S eq 0 then
              takencareof := true;
            end if;
            if a lt n then
              if M[a][a+1] ne 0 then
                bb := a+1;
                if p eq 2 then
                  if (M[a][a] eq M[bb][bb]) and (((M[a][a] eq 0) and (M[bb][bb] eq 0) and (M[a][bb] eq 2^Valuation(M[a][bb],2))) or ((M[a][a] eq M[bb][bb]) and (M[a][a] eq 2^Valuation(M[a][a],2)) and (M[a][bb] eq (1/2)*M[a][a]))) then
                    S1 :=
                    { M[a][c] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S2 :=
                    { M[c][a] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S3 :=
                    { M[b][c] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S4 :=
                    { M[c][b] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S := S1 join S2 join S3 join S4;
                    if S eq {0} or #S eq 0 then
                      takencareof := true;
                    end if;
                  end if;
                end if;
              end if;
            end if;
            if a gt 1 then
              if M[a-1][a] ne 0 then
                bb := a-1;
                if p eq 2 then
                  if (M[a][a] eq M[bb][bb]) and (((M[a][a] eq 0) and (M[bb][bb] eq 0) and M[a][bb] eq 2^Valuation(M[a][bb],2)) or ((M[a][a] eq M[bb][bb]) and (M[a][a] eq 2^Valuation(M[a][a],2)) and (M[a][bb] eq (1/2)*M[a][a]))) then
                    S1 :=
                      { M[a][c] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S2 :=
                      { M[c][a] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S3 :=
                      { M[b][c] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S4 :=
                      { M[c][b] : c in [1..Min(a,bb)-1] cat [Max(a,bb)+1..n] };
                    S := S1 join S2 join S3 join S4;
                    if S eq {0} or #S eq 0 then
                      takencareof := true;
                    end if;
                  end if;
                end if;
              end if;
            end if;
          end if;
          if a ne b then
            if AbsoluteValue(a-b) eq 1 and p eq 2 then
              bb := b;
              if (M[a][a] eq M[bb][bb]) and (((M[a][a] eq 0) and (M[bb][bb] eq 0) and (M[a][bb] eq 2^Valuation(M[a][bb],2))) or ((M[a][a] eq M[bb][bb]) and (M[a][a] eq 2^Valuation(M[a][a],2)) and (M[a][bb] eq (1/2)*M[a][a]))) then
                S1 :=
                { M[a][c] : c in [1..Min(a,b)-1] cat [Max(a,b)+1..n] };
                S2 :=
                { M[c][a] : c in [1..Min(a,b)-1] cat [Max(a,b)+1..n] };
                S3 :=
                { M[b][c] : c in [1..Min(a,b)-1] cat [Max(a,b)+1..n] };
                S4 :=
                { M[c][b] : c in [1..Min(a,b)-1] cat [Max(a,b)+1..n] };
                S := S1 join S2 join S3 join S4;
                if S eq {0} or #S eq 0 then
                  takencareof := true;
                end if;
              end if;
            end if;
          end if;
        else
          takencareof := true;
        end if;
        if takencareof eq false then
          //printf "Found bad spot at %o, %o.\n",a,b;
          //printf "Current matrix = %o.\n",M;
          found := true;
          if Valuation(M[a][b],p) lt minval then
            loc := [a,b];
            minval := Valuation(M[a][b],p);
          end if;
          if Valuation(M[a][b],p) eq minval and a eq b then
            loc := [a,b];
          end if;
        end if;
      end for;
    end for;
    if found eq false then
      done := true;
    else
      a := loc[1];
      b := loc[2];
      //printf "Bad spot at %o, %o.\n",a,b;
      //printf "Current matrix = %o.\n",M;
      if a eq b then
        // Change variables to clear diagonal
        for m in [1..n] do
          if m ne a then
            alpha := M[m][m];
            beta := M[m][a];
            gamma := M[a][a];
            c := -beta/gamma;
            // Make change of variables (x_m, x_a) -> (x_m, x_a + cx_m)
            AddColumn(~M,c,a,m);
            AddRow(~M,c,a,m);
          end if;
        end for;
      end if;
      if a ne b then
        // Permute rows and columns so that |a-b| = 1
        larg := Max({a,b});
        smal := Min({a,b});
        SwapColumns(~M,smal+1,larg);
        SwapRows(~M,smal+1,larg);
        // Clear this 2x2 subblock
        for m in [1..n] do
          if (m ne smal) and (m ne smal+1) then
            alpha := M[smal][smal];
            beta := M[smal][smal+1];
            gamma := M[smal+1][smal+1];
            a := M[smal][m];
            b := M[smal+1][m];
            c := M[m][m];
            // Find d and e so that
            // alpha * d + beta * e = -a
            // beta * d + gamma * e = -b
            d := (b*beta - gamma*a)/(alpha*gamma-beta^2);
            e := (a*beta - alpha*b)/(alpha*gamma-beta^2);
            // Make the substitution (x_smal,x_(smal+1),x_m) ->
            // (x_smal + d*x_m, x_(smal+1) + e*x_m, x_m)
            AddColumn(~M,d,smal,m);
            AddRow(~M,d,smal,m);
            AddColumn(~M,e,smal+1,m);
            AddRow(~M,e,smal+1,m);
          end if;
        end for;
        alpha := M[smal][smal];
        beta := M[smal][smal+1];
        gamma := M[smal+1][smal+1];
        if p gt 2 then
          //printf "Diagonalizing.\n";
          //printf "Starting with %o x^2 + %o xy + %o y^2.\n",
          //alpha,2*beta,gamma;
          if Valuation(beta,p) lt
                    Min(Valuation(alpha,p),Valuation(gamma,p)) then
            temp := MakeInt(alpha+2*beta+gamma,p);
            M[smal][smal] := temp;
            temp := beta + gamma;
            M[smal+1][smal] := temp;
            M[smal][smal+1] := temp;
            alpha := M[smal][smal];
            beta := M[smal+1][smal];
            gamma := M[smal+1][smal+1];
            // Now, ord_p(alpha) = ord_p(beta)
          end if;
          // Diagonalize
          if Valuation(beta,p) ge Valuation(alpha,p) then
            c := -beta/alpha;
            temp := MakeInt(alpha*c^2 + 2*beta*c + gamma,p);
            M[smal+1][smal+1] := temp;
            M[smal][smal+1] := 0;
            M[smal+1][smal] := 0;
            alpha := M[smal][smal];
            beta := M[smal+1][smal];
            gamma := M[smal+1][smal+1];
            else
            c := -beta/gamma;
            M[smal][smal+1] := 0;
            M[smal+1][smal] := 0;
            temp := MakeInt(a + 2*beta*c + gamma*c^2,p);
            M[smal][smal] := temp;
            alpha := M[smal][smal];
            beta := M[smal+1][smal];
            gamma := M[smal+1][smal+1];
          end if;
          //printf "Finally %o x^2 + %o xy + %o y^2.\n",alpha,2*beta,gamma;
        end if;
        if p eq 2 then
          // Use the theorem of BW Jones. Two binary forms
          // ax^2 + bxy + cy^2 with b odd are equivalent over Z_2
          // if and only if the two expressions 4ac - b^2 are congruent
          // modulo 8.
          v := Valuation(beta,p)+1;
          reda := alpha/(2^(v));
          redb := beta/(2^(v));
          redc := gamma/(2^(v));
          disc := Integers()!MakeInt((4*reda*redc-4*redb^2),2);
          if (disc mod 8 eq 3) then
            M[smal][smal] := 2^(v);
            M[smal][smal+1] := 2^(v-1);
            M[smal+1][smal] := 2^(v-1);
            M[smal+1][smal+1] := 2^(v);
          else
            M[smal][smal] := 0;
            M[smal][smal+1] := 2^(v-1);
            M[smal+1][smal] := 2^(v-1);
            M[smal+1][smal+1] := 0;
          end if;
        end if;
      end if;
    end if;
    //printf "Current matrix is %o.\n",M;
  end while;
  // Clear all denominators that are coprime to p
  for i in [1..n] do
    a := M[i][i];
    if (a ne 0) then
      M[i][i] := MakeInt(a,p);
    end if;
  end for;
  ret2 := M;
  return ret2;
end function;

M := Matrix([[1,0,-1/2,0],[0,2,1/2,-1],[-1/2,1/2,5,1/2],[0,-1,1/2,29]]);