attachment:g1reduction-2008.m of days23/CodingProjects

   1 /************************************************************************\
   2 *  Reduction of Genus One Models over Q (of degrees n = 2,3,4 and 5)     *
   3 *  Author: T. Fisher                                                     *
   4 *  Date: 22nd July 2008                                                  *
   5 *  Version for n = 2 calls "QuarticReduce"                               *
   6 *  Versions for n = 3,4 based on joint work with J. Cremona and M. Stoll *
   7 \************************************************************************/
   8 
   9 declare verbose Reduction,4;
  10 
  11 // Find pathname by typing e.g. "Hessian:Maximal;" at the magma prompt
  12 import "/usr/magma/package/Geometry/CrvG1/g1invariants.m" :
  13   SL4Covariants,MultiFact,AuxiliaryQuadrics,ExactPfaffianMatrix;
  14 import "/usr/magma/package/Geometry/CrvG1/3descent-testeq.m" :
  15   RedTorsPts,ReductionInnerProduct;
  16 
  17 MC := MonomialCoefficient;
  18 MD := MonomialsOfDegree;
  19 Id := IdentityMatrix;
  20 Deriv := Derivative;
  21 MySign := func<x|x ge 0 select 1 else -1>;
  22 dyadic_rep := func<x|x eq 0 select 0 
  23   else a*2^b where a,b is MantissaExponent(x)>;
  24 L2Norm := func<model|&+[x^2: x in Eltseq(model)]>;
  25 
  26 function FirstSign(seq)
  27   xx := [x : x in seq | x ne 0];
  28   fnz := #xx eq 0 select 1 else xx[1];
  29   return MySign(fnz);
  30 end function;
  31 
  32 function IntegerMatrix(mat)
  33   n := Nrows(mat);
  34   mat1 := Matrix(n,n,[dyadic_rep(mat[i,j]):i,j in [1..n]]);
  35   denom := LCM([Integers()|Denominator(x) : x in Eltseq(mat1)]);
  36   return MatrixAlgebra(Integers(),n)!(denom*mat1);
  37 end function;  
  38 
  39 function MyEchelonForm(mat,tol)
  40 // Row reduce an m by n matrix over the reals or complexes,
  41 // always using largest available entry as the next pivot.
  42 // Matrix returned is of form ( Id(m) | B )  
  43   m := Nrows(mat);
  44   n := Ncols(mat);
  45   for j in [1..m] do
  46     _,jj := Maximum([Abs(mat[k,j]): k in [j..m]]);
  47     jj +:= (j-1);
  48     temp := mat[j]; mat[j] := mat[jj]; mat[jj] := temp;
  49     if Abs(mat[j,j]) lt tol then return 0; end if;
  50     mat[j] *:= (1/mat[j,j]);
  51     for k in [1..m] do
  52       if k ne j then mat[k] -:= mat[k,j]*mat[j]; end if;
  53     end for;
  54   end for;
  55   return mat;
  56 end function;
  57 
  58 function SizeString(model)
  59   size := Round(L2Norm(model));
  60   logsize := Ilog(10,size);
  61   if logsize lt 10 then 
  62     return "L2Norm = " cat IntegerToString(size);
  63   else 
  64     return "Log10(L2Norm) = " cat IntegerToString(logsize);
  65   end if;
  66 end function;
  67 
  68 function CoefficientSize(model:name:="Model");
  69   cs := Ilog(10,Max([Round(Abs(c)): c in Eltseq(model)])); 
  70   if cs lt 10 then 
  71     vprintf Reduction, 2: "%o has coefficients\n%o\n",name,Eltseq(model);
  72   else
  73     vprintf Reduction, 2: "%o has coefficients of length %o\n",name,cs;
  74   end if;
  75   return cs;
  76 end function;
  77 
  78 //////////////////////////////////////////////////////////////
  79 //                      Case n = 2                          //
  80 //////////////////////////////////////////////////////////////
  81 
  82 function ReduceTwo(model)
  83   assert Degree(model) eq 2;
  84   P<X> := PolynomialRing(Integers());
  85   CrossTerms := (#Eltseq(model) eq 8);
  86   model1 := CrossTerms select CompleteTheSquare(model) else model;
  87   quartic := Evaluate(Equation(model1),[X,1]);
  88   _,M := QuarticReduce(quartic); // from CrvEll/FourDesc/quartred.m
  89   M := ChangeRing(M,Rationals());
  90   tr := <1,Transpose(M)>;
  91   model := tr*model;
  92   if CrossTerms then 
  93     coeffs := Eltseq(model);
  94     r := [ -(Integers()!(coeffs[i]) div 2) : i in [1..3]];
  95     tr1 := <1,r,Id(Rationals(),2)>;
  96     model := tr1*model;
  97     tr := tr1*tr;
  98   end if;
  99   return model,tr;
 100 end function;
 101 
 102 //////////////////////////////////////////////////////////////
 103 //                      Case n = 4                          //
 104 //////////////////////////////////////////////////////////////
 105 
 106 SRCompare := func<b,a| Abs(a[2]) ne Abs(b[2]) 
 107     select Abs(b[2])-Abs(a[2]) else Abs(b[3])-Abs(a[3]) >;
 108 
 109 function WomackTransformation(model)
 110 // Returns a pair of "signed permutations" in GL(2,Z) and GL(4,Z)
 111 // (just to make the final answer look pretty)
 112   Q := Rationals();
 113   assert Degree(model) eq 4;
 114   assert CoefficientRing(model) eq Q;
 115   R := PolynomialRing(model);
 116   q1,q2 := Explode(Equations(model));
 117   d0 := [<i,MC(q1,R.i^2),MC(q2,R.i^2)> : i in [1..4]]; 
 118   Sort(~d0, SRCompare);
 119   N1 := Matrix(Q,4,4,[<i,d0[i][1],1>: i in [1..4]]);
 120   model := <Id(Q,2),N1>*model;
 121   qq := Equations(model);
 122   signs := [Q|FirstSign([MC(q,R.i^2): i in [1..4]]): q in qq];
 123   M := DiagonalMatrix(signs);
 124   model := <M,Id(Q,4)>*model;
 125   q1,q2 := Explode(Equations(model));
 126   s1 := [Sign(MC(q1,R.1*R.j)): j in [1..4]];
 127   s2 := [MySign(MC(q2,R.1*R.j)): j in [1..4]];
 128   s := [Q|s1[i] ne 0 select s1[i] else s2[i]: i in [1..4]];
 129   N2 := DiagonalMatrix(s); 
 130   return <M,N2*N1>;
 131 end function;
 132 
 133 function ReductionCovariantFour(model,prec)
 134   Q := Rationals();
 135   R := RealField(prec);
 136   C := ComplexField(prec);
 137   P<X> := PolynomialRing(Q);
 138   a,bb,c,dd,e := Explode(SL4Invariants(model));
 139   covars := SL4Covariants(model); // from g1invariants.m
 140   g := a*X^4 + 2*bb*X^3 + c*X^2 + 2*dd*X + e;
 141   ff := [f[1]: f in Factorization(g)| f[1] ne X];
 142   quadrics := <>;
 143   if a eq 0 then 
 144     qA,T1,T2,qB := Explode(covars);
 145     Append(~quadrics,T2 - 2*bb*qB);
 146   end if;    
 147   if e eq 0 then 
 148     qA,T1,T2,qB := Explode(covars);
 149     Append(~quadrics,T1 - 2*dd*qA);
 150   end if;    
 151   for f in ff do
 152     if Degree(f) gt 1 then
 153       F<theta> := NumberField(f);
 154 // Since we use "Conjugates" below, we need to set the
 155 // precision of F. This seems quite slow for large precisons, so it 
 156 // might be worth re-writing to avoid the Kant machinery.
 157       SetKantPrecision(F,prec);
 158     else
 159       F := Q; 
 160       theta := Roots(f)[1][1];
 161     end if;   
 162     RR := PolynomialRing(F,4);
 163     covarsF := ChangeUniverse(covars,RR);
 164     qA,T1,T2,qB := Explode(covarsF);
 165     GG := theta^(-1)*e*qA + T1 + theta*T2 + theta^2*a*qB;
 166     Append(~quadrics,GG);
 167   end for;
 168   mats := <>;
 169   diagelts := [];
 170   for GG in quadrics do
 171     RR := Parent(GG);
 172     F := BaseRing(RR);
 173     G := Factorization(GG)[1][1];
 174     alpha := F!ExactQuotient(GG,G^2);
 175     assert alpha ne 0;
 176     assert GG eq alpha*G^2;
 177     M1 := Matrix(C,4,Degree(F),[Conjugates(MC(G,RR.i)): i in [1..4]]);
 178     Append(~mats,M1); 
 179     aa := [C|Abs(x): x in Conjugates(alpha)];
 180     dg := [C|Abs(x): x in Conjugates(Evaluate(Derivative(g),theta))];
 181     ll := [aa[i]/dg[i]^(1/2) : i in [1..Degree(F)]];
 182     diagelts cat:= ll;
 183   end for;
 184   M1 := HorizontalJoin(mats);
 185   M2 := Matrix(C,4,4,[Conjugate(M1[j,i]):i,j in [1..4]]);
 186   D := DiagonalMatrix(C,diagelts);
 187   M := M1*D*M2;
 188   Qmat := Matrix(R,4,4,[Re(M[i,j]): i,j in [1..4]]);
 189   dd := &+[Abs(x)^2: x in Eltseq(Qmat)];
 190   Qmat *:= dd^(-1/2);
 191   return Qmat + Transpose(Qmat); 
 192 end function;
 193 
 194 //////////////////////////////////////////////////////////////
 195 //                      Case n = 5                          //
 196 //////////////////////////////////////////////////////////////
 197 
 198 function InvariantsAndHessian(model)
 199 // Computes the invariants c4,c6 and Hessian H of a genus one model 
 200 // of degree 5. (Calling the existing functions cInvariants and 
 201 // Hessian separately leads to some duplication of effort.)
 202   assert Degree(model) eq 5;
 203   assert BaseRing(model) cmpeq Rationals();
 204   phi := Matrix(model);
 205   p := Equations(model);
 206   q := AuxiliaryQuadrics(p); // from g1invariants.m
 207   R := BaseRing(phi);
 208   RR<u> := PolynomialRing(R);
 209   SS<t> := quo<RR|u^4>; // this trick would speed up cInvariants
 210   U := Matrix(SS,5,5,[&+[MC(phi[j,k],R.i)*R.k: k in [1..5]]: i,j in [1..5]]);
 211   P := Matrix(R,5,5,
 212          [&+[MC(Deriv(p[k],j),R.i)*R.k: k in [1..5]]: i,j in [1..5]]);
 213   Q := Matrix(SS,5,5,[Deriv(q[i],j):i,j in [1..5]]);
 214   det1 := Determinant(P);
 215   det2 := Determinant(U+t*Q);
 216   T := [R|Coefficient(det2,i): i in [1,3]];
 217   cc := [&+[MultiFact(mon)*MC(det1,mon)*MC(d2,mon):mon in MD(R,5)]:d2 in T];
 218   c4 := cc[1]/40;
 219   c6 := -cc[2]/320;
 220   RR<y1,y2,y3,y4,y5> := PolynomialRing(R,5);
 221   M := Determinant(Matrix(RR,5,5,
 222          [&+[MC(Deriv(p[k],j),R.i)*RR.k : k in [1..5]]:i,j in [1..5]]));
 223   for ctr in [1..2] do
 224     M := &+[R.i*MC(q[i],R.j*R.k)*Deriv(Deriv(M,RR.j),RR.k)
 225 	   : j,k in [1..5], i in [1..5] | j le k ];
 226   end for;
 227   q2m := [R|MC(M,RR.i): i in [1..5]];
 228   p22 := Eltseq(4*c4*Vector(p)-(3/16)*Vector(q2m));
 229   hessian := ExactPfaffianMatrix(p22,1,1); // from g1invariants.m
 230   coeffs := func<f|[MC(f,R.k): k in [1..5]]>;
 231   Umat := Matrix(10,5,[coeffs(phi[i,j]): i,j in [1..5] | i lt j]);
 232   Hmat := Matrix(10,5,[coeffs(hessian[i,j]): i,j in [1..5] | i lt j]);  
 233   Dmat := HorizontalJoin(Umat,Hmat);
 234   discr1 := Determinant(Dmat);
 235   discr2 := 12^2*(c4^3-c6^2);
 236   if discr1 eq -discr2 then hessian*:= -1; end if;
 237   return c4,c6,GenusOneModel(hessian);
 238 end function;
 239 
 240 function SyzygeticMatrix(model)
 241 // This is a 3 by 5 matrix of quadratic forms. The first row contains
 242 // the 4 by 4 Pfaffians of the model, and the last row contains the
 243 // the 4 by 4 Pfaffians of the Hessian.
 244 // The Hesse family has 12 singular fibres: the "syzygetic pentagons".
 245 // These pentagons have 30 vertices in total. These make up the locus 
 246 // where the syzygetic matrix has rank 1. (The total space of the Hesse 
 247 // family is the locus where the matrix has rank 2.)
 248   assert Degree(model) eq 5;
 249   assert BaseRing(model) cmpeq Rationals();
 250   P<t> := PolynomialRing(Rationals());
 251   vprintf Reduction,2 : 
 252     "Computing the invariants, Hessian and syzygetic matrix\n";
 253   c4,c6,H := InvariantsAndHessian(model);
 254   vprintf Reduction,3 : "cInvariants = %o\n",[c4,c6]; 
 255   coeffsize := CoefficientSize(H:name:="Hessian");
 256   F := t*Vector(P,50,Eltseq(model))+Vector(P,50,Eltseq(H));
 257   FF := GenusOneModel(5,Eltseq(F));
 258   eqns := Equations(FF);
 259   R := PolynomialRing(model);
 260   RR := Universe(eqns);
 261   mons := MD(R,2);
 262   coeffs := [[MC(f,RR!mon): mon in mons]: f in eqns];
 263   p2 := [&+[Coefficient(c[j],2)*mons[j]: j in [1..#mons]]: c in coeffs];
 264   p12 := [&+[Coefficient(c[j],1)*mons[j]: j in [1..#mons]]: c in coeffs];
 265   p22 := [&+[Coefficient(c[j],0)*mons[j]: j in [1..#mons]]: c in coeffs];
 266   assert p2 eq Equations(model);
 267   assert p22 eq Equations(H);
 268   return Matrix(R,3,5,[p2,p12,p22]),c4,c6;
 269 end function;
 270 
 271 function ApplyToMat(tr,mat)
 272 // Since the Hessian is a covariant, there is no need to to keep 
 273 // recomputing the syzygetic matrix.
 274   R := BaseRing(mat);
 275   S := ChangeRing(Adjoint(tr[1]),R);
 276   T := Transpose(tr[2]);
 277   mat := Matrix(3,5,[mat[i,j]^T:j in [1..5],i in [1..3]]);
 278   return mat*S;
 279 end function;
 280 
 281 function SolveForPentagons(qq,C,tol)
 282   matR := 0;
 283   matC := 0;
 284   for quads in qq do 
 285     assert #quads eq 10;   // 10 homogeneous quadrics 
 286     RR := Universe(quads);
 287     assert Rank(RR) eq 5;  // in 5 variables 
 288     R := BaseRing(RR);     // over the reals
 289     mons := [RR.i*RR.j: i,j in [1..4]|i le j] cat [RR.i*RR.5: i in [1..5]];
 290     coeffmat := Matrix(10,15,[MC(q,mon): mon in mons,q in quads]);
 291     coeffmat := MyEchelonForm(coeffmat,tol);
 292     if coeffmat eq 0 then 
 293       vprintf Reduction,2 : "Pentagon not in general position\n";
 294       return 0,0; 
 295     end if;
 296     charmat := Matrix(4,5,[-coeffmat[i,j]: j in [11..15],i in [1..4]]);
 297     charmat := VerticalJoin(charmat,Vector(R,5,[1,0,0,0,0]));
 298     evals := Eigenvalues(charmat);
 299     s := func<n|n eq 1 select "ex" else "ices">;
 300     vprintf Reduction,2 : "Pentagon has %o real vert%o\n",#evals,s(#evals);
 301     realflag := true;
 302     if #evals ne 5 then
 303       charmat := ChangeRing(charmat,C);
 304       charpoly := CharacteristicPolynomial(charmat);
 305       evals := Roots(charpoly,C);
 306 // Don't use Eigenvalues(charmat), since in "CovariantFromHyperplanes" 
 307 // we assume the pairs of complex conjugate roots are adjacent.
 308       if #evals ne 5 then 
 309         vprintf Reduction,2 : "Unable to find all complex vertices\n";
 310         return 0,0; 
 311       end if;
 312       realflag := false;
 313     end if;
 314     F := realflag select R else C;
 315     evals := [F| x[1]: x in evals];
 316     mm := Minimum([Abs(evals[i]-evals[j]):i,j in [1..5]|i lt j]);
 317     if mm lt tol then 
 318       vprintf Reduction,2 : "Unable to distinguish vertices\n";
 319       return 0,0;
 320     end if;
 321     VdM := Matrix(F,5,5,[x^i : x in evals,i in [0..4]]);
 322     v := Vector(F,5,[0,0,0,0,1]);
 323     vv := [];
 324     for i in [1..5] do vv cat:= [v]; v := v*charmat; end for;
 325     augmat := HorizontalJoin(VdM,VerticalJoin(vv));
 326     augmat := MyEchelonForm(augmat,tol);
 327     if augmat eq 0 then return 0,0; end if;
 328     mat := Matrix(5,5,[augmat[i,j]: j in [6..10],i in [1..5]]);
 329     dd := [1/Sqrt(&+[Abs(x)^2: x in Eltseq(mat[i])]): i in [1..5]];
 330     mat := DiagonalMatrix(F,dd)*mat;
 331     if realflag then matR := mat; else matC := mat; end if;
 332   end for;
 333   return matR,matC;
 334 end function;
 335 
 336 function GramSchmidt(mat)
 337   d := Ncols(mat);
 338   for i := 1 to Nrows(mat) do
 339     s := [&+[mat[j,k]*mat[i,k]: k in [1..d]]: j in [1..i-1]];
 340     if i gt 1 then
 341       mat[i] -:= &+[s[j]*mat[j]: j in [1..i-1]];
 342     end if;
 343     n := Sqrt(&+[mat[i,j]^2 : j in [1..d]]);
 344     mat[i] *:= (1/n);
 345   end for;
 346   return mat;
 347 end function;
 348 
 349 function MyProjection(mat,vec)
 350   m := Nrows(mat);
 351   n := Ncols(mat);
 352   a := [&+[mat[j,k]*vec[k]: k in [1..n]]: j in [1..m]];
 353   vec -:= &+[a[j]*mat[j]: j in [1..m]];
 354   n := Sqrt(&+[vec[i]^2 : i in [1..n]]);
 355   return (1/n)*vec;
 356 end function;
 357 
 358 function RotationMatrix(theta)
 359   c := Cos(theta); 
 360   s := Sin(theta);
 361   return Matrix(2,2,[c,s,-s,c]);
 362 end function;
 363 
 364 function ClosestPointOnCircle(A,B)
 365 // Given matrices A and B whose rows are orthonormal bases
 366 // for subspaces V and W of Euclidean space R^d, finds 
 367 // cos(phi)*B[1] + sin(phi)*B[2] such that distance to V
 368 // is minimal.
 369   assert Ncols(A) eq Ncols(B);
 370   assert Nrows(B) eq 2;
 371   n := Nrows(A);
 372   M := B*Transpose(A); 
 373   x := &+[M[1,i]^2-M[2,i]^2: i in [1..n]];
 374   y := &+[2*M[1,i]*M[2,i]: i in [1..n]];
 375   phi := Arctan2(x,y)/2; // N.B. don't reverse the arguments!
 376   R := RotationMatrix(phi);
 377   S := R*M;  
 378   dd := [1-&+[S[i,j]^2: j in [1..n]]: i in [1..2]];
 379   dist,m := Minimum(dd);
 380   return R[m]*B,dist;
 381 end function;
 382 
 383 function CommonVector(A,B)
 384 // Computes a basis for the intersection of the row spaces 
 385 // of A and B, where the intersection is assumed 1-dimensional 
 386   assert Nrows(A) eq 5;
 387   assert Nrows(B) eq 3;
 388   A := GramSchmidt(A);
 389   B := GramSchmidt(B);
 390   M := A*Transpose(B);
 391   dd := [1-&+[M[i,j]^2: i in [1..5]]: j in [1..3]];
 392   _,m := Maximum(dd);
 393   AA := VerticalJoin(A,MyProjection(A,B[m]));
 394   BB := Matrix([Eltseq(B[i]): i in [1..3]| i ne m]);
 395   vec,d1 := ClosestPointOnCircle(AA,BB);
 396   BB := VerticalJoin(vec,B[m]);
 397   vec,d2 := ClosestPointOnCircle(A,BB); 
 398   return vec,Maximum(d1,d2);
 399 end function;
 400 
 401 function CovariantFromHyperplanes(realforms,complexforms,prec)
 402 // Finds the covariant positive definite Hermitian form from the
 403 // given linear forms associated to a elliptic normal quintic
 404   P := Universe(realforms);
 405   PP := Universe(complexforms);
 406   R := BaseRing(P); // real field
 407   Rprint := RealField(5);
 408   mons := MD(P,2);
 409   A := Matrix(R,5,15,[MC(realforms[i]^2,mon): mon in mons, i in [1..5]]);
 410 //  vprintf Reduction,3 :  
 411 //    "Basis for 5-dimensional vector space :\n%o\n",ChangeRing(A,Rprint);
 412   cc := complexforms;
 413   scores := [&+[Im(MC(c,PP.i))^2: i in [1..5]]: c in complexforms];
 414   _,m := Minimum(scores);
 415   if m notin {1,3,5} then
 416     vprintf Reduction,1 : "Failed to recognise real form\n";
 417     return 0;
 418   end if;
 419   case m :
 420     when 1 : polysC := [cc[1]^2,cc[2]*cc[3],cc[4]*cc[5]];
 421     when 3 : polysC := [cc[1]*cc[2],cc[3]^2,cc[4]*cc[5]];
 422     when 5 : polysC := [cc[1]*cc[2],cc[3]*cc[4],cc[5]^2];
 423   end case;
 424   coeffs := [[MC(f,PP!mon): mon in mons]: f in polysC];
 425   imsize := &+[Im(c[i])^2: i in [1..15],c in coeffs];
 426   vprintf Reduction,2 : "Rounding error = %o\n",Rprint!imsize;
 427   if imsize gt 10^(-3) then 
 428     vprint Reduction,1 : "Rounding error is too large";
 429     return 0;
 430   end if;
 431   B := Matrix(R,3,15,[Re(c[i]): i in [1..15], c in coeffs]);
 432 //  vprintf Reduction,3 : 
 433 //    "Basis for 3-dimensional vector space :\n%o\n",ChangeRing(B,Rprint);
 434   vec,dist := CommonVector(A,B);
 435   vprintf Reduction,2 : "Minimum distance = %o\n",Rprint!dist;
 436   if dist gt 10^(-3) then 
 437     vprint Reduction,1 : "Minimum distance is too large";
 438     return 0;
 439   end if;
 440   Q := &+[vec[i]*mons[i] : i in [1..#mons]];
 441   if MC(Q,P.1^2) lt 0 then Q := -Q; end if;
 442   Qmat := Matrix(R,5,5,[Deriv(Deriv(Q,i),j): i,j in [1..5]]);
 443   if not IsPositiveDefinite(Qmat) then 
 444     vprint Reduction,1 : "Covariant matrix is not positive definite";
 445     return 0;
 446   end if;
 447   dd := &+[Abs(x)^2: x in Eltseq(Qmat)];
 448   Qmat *:= dd^(-1/2);
 449   return Qmat + Transpose(Qmat); 
 450 end function;
 451 
 452 RandomMats := [
 453   [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],
 454   [-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],
 455   [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],
 456   [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],
 457   [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],
 458   [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]
 459 ];
 460 
 461 function RealMultipliers(c4,c6,prec)
 462   R := RealField(prec);
 463   P<X> := PolynomialRing(Rationals());
 464   if c6 ne 0 then 
 465     f1 := X^6 - 15*c4*X^4 + 40*c6*X^3 - 45*c4^2*X^2 + 24*c4*c6*X - 5*c6^2;
 466     f2 := (-2*X^5 + 30*c4*X^3 - 82*c6*X^2 + 90*c4^2*X - 33*c4*c6)/(3*c6);
 467     realrts := Roots(f1,R);
 468     vprintf Reduction,3 : 
 469       "RealRoots = %o\n",[RealField(15)|rt[1]: rt in realrts];
 470     ans := [[2*rt[1],Evaluate(f2,rt[1])]: rt in realrts];
 471     error if #ans ne 2,
 472        "Failed to compute real roots (prec = ",prec,
 473        ") - try increasing precision";
 474   else
 475     rt5 := SquareRoot(R!5);
 476     u := SquareRoot(3*rt5*c4);
 477     ans := [[sgn*u*(1+rt5),-3*rt5*c4]: sgn in [-1,1]];
 478   end if;
 479   return ans;
 480 end function;
 481 
 482 function ReductionCovariantFive(A,nos,prec)
 483   R := RealField(prec);
 484   C := ComplexField(prec);
 485   RR := PolynomialRing(R,5);
 486   CC := PolynomialRing(C,5);
 487   AR := ChangeRing(A,RR);
 488   qq := [ quads1 cat quads2
 489       where quads1 is [AR[2,j] - nn[1]*AR[1,j]: j in [1..5]] 
 490       where quads2 is [AR[3,j] - nn[2]*AR[1,j]: j in [1..5]]
 491 		 : nn in nos ];
 492   ctr := 0;
 493   repeat
 494     if ctr eq 0 then 
 495       matR,matC := SolveForPentagons(qq,C,10^(-prec/2));
 496     elif ctr le #RandomMats then 
 497       T := Matrix(Integers(),5,5,RandomMats[ctr]);
 498       vprintf Reduction,2 : "We make a change of co-ordinates\n";
 499       TT := ChangeRing(T,R);
 500       qq1 := [[q^TT: q in q1]: q1 in qq];
 501       matR,matC := SolveForPentagons(qq1,C,10^(-prec/2));
 502       matR := matR*T^(-1);
 503       matC := matC*T^(-1);
 504     else 
 505       vprint Reduction,1 :
 506         "Problem was not resolved by a change of co-ordinates";
 507       return 0;
 508     end if;
 509     ctr +:= 1;
 510   until matR ne 0 and matC ne 0;
 511   forms1 := [&+[matR[i,j]*RR.j: j in [1..5]]: i in [1..5]];
 512   forms2 := [&+[matC[i,j]*CC.j: j in [1..5]]: i in [1..5]]; 
 513   Rprint := RealField(8);
 514   Cprint<i> := ComplexField(8);
 515   RRR<x1,x2,x3,x4,x5> := PolynomialRing(Rprint,5);
 516   CCC<x1,x2,x3,x4,x5> := PolynomialRing(Cprint,5);
 517   vprintf Reduction,3 : "Real pentagon defined by linear forms : \n%o\n",
 518     [&+[matR[i,j]*RRR.j: j in [1..5]]: i in [1..5]];
 519   vprintf Reduction,3 : "Imaginary pentagon defined by linear forms : \n%o\n",
 520     [&+[matC[i,j]*CCC.j: j in [1..5]]: i in [1..5]]; 
 521   return CovariantFromHyperplanes(forms1,forms2,prec);
 522 end function;
 523 
 524 //////////////////////////////////////////////////////////////
 525 //                   Cases n = 2,3,4,5                      //
 526 //////////////////////////////////////////////////////////////
 527 
 528 function LLLReduceOnEquations(model)
 529   Q := Rationals();
 530   n := Degree(model);
 531   assert n in {4,5};
 532   eqns := Equations(model);
 533   R := Universe(eqns);
 534   mons := MD(R,2);
 535   coeffmat := Matrix(Q,#eqns,#mons,[MC(f,mon): mon in mons,f in eqns]);
 536   _,T := LLL(coeffmat : Delta:= 0.999);
 537   if n eq 5 then T := Transpose(T^(-1)); end if;
 538   tr := <ChangeRing(T,Q),Id(Q,n)>;
 539   return tr*model,tr;
 540 end function;
 541 
 542 function LLLReduceOnAmbient(model)
 543   Q := Rationals();
 544   n := Degree(model);
 545   assert n in {4,5};
 546   d := n eq 4 select 2 else 5;
 547   if n eq 4 then 
 548     mat := ChangeRing(Matrices(model)[1],Integers());
 549     _,T := LLLGram(mat); // This is (usually) indefinite LLL
 550   else 
 551     R := PolynomialRing(model);
 552     phi := Matrix(model);
 553     coeffmat := Matrix(Q,5,10,
 554       [[MC(phi[j,k],R.i):j,k in [1..5]|j lt k]: i in [1..5]]);
 555     _,T := LLL(coeffmat : Delta:= 0.999);
 556   end if;
 557   tr := <Id(Q,d),ChangeRing(T,Q)>;
 558   return tr*model,tr;
 559 end function;
 560 
 561 function AdHocReduce(model)
 562   Q := Rationals();
 563   n := Degree(model);
 564   assert n in {4,5};
 565   tr := n eq 4 select <Id(Q,2),Id(Q,4)> else <Id(Q,5),Id(Q,5)>;
 566   vprintf Reduction,2: "Performing ad hoc reduction \n";
 567   IndentPush();
 568   while true do
 569     model1,tr1 := LLLReduceOnEquations(model); 
 570     vprintf Reduction,2: 
 571       "LLL reduction on equations => %o\n",SizeString(model);
 572     if L2Norm(model1) lt L2Norm(model) then 
 573       model := model1; 
 574       tr := tr1*tr;
 575     end if;
 576     model1,tr1 := LLLReduceOnAmbient(model);
 577     vprintf Reduction,2: 
 578       "LLL reduction on ambient   => %o\n",SizeString(model);
 579     if L2Norm(model1) lt L2Norm(model) then 
 580       model := model1; 
 581       tr := tr1*tr;
 582     else 
 583       break; 
 584     end if; 
 585   end while;
 586   IndentPop();
 587   vprintf Reduction,1: "Ad hoc reduction     => %o\n",SizeString(model);
 588   return model,tr; 
 589 end function;
 590 
 591 function DefaultPrecision(n,coeffsize)
 592 // TO DO : work out sensible default precisions
 593   case n :
 594     when 3 : return Max(10, coeffsize*3 div 2);  
 595        // from previous magma code
 596     when 4 : return 10 + 2*coeffsize; 
 597        // a wild guess 
 598     when 5 : return 50 + 12*coeffsize;
 599        // a guess based on a few examples
 600   end case;
 601 end function;
 602 
 603 intrinsic Reduce2008(model::ModelG1:prec:=0,Minkowski:=false) -> ModelG1, Tup
 604 {A reduced genus one model integrally equivalent to the given one.}
 605   Q := Rationals();
 606   require CoefficientRing(model) cmpeq Q
 607     : "Model must be defined over the rationals";
 608   require Discriminant(model) ne 0 : "Model is singular";
 609   n := Degree(model);
 610   vprintf Reduction, 1 :"Reducing a genus one model of degree %o\n",n;
 611   require n in {2,3,4,5} : "Model must have degree 2,3,4 or 5.";
 612   require forall{ x : x in Eltseq(model) | IsIntegral(x) } 
 613     : "Coefficients of model must be integers";
 614   if n eq 2 then 
 615     return ReduceTwo(model); // see above
 616   end if;
 617   idtrans1 := case<n| 3:1, 4:Id(Q,2), 5:Id(Q,5), default:"">;
 618   tr := <idtrans1,Id(Q,n)>;
 619   vprintf Reduction,1: "Size of model given     %o\n",SizeString(model);
 620   if n in {4,5} then 
 621     model,tr1 := AdHocReduce(model);
 622     tr := tr1*tr;
 623   end if;
 624   coeffsize := CoefficientSize(model);
 625   if n eq 5 then syzmat,c4,c6 := SyzygeticMatrix(model); end if;
 626   if prec eq 0 then prec := DefaultPrecision(n,coeffsize); end if;
 627   vprintf Reduction,1 : "Working to real precision %o\n",prec;
 628   case n :
 629     when 3 : 
 630       torspts := RedTorsPts(Equation(model),prec);
 631     when 5 : 
 632       nos := RealMultipliers(c4,c6,prec);
 633   end case;
 634   history := {};
 635   done := false;
 636   earlyexit := false;
 637   while not done do
 638     vprintf Reduction,2 : 
 639       "Computing the reduction covariant (prec = %o)\n",prec;
 640     IndentPush();
 641     case n :
 642       when 3 : cov := ReductionInnerProduct(Equation(model),torspts,prec);
 643       when 4 : cov := ReductionCovariantFour(model,prec);
 644       when 5 : cov := ReductionCovariantFive(syzmat,nos,prec);
 645     end case;
 646     IndentPop();
 647     error if cov eq 0, 
 648       "Failed to compute reduction covariant ( prec =",prec,
 649       ") - try increasing precision";
 650     if n eq 3 then // scale to make nice for printing
 651       dd := &+[Abs(x)^2: x in Eltseq(cov)];
 652       cov *:= dd^(-1/2);
 653       cov := cov + Transpose(cov);
 654     end if;
 655     vprintf Reduction, 3:  
 656       "Reduction covariant:\n%o\n",ChangeRing(cov,RealField(Round(50/n)));
 657     if n eq 3 and IsDiagonal(cov) then 
 658        vprint Reduction, 2: "Diagonal matrix => early exit";
 659        earlyexit := true;
 660        break;
 661     end if;
 662     vprintf Reduction, 2: "Performing Heisenberg reduction\n";
 663     covZ := IntegerMatrix(cov); // Taken from earlier code for n = 4
 664     _,T := LLLGram(covZ : Delta:=0.999);
 665     vprintf Reduction, 3: "Transformation:\n%o\n", T;
 666     tr1 := <idtrans1,ChangeRing(T,Q)>;
 667     model := tr1*model; 
 668     if n in {4,5} then 
 669       model,tr2 := LLLReduceOnEquations(model);
 670       tr1 := tr2*tr1;
 671     end if;
 672     vprintf Reduction,1: "Heisenberg reduction => %o\n",SizeString(model);
 673     done := (model in history);
 674     Include(~history,model);
 675     if not done and n in {4,5} then 
 676       model,tr2 := AdHocReduce(model);
 677       tr1 := tr2*tr1;
 678     end if;
 679     if not done and n eq 5 then 
 680       vprint Reduction,2: "Updating the syzygetic matrix";
 681       syzmat := ApplyToMat(tr1,syzmat);
 682     end if;
 683 // Perhaps could decrease precision for next iteration? 
 684     tr := tr1*tr;
 685   end while;
 686   if Minkowski and (not earlyexit) then 
 687     coeffsize := CoefficientSize(model);
 688     vprint Reduction,1: "Final stage - Minkowski reduction";
 689 // Taken from earlier code for n = 3
 690     cov := tr1[2]*cov*Transpose(tr1[2]);
 691     c := 16;
 692     repeat
 693       c *:= 2;
 694       M := Matrix(n,n,[Round(c*e) : e in Eltseq(cov)]);
 695     until IsPositiveDefinite(M);
 696 // Could we just use "IntegerMatrix" here?
 697     M := Matrix(n,n,[Round(1000*c*e) : e in Eltseq(cov)]);
 698 // TO DO: Sort out why SuccessiveMinima takes so long, 
 699 //        in some examples, especially for diagonal matrices!!!
 700     vprint Reduction, 3: "Calling SuccessiveMinima"; 
 701     _,mins1 := SuccessiveMinima(LatticeWithGram(M));
 702     T := Matrix([Eltseq(v) : v in mins1]);
 703     vprintf Reduction, 3: "Transformation:\n%o\n", T;
 704     tr1 := <idtrans1,ChangeRing(T,Q)>;
 705     model := tr1*model; 
 706     if n in {4,5} then 
 707       model,tr2 := LLLReduceOnEquations(model);
 708       tr1 := tr2*tr1;
 709     end if;
 710     tr := tr1*tr;
 711   end if;
 712   if n eq 4 then 
 713     tr1 := WomackTransformation(model);
 714 // This is just a signed permutation to make the answer look pretty.
 715 // Perhaps do something similar for n=3,5?
 716     model := tr1*model;
 717     tr := tr1*tr;
 718   end if;
 719   return model,tr;
 720 end intrinsic;