// redqf.m -- minimisation and reduction of quadratic forms
// M. Stoll, started 2002-05-06

declare verbose QFReduce, 3;

intrinsic pMinimise(Q::Mtrx, p::RngIntElt) -> Mtrx, Mtrx
{Minimises the square symmetric nonsingular matrix Q with integral entries
 at p, where Q is interpreted as a quadratic form.}
  require NumberOfRows(Q) eq NumberOfColumns(Q): "Q must be square.\n";
  require Q eq Transpose(Q): "Q must be symmetric.";
  Q := ChangeRing(Q, Rationals());
  require forall{c : c in Eltseq(Q) | Denominator(c) eq 1}:
          "Q must have integral entries.";
  require Determinant(Q) ne 0: "Q must be nonsingular.";
  n := NumberOfRows(Q);
  F := GF(p);
  V := KSpace(F, n);
  tr := IdentityMatrix(Rationals(), n);
  ind := function(b)
           i := 0;
           repeat i +:= 1; until b[i] ne 0;
           return i;
         end function; 
  repeat
    Ker := Kernel(ChangeRing(Q, F));
    dim := Dimension(Ker);
    if dim eq 0 then break; end if;
    inds := [ind(b) : b in Basis(Ker)];
    assert #inds eq #Seqset(inds);
    assert forall{i : i in [1..dim] | Basis(Ker)[i,inds[i]] eq 1};
    co := [V.i : i in [1..n] | i notin inds];
    mat := Matrix(Rationals(), [ChangeUniverse(Eltseq(b), Integers())
                                 : b in Basis(Ker) cat co]);
    assert Abs(Determinant(mat)) eq 1;
    tr := mat*tr;
    Q := mat*Q*Transpose(mat);
    Qsub := ChangeRing(Submatrix(Q, 1, 1, dim, dim)/p, F);
    Ker := Kernel(Qsub);
    dims := Dimension(Ker);
    if dims eq 0 then
      // find isotropic subspace
      Ker := IsotropicSubspace(Qsub);
      dims := Dimension(Ker);
      if dims eq 0 then break; end if;
    end if;
    V1 := KSpace(F, dim);
    inds := [ind(b) : b in Basis(Ker)];
    assert #inds eq #Seqset(inds);
    assert forall{i : i in [1..dims] | Basis(Ker)[i,inds[i]] eq 1};
    co := [V1.i : i in [1..dim] | i notin inds];
    mat := Matrix(Rationals(), [ChangeUniverse(Eltseq(b), Integers())
                                 cat [0 : j in [dim+1..n]]
                                : b in Basis(Ker) cat co]
                                cat [[i eq j select 1 else 0 : i in [1..n]]
                                       : j in [dim+1..n]]);
    assert Abs(Determinant(mat)) eq 1;
    tr := mat*tr;
    Q := mat*Q*Transpose(mat);
    diag := DiagonalMatrix([1/p : i in [1..dims]]
                             cat [1 : i in [dims+1..n]]);
    tr := diag*tr;
    Q := diag*Q*diag;
  until Valuation(Determinant(Q), p) le 1;
  return Q, tr;
end intrinsic;

intrinsic pMinimise(Q::Mtrx, p::RngOrdIdl) -> Mtrx, Mtrx
{Minimises the square symmetric nonsingular matrix Q with integral entries
 at p, where Q is interpreted as a quadratic form.}
  require NumberOfRows(Q) eq NumberOfColumns(Q): "Q must be square.\n";
  require Q eq Transpose(Q): "Q must be symmetric.";
  O := Order(p);
  K := FieldOfFractions(O);
  Q := ChangeRing(Q, K);
  require forall{c : c in Eltseq(Q) | IsIntegral(c)}:
          "Q must have integral entries.";
  require Determinant(Q) ne 0: "Q must be nonsingular.";
  flag, gen := IsPrincipal(p);
  require flag: "p must be a principal ideal.";
  n := NumberOfRows(Q);
  F, m := ResidueClassField(p);
  V := KSpace(F, n);
  tr := IdentityMatrix(K, n);
  ind := function(b)
           i := 0;
           repeat i +:= 1; until b[i] ne 0;
           return i;
         end function; 
  repeat
    Ker := Kernel(Matrix(F, n, [m(c) : c in Eltseq(Q)]));
    dim := Dimension(Ker);
    if dim eq 0 then break; end if;
    inds := [ind(b) : b in Basis(Ker)];
    assert #inds eq #Seqset(inds);
    assert forall{i : i in [1..dim] | Basis(Ker)[i,inds[i]] eq 1};
    co := [V.i : i in [1..n] | i notin inds];
    mat := Matrix(K, [[c @@ m : c in Eltseq(b)] : b in Basis(Ker) cat co]);
    assert IsUnit(O!Determinant(mat));
    tr := mat*tr;
    Q := mat*Q*Transpose(mat);
    Qsub := Matrix(F, dim, 
                      [m(c/gen) : c in Eltseq(Submatrix(Q, 1, 1, dim, dim))]);
    Ker := Kernel(Qsub);
    dims := Dimension(Ker);
    if dims eq 0 then
      // find isotropic subspace
      Ker := IsotropicSubspace(Qsub);
      dims := Dimension(Ker);
      if dims eq 0 then break; end if;
    end if;
    V1 := KSpace(F, dim);
    inds := [ind(b) : b in Basis(Ker)];
    assert #inds eq #Seqset(inds);
    assert forall{i : i in [1..dims] | Basis(Ker)[i,inds[i]] eq 1};
    co := [V1.i : i in [1..dim] | i notin inds];
    mat := Matrix(K, [[c @@ m : c in Eltseq(b)] cat [0 : j in [dim+1..n]]
                        : b in Basis(Ker) cat co]
                      cat [[i eq j select 1 else 0 : i in [1..n]]
                            : j in [dim+1..n]]);
    assert IsUnit(O!Determinant(mat));
    tr := mat*tr;
    Q := mat*Q*Transpose(mat);
    diag := DiagonalMatrix([1/gen : i in [1..dims]]
                             cat [1 : i in [dims+1..n]]);
    tr := diag*tr;
    Q := diag*Q*diag;
  until Valuation(Determinant(Q), p) le 1;
  return Q, tr;
end intrinsic;

intrinsic IsotropicSubspace(Q::Mtrx) -> ModTupFld
{Finds a maximal isotropic subspace for Q over a finite field.}
  require NumberOfRows(Q) eq NumberOfColumns(Q): "Q must be square.\n";
  require Q eq Transpose(Q): "Q must be symmetric.";
  require Determinant(Q) ne 0: "Q must be nonsingular.";
  F := BaseRing(Q);
  require Type(F) eq FldFin: "Base ring must be a finite field.";
  n := NumberOfRows(Q);
  V := KSpace(F, n);
  if n le 1 then return sub<V | >; end if;
  Q0 := Q;
  tr := IdentityMatrix(F, n);
  // first find one isotropic element
  if Characteristic(F) eq 2 then
    V1 := KSpace(F, 1);
    Ker := Kernel(hom<V -> V1 | [V1 | [Q[i,i]] : i in [1..n]]>);
    if Dimension(Ker) eq 0 then return sub<V | >; end if;
    v := Basis(Ker)[1];
  elif Q[1,1] eq 0 then
    // take v = (1 0 0 ...)
    v := V.1;
  else
    mat := Matrix([[1] cat [-Q[1,i]/Q[1,1] : i in [2..n]]]
                   cat [[j eq i select 1 else 0 : j in [1..n]] : i in [2..n]]);
    tr := Transpose(mat)*tr;
    Q := Transpose(mat)*Q*mat;
    if Q[2,2] eq 0 then
      v := V.2;
    else
      mat := Matrix(F, [[j eq 1 select 1 else 0 : j in [1..n]]]
                        cat [[0, 1] cat [-Q[2,i]/Q[2,2] : i in [3..n]]]
                        cat [[j eq i select 1 else 0 : j in [1..n]]
                              : i in [3..n]]);
      tr := Transpose(mat)*tr;
      Q := Transpose(mat)*Q*mat;
      flag, sqrt := IsSquare(-Q[2,2]/Q[1,1]);
      if flag then
        v := sqrt*V.1 + V.2;
      elif n eq 2 then // form is anisotropic
        return sub<V | >;
      elif Q[3,3] eq 0 then
        v := V.3;
      else
        // one step further
        mat := Matrix([[j eq i select 1 else 0 : j in [1..n]] : i in [1..2]]
                       cat [[0, 0, 1] cat [-Q[3,i]/Q[3,3] : i in [4..n]]]
                       cat [[j eq i select 1 else 0 : j in [1..n]]
                              : i in [4..n]]);
        tr := Transpose(mat)*tr;
        Q := Transpose(mat)*Q*mat;
        repeat
          y := Random(F);
          flag, sqrt := IsSquare((-Q[3,3]-Q[2,2]*y^2)/Q[1,1]);
        until flag;
        v := sqrt*V.1 + y*V.2 + V.3;
      end if;
    end if;
  end if;
  sp := func<v1, v2, Q | (v1*Q*Matrix([[v2[i]] : i in [1..n]]))[1]>;
  // now v is an isotropic vector
  assert sp(v, v, Q) eq 0;
  assert sp(v0, v0, Q0) eq 0 where v0 := v*tr;
  // find orthogonal complement
  Orth := Kernel(Q*Matrix([[v[i]] : i in [1..n]]));
  comp := Basis(qu)[1] @@ pr where qu, pr := quo<V | Orth>;
  // cc := sp(comp, comp, Q);
  // cv := sp(comp, v, Q);
  // comp -:= cc/(2*cv)*v;
  Sub1 := sub<V | v, comp>; // a hyperbolic plane
  Orth := Kernel(Q*Matrix([[v[i], comp[i]] : i in [1..n]]));
  assert forall{0 : b in Basis(Orth) | sp(v, b, Q) eq 0};
  mat := Matrix([Eltseq(v), Eltseq(comp)]
                 cat [Eltseq(b) : b in Basis(Orth)]);
  tr := mat*tr;
  Q := mat*Q*Transpose(mat);
  assert forall{0 : i in [3..n] | sp(V.1, V.i, Q) eq 0};
  // recurse
  Iso1 := IsotropicSubspace(Submatrix(Q, 3, 3, n-2, n-2));
  // put together
  IsoBas := [V.1] cat [V | [0,0] cat Eltseq(b) : b in Basis(Iso1) ];
  assert forall{0 : bi in IsoBas | sp(bi, V.1, Q) eq 0};
  // undo transformations
  Iso := sub<V | [b*tr : b in IsoBas]>;
  assert forall{0 : bi, bj in Basis(Iso) | sp(bi, bj, Q0) eq 0};
  return Iso;
end intrinsic;

function red1(Q)
  n := NumberOfRows(Q);
  tr := IdentityMatrix(Rationals(), n);
  for i := 1 to n do
    if Q[i,i] ne 0 then
      coli := [j eq i select 1 else -Round(Q[i,j]/Q[i,i]) : j in [1..n]];
      mat := Matrix(Rationals(),
                    [j eq i select coli 
                            else [k eq j select 1 else 0 : k in [1..n]]
                      : j in [1..n]]);
      tr := Transpose(mat)*tr;
      Q := Transpose(mat)*Q*mat;
    end if;
  end for;
  return Q, tr;
end function;

function em(i,j,k,N) // elementary matrix
  return IdentityMatrix(Rationals(),N)
          + Matrix(N, [(ii eq i and jj eq j) select k else 0
                          : ii,jj in [1..N]]);
end function;

function red2(Q)
  n := NumberOfRows(Q);
  tr := IdentityMatrix(Rationals(), n);
  inds0 := [<i, Abs(Q[i,i])> : i in [1..n] | Abs(Q[i,i]) ne 1];
  inds := [a : a in inds0 | a[2] ne 0];
  if #inds eq 0 then return Q, tr; end if;
  // sort ascendingly
  Sort(~inds, func<a,b | b[2]-a[2]>);
  for a in inds do
    i := a[1];
    discs := [<j, -d> : b in inds0 | d lt 0 where d := Q[i,i]*Q[j,j]-Q[i,j]^2
                                            where j := b[1]];
    if #discs gt 0 then
      _, pos := Min([b[2] : b in discs]);
      j := discs[pos, 1];
      if Q[j,j] eq 0 then
        l := Round(-Q[i,i]/(2*Q[i,j]));
      else
        sd := Sqrt(discs[pos, 2]);
        ls := [Round((-Q[i,j] + sd)/Q[j,j]), Round((-Q[i,j]-sd)/Q[j,j])];
        as := [Q[i,i] + 2*l*Q[i,j] + l^2*Q[j,j] : l in ls];
        l := as[1] lt as[2] select ls[1] else ls[2];
      end if;
      tr := em(i,j,l,n);
      Q := tr*Q*Transpose(tr);
      break;
    end if;
  end for;
  return Q, tr;
end function;

intrinsic Reduce(Q::Mtrx) -> Mtrx, Mtrx
{Reduce size of matrix entries of Q using unimodular transformations.}
  require NumberOfRows(Q) eq NumberOfColumns(Q): "Q must be square.\n";
  require Q eq Transpose(Q): "Q must be symmetric.";
  Q := ChangeRing(Q, Rationals());
  require forall{c : c in Eltseq(Q) | Denominator(c) eq 1}:
          "Q must have integral entries.";
  require Determinant(Q) ne 0: "Q must be nonsingular.";
  N := NumberOfRows(Q);
  tr1 := IdentityMatrix(Rationals(), N);
  // set up a list of (monoid) generators
  gs1 := [em(i,j,1,N) : i,j in [1..N] | i ne j];
  gs2 := [em(i,j,-1,N) : i,j in [1..N] | i ne j];
  gens1 := gs1 cat gs2;
  size := func<Q | &+[Abs(c) : c in Eltseq(Q)]>;
  s1 := size(Q);
  Q1 := Q;
  repeat
    s0 := s1;
    // first use red1
    repeat
      s := s1;
      Q := Q1;
      tr := tr1;
      Q1, mat := red1(Q);
      tr1 := mat*tr;
      for k := 1 to 3 do
        Q1, mat := red1(Q1);
        tr1 := mat*tr1;
      end for;
      s1 := size(Q1);
    until s1 ge s;
    // then use red2
    s1 := s;
    Q1 := Q;
    tr1 := tr;
    repeat
      s := s1;
      Q := Q1;
      tr := tr1;
      Q1, mat := red2(Q);
      tr1 := mat*tr;
      s1 := size(Q1);
    until s1 ge s;
    // then do a step with the generators
    Qs := [g*Q*Transpose(g) : g in gens1];
    s1, pos := Min([size(Q1) : Q1 in Qs]);
    Q1 := Qs[pos];
    tr1 := gens1[pos]*tr;
  until s1 ge s0;
  return Q, tr;
end intrinsic;

intrinsic GramSchmidt(seq::SeqEnum) -> SeqEnum
{Perform Gram-Schmidt orthonormalization.}
  for i := 1 to #seq do
    s := Sqrt((seq[i], seq[i]));
    seq[i] *:= 1/s;
    for j := i+1 to #seq do
      seq[j] := seq[j] - (seq[i], seq[j])*seq[i];
    end for;
  end for;
  return seq;
end intrinsic;

/*
function MyConjugates(z, prec)
  cs := Conjugates(z);
  rs := [r[1] : r in Roots(MinimalPolynomial(z), ComplexField())];
  return [rs[pos] where _, pos := Min([Norm(c-r) : r in rs]) : c in cs];
end function;
*/

intrinsic Upper(M::Mtrx) -> Mtrx
{Given a symmetric nondegenerate matrix, this returns a positive definite
 matrix over R bounding the quadratic form.}
  n := NumberOfRows(M);
  require n eq NumberOfColumns(M): "M must be a square matrix.";
  require M eq Transpose(M): "M must be symmetric.";
  prec := Ceiling(Log(10, &+[c^2 : c in Eltseq(M)])) + 20;
  pol := CharacteristicPolynomial(M);
  fact := Factorisation(pol);
  evs := [];
  lambdas := [];
  for f in fact do
    polf := f[1];
    mult := f[2];
    if Degree(polf) eq 1 then
      a := -Coefficient(polf, 0);
      E := Eigenspace(ChangeRing(M, RealField()), a);
      evs cat:= GramSchmidt(Basis(E));
      lambdas cat:= [a : j in [1..Dimension(E)]];
    else
      K := NumberField(polf);
      SetKantPrecision(K, prec);
      EK := Eigenspace(ChangeRing(M, K), K.1);
      BEK := Basis(EK);
      d := Dimension(EK);
      ls := ChangeUniverse(Conjugates(K.1), RealField());
      conjs := [ChangeUniverse(Conjugates(b[i]), RealField())
                 : i in [1..n], b in BEK];
      for k := 1 to Degree(polf) do
        vecs := [Vector([conjs[(j-1)*n+i, k] : i in [1..n]]) : j in [1..d]];
        evs cat:= GramSchmidt(vecs);
        lambdas cat:= [ls[k] : j in [1..d]];
      end for;
    end if;
  end for;
  Emat := Matrix(evs);
  D := DiagonalMatrix([Abs(l) : l in lambdas]);
  return Transpose(Emat)*D*Emat;
end intrinsic;

function Red2(Fseq, i, j)
  PP := Universe(Fseq);
  n := Rank(PP);
  P := PolynomialRing(CoefficientRing(PP));
  subst := [k eq i select P.1 else k eq j select 1 else 0 : k in [1..n]];
  F1 := &+[Evaluate(F, subst)^2 : F in Fseq];
  tr := IdentityMatrix(Integers(), n);
  deg := 2*TotalDegree(Fseq[1]);
  if F1 eq 0 or Degree(F1) lt deg-1 or Discriminant(F1) eq 0 then
    return Fseq, tr;
  end if;
  F1, mat := Reduce(F1, deg);
  tr[i,i] := mat[1,1];
  tr[i,j] := mat[1,2];
  tr[j,i] := mat[2,1];
  tr[j,j] := mat[2,2];
  Fseq := [Evaluate(F, [k eq i select mat[1,1]*PP.i + mat[1,2]*PP.j else
                        k eq j select mat[2,1]*PP.i + mat[2,2]*PP.j
                               else PP.k
                         : k in [1..n]])
            : F in Fseq];
  return Fseq, tr;
end function;

intrinsic RedForms(Fseq::[RngMPolElt]) -> SeqEnum
{}
  vprintf QFReduce, 1: "RedForms ...\n";
  s1 := &+[&+[c^2 : c in Coefficients(F)] : F in Fseq];
  PP := Universe(Fseq);
  n := Rank(PP);
  deg := TotalDegree(Fseq[1]);
  mons := MonomialsOfDegree(PP, deg);
  function LLLR(Fseq)
    fmat := LLL(Matrix([[MonomialCoefficient(f, m) : m in mons] : f in Fseq]));
    return [&+[fmat[i,j]*mons[j] : j in [1..#mons]] : i in [1..#Fseq]];
  end function;
  Fseq1 := Fseq;
  tr1 := IdentityMatrix(Integers(), n);
  repeat
    Fseq := Fseq1;
    s := s1;
    tr := tr1;
    Fs := [<f, t> where f, t := Red2(Fseq, i, j) : j in [i+1..n], i in [1..n]];
    Fs cat:= [<LLLR(a[1]), a[2]> : a in Fs];
    s1, pos := Min([&+[&+[c^2 : c in Coefficients(p)] : p in F[1]] : F in Fs]);
    Fseq1 := Fs[pos, 1];
    tr1 := tr*Fs[pos, 2];
    vprintf QFReduce, 3: "    size = %o\n", s1;
  until s1 ge s;
  return Fseq, tr;
end intrinsic;

function Red2F(F, i, j)
  PP := Parent(F);
  n := Rank(PP);
  P := PolynomialRing(CoefficientRing(PP));
  subst := [k eq i select P.1 else k eq j select 1 else 0 : k in [1..n]];
  F1 := Evaluate(F, subst);
  tr := IdentityMatrix(Integers(), n);
  deg := TotalDegree(F);
  if F1 eq 0 or Degree(F1) lt deg-1 or Discriminant(F1) eq 0 then
    return F, tr;
  end if;
  F1, mat := Reduce(F1, deg);
  tr[i,i] := mat[1,1];
  tr[i,j] := mat[1,2];
  tr[j,i] := mat[2,1];
  tr[j,j] := mat[2,2];
  F := Evaluate(F, [k eq i select mat[1,1]*PP.i + mat[1,2]*PP.j else
                    k eq j select mat[2,1]*PP.i + mat[2,2]*PP.j
                           else PP.k
                     : k in [1..n]]);
  return F, tr;
end function;

function Red3F(F, i, j, k, s)
  PP := Parent(F);
  n := Rank(PP);
  P3 := PolynomialRing(CoefficientRing(PP), 3);
  subst := [l eq i select P3.1 else 
            l eq j select P3.2 else 
            l eq k select P3.3 else 0 : l in [1..n]];
  F1 := Evaluate(F, subst);
  tr := IdentityMatrix(Integers(), n);
  if F1 eq 0 then
    return F, tr;
  end if;
  _, _, D := T3Invariants(F1);
  if D eq 0 then
    return F, tr;
  end if;
  F1, mat := Reduce(F1 : Steps := 1, AdHoc := s ge 10^20);
  inds := [i,j,k];
  for r,s in [1..3] do tr[inds[r], inds[s]] := mat[r, s]; end for;
  F := Evaluate(F, [l eq i select mat[1,1]*PP.i + mat[1,2]*PP.j + mat[1,3]*PP.k
               else l eq j select mat[2,1]*PP.i + mat[2,2]*PP.j + mat[2,3]*PP.k
               else l eq k select mat[3,1]*PP.i + mat[3,2]*PP.j + mat[3,3]*PP.k
                           else PP.l
                     : l in [1..n]]);
  return F, tr;
end function;

function RedFormHelp(F)
  s1 := &+[c^2 : c in Coefficients(F)];
  PP := Parent(F);
  n := Rank(PP);
  deg := TotalDegree(F);
  F1 := F;
  tr1 := IdentityMatrix(Integers(), n);
  repeat
    F := F1;
    s := s1;
    tr := tr1;
    Fs := [<f, t> where f, t := Red2F(F, i, j) : j in [i+1..n], i in [1..n]];
    s1, pos := Min([&+[c^2 : c in Coefficients(p[1])] : p in Fs]);
    F1 := Fs[pos, 1];
    tr1 := tr*Fs[pos, 2];
    vprintf QFReduce, 3: "    size = %o\n", s1;
  until s1 ge s;
  return F, tr;
end function;

intrinsic RedForm(F::RngMPolElt) -> RngMPolElt, Mtrx
{}
  vprintf QFReduce, 1: "RedForm ...\n";
  s1 := &+[c^2 : c in Coefficients(F)];
  PP := Parent(F);
  n := Rank(PP);
  deg := TotalDegree(F);
  if deg ge 4 then return RedFormHelp(F); end if;
  F1 := F;
  tr1 := IdentityMatrix(Integers(), n);
  repeat
    F, T := RedFormHelp(F1);
    s := s1;
    tr := tr1*T;
    vprintf QFReduce, 3: "    doing ternary reduction...\n";
    Fs := [<f, t> where f, t := Red3F(F, i, j, k, s)
            : k in [j+1..n], j in [i+1..n], i in [1..n]];
    s1, pos := Min([&+[c^2 : c in Coefficients(p[1])] : p in Fs]);
    F1 := Fs[pos, 1];
    tr1 := tr*Fs[pos, 2];
    vprintf QFReduce, 3: "    size = %o\n", s1;
  until s1 ge s;
  return F, tr;
end intrinsic;