/************************************************************************\
*  Reduction of Genus One Models over Q (of degrees n = 2,3,4 and 5)     *
*  Author: T. Fisher                                                     *
*  Date: 22nd July 2008                                                  *
*  Version for n = 2 calls "QuarticReduce"                               *
*  Versions for n = 3,4 based on joint work with J. Cremona and M. Stoll *
\************************************************************************/

declare verbose Reduction,4;

// Find pathname by typing e.g. "Hessian:Maximal;" at the magma prompt
import "/usr/magma/package/Geometry/CrvG1/g1invariants.m" :
  SL4Covariants,MultiFact,AuxiliaryQuadrics,ExactPfaffianMatrix;
import "/usr/magma/package/Geometry/CrvG1/3descent-testeq.m" :
  RedTorsPts,ReductionInnerProduct;

MC := MonomialCoefficient;
MD := MonomialsOfDegree;
Id := IdentityMatrix;
Deriv := Derivative;
MySign := func<x|x ge 0 select 1 else -1>;
dyadic_rep := func<x|x eq 0 select 0 
  else a*2^b where a,b is MantissaExponent(x)>;
L2Norm := func<model|&+[x^2: x in Eltseq(model)]>;

function FirstSign(seq)
  xx := [x : x in seq | x ne 0];
  fnz := #xx eq 0 select 1 else xx[1];
  return MySign(fnz);
end function;

function IntegerMatrix(mat)
  n := Nrows(mat);
  mat1 := Matrix(n,n,[dyadic_rep(mat[i,j]):i,j in [1..n]]);
  denom := LCM([Integers()|Denominator(x) : x in Eltseq(mat1)]);
  return MatrixAlgebra(Integers(),n)!(denom*mat1);
end function;  

function MyEchelonForm(mat,tol)
// Row reduce an m by n matrix over the reals or complexes,
// always using largest available entry as the next pivot.
// Matrix returned is of form ( Id(m) | B )  
  m := Nrows(mat);
  n := Ncols(mat);
  for j in [1..m] do
    _,jj := Maximum([Abs(mat[k,j]): k in [j..m]]);
    jj +:= (j-1);
    temp := mat[j]; mat[j] := mat[jj]; mat[jj] := temp;
    if Abs(mat[j,j]) lt tol then return 0; end if;
    mat[j] *:= (1/mat[j,j]);
    for k in [1..m] do
      if k ne j then mat[k] -:= mat[k,j]*mat[j]; end if;
    end for;
  end for;
  return mat;
end function;

function SizeString(model)
  size := Round(L2Norm(model));
  logsize := Ilog(10,size);
  if logsize lt 10 then 
    return "L2Norm = " cat IntegerToString(size);
  else 
    return "Log10(L2Norm) = " cat IntegerToString(logsize);
  end if;
end function;

function CoefficientSize(model:name:="Model");
  cs := Ilog(10,Max([Round(Abs(c)): c in Eltseq(model)])); 
  if cs lt 10 then 
    vprintf Reduction, 2: "%o has coefficients\n%o\n",name,Eltseq(model);
  else
    vprintf Reduction, 2: "%o has coefficients of length %o\n",name,cs;
  end if;
  return cs;
end function;

//////////////////////////////////////////////////////////////
//                      Case n = 2                          //
//////////////////////////////////////////////////////////////

function ReduceTwo(model)
  assert Degree(model) eq 2;
  P<X> := PolynomialRing(Integers());
  CrossTerms := (#Eltseq(model) eq 8);
  model1 := CrossTerms select CompleteTheSquare(model) else model;
  quartic := Evaluate(Equation(model1),[X,1]);
  _,M := QuarticReduce(quartic); // from CrvEll/FourDesc/quartred.m
  M := ChangeRing(M,Rationals());
  tr := <1,Transpose(M)>;
  model := tr*model;
  if CrossTerms then 
    coeffs := Eltseq(model);
    r := [ -(Integers()!(coeffs[i]) div 2) : i in [1..3]];
    tr1 := <1,r,Id(Rationals(),2)>;
    model := tr1*model;
    tr := tr1*tr;
  end if;
  return model,tr;
end function;

//////////////////////////////////////////////////////////////
//                      Case n = 4                          //
//////////////////////////////////////////////////////////////

SRCompare := func<b,a| Abs(a[2]) ne Abs(b[2]) 
    select Abs(b[2])-Abs(a[2]) else Abs(b[3])-Abs(a[3]) >;

function WomackTransformation(model)
// Returns a pair of "signed permutations" in GL(2,Z) and GL(4,Z)
// (just to make the final answer look pretty)
  Q := Rationals();
  assert Degree(model) eq 4;
  assert CoefficientRing(model) eq Q;
  R := PolynomialRing(model);
  q1,q2 := Explode(Equations(model));
  d0 := [<i,MC(q1,R.i^2),MC(q2,R.i^2)> : i in [1..4]]; 
  Sort(~d0, SRCompare);
  N1 := Matrix(Q,4,4,[<i,d0[i][1],1>: i in [1..4]]);
  model := <Id(Q,2),N1>*model;
  qq := Equations(model);
  signs := [Q|FirstSign([MC(q,R.i^2): i in [1..4]]): q in qq];
  M := DiagonalMatrix(signs);
  model := <M,Id(Q,4)>*model;
  q1,q2 := Explode(Equations(model));
  s1 := [Sign(MC(q1,R.1*R.j)): j in [1..4]];
  s2 := [MySign(MC(q2,R.1*R.j)): j in [1..4]];
  s := [Q|s1[i] ne 0 select s1[i] else s2[i]: i in [1..4]];
  N2 := DiagonalMatrix(s); 
  return <M,N2*N1>;
end function;

function ReductionCovariantFour(model,prec)
  Q := Rationals();
  R := RealField(prec);
  C := ComplexField(prec);
  P<X> := PolynomialRing(Q);
  a,bb,c,dd,e := Explode(SL4Invariants(model));
  covars := SL4Covariants(model); // from g1invariants.m
  g := a*X^4 + 2*bb*X^3 + c*X^2 + 2*dd*X + e;
  ff := [f[1]: f in Factorization(g)| f[1] ne X];
  quadrics := <>;
  if a eq 0 then 
    qA,T1,T2,qB := Explode(covars);
    Append(~quadrics,T2 - 2*bb*qB);
  end if;    
  if e eq 0 then 
    qA,T1,T2,qB := Explode(covars);
    Append(~quadrics,T1 - 2*dd*qA);
  end if;    
  for f in ff do
    if Degree(f) gt 1 then
      F<theta> := NumberField(f);
// Since we use "Conjugates" below, we need to set the
// precision of F. This seems quite slow for large precisons, so it 
// might be worth re-writing to avoid the Kant machinery.
      SetKantPrecision(F,prec);
    else
      F := Q; 
      theta := Roots(f)[1][1];
    end if;   
    RR := PolynomialRing(F,4);
    covarsF := ChangeUniverse(covars,RR);
    qA,T1,T2,qB := Explode(covarsF);
    GG := theta^(-1)*e*qA + T1 + theta*T2 + theta^2*a*qB;
    Append(~quadrics,GG);
  end for;
  mats := <>;
  diagelts := [];
  for GG in quadrics do
    RR := Parent(GG);
    F := BaseRing(RR);
    G := Factorization(GG)[1][1];
    alpha := F!ExactQuotient(GG,G^2);
    assert alpha ne 0;
    assert GG eq alpha*G^2;
    M1 := Matrix(C,4,Degree(F),[Conjugates(MC(G,RR.i)): i in [1..4]]);
    Append(~mats,M1); 
    aa := [C|Abs(x): x in Conjugates(alpha)];
    dg := [C|Abs(x): x in Conjugates(Evaluate(Derivative(g),theta))];
    ll := [aa[i]/dg[i]^(1/2) : i in [1..Degree(F)]];
    diagelts cat:= ll;
  end for;
  M1 := HorizontalJoin(mats);
  M2 := Matrix(C,4,4,[Conjugate(M1[j,i]):i,j in [1..4]]);
  D := DiagonalMatrix(C,diagelts);
  M := M1*D*M2;
  Qmat := Matrix(R,4,4,[Re(M[i,j]): i,j in [1..4]]);
  dd := &+[Abs(x)^2: x in Eltseq(Qmat)];
  Qmat *:= dd^(-1/2);
  return Qmat + Transpose(Qmat); 
end function;

//////////////////////////////////////////////////////////////
//                      Case n = 5                          //
//////////////////////////////////////////////////////////////

function InvariantsAndHessian(model)
// Computes the invariants c4,c6 and Hessian H of a genus one model 
// of degree 5. (Calling the existing functions cInvariants and 
// Hessian separately leads to some duplication of effort.)
  assert Degree(model) eq 5;
  assert BaseRing(model) cmpeq Rationals();
  phi := Matrix(model);
  p := Equations(model);
  q := AuxiliaryQuadrics(p); // from g1invariants.m
  R := BaseRing(phi);
  RR<u> := PolynomialRing(R);
  SS<t> := quo<RR|u^4>; // this trick would speed up cInvariants
  U := Matrix(SS,5,5,[&+[MC(phi[j,k],R.i)*R.k: k in [1..5]]: i,j in [1..5]]);
  P := Matrix(R,5,5,
         [&+[MC(Deriv(p[k],j),R.i)*R.k: k in [1..5]]: i,j in [1..5]]);
  Q := Matrix(SS,5,5,[Deriv(q[i],j):i,j in [1..5]]);
  det1 := Determinant(P);
  det2 := Determinant(U+t*Q);
  T := [R|Coefficient(det2,i): i in [1,3]];
  cc := [&+[MultiFact(mon)*MC(det1,mon)*MC(d2,mon):mon in MD(R,5)]:d2 in T];
  c4 := cc[1]/40;
  c6 := -cc[2]/320;
  RR<y1,y2,y3,y4,y5> := PolynomialRing(R,5);
  M := Determinant(Matrix(RR,5,5,
         [&+[MC(Deriv(p[k],j),R.i)*RR.k : k in [1..5]]:i,j in [1..5]]));
  for ctr in [1..2] do
    M := &+[R.i*MC(q[i],R.j*R.k)*Deriv(Deriv(M,RR.j),RR.k)
	   : j,k in [1..5], i in [1..5] | j le k ];
  end for;
  q2m := [R|MC(M,RR.i): i in [1..5]];
  p22 := Eltseq(4*c4*Vector(p)-(3/16)*Vector(q2m));
  hessian := ExactPfaffianMatrix(p22,1,1); // from g1invariants.m
  coeffs := func<f|[MC(f,R.k): k in [1..5]]>;
  Umat := Matrix(10,5,[coeffs(phi[i,j]): i,j in [1..5] | i lt j]);
  Hmat := Matrix(10,5,[coeffs(hessian[i,j]): i,j in [1..5] | i lt j]);  
  Dmat := HorizontalJoin(Umat,Hmat);
  discr1 := Determinant(Dmat);
  discr2 := 12^2*(c4^3-c6^2);
  if discr1 eq -discr2 then hessian*:= -1; end if;
  return c4,c6,GenusOneModel(hessian);
end function;

function SyzygeticMatrix(model)
// This is a 3 by 5 matrix of quadratic forms. The first row contains
// the 4 by 4 Pfaffians of the model, and the last row contains the
// the 4 by 4 Pfaffians of the Hessian.
// The Hesse family has 12 singular fibres: the "syzygetic pentagons".
// These pentagons have 30 vertices in total. These make up the locus 
// where the syzygetic matrix has rank 1. (The total space of the Hesse 
// family is the locus where the matrix has rank 2.)
  assert Degree(model) eq 5;
  assert BaseRing(model) cmpeq Rationals();
  P<t> := PolynomialRing(Rationals());
  vprintf Reduction,2 : 
    "Computing the invariants, Hessian and syzygetic matrix\n";
  c4,c6,H := InvariantsAndHessian(model);
  vprintf Reduction,3 : "cInvariants = %o\n",[c4,c6]; 
  coeffsize := CoefficientSize(H:name:="Hessian");
  F := t*Vector(P,50,Eltseq(model))+Vector(P,50,Eltseq(H));
  FF := GenusOneModel(5,Eltseq(F));
  eqns := Equations(FF);
  R := PolynomialRing(model);
  RR := Universe(eqns);
  mons := MD(R,2);
  coeffs := [[MC(f,RR!mon): mon in mons]: f in eqns];
  p2 := [&+[Coefficient(c[j],2)*mons[j]: j in [1..#mons]]: c in coeffs];
  p12 := [&+[Coefficient(c[j],1)*mons[j]: j in [1..#mons]]: c in coeffs];
  p22 := [&+[Coefficient(c[j],0)*mons[j]: j in [1..#mons]]: c in coeffs];
  assert p2 eq Equations(model);
  assert p22 eq Equations(H);
  return Matrix(R,3,5,[p2,p12,p22]),c4,c6;
end function;

function ApplyToMat(tr,mat)
// Since the Hessian is a covariant, there is no need to to keep 
// recomputing the syzygetic matrix.
  R := BaseRing(mat);
  S := ChangeRing(Adjoint(tr[1]),R);
  T := Transpose(tr[2]);
  mat := Matrix(3,5,[mat[i,j]^T:j in [1..5],i in [1..3]]);
  return mat*S;
end function;

function SolveForPentagons(qq,C,tol)
  matR := 0;
  matC := 0;
  for quads in qq do 
    assert #quads eq 10;   // 10 homogeneous quadrics 
    RR := Universe(quads);
    assert Rank(RR) eq 5;  // in 5 variables 
    R := BaseRing(RR);     // over the reals
    mons := [RR.i*RR.j: i,j in [1..4]|i le j] cat [RR.i*RR.5: i in [1..5]];
    coeffmat := Matrix(10,15,[MC(q,mon): mon in mons,q in quads]);
    coeffmat := MyEchelonForm(coeffmat,tol);
    if coeffmat eq 0 then 
      vprintf Reduction,2 : "Pentagon not in general position\n";
      return 0,0; 
    end if;
    charmat := Matrix(4,5,[-coeffmat[i,j]: j in [11..15],i in [1..4]]);
    charmat := VerticalJoin(charmat,Vector(R,5,[1,0,0,0,0]));
    evals := Eigenvalues(charmat);
    s := func<n|n eq 1 select "ex" else "ices">;
    vprintf Reduction,2 : "Pentagon has %o real vert%o\n",#evals,s(#evals);
    realflag := true;
    if #evals ne 5 then
      charmat := ChangeRing(charmat,C);
      charpoly := CharacteristicPolynomial(charmat);
      evals := Roots(charpoly,C);
// Don't use Eigenvalues(charmat), since in "CovariantFromHyperplanes" 
// we assume the pairs of complex conjugate roots are adjacent.
      if #evals ne 5 then 
        vprintf Reduction,2 : "Unable to find all complex vertices\n";
        return 0,0; 
      end if;
      realflag := false;
    end if;
    F := realflag select R else C;
    evals := [F| x[1]: x in evals];
    mm := Minimum([Abs(evals[i]-evals[j]):i,j in [1..5]|i lt j]);
    if mm lt tol then 
      vprintf Reduction,2 : "Unable to distinguish vertices\n";
      return 0,0;
    end if;
    VdM := Matrix(F,5,5,[x^i : x in evals,i in [0..4]]);
    v := Vector(F,5,[0,0,0,0,1]);
    vv := [];
    for i in [1..5] do vv cat:= [v]; v := v*charmat; end for;
    augmat := HorizontalJoin(VdM,VerticalJoin(vv));
    augmat := MyEchelonForm(augmat,tol);
    if augmat eq 0 then return 0,0; end if;
    mat := Matrix(5,5,[augmat[i,j]: j in [6..10],i in [1..5]]);
    dd := [1/Sqrt(&+[Abs(x)^2: x in Eltseq(mat[i])]): i in [1..5]];
    mat := DiagonalMatrix(F,dd)*mat;
    if realflag then matR := mat; else matC := mat; end if;
  end for;
  return matR,matC;
end function;

function GramSchmidt(mat)
  d := Ncols(mat);
  for i := 1 to Nrows(mat) do
    s := [&+[mat[j,k]*mat[i,k]: k in [1..d]]: j in [1..i-1]];
    if i gt 1 then
      mat[i] -:= &+[s[j]*mat[j]: j in [1..i-1]];
    end if;
    n := Sqrt(&+[mat[i,j]^2 : j in [1..d]]);
    mat[i] *:= (1/n);
  end for;
  return mat;
end function;

function MyProjection(mat,vec)
  m := Nrows(mat);
  n := Ncols(mat);
  a := [&+[mat[j,k]*vec[k]: k in [1..n]]: j in [1..m]];
  vec -:= &+[a[j]*mat[j]: j in [1..m]];
  n := Sqrt(&+[vec[i]^2 : i in [1..n]]);
  return (1/n)*vec;
end function;

function RotationMatrix(theta)
  c := Cos(theta); 
  s := Sin(theta);
  return Matrix(2,2,[c,s,-s,c]);
end function;

function ClosestPointOnCircle(A,B)
// Given matrices A and B whose rows are orthonormal bases
// for subspaces V and W of Euclidean space R^d, finds 
// cos(phi)*B[1] + sin(phi)*B[2] such that distance to V
// is minimal.
  assert Ncols(A) eq Ncols(B);
  assert Nrows(B) eq 2;
  n := Nrows(A);
  M := B*Transpose(A); 
  x := &+[M[1,i]^2-M[2,i]^2: i in [1..n]];
  y := &+[2*M[1,i]*M[2,i]: i in [1..n]];
  phi := Arctan2(x,y)/2; // N.B. don't reverse the arguments!
  R := RotationMatrix(phi);
  S := R*M;  
  dd := [1-&+[S[i,j]^2: j in [1..n]]: i in [1..2]];
  dist,m := Minimum(dd);
  return R[m]*B,dist;
end function;

function CommonVector(A,B)
// Computes a basis for the intersection of the row spaces 
// of A and B, where the intersection is assumed 1-dimensional 
  assert Nrows(A) eq 5;
  assert Nrows(B) eq 3;
  A := GramSchmidt(A);
  B := GramSchmidt(B);
  M := A*Transpose(B);
  dd := [1-&+[M[i,j]^2: i in [1..5]]: j in [1..3]];
  _,m := Maximum(dd);
  AA := VerticalJoin(A,MyProjection(A,B[m]));
  BB := Matrix([Eltseq(B[i]): i in [1..3]| i ne m]);
  vec,d1 := ClosestPointOnCircle(AA,BB);
  BB := VerticalJoin(vec,B[m]);
  vec,d2 := ClosestPointOnCircle(A,BB); 
  return vec,Maximum(d1,d2);
end function;

function CovariantFromHyperplanes(realforms,complexforms,prec)
// Finds the covariant positive definite Hermitian form from the
// given linear forms associated to a elliptic normal quintic
  P := Universe(realforms);
  PP := Universe(complexforms);
  R := BaseRing(P); // real field
  Rprint := RealField(5);
  mons := MD(P,2);
  A := Matrix(R,5,15,[MC(realforms[i]^2,mon): mon in mons, i in [1..5]]);
//  vprintf Reduction,3 :  
//    "Basis for 5-dimensional vector space :\n%o\n",ChangeRing(A,Rprint);
  cc := complexforms;
  scores := [&+[Im(MC(c,PP.i))^2: i in [1..5]]: c in complexforms];
  _,m := Minimum(scores);
  if m notin {1,3,5} then
    vprintf Reduction,1 : "Failed to recognise real form\n";
    return 0;
  end if;
  case m :
    when 1 : polysC := [cc[1]^2,cc[2]*cc[3],cc[4]*cc[5]];
    when 3 : polysC := [cc[1]*cc[2],cc[3]^2,cc[4]*cc[5]];
    when 5 : polysC := [cc[1]*cc[2],cc[3]*cc[4],cc[5]^2];
  end case;
  coeffs := [[MC(f,PP!mon): mon in mons]: f in polysC];
  imsize := &+[Im(c[i])^2: i in [1..15],c in coeffs];
  vprintf Reduction,2 : "Rounding error = %o\n",Rprint!imsize;
  if imsize gt 10^(-3) then 
    vprint Reduction,1 : "Rounding error is too large";
    return 0;
  end if;
  B := Matrix(R,3,15,[Re(c[i]): i in [1..15], c in coeffs]);
//  vprintf Reduction,3 : 
//    "Basis for 3-dimensional vector space :\n%o\n",ChangeRing(B,Rprint);
  vec,dist := CommonVector(A,B);
  vprintf Reduction,2 : "Minimum distance = %o\n",Rprint!dist;
  if dist gt 10^(-3) then 
    vprint Reduction,1 : "Minimum distance is too large";
    return 0;
  end if;
  Q := &+[vec[i]*mons[i] : i in [1..#mons]];
  if MC(Q,P.1^2) lt 0 then Q := -Q; end if;
  Qmat := Matrix(R,5,5,[Deriv(Deriv(Q,i),j): i,j in [1..5]]);
  if not IsPositiveDefinite(Qmat) then 
    vprint Reduction,1 : "Covariant matrix is not positive definite";
    return 0;
  end if;
  dd := &+[Abs(x)^2: x in Eltseq(Qmat)];
  Qmat *:= dd^(-1/2);
  return Qmat + Transpose(Qmat); 
end function;

RandomMats := [
  [1,0,0,2,-2,0,-7,4,-8,6,3,0,1,6,-6,0,5,-6,9,-6,0,-2,0,-1,1],
  [-1,4,2,1,0,-1,2,1,1,0,-1,5,3,0,0,-2,9,5,2,0,1,-3,-1,-1,1],
  [3,0,-2,4,6,1,1,-2,4,1,-1,0,1,-2,-2,3,0,0,1,7,5,0,0,1,12],
  [1,1,-1,0,-1,0,1,3,-2,0,-1,-1,1,1,1,1,2,1,-2,-1,-1,-3,-8,8,2],
  [1,-1,-4,-7,-2,-2,2,3,6,0,3,-2,-5,-11,-1,1,-1,-3,-5,-1,1,-1,-2,-4,0],
  [0,-2,1,2,3,1,1,0,0,-1,0,0,0,-1,-1,0,-1,1,2,3,-2,-1,-1,-1,1]
];

function RealMultipliers(c4,c6,prec)
  R := RealField(prec);
  P<X> := PolynomialRing(Rationals());
  if c6 ne 0 then 
    f1 := X^6 - 15*c4*X^4 + 40*c6*X^3 - 45*c4^2*X^2 + 24*c4*c6*X - 5*c6^2;
    f2 := (-2*X^5 + 30*c4*X^3 - 82*c6*X^2 + 90*c4^2*X - 33*c4*c6)/(3*c6);
    realrts := Roots(f1,R);
    vprintf Reduction,3 : 
      "RealRoots = %o\n",[RealField(15)|rt[1]: rt in realrts];
    ans := [[2*rt[1],Evaluate(f2,rt[1])]: rt in realrts];
    error if #ans ne 2,
       "Failed to compute real roots (prec = ",prec,
       ") - try increasing precision";
  else
    rt5 := SquareRoot(R!5);
    u := SquareRoot(3*rt5*c4);
    ans := [[sgn*u*(1+rt5),-3*rt5*c4]: sgn in [-1,1]];
  end if;
  return ans;
end function;

function ReductionCovariantFive(A,nos,prec)
  R := RealField(prec);
  C := ComplexField(prec);
  RR := PolynomialRing(R,5);
  CC := PolynomialRing(C,5);
  AR := ChangeRing(A,RR);
  qq := [ quads1 cat quads2
      where quads1 is [AR[2,j] - nn[1]*AR[1,j]: j in [1..5]] 
      where quads2 is [AR[3,j] - nn[2]*AR[1,j]: j in [1..5]]
		 : nn in nos ];
  ctr := 0;
  repeat
    if ctr eq 0 then 
      matR,matC := SolveForPentagons(qq,C,10^(-prec/2));
    elif ctr le #RandomMats then 
      T := Matrix(Integers(),5,5,RandomMats[ctr]);
      vprintf Reduction,2 : "We make a change of co-ordinates\n";
      TT := ChangeRing(T,R);
      qq1 := [[q^TT: q in q1]: q1 in qq];
      matR,matC := SolveForPentagons(qq1,C,10^(-prec/2));
      matR := matR*T^(-1);
      matC := matC*T^(-1);
    else 
      vprint Reduction,1 :
        "Problem was not resolved by a change of co-ordinates";
      return 0;
    end if;
    ctr +:= 1;
  until matR ne 0 and matC ne 0;
  forms1 := [&+[matR[i,j]*RR.j: j in [1..5]]: i in [1..5]];
  forms2 := [&+[matC[i,j]*CC.j: j in [1..5]]: i in [1..5]]; 
  Rprint := RealField(8);
  Cprint<i> := ComplexField(8);
  RRR<x1,x2,x3,x4,x5> := PolynomialRing(Rprint,5);
  CCC<x1,x2,x3,x4,x5> := PolynomialRing(Cprint,5);
  vprintf Reduction,3 : "Real pentagon defined by linear forms : \n%o\n",
    [&+[matR[i,j]*RRR.j: j in [1..5]]: i in [1..5]];
  vprintf Reduction,3 : "Imaginary pentagon defined by linear forms : \n%o\n",
    [&+[matC[i,j]*CCC.j: j in [1..5]]: i in [1..5]]; 
  return CovariantFromHyperplanes(forms1,forms2,prec);
end function;

//////////////////////////////////////////////////////////////
//                   Cases n = 2,3,4,5                      //
//////////////////////////////////////////////////////////////

function LLLReduceOnEquations(model)
  Q := Rationals();
  n := Degree(model);
  assert n in {4,5};
  eqns := Equations(model);
  R := Universe(eqns);
  mons := MD(R,2);
  coeffmat := Matrix(Q,#eqns,#mons,[MC(f,mon): mon in mons,f in eqns]);
  _,T := LLL(coeffmat : Delta:= 0.999);
  if n eq 5 then T := Transpose(T^(-1)); end if;
  tr := <ChangeRing(T,Q),Id(Q,n)>;
  return tr*model,tr;
end function;

function LLLReduceOnAmbient(model)
  Q := Rationals();
  n := Degree(model);
  assert n in {4,5};
  d := n eq 4 select 2 else 5;
  if n eq 4 then 
    mat := ChangeRing(Matrices(model)[1],Integers());
    _,T := LLLGram(mat); // This is (usually) indefinite LLL
  else 
    R := PolynomialRing(model);
    phi := Matrix(model);
    coeffmat := Matrix(Q,5,10,
      [[MC(phi[j,k],R.i):j,k in [1..5]|j lt k]: i in [1..5]]);
    _,T := LLL(coeffmat : Delta:= 0.999);
  end if;
  tr := <Id(Q,d),ChangeRing(T,Q)>;
  return tr*model,tr;
end function;

function AdHocReduce(model)
  Q := Rationals();
  n := Degree(model);
  assert n in {4,5};
  tr := n eq 4 select <Id(Q,2),Id(Q,4)> else <Id(Q,5),Id(Q,5)>;
  vprintf Reduction,2: "Performing ad hoc reduction \n";
  IndentPush();
  while true do
    model1,tr1 := LLLReduceOnEquations(model); 
    vprintf Reduction,2: 
      "LLL reduction on equations => %o\n",SizeString(model);
    if L2Norm(model1) lt L2Norm(model) then 
      model := model1; 
      tr := tr1*tr;
    end if;
    model1,tr1 := LLLReduceOnAmbient(model);
    vprintf Reduction,2: 
      "LLL reduction on ambient   => %o\n",SizeString(model);
    if L2Norm(model1) lt L2Norm(model) then 
      model := model1; 
      tr := tr1*tr;
    else 
      break; 
    end if; 
  end while;
  IndentPop();
  vprintf Reduction,1: "Ad hoc reduction     => %o\n",SizeString(model);
  return model,tr; 
end function;

function DefaultPrecision(n,coeffsize)
// TO DO : work out sensible default precisions
  case n :
    when 3 : return Max(10, coeffsize*3 div 2);  
       // from previous magma code
    when 4 : return 10 + 2*coeffsize; 
       // a wild guess 
    when 5 : return 50 + 12*coeffsize;
       // a guess based on a few examples
  end case;
end function;

intrinsic Reduce2008(model::ModelG1:prec:=0,Minkowski:=false) -> ModelG1, Tup
{A reduced genus one model integrally equivalent to the given one.}
  Q := Rationals();
  require CoefficientRing(model) cmpeq Q
    : "Model must be defined over the rationals";
  require Discriminant(model) ne 0 : "Model is singular";
  n := Degree(model);
  vprintf Reduction, 1 :"Reducing a genus one model of degree %o\n",n;
  require n in {2,3,4,5} : "Model must have degree 2,3,4 or 5.";
  require forall{ x : x in Eltseq(model) | IsIntegral(x) } 
    : "Coefficients of model must be integers";
  if n eq 2 then 
    return ReduceTwo(model); // see above
  end if;
  idtrans1 := case<n| 3:1, 4:Id(Q,2), 5:Id(Q,5), default:"">;
  tr := <idtrans1,Id(Q,n)>;
  vprintf Reduction,1: "Size of model given     %o\n",SizeString(model);
  if n in {4,5} then 
    model,tr1 := AdHocReduce(model);
    tr := tr1*tr;
  end if;
  coeffsize := CoefficientSize(model);
  if n eq 5 then syzmat,c4,c6 := SyzygeticMatrix(model); end if;
  if prec eq 0 then prec := DefaultPrecision(n,coeffsize); end if;
  vprintf Reduction,1 : "Working to real precision %o\n",prec;
  case n :
    when 3 : 
      torspts := RedTorsPts(Equation(model),prec);
    when 5 : 
      nos := RealMultipliers(c4,c6,prec);
  end case;
  history := {};
  done := false;
  earlyexit := false;
  while not done do
    vprintf Reduction,2 : 
      "Computing the reduction covariant (prec = %o)\n",prec;
    IndentPush();
    case n :
      when 3 : cov := ReductionInnerProduct(Equation(model),torspts,prec);
      when 4 : cov := ReductionCovariantFour(model,prec);
      when 5 : cov := ReductionCovariantFive(syzmat,nos,prec);
    end case;
    IndentPop();
    error if cov eq 0, 
      "Failed to compute reduction covariant ( prec =",prec,
      ") - try increasing precision";
    if n eq 3 then // scale to make nice for printing
      dd := &+[Abs(x)^2: x in Eltseq(cov)];
      cov *:= dd^(-1/2);
      cov := cov + Transpose(cov);
    end if;
    vprintf Reduction, 3:  
      "Reduction covariant:\n%o\n",ChangeRing(cov,RealField(Round(50/n)));
    if n eq 3 and IsDiagonal(cov) then 
       vprint Reduction, 2: "Diagonal matrix => early exit";
       earlyexit := true;
       break;
    end if;
    vprintf Reduction, 2: "Performing Heisenberg reduction\n";
    covZ := IntegerMatrix(cov); // Taken from earlier code for n = 4
    _,T := LLLGram(covZ : Delta:=0.999);
    vprintf Reduction, 3: "Transformation:\n%o\n", T;
    tr1 := <idtrans1,ChangeRing(T,Q)>;
    model := tr1*model; 
    if n in {4,5} then 
      model,tr2 := LLLReduceOnEquations(model);
      tr1 := tr2*tr1;
    end if;
    vprintf Reduction,1: "Heisenberg reduction => %o\n",SizeString(model);
    done := (model in history);
    Include(~history,model);
    if not done and n in {4,5} then 
      model,tr2 := AdHocReduce(model);
      tr1 := tr2*tr1;
    end if;
    if not done and n eq 5 then 
      vprint Reduction,2: "Updating the syzygetic matrix";
      syzmat := ApplyToMat(tr1,syzmat);
    end if;
// Perhaps could decrease precision for next iteration? 
    tr := tr1*tr;
  end while;
  if Minkowski and (not earlyexit) then 
    coeffsize := CoefficientSize(model);
    vprint Reduction,1: "Final stage - Minkowski reduction";
// Taken from earlier code for n = 3
    cov := tr1[2]*cov*Transpose(tr1[2]);
    c := 16;
    repeat
      c *:= 2;
      M := Matrix(n,n,[Round(c*e) : e in Eltseq(cov)]);
    until IsPositiveDefinite(M);
// Could we just use "IntegerMatrix" here?
    M := Matrix(n,n,[Round(1000*c*e) : e in Eltseq(cov)]);
// TO DO: Sort out why SuccessiveMinima takes so long, 
//        in some examples, especially for diagonal matrices!!!
    vprint Reduction, 3: "Calling SuccessiveMinima"; 
    _,mins1 := SuccessiveMinima(LatticeWithGram(M));
    T := Matrix([Eltseq(v) : v in mins1]);
    vprintf Reduction, 3: "Transformation:\n%o\n", T;
    tr1 := <idtrans1,ChangeRing(T,Q)>;
    model := tr1*model; 
    if n in {4,5} then 
      model,tr2 := LLLReduceOnEquations(model);
      tr1 := tr2*tr1;
    end if;
    tr := tr1*tr;
  end if;
  if n eq 4 then 
    tr1 := WomackTransformation(model);
// This is just a signed permutation to make the answer look pretty.
// Perhaps do something similar for n=3,5?
    model := tr1*model;
    tr := tr1*tr;
  end if;
  return model,tr;
end intrinsic;