attachment:g1minimisation-2008.m of days23/CodingProjects

   1 /************************************************************************\
   2 *  Minimisation of Genus One Models over Q (of degrees n = 2,3,4 and 5)  *
   3 *  Author: T. Fisher (based on joint work with J. Cremona and M. Stoll)  *
   4 *  Date: 22nd July 2008                                                  *
   5 *  Version for n = 5 is not yet proven to work in all cases              *
   6 *  This file also contains an intrinsic "GenusOneModel" that takes as    *
   7 *    input n and P, and computes a minimal genus one model of degree n   * 
   8 *    for (E,[(n-1)0+P]) where P in E(Q).                                 *
   9 \************************************************************************/
  10 
  11 declare verbose Minimisation,3;
  12 
  13 MC := MonomialCoefficient;
  14 MD := MonomialsOfDegree;
  15 Diag := DiagonalMatrix;
  16 Id := IdentityMatrix;
  17 Char := Characteristic;
  18 Deriv := Derivative;
  19 MAT := func<R,n,Q|Matrix(R,n,n,[Deriv(Deriv(Q,i),j): i,j in [1..n]])>;
  20 
  21 function RationalGCD(S)
  22    Z := Integers();
  23    Q := Rationals();
  24    d := LCM([Denominator(Q!x):x in S| x ne 0] cat [1]);
  25    return Universe(S)!(GCD([Z!(x*d): x in S])/d);
  26 end function;
  27 
  28 PRIM := func<seq|[Integers()|x/d: x in seq] where d is RationalGCD(seq)>;
  29 
  30 function MySaturate(mat)
  31   n := Nrows(mat);
  32   D := DiagonalMatrix([RationalGCD(Eltseq(mat[i])): i in [1..n]]);
  33   mat := D^(-1)*mat;
  34   L := Lattice(mat);
  35   L1 := PureLattice(L);
  36   return D*Matrix(Rationals(),n,n,[Coordinates(L1,L!mat[i]): i in [1..n]]);
  37 end function;
  38 
  39 function CoeffMatrix(polys,d)
  40   R := Universe(polys);
  41   mons := MD(R,d);
  42   return Matrix(BaseRing(R),#polys,#mons,[MC(f,mon): mon in mons,f in polys]);
  43 end function;  
  44 
  45 function GlobalLevel(phi)
  46   vprintf Minimisation,1 : "Computing the invariants\n";
  47   E := Jacobian(phi);
  48   _,iso := MinimalModel(E);
  49   u := IsomorphismData(iso)[4];
  50 // When n = 5, we compensate for the fact "Jacobian" uses the 
  51 // cInvariants (instead of the aInvariants).
  52   if Degree(phi) eq 5 then u := 6*u; end if;
  53   return Integers()!(1/u);
  54 end function;
  55 
  56 function MultipleRoots(F,m) 
  57 // Find roots of F of multiplicity at least m
  58   k := BaseRing(F);
  59   if Char(k) gt m then
  60     for i in [2..m] do 
  61       F := GCD([F,Deriv(F)]); 
  62     end for;
  63     rts := [rt[1]: rt in Roots(F)];
  64   else
  65     rts := [rt[1]: rt in Roots(F)| rt[2] ge m];
  66   end if;
  67   return rts;
  68 end function;
  69 
  70 function PrettyFactorization(F) 
  71 // Only used for verbose printing 
  72   P<t> := Parent(F);
  73   if F eq 0 or Degree(F) eq 0 then return Sprint(F); end if;
  74   ff := Factorization(F);
  75   str := "";
  76   for i in [1..#ff] do
  77     f := ff[i];
  78     if (f[1] eq t) or (#ff eq 1 and f[2] eq 1) then 
  79       str cat:= Sprint(f[1]);
  80     else
  81       str cat:= ("(" cat Sprint(f[1]) cat ")");
  82     end if;
  83     if f[2] ne 1 then 
  84       str cat:= ("^" cat Sprint(f[2]));
  85     end if;
  86     if i lt #ff then str cat:= "*";end if;
  87   end for;
  88   return str;
  89 end function;
  90 
  91 function PrettySubstitution(tr)
  92 // Only used for verbose printing 
  93   R<x,z> := PolynomialRing(Rationals(),2);
  94   S<y> := PolynomialRing(R);
  95   mons := [x^2,x*z,x^2];
  96   rr := &+[tr[2][i]*mons[i]: i in [1..3]];
  97   return "y <- " cat Sprint((1/tr[1])*y + rr);
  98 end function;
  99 
 100 //////////////////////////////////////////////////////////////
 101 //                      Case n = 2                          //
 102 //////////////////////////////////////////////////////////////
 103 
 104 function RemoveCrossTerms(phi)
 105 // Same as "CompleteTheSquare", but also returns the transformation. 
 106   Q := Rationals();
 107   seq := Eltseq(phi);
 108   if #seq eq 5 then 
 109     tr := <1,Id(Q,2)>;
 110   else
 111     a0,a1,a2 := Explode(seq);
 112     tr := <2,[-a0/2,-a1/2,-a2/2],Id(Q,2)>;
 113     seq := Eltseq(tr*phi);
 114     phi := GenusOneModel(Q,2,[seq[i]: i in [4..8]]);
 115   end if;
 116   return phi,tr;
 117 end function;
 118 
 119 function RepeatedRoot(F,d);
 120   assert Degree(F) le d;
 121   m := Floor(d/2) + 1;
 122 // Finds the m-fold root of F, where the latter is given as
 123 // a univariate polynomial, but viewed as a binary form 
 124 // of degree d. 
 125   error if F eq 0, "Error: Reduced polynomial should be non-zero.";
 126   vprintf Minimisation,3 : "F = %o\n",PrettyFactorization(F);
 127   if Degree(F) le d - m then 
 128     ans := [1,0];
 129   else 
 130     rts := MultipleRoots(F,m);
 131     if #rts eq 0 then return false,_; end if;
 132     ans := [rts[1],1];
 133   end if;
 134   return true,ans;
 135 end function;
 136 
 137 function PolynomialQ1(phi,p)
 138   k := GF(p);
 139   Z := Integers();
 140   P<t> := PolynomialRing(k);
 141   seq := Eltseq(phi);
 142   if #seq eq 5 then 
 143     if forall{x : x in seq | (Z!x) mod p eq 0} then 
 144       seq := [x/p: x in seq];
 145     end if;
 146     F := &+[(k!seq[i+1])*t^(4-i): i in [0..4]];
 147     dd := 4;
 148   else   
 149     assert p eq 2;
 150     l,m,n,a,b,c,d,e := Explode(ChangeUniverse(seq,k));
 151     dd := 2;
 152     F := l*t^2 + m*t + n;
 153     if F eq 0 then F := b*t^2 + d; end if;
 154     if F eq 0 then 
 155       F := &+[k!(seq[i+4]/2)*t^(4-i): i in [0..4]];
 156       dd := 4;
 157     end if;
 158   end if;
 159   return F,dd; 
 160 end function;
 161 
 162 function SaturateWithCrossTerms(phi)
 163 // Apply y-substitutions to a genus one model of degree 2,
 164 // with the aim of decreasing the level at p = 2.
 165 // Returns (P,Q) with either v(P,Q) = 0, or v(P) > 0 and v(Q) = 1.
 166   Q := Rationals();  
 167   tr := <1,[0,0,0],Id(Q,2)>;
 168   while true do
 169     seq := ChangeUniverse(Eltseq(phi),Integers());
 170     vP := Minimum([Valuation(seq[i],2): i in [1..3]]);
 171     vQ := Minimum([Valuation(seq[i],2): i in [4..8]]);
 172     vPQ := Minimum(2*vP,vQ);
 173     if vPQ eq 0 and forall{i : i in [1,2,3,5,7] | seq[i] mod 2 eq 0} then 
 174       tr1 := <1,[seq[i] mod 2: i in [4,6,8]],Id(Q,2)>;
 175     else
 176       if vPQ lt 2 then break; end if;
 177       tr1 := <1/(2^Floor(vPQ/2)),[0,0,0],Id(Q,2)>; 
 178     end if;
 179     phi := tr1*phi;
 180     tr := tr1*tr;
 181   end while;
 182   vprintf Minimisation,3 : "y-substitution: %o\n",PrettySubstitution(tr);
 183   return phi,tr,1;
 184 end function;
 185 
 186 //////////////////////////////////////////////////////////////
 187 //                      Case n = 3                          //
 188 //////////////////////////////////////////////////////////////
 189 
 190 function SingularLocusThree(V,phi);
 191 // Computes the singular locus of a ternary cubic.
 192 // The answer is returned as a subspace of V = k^3.
 193   F := Equation(phi);
 194   k := CoefficientRing(phi);
 195   R<x,y,z> := PolynomialRing(phi);
 196   error if Discriminant(phi) ne 0,
 197     "Error: Reduced cubic should be singular.";
 198   error if cInvariants(phi) ne [0,0],
 199     "Error: Reduced cubic should be a null-form.";
 200   error if F eq 0, "Error: Reduced cubic should be non-zero.";
 201   eqns := [F] cat [Derivative(F,i): i in [1..3]];
 202   gcd := GCD(eqns);
 203   if TotalDegree(gcd) gt 0 then
 204     f := Factorization(gcd)[1][1];
 205     assert TotalDegree(f) eq 1;
 206     vprintf Minimisation,3 : "Singular locus is { %o = 0 }.\n",f;
 207     U := sub<V|V![MC(f,R.i): i in [1..3]]>;
 208   else
 209     PP := Points(Scheme(Proj(R),eqns));
 210     assert #PP eq 1;
 211     vprintf Minimisation,3 : "Singular locus is { %o }.\n",PP[1];
 212     mat := Matrix(k,3,1,Eltseq(PP[1]));
 213     km := KernelMatrix(mat);
 214     assert Nrows(km) eq 2;
 215     U := sub<V|[V!km[i]: i in [1..2]]>;
 216   end if;
 217   return U;
 218 end function;
 219   
 220 //////////////////////////////////////////////////////////////
 221 //                      Case n = 4                          //
 222 //////////////////////////////////////////////////////////////
 223 
 224 function KernelOfQuadric(V,Q)
 225 // V is a copy of the k-vector space k^d.
 226 // Q is a homogeneous form of degree 2 in k[x_1,..,x_d].
 227 // We compute ker(Q), i.e. the subspace of V where Q and all its
 228 // partial derivatives vanish.
 229   k := BaseField(V);
 230   d := Dimension(V);
 231   mat := MAT(k,Dimension(V),Q);
 232   km := KernelMatrix(mat);
 233   n := Nrows(km);
 234   if Char(k) eq 2 and n ge 1 then
 235     R := PolynomialRing(k,n);
 236     subst := [&+[km[i,j]*R.i: i in [1..n]]: j in [1..d]];
 237     Q1 := Evaluate(Q,subst);
 238     if Q1 ne 0 then  
 239       f := SquareRoot(Q1);
 240       mat1 := Matrix(k,n,1,[MC(f,R.i): i in [1..n]]);
 241       km := KernelMatrix(mat1)*km;
 242     end if;
 243   end if;
 244   return sub<V|[V!km[i]: i in [1..Nrows(km)]]>;
 245 end function;
 246 
 247 function DeterminantPolynomial(phi4,t)
 248 // The polynomial whose roots correspond to singular quadrics
 249 // (with suitable modifications in characteristic 2)
 250   k := CoefficientRing(phi4);
 251   P := Parent(t);
 252   if Char(k) ne 2 then 
 253     M := [MAT(P,4,Q): Q in Equations(phi4)];
 254     return Determinant(t*M[1]+M[2]),4;
 255   else 
 256     phi2 := DoubleGenusOneModel(phi4);
 257     l,m,n,a,b,c,d,e := Explode(Eltseq(phi2));
 258     f := l*t^2 + m*t + n;
 259     return (f ne 0 select f else b*t^2 + d),2;
 260   end if;
 261 end function;
 262 
 263 function SubDiscriminantThree(Q)
 264 // Find a non-zero 3 by 3 "sub-discriminant" of the quadric Q.
 265   R := Parent(Q);
 266   for i in [1..4] do
 267     x,y,z := Explode([R.j : j in [1..4]| j ne i]);
 268     a,b,c := Explode([MC(Q,mon): mon in [x^2,y^2,z^2]]);
 269     d,e,f := Explode([MC(Q,mon): mon in [y*z,z*x,x*y]]);
 270     Delta3 := 4*a*b*c - a*d^2 - b*e^2 - c*f^2 + d*e*f;
 271     if Delta3 ne 0 then return Delta3; end if;
 272   end for;
 273   return 0;
 274 end function;
 275 
 276 function SpanOfSingularLocusInternal(V,quads,F,d) 
 277   assert Degree(F) le d;
 278   m := Floor(d/2) + 1;
 279 // Finds the m-fold root of F, where the latter is given as
 280 // a univariate polynomial, but viewed as a binary form 
 281 // of degree d. Then lets Q1,Q2 be a basis for the pencil
 282 // of quadrics, where Q2 corresponds to this root, and
 283 // returns the linear span of (ker(Q2) meet {Q1 = 0}),
 284 // as a subspace of V = k^4.
 285   if Degree(F) le (d - m) then 
 286     Q2,Q1 := Explode(quads);
 287   else 
 288     rts := MultipleRoots(F,m);
 289     if #rts eq 0 then return false,_; end if;
 290     Q1,Q2 := Explode(quads);   
 291     Q2 := rts[1]*Q1 + Q2;
 292   end if;
 293   k := BaseRing(Universe(quads));
 294   U := KernelOfQuadric(V,Q2);
 295   mat := Matrix(Basis(U));
 296   d := Nrows(mat);
 297   R := PolynomialRing(k,d);
 298   subst := [&+[mat[i,j]*R.i: i in [1..d]]: j in [1..4]];
 299   Q1 := Evaluate(Q1,subst);
 300   if Q1 ne 0 then
 301     ff := Factorization(Q1);
 302     if ff[1,2] eq 2 then // Q1 restricted to ker(Q2) has rank 1
 303       mat1 := Matrix(k,d,1,[MC(ff[1,1],R.i): i in [1..d]]);
 304       mat := KernelMatrix(mat1)*mat;
 305       U := sub<V|[V!mat[i]: i in [1..d-1]]>;
 306     end if;
 307   end if;
 308   return true,U;
 309 end function;
 310 
 311 function SpanOfSingularLocus(V,phi);
 312 // Computes the span of the singular locus of a quadric intersection.
 313 // The answer is returned as a subspace of V = k^4.
 314   quads := Equations(phi);
 315   k := CoefficientRing(phi);
 316   R := PolynomialRing(phi);
 317   error if Discriminant(phi) ne 0,
 318     "Error: Reduced quadric intersection should be singular.";
 319   error if cInvariants(phi) ne [0,0],
 320     "Error: Reduced quadric intersection should be a null-form.";
 321   error if GCD(quads) ne 1,
 322     "Error: Reduced quadrics should be coprime.";
 323   U := &meet[KernelOfQuadric(V,q): q in quads];
 324   P<t> := PolynomialRing(k);
 325   common_nullity := Dimension(U);
 326   vprintf Minimisation,3 : "Common Nullity = %o\n",common_nullity;
 327   error if common_nullity gt 2, 
 328     "Error : Common nullity =",common_nullity," Please report";
 329   case common_nullity :
 330     when 0:
 331       F,d := DeterminantPolynomial(phi,t);
 332       vprintf Minimisation,3 : "F = %o\n",PrettyFactorization(F);
 333       if F eq 0 then 
 334         Q1,Q2 := Explode(quads);
 335         UU := [KernelOfQuadric(V,q): q in [Q1,Q2,Q1+Q2]];
 336         U := &+[U : U in UU | Dimension(U) eq 1];
 337       else
 338         flag,U := SpanOfSingularLocusInternal(V,quads,F,d);
 339         assert flag;
 340       end if;
 341     when 1 :
 342       Q1,Q2 := Explode(Equations(ChangeRing(phi,P)));
 343       F := SubDiscriminantThree(t*Q1+Q2);
 344       error if F eq 0, 
 345         "Error: Generic rank is less than 3. Please report.";
 346       vprintf Minimisation,3 : "F = %o\n",PrettyFactorization(F);
 347       flag,Unew := SpanOfSingularLocusInternal(V,quads,F,3);
 348       if flag then U := Unew; end if;
 349   end case;
 350   vprintf Minimisation,3 : 
 351     "Span of singular locus has basis \n%o\n",Matrix(Basis(U));  
 352   return U;
 353 end function;
 354 
 355 //////////////////////////////////////////////////////////////
 356 //                      Case n = 5                          //
 357 //////////////////////////////////////////////////////////////
 358 
 359 function GetQuadrics(phi,p,irreg)
 360   cc := p^irreg;
 361   quads := [(1/cc)*q : q in Equations(phi)];
 362   if irreg gt 0 then 
 363     R := PolynomialRing(phi);
 364     Q := Rationals();
 365     M := Matrix(phi);
 366     mat := Matrix(Q,5,25,[[MC(M[i,r],R.s): r,s in [1..5]]: i in [1..5]]);
 367     mat := MySaturate(mat);
 368     quads := [&+[mat[i,j]*quads[i]: i in [1..5]]: j in [1..5]];
 369   end if;
 370   return quads;
 371 end function;
 372 
 373 function SpanOfSingularLocusFive(quads)
 374   R := Universe(quads);
 375   k := CoefficientRing(R);
 376   J := Matrix(R,5,5,[Deriv(q,i): i in [1..5],q in quads]);
 377   I := ideal<R|quads cat Minors(J,3)>;
 378   ff := MinimalBasis(Radical(I));
 379   ff := [f : f in ff | TotalDegree(f) eq 1];
 380   mat := Matrix(k,#ff,5,[MC(f,R.i): i in [1..5],f in ff]);
 381   vprintf Minimisation,3 :
 382     "Singular locus defined by linear forms\n%o\n",mat;
 383   return mat;
 384 end function;
 385 
 386 //////////////////////////////////////////////////////////////
 387 //                   Cases n = 2,3,4,5                      //
 388 //////////////////////////////////////////////////////////////
 389 
 390 function SaturateModel(phi)
 391 // n = 2  : makes GCD(coeffs) square-free
 392 // n = 3  : makes GCD(coeffs) = 1
 393 // n = 4  : makes quadrics linearly indep. mod p for all p.
 394 // n = 5  : makes Pfaffians linearly indep. mod p for all p.
 395   Q := Rationals();  
 396   case Degree(phi) :
 397     when 2:
 398       if #Eltseq(phi) eq 5 then 
 399         d := RationalGCD(Eltseq(phi));
 400         a1,a2 := SquareFree(Numerator(d));
 401         b1,b2 := SquareFree(Denominator(d));
 402         tr := <b1*b2/a2,Id(Q,2)>;
 403       else
 404         return SaturateWithCrossTerms(phi);
 405       end if;
 406     when 3:
 407       tr := <1/RationalGCD(Eltseq(phi)),Id(Q,3)>;
 408     when 4:  
 409       tr := <MySaturate(CoeffMatrix(Equations(phi),2))^(-1),Id(Q,4)>;
 410       vprintf Minimisation,3 : 
 411         "Change of basis for space of quadrics \n%o\n",tr[1];
 412     when 5:
 413       mat := MySaturate(CoeffMatrix(Equations(phi),2));
 414       d1 := RationalGCD(Eltseq(mat));
 415       mat := (1/d1)*mat;
 416       tr := <Transpose(mat),Id(Q,5)>;
 417       phi := tr*phi;
 418       d2 := RationalGCD(Eltseq(phi));
 419       irreg := Integers()!(d1*Determinant(mat)/d2^2);
 420       tr1 := <Id(Q,5),(1/d2)*Id(Q,5)>;
 421       vprintf Minimisation,3 : 
 422         "Change of basis for space of quadrics \n%o\n",tr[1];
 423       vprintf Minimisation,3 : 
 424         "Irregularity = %o\n",Factorization(irreg);
 425       return tr1*phi,tr1*tr,irreg;
 426   end case;
 427   return tr*phi,tr,1;
 428 end function;
 429 
 430 function MinimiseInternal(phi,tr,leveldata,tflag)
 431   Z := Integers();
 432   Q := Rationals();
 433   n := Degree(phi);
 434   assert n in {2,3,4,5};
 435   maxsteps := case< n | 2:2, 3:4, 4:5, 5:6, default:0 >;// a guess when n = 5
 436   idqtrans := case< n | 4:Id(Q,2), 5:Id(Q,5), default:1 >;
 437   vprintf Minimisation,1 : "Level = %o\n",leveldata;
 438   vprintf Minimisation,2 : 
 439      "Notation: \"/\" = level decreases  \".\" = level preserved\n";
 440   failedprimes := [];
 441   for pr in leveldata do 
 442     p,level := Explode(pr);
 443     vprintf Minimisation,1 : "Minimising at p = %o (level = %o)\n",p,level;
 444     k := GF(p);
 445     V := VectorSpace(k,n);
 446     R := PolynomialRing(k,n);
 447     if n eq 5 then irreg := pr[3]; end if;
 448     nsteps := 0;
 449     oldphi := phi;
 450     tr0 := <idqtrans,Id(Q,n)>;
 451     IndentPush();
 452     while level gt 0 do
 453       if n in {3,4} then 
 454         seq := ChangeUniverse(Eltseq(phi),k);
 455         phibar := GenusOneModel(n,seq:PolyRing:=R);
 456       end if;
 457       case n :
 458         when 2 :
 459           F,d := PolynomialQ1(phi,p);
 460           flag,rr := RepeatedRoot(F,d);
 461           if not flag then 
 462             vprintf Minimisation,2 : 
 463               " no triple root => minimal at level %o\n",level; 
 464             error if p ne 2, "Error:  Null-form without triple root.";
 465             failedprimes cat:= [<p,level>]; 
 466             break;
 467           end if;
 468           mat := Matrix(k,1,2,[-rr[2],rr[1]]);
 469         when 3 :         
 470           U := SingularLocusThree(V,phibar);
 471           mat := Matrix(Basis(U));
 472         when 4 :
 473           f := GCD(Equations(phibar));
 474           if TotalDegree(f) eq 1 then
 475             vprintf Minimisation,3 : "GCD = %o\n",f;
 476             mat := Matrix(k,1,4,[MC(f,R.i): i in [1..4]]);
 477           else
 478             U1 := SpanOfSingularLocus(V,phibar);
 479             mat := Matrix(Basis(U1));
 480             mat := KernelMatrix(Transpose(mat));
 481           end if;
 482         when 5 :
 483 	  quads := ChangeUniverse(GetQuadrics(phi,p,irreg),R);
 484           vprintf Minimisation,3 :
 485             "Space of quadrics has dimension %o\n",Rank(CoeffMatrix(quads,2));
 486           mat := SpanOfSingularLocusFive(quads);
 487       end case;
 488       q := Nrows(mat);
 489       assert q lt n;
 490       _,_,mat := SmithForm(ChangeRing(mat,Z));
 491       dd := [i le q select p else 1 : i in [1..n]];
 492       M := Diag(Q,dd)*Transpose(mat); 
 493       vprintf Minimisation,3 : 
 494         "Change of co-ordinates on P^%o\n%o\n",n-1,M;
 495       tr1 := <idqtrans,M>;
 496       phi := tr1*phi; 
 497       if tflag then tr0 := tr1*tr0; end if;
 498       phi,tr1,globirreg := SaturateModel(phi);
 499       if tflag then tr0 := tr1*tr0; end if;
 500       ch := q + Valuation(ScalingFactor(n,tr1),p:Check:=false);
 501       if n eq 5 then 
 502         newirreg := Valuation(globirreg,p:Check:=false);
 503         ch := ch - 2*(newirreg - irreg);
 504         irreg := newirreg;
 505       end if;
 506       level +:= ch;
 507       error if ch gt 0, "Error : The level increased. Please report.";
 508       vprintf Minimisation,3 : "Level decreases by %o   (",-ch;
 509       if ch lt 0 then
 510         nsteps := 0;
 511         oldphi := phi;
 512         if tflag then tr := tr0*tr; end if;
 513         tr0 := <idqtrans,Id(Q,n)>;
 514         vprintf Minimisation,2 :"/";
 515       else 
 516         nsteps +:= 1;
 517         vprintf Minimisation,2 :".";
 518       end if;
 519       vprintf Minimisation,3 : ")\n";
 520       if nsteps ge maxsteps then 
 521         vprintf Minimisation,2 : " => minimal at level %o\n",level;
 522         phi := oldphi;
 523         failedprimes cat:= [<p,level>];
 524         break;
 525       end if;
 526     end while;
 527     if level eq 0 then 
 528       vprintf Minimisation,2 : "  => level 0 \n";
 529     end if;
 530     IndentPop();
 531   end for;
 532   return phi,tr,failedprimes;
 533 end function;
 534 
 535 function MyFactorization(n,plist)
 536   n := Integers()!n;
 537   if #plist ne 0 then 
 538     plist := Sort(SetToSequence(SequenceToSet(plist)));
 539     return [<p,Valuation(n,p:Check:=false)>: p in plist];
 540   else 
 541     vprintf Minimisation,1 : "Factoring the level\n";
 542     return Factorization(n);
 543   end if;
 544 end function; 
 545 
 546 intrinsic Minimise2008(model::ModelG1:
 547   Transformation := true,
 548   CrossTerms := false,
 549   UsePrimes := [] ) 
 550   -> ModelG1, Tup, SetEnum
 551 {Minimises a genus one model of degree 2,3,4 or 5 over Q.
 552 Also returns the transformation (unless Transformation = false)
 553 and the set of primes p where the minimal model has positive level. 
 554 Except possibly when p = 2 (for models of degree 2 with CrossTerms = false), 
 555 these are the primes where the model is not Q_p^nr-soluble. The current
 556 version for models of degree 5 is not yet proven to work in all cases.}
 557   R := CoefficientRing(model);
 558   require (R cmpeq Rationals()) or (R cmpeq Integers()) 
 559     : "Model must be defined over the rationals";
 560   require Discriminant(model) ne 0 : "Model is singular";
 561   n := Degree(model);
 562   phi := ChangeRing(model,Rationals());
 563   vprintf Minimisation, 1 :"Minimising a genus one model of degree %o\n",n;
 564   require n in {2,3,4,5} : "Model must have degree 2,3,4 or 5.";
 565   if n eq 2 then 
 566     phi,tr1 := RemoveCrossTerms(phi); 
 567     phi,tr2 := SaturateModel(phi);
 568     tr := tr2*tr1;
 569   else
 570     phi,tr,irreg := SaturateModel(phi);
 571   end if;
 572   if n eq 5 then 
 573     globlev := GlobalLevel(phi)/irreg^2;
 574     ff := MyFactorization(globlev,UsePrimes);
 575     leveldata := [<f[1],f[2],Valuation(irreg,f[1]:Check:=false)>: f in ff];
 576   else
 577     leveldata := MyFactorization(GlobalLevel(phi),UsePrimes);
 578   end if;
 579   phi,tr,leveldata := MinimiseInternal(phi,tr,leveldata,Transformation);
 580   failedprimes := {f[1]: f in leveldata};
 581   if n eq 2 and CrossTerms then 
 582     vprintf Minimisation,1 : "Final stage : introducing cross terms\n";
 583     phi := <1,[0,0,0],Id(Rationals(),2)>*phi; 
 584     ld2 := [pr : pr in leveldata | pr[1] eq 2];
 585     level := #ld2 ne 0 select ld2[1][2] else 0;
 586     phi,tr1 := SaturateModel(phi);
 587     tr := tr1*tr;
 588     level +:= Valuation(ScalingFactor(2,tr1),2);
 589     vprintf Minimisation,1 : 
 590       "Making a y-substutition => level = %o\n",level;
 591     if level gt 0 then 
 592       phi,tr,ld2 := MinimiseInternal(phi,tr,[<2,level>],Transformation);
 593       level := #ld2 ne 0 select ld2[1][2] else 0;
 594     end if;
 595     if level eq 0 then Exclude(~failedprimes,2); end if;
 596   end if;
 597   if #failedprimes gt 0 then
 598     vprintf Minimisation,1 : 
 599       "Minimal model has positive level at primes %o.\n",failedprimes;
 600   else
 601     vprintf Minimisation,1 : "Model has level 0 at all primes.\n";
 602   end if;
 603   if Transformation then
 604     return phi,tr,failedprimes;
 605   else
 606     return phi,_,failedprimes;
 607   end if;
 608 end intrinsic;
 609 
 610 /////////////////////////////////////////////////////////////////////////
 611 //       Alternative treatment of (globally) soluble models            //
 612 /////////////////////////////////////////////////////////////////////////
 613 
 614 // The following is a bug-fix for ApplyTransformation (see g1invariants.m)
 615 // The modification required is 
 616 //      R0 := BaseRing(R);
 617 //      subst := subst cat [R0!(1/mu)*y + r[1]*x^2 + r[2]*x*z + r[3]*z^2]; 
 618 // This is to cope with the fact Magma doesn't recognise 1/1 as an integer!
 619 
 620 function ApplyTransformation2(g,model)
 621 // {Apply the transformation g to the given genus one model.}
 622   n := model`Degree;
 623   flag,lambda := IsTransformation(n,g);
 624 //  require flag : "g must be a transformation of genus one curves of degree n (matching the given model).";
 625   if n in {2,3} then
 626     f := model`Equation;
 627     R := Parent(f);
 628     if #g eq 2 then 
 629       mu,B := Explode(g);
 630       r := [0,0,0];
 631     else
 632       mu,r,B := Explode(g);
 633     end if;
 634     if VariableWeights(R) eq [1,1] and #g eq 3 then
 635       f := Equations(model)[1];
 636       R := Parent(f);
 637     end if;
 638     subst := [&+[B[i,j]*R.i: i in [1..n]]: j in [1..n]];
 639     if VariableWeights(R) eq [1,1,2] then 
 640       R<x,z,y> := Parent(f);
 641 // the modification is here
 642       R0 := BaseRing(R);
 643       subst := subst cat [R0!(1/mu)*y + r[1]*x^2 + r[2]*x*z + r[3]*z^2]; 
 644     end if;
 645     f_new := (n eq 2 select mu^2 else mu)*Evaluate(f,subst);
 646     return GenusOneModel(f_new);
 647   elif n eq 4 then
 648     R := PolynomialRing(model);
 649     phi := model`Equations;
 650     A,B := Explode(g);
 651     subst:=[&+[B[i,j]*R.i :i in [1..4]]: j in [1..4]];
 652     phi_B := [ Evaluate(phi[j],subst) : j in [1..2] ];  
 653     return GenusOneModel([ A[1,1]*phi_B[1] + A[1,2]*phi_B[2], A[2,1]*phi_B[1] + A[2,2]*phi_B[2] ]);
 654   else 
 655     R := PolynomialRing(model);
 656     phi := model`DefiningMatrix;
 657     A,B := Explode(g);
 658     A := MatrixRing(R,5)!A; 
 659     subst:=[&+[B[i,j]*R.i :i in [1..5]]: j in [1..5]];
 660     phi_B := Parent(phi)![Evaluate(phi[i,j],subst):j in [1..5],i in [1..5]];
 661     phi1 := Parent(phi)!(A*phi_B*Transpose(A));
 662     for i in [1..5] do phi1[i,i] := 0; end for;
 663 for i in [1..5] do for j in [1..i-1] do phi1[i,j] := -phi1[j,i]; 
 664 end for; end for;
 665     return GenusOneModel(phi1);
 666 
 667 //    _,ans := IsGenusOneModel( Parent(phi)!(A*phi_B*Transpose(A)) );
 668 //    return ans;
 669   end if;
 670 end function;
 671 
 672 function IntegralGenusOneModel(eqns)
 673   RZ := Universe(eqns);
 674   case Rank(RZ) :
 675     when 2 : 
 676       P,Q := Explode(eqns);
 677       coeffs1 := [MC(P,RZ.1^(2-i)*RZ.2^i): i in [0..2]];
 678       coeffs2 := [MC(Q,RZ.1^(4-i)*RZ.2^i): i in [0..4]];
 679       model := GenusOneModel(2,coeffs1 cat coeffs2);
 680     when 3 : 
 681       model := GenusOneModel(eqns[1]);
 682     when 4 : 
 683       model := GenusOneModel(eqns);
 684     when 5 :
 685       RQ := PolynomialRing(Rationals(),5);
 686       eqns := ChangeUniverse(eqns,RQ);
 687       model := GenusOneModel(eqns);
 688       coeffs := PRIM(Eltseq(model));
 689       model := GenusOneModel(5,coeffs:PolyRing:=RZ);
 690   end case;
 691   return model;
 692 end function;
 693 
 694 function WeierstrassEquations(E,n)
 695   assert n in {3,4,5,6};
 696   R<x,y,z> := PolynomialRing(Integers(),3); 
 697   if n eq 3 then 
 698     return [Evaluate(Equation(E),[y,z,x])],[z,x,y];
 699   end if;
 700   ainv := ChangeUniverse(aInvariants(E),Integers());
 701   a1,a2,a3,a4,a6 := Explode(ainv);
 702   case n :
 703     when 4 :
 704       RR<x0,x2,x3,x4> := PolynomialRing(Integers(),4);
 705       ff := [z^2,x*z,y*z,x^2];
 706       quads := [ x0*x4 - x2^2,
 707 	x3^2 + a1*x2*x3 + a3*x0*x3 - (x2*x4 + a2*x0*x4 + a4*x0*x2 + a6*x0^2)];
 708     when 5 :
 709       RR<x0,x2,x3,x4,x5> := PolynomialRing(Integers(),5);
 710       ff := [z^2,x*z,y*z,x^2,x*y];
 711       quads := [ x0*x4 - x2^2,x0*x5 - x2*x3,x2*x5 - x3*x4,
 712 	x3^2 + a1*x0*x5 + a3*x0*x3 - (x2*x4 + a2*x0*x4 + a4*x0*x2 + a6*x0^2),
 713         x3*x5 + a1*x2*x5 + a3*x2*x3 - (x4^2 + a2*x2*x4 + a4*x2^2 + a6*x0*x2)];
 714     when 6 :
 715       RR<x0,x2,x3,x4,x5,x6> := PolynomialRing(Integers(),6);
 716       ff := [z^3,x*z^2,y*z^2,x^2*z,x*y*z,x^3];
 717       quads := [ x0*x4 - x2^2,x0*x5 - x2*x3,x0*x6 - x2*x4,
 718         x2*x5 - x3*x4,x2*x6 - x4^2,x3*x6 - x4*x5,
 719         x3^2 + a1*x0*x5 + a3*x0*x3 - x0*(x6 + a2*x4 + a4*x2 + a6*x0),
 720         x3*x5 + a1*x2*x5 + a3*x2*x3 - x2*(x6 + a2*x4 + a4*x2 + a6*x0),
 721 	x5^2 + a1*x3*x6 + a3*x2*x5 - x4*(x6 + a2*x4 + a4*x2 + a6*x0)];
 722   end case;
 723   return quads,ff;
 724 end function;
 725 
 726 function GenusOneModelFromPoint(n,eqns,pt)
 727   assert forall{f: f in eqns| Evaluate(f,pt) eq 0};
 728   target := Proj(Universe(eqns))!([1] cat [0: i in [1..n]]);
 729   vprintf Minimisation,1 : "Moving P to %o\n",target;
 730   _,_,mat := SmithForm(Matrix(1,n+1,pt));
 731   matdual := Transpose(mat^(-1));
 732   eqns := [q^matdual : q in eqns];
 733   pt := [mat[n+1,j]: j in [2..n+1]]; 
 734   vprintf Minimisation,1 : "Projecting away from %o\n",target;
 735   R := PolynomialRing(Integers(),n);
 736   S<t> := PolynomialRing(R);
 737   subst := [t] cat [R.i : i in [1..n]];
 738   eqns := [Evaluate(q,subst): q in eqns];
 739   if n eq 2 then 
 740     f1,f2,f3 := Explode([Coefficient(eqns[1],i): i in [2,1,0]]);
 741     eqns := [f2,-f1*f3];
 742     yy := mat[3,1]*Evaluate(f1,pt);
 743   else
 744     gg := [Coefficient(f,1): f in eqns];
 745     hh := [Coefficient(f,0): f in eqns];
 746     if n eq 3 then
 747       eqns := [gg[1]*hh[2]-gg[2]*hh[1]];  
 748     else
 749       km := KernelMatrix(CoeffMatrix(gg,1)); 
 750       eqns := [&+[km[i,j]*hh[j]: j in [1..#hh]]: i in [1..Nrows(km)]];
 751     end if;
 752   end if;
 753   vprint Minimisation,1 : "Performing some ad hoc reduction";
 754   submat := Matrix(n,n+1,[mat[i,j]: i in [1..n+1],j in [2..n+1]]);
 755   _,T := LLL(submat);
 756   Tinv := T^(-1);
 757   eqns := [q^Tinv : q in eqns];
 758   pt := [&+[T[i,j]*pt[j]: j in [1..n]]: i in [1..n]];
 759   if n in {4,5} then 
 760     _,U := LLL(CoeffMatrix(eqns,2));
 761     eqns := [&+[U[i,j]*eqns[j]: j in [1..#eqns]]: i in [1..#eqns]];
 762   end if;
 763   vprintf Minimisation,1 : "Re-writing as a genus one model\n";
 764   model := IntegralGenusOneModel(eqns);
 765   if n eq 2 then 
 766     rr := [-Floor(seq[i]/2): i in [1..3]] where seq is Eltseq(model);
 767 // a bug fix to ApplyTransformation was required (see above)
 768     model := ApplyTransformation2(<1,rr,Id(Integers(),2)>,model);
 769     yy := yy - rr[1]*pt[1]^2 - rr[2]*pt[1]*pt[2] - rr[3]*pt[2]^2;
 770     pt := pt cat [yy];
 771   end if;
 772   if n eq 5 then
 773     quads := Equations(model);
 774     assert (eqns eq quads) or (eqns eq [-q: q in quads]);
 775   end if;
 776   return model,pt;
 777 end function;
 778 
 779 intrinsic GenusOneModel(n::RngIntElt,P::PtEll:NoReduction := false)
 780   -> ModelG1, Pt
 781 {Given a point P on an elliptic curve E over Q, and n = 2,3,4 or 5, maps E to P^\{n-1\} via the complete linear system |(n-1).O+P| and returns the corresponding genus one model of degree n. The formulae used give a genus one model that is minimised and is close to being reduced. Finally the function Reduce is called (unless NoReduction = true). The second returned value is a point on the covering curve that maps down to P (or -P).}
 782   E := Curve(P);
 783   K := BaseRing(E);
 784   require K cmpeq Rationals() 
 785     : "Elliptic curve must be defined over the rationals";
 786   K := Ring(Parent(P));
 787   require K cmpeq Rationals() 
 788     : "Point must be defined over the rationals";
 789   vprintf Minimisation, 1 :"Computing a minimal genus one model of degree %o\n",n;
 790   require n in {2,3,4,5} : "Model must have degree 2,3,4 or 5.";
 791   if P eq E!0 then 
 792     error "Point given is the identity - use GenusOneModel(n,E) instead";
 793   end if;
 794   E,iso := MinimalModel(E);
 795   P := iso(P);
 796   eqns,ff := WeierstrassEquations(E,n+1);
 797   pt := PRIM([Evaluate(f,Eltseq(P)): f in ff]);
 798   vprintf Minimisation, 1 :
 799     "We embed E in P^%o via the linear system |%o.0|\n",n,n+1;
 800   model,pt := GenusOneModelFromPoint(n,eqns,pt);
 801   model := ChangeRing(model,Rationals());  
 802   vprintf Minimisation, 1 : "Model has coefficients\n%o\n",Eltseq(model);
 803   if not NoReduction then 
 804     model,tr := Reduce2008(model);
 805     if n eq 2 then yy := pt[3]; end if;
 806     T := Transpose(tr[#tr])^(-1);
 807     pt := [&+[T[i,j]*pt[j]: j in [1..n]]: i in [1..n]];
 808     cc := RationalGCD(pt);
 809     pt := [x/cc: x in pt];
 810     if n eq 2 then
 811       rr := tr[2];   
 812       yy := yy/cc^2 - rr[1]*pt[1]^2 - rr[2]*pt[1]*pt[2] - rr[3]*pt[2]^2;
 813       pt := pt cat [yy];
 814     end if;
 815   end if;
 816   vprintf Minimisation, 1 :"Checking the invariants\n";
 817   error if cInvariants(model) ne cInvariants(E),
 818     "Model is not minimal. Please report";
 819   P := Curve(model)!pt; 
 820   return model,P;
 821 end intrinsic;